原子范数最小化(Atomic Norm Minimization)

原子范数最小化(Atomic Norm Minimization)

$l_0$-原子范数

原子集
A = { a ( θ , ϕ ) = a ( θ ) ϕ : θ ∈ [ − 9 0 ∘ , 9 0 ∘ ] , ϕ ∈ C , ∣ ϕ ∣ = 1 } \mathcal{A} = \{\boldsymbol{a}(\theta,\phi) = \boldsymbol{a}(\theta)\phi : \theta\in [-90^\circ,90^\circ], \phi \in \mathbb{C}, |\phi| = 1 \} A={a(θ,ϕ)=a(θ)ϕ:θ∈[−90∘,90∘],ϕ∈C,∣ϕ∣=1}
对于无噪声条件下的单快拍接收信号。
$$z=\mathbf{A}(\theta )\mathbf{s}$$
$l_0$范数的定义
∥ z ∥ A , 0 = inf ⁡ s k , θ k { K : z = ∑ k = 1 K a ( θ k ) s k , θ k ∈ [ − 9 0 ∘ , 9 0 ∘ ] } \|\mathit{z}\|_{\mathcal{A},0} = \inf_{s_k,\theta_k}\{\mathcal{K}:z = \sum_{k = 1}^{\mathcal{K}} \boldsymbol{a} \left(\theta_k\right)s_k, \theta_k \in [-90^\circ,90^\circ] \} ∥z∥A,0​=sk​,θk​inf​{K:z=k=1∑K​a(θk​)sk​,θk​∈[−90∘,90∘]}
$l_0$范数指的是，组成向量$z$的最小原子的个数。DOA估计反映在该问题中就是找出原子个数最小的情况下的方向向量的组合。

而对于稀疏阵列而言，信号$z$往往是残缺的，不完整的。信号恢复的过程，可以理解为：填充向量$z$，使其在原子集$\mathcal{A}$下，原子范数最小，即：
$$\underset{z}{\min }\Vert z{\Vert }_{\mathcal{A},0}$$

最小化Toeplitz矩阵的秩

引理1

对于任意的半正定矩阵$\mathcal{T}\left(z\right)\in {\mathbb{C}}^{L\times L}$，如果$r=rank(\mathcal{T}(z))\le L-1$，则Hermitian Toeplitz矩阵能够被唯一范德蒙分解为：
$$\mathcal{T}\left(z\right)=\sum _{k=1}^r{p}_k\mathbf{a}({\theta }_k){\mathbf{a}}^H({\theta }_k)=\mathbf{A}(\theta )\mathbf{P}{\mathbf{A}}^H(\theta )$$
其中，$\mathbf{A}\in {\mathbb{C}}^{L\times r}$是范德蒙矩阵，它的列视为来自不同方向的信号的方向向量，$\mathbf{P}\in {\mathbb{C}}^{r\times r}$是对角矩阵。

引理2

考虑以下的式子：
[ x z H z T ( z ) ] ⪰ 0

$$\begin{bmatrix} x & \boldsymbol{z}^H \\ \boldsymbol{z} & \mathcal{T}\left(z\right) \end{bmatrix}$$ \succeq 0 [xz​zHT(z)​]⪰0
其中，$x$是一个待优化的变量。这个式子对应了以下的结论：

1. $\mathcal{T}\left(z\right)⪰0$；
2. $\mathbf{z}$ 一定位于Toeplitz矩阵 $\mathcal{T}\left(z\right)$ 的列向量中。

问题的转化

由于上面的结论2，所以$\mathbf{z}$一定可以被$\mathbf{A}(\theta )$中的$r$个子向量线性表出。而这$r$个子向量就是矩阵$\mathcal{T}\left(z\right)$范德蒙分解后的$r$个特征向量。因此，向量$\mathbf{z}$的$l_0$范数最小化就是$\mathcal{T}\left(z\right)$矩阵秩的最小化。问题即可转换为：
min ⁡ x , z r a n k ( T ( z ) ) s . t [ x z H z T ( z ) ] ⪰ 0 \min_{x,z} rank(\mathcal{T}\left(z\right)) \\ s.t \quad

$$\begin{bmatrix} x & \boldsymbol{z}^H \\ \boldsymbol{z} & \mathcal{T}\left(z\right) \end{bmatrix}$$ \succeq 0 x,zmin​rank(T(z))s.t[xz​zHT(z)​]⪰0
然而，该优化问题并不是一个凸优化，无法直接求解。

转换为凸优化问题

凸松弛(Relaxation)

解释：通过对问题限制条件的松弛，将原非凸问题转换为凸优化问题。注：松弛原问题，只能得到一个可行域更大的问题，如果原问题是求最小，则松弛后的问题的最优值一定小于等于原问题的最优值，这也是一种给出下界的方法。松弛不仅仅用于整数约束，只要利于将定义域非凸变为凸集即可。

例如，$l_0$范数（非凸）的最紧凸松弛为$l_1$范数（凸），即：将向量中非零元素的个数转为元素绝对值的和。

那么，类似的，可以将$l_0$-原子范数松弛为$l_1$-原子范数。或者，将$rank(\mathcal{T}\left(z\right))$松弛为核范数$\Vert \mathcal{T}\left(z\right){\Vert }_{\ast }$。

$l_1$-原子范数最小化

$l_1$-原子范数的定义如下：
∥ z ∥ A , 1 = inf ⁡ s k , θ { ∑ k ∣ s k ∣ : z = ∑ k = 1 K a ( θ k ) s k , θ k ∈ T } \|z\|_{\mathcal{A},1} = \inf_{s_k,\theta}\{\sum_{k}{|s_k|}:z=\sum_{k=1}^{\mathcal{K}}\boldsymbol{a}(\theta_k)s_k,\theta_k\in\mathbb{T}\} ∥z∥A,1​=sk​,θinf​{k∑​∣sk​∣:z=k=1∑K​a(θk​)sk​,θk​∈T}
从形式上看，问题已经转换成：恢复一个单快拍信号向量$z$，使得其$l_1$-原子范数最小，即：
$$\underset{z}{\min }\Vert z{\Vert }_{\mathcal{A},1}$$
此时，最小化$l_1$范数依然看起来十分困难，所以再次利用范德蒙分解，将该问题转化为：
min ⁡ x , z 1 1 2 ( x + z 1 ) s . t [ x z H z T ( z ) ] ⪰ 0 \min_{x,z_1} \frac{1}{2}(x+z_1) \\ s.t \quad

$$\begin{bmatrix} x & \boldsymbol{z}^H \\ \boldsymbol{z} & \mathcal{T}\left(z\right) \end{bmatrix}$$ \succeq 0 x,z1​min​21​(x+z1​)s.t[xz​zHT(z)​]⪰0
其中，$z_1$表示单快拍向量$z$中的第一个元素。至此，该问题可以用CVX工具箱来求解了。

核范数最小化

核范数（矩阵奇异值的和）的定义如下：
$$\Vert \mathcal{T}(z){\Vert }_{\ast }=Tr\left(\sqrt{{\mathcal{T}}^H(z)\mathcal{T}(z)}\right)$$
核范数$\Vert \mathcal{T}(z){\Vert }_{\ast }$是一个凸函数，证明如下：

由于$\Vert \mathcal{T}(z){\Vert }_{\ast }$和$\Vert \mathcal{T}(z){\Vert }_2$为一对对偶范数，因此，若$\Vert \mathcal{T}(z){\Vert }_2$为凸函数，那么$\Vert \mathcal{T}(z){\Vert }_{\ast }$即为凸函数。
∥ T ∥ 2 = sup ⁡ u , v { u H T v : ∥ u ∥ 2 = 1 , ∥ v ∥ 2 = 1 } \|\mathcal{T}\|_2 = \sup_{\boldsymbol{u},\boldsymbol{v}}\{\boldsymbol{u}^H\mathcal{T}\boldsymbol{v} : \|\mathcal{u}\|_2=1,\|\mathcal{v}\|_2=1\} ∥T∥2​=u,vsup​{uHTv:∥u∥2​=1,∥v∥2​=1}
向量$\mathbf{u}$和$\mathbf{v}$为各个奇异值对应的向量，${\mathbf{u}}^H\mathcal{T}\mathbf{v}$是对$\mathcal{T}$的仿射变换，因此凸。由下面”逐点上确界“的性质，所以$\Vert \mathcal{T}{\Vert }_2$凸。
f ( x ) = sup ⁡ y , z { g ( x , y , z ) : y ∈ A y , z ∈ A z } f(x) = \sup_{y,z}\{g(x,y,z):y\in\mathcal{A}_y,z\in\mathcal{A}_z\} f(x)=y,zsup​{g(x,y,z):y∈Ay​,z∈Az​}
$g$为凸函数，那么$f$亦为凸函数，因而$\Vert \mathcal{T}(z){\Vert }_{\ast }$的凸性得证。

而取矩阵秩的操作的最紧凸松弛为核范数，所以，可以将之前的非凸问题转换成：
min ⁡ x , z ∥ T ( z ) ∥ ∗ s . t [ x z H z T ( z ) ] ⪰ 0 \min_{x,z}\|\mathcal{T}(z)\|_* \\ s.t \quad

$$\begin{bmatrix} x & \boldsymbol{z}^H \\ \boldsymbol{z} & \mathcal{T}\left(z\right) \end{bmatrix}$$ \succeq 0 x,zmin​∥T(z)∥∗​s.t[xz​zHT(z)​]⪰0
该问题也可以通过CVX工具箱来求解。

变换为矩阵的迹

由于半正定的Hermitian Toeplitz矩阵$\Vert \mathcal{T}(z)\Vert$，所以其存在唯一的范德蒙分解：
$$\mathcal{T}\left(z\right)=\sum _{k=1}^r{p}_k\mathbf{a}({\theta }_k){\mathbf{a}}^H({\theta }_k)$$
取该矩阵的迹：
$$\begin{gathered} Tr(\mathcal{T}\left(z\right))=L\sum _{k=1}^r{p}_k \\ \Vert z{\Vert }_{\mathcal{A},1}=\sum _{k=1}^r{p}_k=Tr(\mathcal{T}\left(z\right))/L \end{gathered}$$
上式中，第一个等号可以理解为：由于方向向量$\mathbf{a}({\theta }_k),k=1,2,\dots$ 彼此相互独立，所以只要确定了信号向量$\mathbf{z}$存在由$\mathbf{a}({\theta }_k)$张成的向量空间中，那么这一组系数$p_k,k=1,2,\dots ,r$便可唯一确定。因此，取下界操作$\inf$就可无视。

因此，该凸问题就转换为：
$$\begin{gathered} \underset{z}{\min }Tr(\mathcal{T}\left(z\right)) \\ s.t\mathcal{T}\left(z\right)⪰0; \end{gathered}$$
该问题也可以通过CVX工具箱来求解。

稀疏阵列DOA估计

step1 通过图(a)中实际布置的天线阵列接收数据$\mathbf{X}$，计算接收阵列的协方差矩阵${\mathbf{R}}_x$;

step2 向量化${\mathbf{R}}_x$，并排序去冗余，得到图(b)中等效的单快拍信号向量${\mathbf{z}}_{ex}$，包括的位置-12~12，空洞处用0元素代替;

step3 用${\mathbf{z}}_{ex}$构造Toeplitz矩阵${\mathbf{R}}_v$，另外，定义与${\mathbf{R}}_v$同等大小的二进制矩阵$\mathbf{G}$，用来反映${\mathbf{R}}_v$有值的位置，$\mathbf{z}$代表待恢复的单快拍信号(对应的位置0~12);

step4 用CVX，求解以下任一优化问题，获取重构的Toeplitz矩阵$\mathcal{T}\left(z\right)$;

仿真结果

信噪比$20dB$，快拍数$200$，来波方向$-50^{\circ }$至$50^{\circ }$，间隔$5^{\circ }$，共21个信源。

仿真使用的是增广展开的互质阵列（Augmented-Expanding Co-Prime Array），共10个阵元（$M=3,N=5$）。

原子范数最小化 ANM

核范数最小化 NNM

最小化迹的方法

相比之下，ANM算法的效果更好一些。

自己写的代码

如有不足，还望各位在评论区不吝赐教。

clc;clear all;close all;
%%%%%增广展开互质阵列（Augmented-Expanding Co-Prime Array）：填充内插阵元算法%%%%%
twpi = 2*pi;
derad = pi/180;
j = sqrt(-1);
c = 1500;%声速
f0 = 1e4;%载频
lambda = c/f0;%波长
dd = lambda/2;%阵元间距
M = 3;%子阵1的阵元数
N = 5;%子阵2的阵元数
kelm = 2*M+N-1;%实际阵元总数
d = [0:-N:(1-2*M)*N,M:M:(N-1)*M];%展开扩展互质阵
D = zeros(1,kelm*kelm);
for m=1:kelm
    D(1,(m-1)*kelm+1:m*kelm) = d-d(m);
end
[loc,index] = sort(D);
%%%%%构造排序去冗余矩阵FD%%%%%
indD = unique(D);%去除重复项
dof = length(indD);%阵列自由度
FD = zeros(dof,length(D));
for ii=1:length(D)
    count = length(find(D(ii)==D));
    temp = find(D(ii)==indD);
    FD(temp,ii) = 1/count;
end
%%%%%接收端%%%%%
d = d*dd;
% theta = [-0.1 0.1];%高精度到达角方向
theta = -45:6:45;%多到达角
iwave = length(theta);%信源数
snr = 30;%信噪比
n = 500;%采样数
fs = 2*f0;%采样频率
t = [1:n]/fs;
Amp = randn(iwave,n);
S = Amp.*(ones(iwave,1)*exp(j*twpi*(f0*t)));%独立信源
Ar = exp(j*twpi*d.'*sin(theta*derad)/lambda);
X = awgn(Ar*S,snr,'measured');
Rx = X*X'/n;
z = reshape(Rx,kelm*kelm,1);%转换为向量形式，成单快拍形式
%%%%%构造虚拟单快拍接收矢量：排序去冗余%%%%%
dm = unique(D);
zD = FD*z;
LD = (dof+1)/2;

%%%%%构造RD矩阵%%%%%
RD = zeros(LD,LD);%重构的协方差矩阵
for mr = 1:LD
    temp = LD+1-mr;
    RD(:,mr) = zD(temp:temp+LD-1);
end
vec_RD = reshape(RD,LD*LD,1);
%%%%%增加阵列自由度%%%%%
dm_ex = min(dm):1:max(dm);
zV = zeros(max(dm)-min(dm)+1,1);
dmax = max(dm);
for ii=1:dof
    zV(dm(ii)+dmax+1,1) = zD(ii,1);
end
%%%%%构造RV矩阵%%%%%%
LV = (length(zV)+1)/2;%Teoplitz矩阵大小
RV = zeros(LV,LV);%重构的协方差矩阵
for mr = 1:LV
    temp = LV+1-mr;
    RV(:,mr) = zV(temp:temp+LV-1);
end

%%%%%构造选择矩阵G%%%%%
G = ones(LV,LV);
posi0V = find(RV==0);
G(posi0V) = 0;

%%%%%利用CVX工具箱求解凸优化问题%%%%%
mu = 0.25;
cvx_begin sdp quiet
% cvx_precision high
cvx_solver sdpt3
    variable T(LV,LV) hermitian toeplitz semidefinite;
    minimize( square_pos(norm(T .* G - RV,'fro')) + mu * trace(T) )	%目标函数
cvx_end

%%%%%%%%%%构造MUSIC的谱函数%%%%%%%%%%%
[EV,Dv] = eig(T);%特征值分解
DD = diag(Dv);%将特征值变为向量形式
[DD,I] = sort(DD);%从小到大
DD = fliplr(DD');%翻转函数，从大到小
EV = fliplr(EV(:,I));
En = EV(:,iwave+1:end);%噪声子空间
dm_ss = dm_ex(1:LV)*dd;
for ii = 1:2001
    angle(ii) = (ii-1001)*90/1000;
    phim = derad*angle(ii);
    a = exp(j*twpi*dm_ss*sin(phim)/lambda).';
    Pmusic(ii) = 1/(a'*En*En'*a);
end
Pmusic = abs(Pmusic);
Pmax = max(Pmusic);
Pmusic_db = 10*log10(Pmusic/Pmax);
% [maxv_music,maxl_music]=findpeaks(Pmusic_db,'minpeakdistance',0.001,'minpeakheight',mean(Pmusic_db));
% [vl,ind] = sort(maxv_music,'descend');
% doa = angle(maxl_music(ind(1:iwave)));
h = plot(angle,Pmusic_db);
set(h,'Linewidth',0.2)
xlabel('angle(degree)')
ylabel('magnitude(dB)')
xlim([-90 90])
set(gca,'XTick',[-90:10:90])
grid on
hold on
m = ylim;
for i = 1:iwave
    line([theta(i) theta(i)],[m(1) 0],'Color','k','LineWidth',0.1,'LineStyle','--');
    hold on
end
legend('ANM','Direction of Arrival');

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118

相关文献与博客

下载文献ANM
提取码：6666

一个大佬的博客，数学证明很详细