压缩感知-Gradient Projection for Sparse Reconstruction(梯度投影稀疏重建)

压缩感知-Gradient Projection for Sparse Reconstruction

GPSR的作用主要是为了解决无约束的凸优化问题：

minx12∥y−Ax∥22+τ∥x∥1

而GPSR的第一个关键步骤就是将公式(1)转化为一个二次规划问题。转化的方法是将x分离成正负两部分：

x=u−v,u≥0,v≥0

则约束项变为：

∥x∥1=1Tnu+1Tnv

这样式(1)就可以重写成带边界约束的二次规划问题(bound-constrained quadratic program, BCQP)：

minx12∥y−A(u−v)∥22+τ1Tnu+τ1Tnv,s.t.u≥0,v≥0

公式(3)可以进一步写成一个标准的BCQP形式：

minzcTz+12zTBz≡F(z),s.t.z≥0

其中F(z)表示目标函数，

z=[uv],b=ATy,c=τ1Tnu+[−bb]

B=[ATA−ATA−ATAATA]

后面有两个经常用到的公式为：

∇F(z)=c+Bz

而：

Bz=B[uv]=[ATA(u−v)−ATA(u−v)]

以及：

zTBz=(u−v)TATA(u−v)=∥A(u−v)∥22

为了求解公式(4),梯度投影GP的核心思想是，首先找到一个标量参数序列α(k)>0且计算：

w(k)=(z(k)−α(k)∇F(z(k)))+

然后选择第二个标量参数序列λ(k)∈[0,1]并更新：

z(k+1)=z(k)+λ(k)(w(k)−z(k))

而根据αk+1和λk+1的产生测量不同产生了两种方法：Basic-Gradient Projection for Sparse Reconstruction (GPSR-Basic)和Barzilai-Borwein Projection for Sparse Reconstruction (GPSR-BB)。

1. GPSR-Basic

GPSR-Basic的思想是让zk+1每次都沿着负梯度方向−∇F(zk+1)进行搜索，并在非负象限上进行投影，执行一个回朔线性搜索直至迭代停止(非常有名的“Armijo rule”)。定义：

g(k)i={(∇F(z(k)))i,ifz(k)>0or(∇F(z(k)))i<00,otherwise

初始化α0的方法是求：

α0=argminαF(z(k)−αg(k))

且可以直接解得：

α0=(g(k))Tg(k)(g(k))TBg(k)

为了防止α0太小或者太大，提前定义[αmin,αmax]，且0≤αmin≤αmax。每次求得的α都要取三者的中间值mid[αmin,α,αmax]。

具体算法步骤如下：

1. 初始化。给定z(0)，设置β∈(0,1)，μ∈(0,1/2)，迭代次数k=0；
2. 使用公式(11)计算α0，并替换α0=mid[αmin,α0,αmax];
3. 回朔线性搜索：从序列α0,βα0,β2α0,...中选取满足：F((z(k)−α(k)∇F(z(k)))+)<F(z(k))−μ∇F(z(k))T(z(k)−(z(k)−α(k)∇F(z(k)))+)的第一个值为αk，并更新
   z(k+1)=(z(k)−α(k)∇F(z(k)))+；
4. 判断是否收敛，若未收敛则返回步骤2继续迭代。

2. GPSR-BB

GPSR-Basic算法保证了目标函数F每一次迭代都会下降，而BB策略则没有这一特性。BB算法最早是用在求解无约束的平滑非线性最小化问题，其每一步都计算一个更新步长为δ(k)=−H−1k∇F(z(k))，其中Hk表示F在z(k)出Hessian矩阵的近似。BB策略给出了一种非常简单的Hk近似方法Hk=η(k)I。其中，η(k)要和真实的Hessian矩阵的特性非常近似才可以，即：

F(z(k))−∇F(z(k−1))≈η(k)[z(k)−z(k−1)]

一般取η(k)为满足上述关系的最小二乘值。在无约束情况下，此时更新公式就变成了：

z(k+1)=z(k)−(η(k))−1∇F(z(k))

且即使目标函数值F增大，也仍然按照上式更新。
GPSR-BB算法是将BB策略扩展到BCQP问题，选择公式(6)中的λk∈[0,1]作为精确的最小化因子，η(k+1)的选择按照上面描述的方法进行。具体算法步骤如下：

1. 初始化。给定给定z(0)，设置αminαmax，,α0∈[αminαmax]，迭代次数k=0；
2. 计算更新步长：δ(k)=(z(k)−α(k)∇F(z(k)))+−z(k)；
3. 线性搜索。在[0,1]区间内找到最小化F(z(k)+λ(k)δ(k))的值为λ(k)，且更新z(k+1)=z(k)+λ(k)δ(k)；
4. 更新α。计算
   
   γ(k)=(δ(k))TBδ(k)
   
   如果γ(k)=0，则令α(k+1)=αmax，否则：
   α(k+1)=mid⎧⎩⎨⎪⎪αmin,∥∥δ(k)∥∥22γ(k),αmax⎫⎭⎬⎪⎪
5. 判断是否收敛，若未收敛则返回步骤2继续迭代。

上述步骤(3)中λ(k)的计算方式为：

λk+1=mid⎧⎩⎨0,(δ(k))T∇F(z(k))(δ(k))TBδ(k),1⎫⎭⎬

参考文献：
1. Figueiredo M A T, Nowak R D, Wright S J. Gradient Projection for Sparse Reconstruction: Application to Compressed Sensing and Other Inverse Problems[J]. IEEE Journal of Selected Topics in Signal Processing, 2008, 1(4):586-597.