原子范数 Atomic norm最小化：简单的Matlab例程

最新推荐文章于 2024-08-02 04:01:52 发布

B417科研笔记

最新推荐文章于 2024-08-02 04:01:52 发布

阅读量7.1k

点赞数 20

分类专栏：通信中的常用数学文章标签： matlab 线性代数矩阵

本文链接：https://blog.csdn.net/weixin_39274659/article/details/121982277

版权

通信中的常用数学专栏收录该内容

48 篇文章 517 订阅

订阅专栏

前言

基于压缩感知的尽头: 原子范数最小化中的原子范数最小化算法，笔者做了一些matlab的仿真，作为简单的例程，希望帮助大家进一步理解算法和自定义的拓展。

由于凸问题的求解需要使用 CVX，因此需要读者先自行安装好 matlab 的 CVX包。

1-D 无噪场景

假设接收天线有 $64$ 根，有 $3$ 个单天线的目标源同时发射信号，接收信号可以表示为：
$z=[a(f_1), a(f_2), a(f_3)][s_1\; s_2\; s_3]^T=A(f)s$
其中， $\boldsymbol{a}(f)=\left[1, e^{i 2 \pi f}, \ldots, e^{i 2 \pi(M-1) f}\right]^{T}$ 代表天线响应矢量。这里假设没有噪声。那么从 $z$ 中估计 $f$ ，我们使用原子范数最小化算法。代码如下：

N = 64;  % number of antennas
f = [0.1, 0.4, 0.3];  % f = cos(theta)
A = exp(1j * 2 * pi * (0 : N-1)' * f);  % A = [a(f1), ..., a(fK)]
s = ones(3, 1);  % 3 sources
z = A * s;


cvx_begin sdp 
cvx_solver sdpt3
variable T(N, N) hermitian toeplitz
variable x 
minimize(0.5 * x + 0.5 * T(1,1))
[ x z'; z T] >= 0;
cvx_end

[Phi,P] = rootmusic(T, 3, 'corr');
Phi / 2 / pi

中间 cvx_begin 和 cvx_end 包围的部分，我们就是在求解和原子范数最小化等价的 SDP问题。 CVX可以返回解得的 T。最后，我们使用 rootmusic函数，可以将 $T$ 对应的 $f$ 得到，即为 Phi。（本来是做范德蒙德分解，然后得到，但通过rootmusic提取出 $f$ 更加便捷。）

输出结果如下：

>> ANM_1D
 
Calling SDPT3 4.0: 4225 variables, 128 equality constraints
   For improved efficiency, SDPT3 is solving the dual problem.
------------------------------------------------------------

 num. of constraints = 128
 dim. of sdp    var  = 130,   num. of sdp  blk  =  1
*******************************************************************
   SDPT3: Infeasible path-following algorithms
*******************************************************************
 version  predcorr  gam  expon  scale_data
   HKM      1      0.000   1        0    
it pstep dstep pinfeas dinfeas  gap      prim-obj      dual-obj    cputime
-------------------------------------------------------------------
 0|0.000|0.000|4.3e+03|1.1e+01|2.2e+05| 0.000000e+00  0.000000e+00| 0:0:00| chol  1  1 
 1|0.893|1.000|4.6e+02|2.3e-01|2.9e+04|-5.933532e+00 -3.889838e+01| 0:0:00| chol  1  1 
 2|0.988|1.000|5.4e+00|2.3e-02|3.8e+02|-2.814220e+00 -3.796329e+01| 0:0:00| chol  1  1 
 3|0.943|1.000|3.1e-01|2.3e-03|3.3e+01|-3.966408e+00 -1.985735e+01| 0:0:00| chol  1  1 
 4|0.888|0.790|3.5e-02|6.6e-04|1.3e+01|-2.602549e+00 -1.421772e+01| 0:0:00| chol  1  1 
 5|1.000|1.000|1.2e-09|7.0e-03|2.8e+00|-2.914921e+00 -5.690264e+00| 0:0:00| chol  1  1 
 6|0.980|0.987|1.4e-10|9.2e-05|3.6e-02|-2.998482e+00 -3.034630e+00| 0:0:00| chol  1  1 
 7|0.979|0.987|4.9e-11|1.4e-06|5.0e-04|-2.999973e+00 -3.000466e+00| 0:0:00| chol  1  1 
 8|0.978|0.986|9.6e-11|2.1e-08|7.3e-06|-2.999999e+00 -3.000007e+00| 0:0:00| chol  1  1 
 9|0.976|0.990|2.1e-10|2.2e-10|2.1e-07|-3.000000e+00 -3.000000e+00| 0:0:00| chol  1  1 
10|0.950|0.988|1.3e-11|2.5e-11|9.0e-09|-3.000000e+00 -3.000000e+00| 0:0:01|
  stop: max(relative gap, infeasibilities) < 1.49e-08
-------------------------------------------------------------------
 number of iterations   = 10
 primal objective value = -3.00000000e+00
 dual   objective value = -3.00000001e+00
 gap := trace(XZ)       = 8.95e-09
 relative gap           = 1.28e-09
 actual relative gap    = 1.31e-09
 rel. primal infeas (scaled problem)   = 1.26e-11
 rel. dual     "        "       "      = 2.48e-11
 rel. primal infeas (unscaled problem) = 0.00e+00
 rel. dual     "        "       "      = 0.00e+00
 norm(X), norm(y), norm(Z) = 7.4e-01, 1.4e+01, 8.0e+01
 norm(A), norm(b), norm(C) = 6.5e+01, 1.7e+00, 1.5e+01
 Total CPU time (secs)  = 0.50  
 CPU time per iteration = 0.05  
 termination code       =  0
 DIMACS: 1.4e-11  0.0e+00  1.5e-10  0.0e+00  1.3e-09  1.3e-09
-------------------------------------------------------------------
 
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +3
 

ans =

    0.4000
    0.1000
    0.3000

中间是CVX的求解过程。可以通过在 cvx_begin 后面输入 quiet参数使其不输出。最后的结果就是根据原子范数最小化求得的 $f$ 的值，可以看到。分毫不差。

把 z= A * s替换为 z = A * s + 0.5 * (randn(N, 1) + 1j * randn(N, 1));，可以人为的加入噪声，当噪声大小增大时，可以看到估计结果会出现一些误差。

s=ones(3,1)可以替换为 s=randn(3,1)之类的，加入随机性。

2-D 无噪声场景

能看到这里，说明已经对原子范数有了一定的了解。那这里就不再多废话，直接上代码：

N = 64;  % number of antennas
L = 5;
f = [0.1, 0.4, 0.3];  % f = cos(theta)
A = exp(1j * 2 * pi * (0 : N-1)' * f);  % A = [a(f1), ..., a(fK)]
S = randn(3, L) + 1j * randn(3, L);  % 3 sources
Z = A * s ;


cvx_begin sdp quiet
cvx_solver sdpt3
variable T(N, N) hermitian toeplitz
variable X(L, L) hermitian 

minimize(trace(X) + trace(T))
[X Z'; Z T] >= 0;
cvx_end

[Phi,P] = rootmusic(T, 3, 'corr');
Phi / 2 / pi

2-D 有噪声场景

N = 64;  % number of antennas
L = 5;
sigma = 0.5; 
f = [0.1, 0.4, 0.3];  % f = cos(theta)
A = exp(1j * 2 * pi * (0 : N-1)' * f);  % A = [a(f1), ..., a(fK)]
S = randn(3, L) + 1j * randn(3, L);  % 3 sources
Y = A * s + sigma * (randn(N, 1) + 1j * randn(N, 1));
lambda = sqrt(N * (L + log(N) + sqrt( 2 * L * log(N)))*sigma);

cvx_begin sdp quiet
cvx_solver sdpt3
variable T(N, N) hermitian toeplitz
variable X(L, L) hermitian 
variable Z(N, L) complex

minimize(lambda * (trace(X) + trace(T)) + 0.5*sum_square_abs(vec(Y - Z)))
[X Z'; Z T] >= 0;
cvx_end

[Phi,P] = rootmusic(T, 3, 'corr');
Phi / 2 / pi