已知$ \bar{u}\in\mathbf{R}^n $是稀疏信号, $ A\in\mathbf{R}^{m\times n}(m\ll n) $是线性算子, $ b\in\mathbf{R}^m $是观测向量, 它们满足下面的关系式$ b=A\bar{u}. $原信号$ \bar{u} $从线性系统中重构是压缩感知问题中稀疏信号重构的核心问题, 然而该线性方程是欠定的, 它有无穷多组解.

在压缩感知问题中, 若先验知道原信号在某个变换下是可压缩的(或是稀疏的), 则可以通过如下一个最小$ \ell_0 $范数问题恢复信号$ \bar{u} $

$ \min\limits_{u\in\mathbf{R}^n}\|u\|_0~~~{\rm s.t.}~Au=b, $

(1.1)

这里$ \ell_0 $范数表示向量非零元素的个数.但是求解最小$ \ell_0 $范数问题是一个NP难问题, 即最小$ \ell_0 $范数问题的求解需要列举原信号中非零值的所有排列组合情况, 计算量巨大, 相关研究见文献[1].

此后, 人们转而求解最小$ \ell_1 $范数问题, 而且证明了在一定的条件下, 最小$\ell_1$范数问题和最小$ \ell_0 $范数问题的解是等价的, 见文献[2].因此问题(1.1) 转化为求解下面的最小$ \ell_1 $范数问题

$ \min\limits_{u\in\mathbf{R}^n}\|u\|_1~~~{\rm s.t.}~Au=b. $

(1.2)

最小$ \ell_1 $范数问题的求解是一个凸优化问题, 这个问题能够转化为线性规划问题, 见文献[3, 4].然而在实际中求解线性规划问题(1.2)是很难的.在压缩感知问题中, 当矩阵$ A $是大且稠密的, 线性规划问题(1.2)常常是病态的.设矩阵A是行满秩的, 问题(1.2)可转化为$\ell_1$正则化最小二乘问题(见文献[5])

$ \min\limits_{u\in\mathbf{R}^n}\lambda\|u\|_1+\frac{1}{2}\|Au-b\|_2^2, $

(1.3)

这里$ \lambda $是正则化参数.

因为$ \|u\|_1 $是关于$u$的非光滑凸函数, 所以对于问题(1.3), 次梯度基本方法的效果可能很差.本文通过求解问题(1.3)的一个适当的光滑形式, 来求解$ \ell_1 $最小化问题(1.2).

由文献[6], 利用Nesterov光滑技术, 光滑化$ \|u\|_1 $为如下光滑函数

$ f_{\tau}(u)=\sum\limits_{i=1}^nH(u(i)), $

其中

$ H(y)=\left\{ \begin{array}{ll} \frac{y^2}{2\tau},&\hbox{$|y|\leq\tau$;} \\ |y|-\frac{\tau}{2},&\hbox{$|y|>\tau$, } \end{array} \right. $

表示Hüber惩罚函数, $ \tau $为光滑参数, 相关研究见文献[3].

光滑函数$ f_{\tau}(u) $是凸的, 梯度$ \nabla f_{\tau}(u) $是Lipchitz连续的, 且其分量为

$ \nabla f_{\tau}(u(i))=\left\{ \begin{array}{ll} \frac{u(i)}{\tau},&\hbox{$|u(i)|\leq\tau$, } \\ {\rm sign}(u(i)),&\hbox{$|u(i)|>\tau$} \end{array} \right.~~(i=1, 2, \cdots, n), $

其中

$ {\rm sign}(x)=\left\{ \begin{array}{ll} 1,&\hbox{$x>0$;} \\ 0,&\hbox{$x=0$;} \\ -1,&\hbox{$x<0$.} \end{array} \right. $

令$ F(u)=\lambda f_{\tau}(u)+\frac{1}{2}\|Au-b\|_2^2 $, 则问题(1.3)转化为如下无约束光滑凸规划问题

$ \min\limits_{u\in\mathbf{R}^n}F(u), $

(2.1)

$ F(u) $是可微的, 其梯度$ \nabla F(u)=\lambda\nabla f_{\tau}(u)+A^T(Au-b) $是Lipchitz连续的.为方便, 记$ g_k=\nabla F(u_k)$, $ y_k=\nabla F(u_{k+1})-\nabla F(u_k)=g_{k+1}-g_k $.

求解问题(2.1), 共轭梯度法的迭代公式如下

$ u_{k+1}=u_k+\alpha_kd_k, $

(2.2)

$ d_k=\left\{ \begin{array}{ll} -g_k,&\hbox{$k=0$;} \\ -g_k+\beta_kd_{k-1},&\hbox{$k\geq1$, } \end{array} \right. $

(2.3)

其中$ \alpha_k $为搜索步长, $ d_k $为搜索方向, $ \beta_k $为参数.

共轭梯度法的关键是参数$ \beta_k $的选取, 常用的$ \beta_k $公式有如下几种

$ \beta_k^{HS}=\frac{g_k^Ty_{k-1}}{d_{k-1}^Ty_{k-1}}, \beta_k^{FR}=\frac{\|g_k\|^2}{\|g_{k-1}\|^2}, \beta_k^{DY}=\frac{\|g_k\|^2}{d_{k-1}^Ty_{k-1}}, \beta_k^{PRP}=\frac{g_k^Ty_{k-1}}{\|g_{k-1}\|^2}. $

一些文献对$ \beta_k $的构造及算法的收敛性做了大量的研究, 见文献[7, 8].文献[8]指出, PRP, HS方法在数值计算中有很好的表现, 但是即使使用精确线搜索PRP方法也不会有全局收敛性. FR, DY方法则具有较好的收敛性质.为寻求兼具良好收敛性和优秀数值结果的算法, 许多研究者在经典共轭梯度法的基础上推出新的$ \beta_k $公式.

本文我们提出用修正HS共轭梯度法解决大规模信号恢复问题(2.1), 且在适当的假设条件下证明了本文提出的算法具有全局收敛性, 最后数值实验结果表明用修正HS共轭梯度法解决大规模信号恢复问题是可行的.

本文构造如下搜索方向

$ d_k=\left\{ \begin{array}{ll} -g_k,&\hbox{$k=0$;} \\ -g_k+\beta_kd_{k-1}-\theta_kg_k,&\hbox{$k\geq1$, } \end{array} \right. $

(2.4)

其中

$ \beta_k=\frac{g_k^Ty_{k-1}}{d_{k-1}^Ty_{k-1}}, \theta_k=\frac{g_k^Ty_{k-1}}{d_{k-1}^Ty_{k-1}}\frac{g_k^Td_{k-1}}{\|g_{k}\|^2}, y_{k-1}=g_{k}-g_{k-1}=\nabla F(u_{k})-\nabla F(u_{k-1}). $

引理 2.1 对任一线搜索, 搜索方向由(2.4)式定义, 则有

$ g_k^Td_{k}=-\|g_{k}\|^2. $

(2.5)

证  由(2.4) 式, 有

$ g_k^Td_{k}=-g_k^Tg_{k} +\frac{g_k^Ty_{k-1}}{d_{k-1}^Ty_{k-1}}g_k^Td_{k-1} -\frac{g_k^Ty_{k-1}}{d_{k-1}^Ty_{k-1}}\frac{g_k^Td_{k-1}}{\|g_{k}\|^2}g_k^Tg_{k} =-g_k^Tg_{k}. $

即$ g_k^Td_{k}=-\|g_{k}\|^2. $这表明搜索方向$ d_k $是目标函数$ F(u) $在$ u_k $处的下降方向, 若使用精确线搜索, 有$ g_k^Td_{k-1}=0 $, 即$ \theta_k=0 $, 则本文提出的修正HS共轭梯度法退化成标准的HS共轭梯度法.

算法 2.1  问题(2.1)的修正HS共轭梯度算法:

步骤 0  给定初始点$ u_0\in\mathbf{R}^n $, 终止误差$ 0\leq\epsilon\ll1 $, $ \tau_0=0.6 $, $ 0<\rho_1<\rho_2<1 $, 令$ k:=0 $;

步骤 1  计算$ g_k=\nabla F(u_k) $, 若$ \|g_k\|<\epsilon $, 且光滑参数$ \tau_k\leq\epsilon $, 算法停止; 否则, 转下一步;

步骤 2  计算搜索方向$ d_k $, 搜索方向由(2.4)式定义;

步骤 3  计算步长$ \alpha_k $:

$ F(u_k+\alpha_kd_k)\leq F(u_k)+\rho_1\alpha_k\nabla F(u_k)^Td_k; $

(2.6)

$ \nabla F(u_k+\alpha_kd_k)^Td_k\geq\rho_2\nabla F(u_k)^Td_k; $

(2.7)

步骤4  更新:令$ u_{k+1}=u_k+\alpha_kd_k $, 更新光滑参数$ \tau_{k+1}=\frac{1}{4}\tau_k $;

步骤5  循环:置$ k:=k+1 $, 转步骤1.

本文提出用修正HS共轭梯度法解决大规模信号恢复问题(2.1), 且在适当的假设条件下证明了本文提出的算法具有全局收敛性.

为证明算法2.1具有全局收敛性, 先作如下基本假设.

H1  水平集$ \Omega=\{x\in R^n|F(x)\leq F(x_0)\} $有界.

H2  在$ \Omega $的某邻域N内, F(x)连续可微, 且其梯度$ \nabla F(x)=g(x) $是Lipschitz连续的, 即存在一常数L > 0使得

$ \|g(x)-g(y)\|\leq L\|x-y\|, ~\forall~x, y\in N. $

(3.1)

由假设H1、H2知, 存在一常数M > 0满足

$ \|g(x)\|\leq M, ~\forall~x\in \Omega. $

(3.2)

参考文献[8, 9], 有如下两个引理.

引理 3.1^[9]  若假设条件H1、H2成立, 则算法是良定的.

注 由(2.6)式知目标函数值序列$ \{F(x_k)\} $是一下降序列, 有$ \sum\limits_{k\geq0}\alpha_k|g_k^Td_k|<+\infty $, 结合(3.2)式有

$ \lim\limits_{k\rightarrow\infty}\alpha_k\|d_k\|=0. $

(3.3)

引理 3.2^[8] 若假设条件H1、H2成立, 对任一共轭梯度法的迭代形式(2.2), 其中d_k满足$ g_k^Td_k<0 $且$ \alpha_k $满足条件(2.6)和(2.7), 则

$ \sum\limits_{k\geq0}\frac{(g_k^Td_k)^2}{\|d_k\|^2}<+\infty. $

(3.4)

式(3.4)结合式(2.5)有

$ \sum\limits_{k\geq0}\frac{\|g_k\|^4}{\|d_k\|^2}<+\infty. $

(3.5)

定理 3.3  如果假设条件H1、H2成立, 算法2.1产生的迭代序列为{u_k}, 若存在常数r > 0满足

$ d_k^Ty_k\geq r, $

(3.6)

则$ \lim\limits_{k\rightarrow\infty}\inf\|g_k\|=0. $

证  由(2.4) 式定义的d_k及式(3.1)、(3.2)、(3.6) 有

$ \|d_k\|\leq\|g_k\|+\frac{2\|g_k\|\|y_{k-1}\|d_{k-1}\|}{d_{k-1}^Ty_{k-1}} \leq M+\frac{2ML\alpha_{k-1}\|d_{k-1}\|}{r}\|d_{k-1}\|. $

由(3.3)式知存在一常数$ t\in(0, 1) $, 存在一正整数k, 使得当$ k\geq k_0 $有

$ \frac{2ML}{r}\alpha_{k-1}\|d_{k-1}\|\leq t. $

因此$ \forall~ k>k_0 $有

$ \|d_k\|\leq M+t\|d_{k-1}\|\leq M+t(M+t\|d_{k-2}\|)\\ \;\;\;\;\;\leq M(1+t+t^2+\cdots+t^{k-k_0-1})+t^{k-k_0}\|d_{k-k_0}\|\\ \;\;\;\;\;\leq\frac{M}{t-1}+\|d_{k_0}\|. $

令$ m=\max\Big\{\|d_1\|, \|d_2\|, \cdots, \|d_{k_0}\|, \frac{M}{t-1}+\|d_{k_0}\|\Big\} $, 则

$ \|d_k\|\leq m, \forall~ k. $

(3.7)

式(3.5)结合式(3.7)有

$ \sum\limits_{k\geq0}\|g_k\|^4<+\infty, $

故$ \lim\limits_{k\rightarrow\infty}\inf\|g_k\|=0. $

在压缩感知理论中, 当$ \{\min\|u\|_1:Au=b\} $有唯一解时, 精确恢复才会实现.参考文献[5], 有下面的定理.

定理 3.4  设$ \ell_1 $范数最小化问题$ \min\{\|u\|_1:Au=b\} $有唯一最优解, 算法2.1产生的迭代序列为$ \{u_k\} $, 当$ \|u\|_1 $的光滑参数$ \tau_k\rightarrow0 $时, 则$ \lim\limits_{k\rightarrow\infty}u_k=u^* $, 这里$ u^*=\arg\min\{\|u\|_1:Au=b\} $.

证  记第k步的目标函数$ F_k(u) $的最小值点为$ u_k^* $, 即$ u_k^*=\arg\min\limits_{u\in\mathbf{R}^n}F_k(u) $.

由Nesterov光滑技术得到的函数$ f_{\tau_k}(u_k) $的光滑性有$ f_{\tau_k}(u_k)\geq \|u_k\|_1-\frac{n\tau_k}{2}, $故

$ \lambda\|u_k\|_1\leq \lambda f_{\tau_k}(u_k)+\frac{\lambda n\tau_k}{2} \leq\lambda f_{\tau_k}(u_k)+\frac{1}{2}\|Au_k-b\|_2^2+\frac{\lambda n\tau_k}{2} =f_{k}(u_k)+\frac{\lambda n\tau_k}{2}. $

因为F_k(u)是凸函数, 有Lipschitz连续梯度, 由文献[5]引理2.3有, Lipschitz常数满足

$ F_k(u_k)\leq F_k(u_k^*) +\frac{\tau_k^2}{2L_k}, $

所以

$ \lambda\|u_k\|_1\leq F_k(u_k^*) +\frac{\tau_k^2}{2L_k} +\frac{\lambda n\tau_k}{2}. $

(3.8)

故当光滑参数$ \tau_k\rightarrow0 $时, (3.8)式表明序列{u_k}有极限点.令u^*表示这个序列的任意极限点, 由文献[5]结合定理3.3, 可证u^*是可行点, 且是问题(1.2)的KKT点.又问题(1.2)是凸规划问题, 故u^*是问题(1.2)的最优解.

实验中, 参数选取$ \rho_1=0.2 $, $ \rho_2=0.7 $, 误差精度$ \epsilon=10^{-8} $. $A$为$m\times n$维随机矩阵, $u_{\rm exact}$表示$n$维$k$稀疏原信号, 满足$Au_{\rm exact}=b$, $u$表示本文算法恢复出来的信号, 相对残差$\frac{\|Au-b\|_2}{\|b\|_2}$, 相对误差$\frac{\|u-u_{\rm exact}\|_2}{\|u_{\rm exact}\|_2}$. 图 1和图 2中的圈和点分别表示原信号和经由我们的算法恢复出来的信号.

表 1、图 1选取$A$的维数都是$512\times1024$, 表 2对应的正则化参数$\lambda=0.1$, 图 2选取$A$的维数分别是$256\times512$, $512\times1024$, $1024\times2048$, $2048\times4096$.

本文提出用基于Nesterov光滑技术和修正HS共轭梯度法的压缩感知重构算法解决大规模信号恢复问题, 最后做了两类数值实验, 第一类反映了随着正则化参数$\lambda$的变化, 相对残差、相对误差、运行时间的变化, 以及随着原信号的维数的变化, 相对残差、相对误差、运行时间的变化; 第二类实验用本文提出的修正HS共轭梯度法与HS共轭梯度法、PRP共轭梯度法解决大规模信号恢复问题作对比, 实验结果表明用修正HS共轭梯度法解决大规模信号恢复问题是可行的.该算法在图像处理中有直接的应用价值.相关研究见文献[10, 11].