已知稀疏信号$\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$范数问题恢复信号:

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

(1.1)

这里$\ell_0$范数表示向量非零元素的个数.

Candés等人证明了可以通过求解线性优化最小$\ell_0$范数问题精确重构信号, 因为最小$\ell_0$范数问题的解是最稀疏的解, 能够得到原信号的最优稀疏逼近.但是求解最小$\ell_0$范数问题是一个NP难问题, 即最小$\ell_0$范数问题的求解需要列举原信号中非零值的所有排列组合情况, 计算量巨大, 相关研究见文献[1-3].

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

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

(1.2)

最小$\ell_1$范数问题的求解是一个凸优化问题, 这个问题能够转化为线性规划问题.Candés, Romberg, Donoho已经证明了当目标信号$\bar{u}$是$\mathcal{K}$-稀疏的, 即$n$个分量只有$k$个分量是非零的, 且测量矩阵满足一些规则性条件, 理论上稀疏信号能够通过解问题(1.2) 得到有效恢复, 见文献[5-7].然而在实际中求解线性规划问题(1.2) 是很难的.在压缩感知问题中, 问题的维数可高达$n=10^6$, 当矩阵$A$是大且稠密的, 线性规划问题(1.2) 常常是病态的.

设矩阵$A$是行满秩的, 问题(1.2) 可转化为$\ell_1$正则化最小二乘问题(见文献[8]):

$ \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).

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

$ f_{\tau}(u)=\sum\limits_{i=1}^{n}{H(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$为光滑参数, 相关研究见文献[6].

光滑函数$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连续的.为方便, 记

$ s_k=u_{k+1}-u_k, 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+d_k, $

其中搜索方向$d_k$通过求解如下线性方程而得:

$ \nabla^2F(u_k)d_k=-g_k. $

(2.2)

牛顿法的优点是具有二阶收敛速度, 但当Hesse矩阵$\nabla^2F(u_k)$不正定时, 不能保证所产生的方向是目标函数$F(u)$在$u$处的下降方向.特别地, 当Hesse矩阵$\nabla^2F(u_k)$奇异时, 算法就无法进行下去.

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

拟牛顿法的基本思想是用某个矩阵$B_k$近似取代(2.2) 式中的Hesse阵$\nabla^2F(u_k)$, 通常$B_k$应具有的特点:

1.对所有的$k$, $B_k$是对称正定的, 从而使得算法所产生的方向是目标函数$F(u)$在$u_k$处的下降方向;

2.矩阵$B_k$更新规则相对较简单.

本文搜索方向$d_k$通过求解如下线性方程而得

$ B_kd_k=-g_k. $

其中矩阵$B_k$更新规则如下:

$ B_{k+1}=B_k-\frac{B_ks_ks_k^TB_k}{s_k^TB_ks_k}+\frac{y_ky_k^T}{y_k^Ts_k}, $

(2.3)

该公式称为BFGS校正公式.

算法2.1  问题(2.1) 的BFGS校正拟牛顿算法:

步骤0:  初始点$u_0\in\mathbf{R}^n$, 终止误差$0\leq\epsilon\ll1$, 初始对称正定矩阵$B_0=I_n$($n$阶单位阵), $\tau_0=0.8$, 且$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$.

$ B_kd_k=-g_k. $

步骤3:  由Wolfe线搜索计算步长$\alpha_k$.

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

(2.4)

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

(2.5)

步骤4:更新:  令$u_{k+1}=u_k+\alpha_kd_k$, 由BFGS校正公式(2.3) 更新$B_k$为$B_{k+1}$, 更新光滑参数$\tau_{k+1}=\frac{1}{2}\tau_k$.

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

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

引理3.1  设$B_k$对称正定, $B_{k+1}$由BFGS校正公式(2.3) 确定, 则$B_{k+1}$对称正定的充要条件是$y_k^Ts_k>0$.

证  详细证明参考文献[10].

引理3.1表明, 若初始矩阵$B_0$对称正定且在迭代过程中保持$y_k^Ts_k>0~(∀k\geq0)$, 则由BFGS校正公式产生的矩阵序列$\{B_k\}$是对称正定的, 从而方程$B_kd_k=-g_k$有唯一解$d_k$, 故有

$ d_k.TB_kd_k=-g_k^Td_k, $

因为$d_k^TB_kd_k>0$, 所以$g_k^Td_k < 0$, 即$d_k$是目标函数$F(u)$在$u_k$处的下降方向.

引理3.2  若在BFGS校正拟牛顿算法中采用Wolfe搜索准则, 则有$y_k^Ts_k>0$.

证  因为

$ y_k^Ts_k=(g_{k+1}-g_k)^T\alpha_kd_k=\alpha_k(g_{k+1}^Td_k-g_k^Td_k), $

由算法步骤$3$中(2.5) 式有

$ g_{k+1}^Td_k\geq\rho_2g_k^Td_k, $

故

$ y_k^Ts_k\geq\alpha_k(\rho_2-1)g_k^Td_k, $

又因为$\rho_2 < 1$, $g_k^Td_k < 0$, 故$y_k^Ts_k>0$.

引理3.3  若$\ell_1$范数$\|u\|_1$按光滑化技术得到的光滑函数为

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

而

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

则有$f_{\tau}(u)\leq\|u\|_1$.

证  当$|u(i)|\leq\tau$时, $f_{\tau}(u)=\sum\limits_{i=1}^n\frac{u^2(i)}{2\tau}=\frac{1}{2\tau}\sum\limits_{i=1}^nu^2(i)\leq\frac{1}{2\tau}\cdot\tau\sum\limits_{i=1}^n|u(i)|=\frac{1}{2}\|u\|_1\leq\|u\|_1;$

当$|u(i)|>\tau$时, $f_{\tau}(u)=\sum\limits_{i=1}^n(|u(i)|-\frac{\tau}{2})=\|u\|_1-\frac{n\tau}{2}\leq\|u\|_1.$

综上有$f_{\tau}(u)\leq\|u\|_1$.

定理3.4  设$\{u_k\}$表示由BFGS校正拟牛顿法迭代产生的序列, 当$\|u\|_1$的光滑参数$\tau_k\rightarrow0$时, 则$\{u_k\}$是有界序列.令$u^*$表示$\{u_k\}$的任意极限点, 则$u^*$是$\ell_1$范数最小化问题(1.2) 的最优解.

证  对于第$k$步的光滑参数$\tau_k$, 令$u_k^*=\arg\min\limits_{u\in\mathbf{R}^n}F_k(u)$表示第$k$步的目标函数$F_k(u)$的最小值点, $v$表示满足$Av=b$的任意向量.

因为$u_k^*=\arg\min\limits_{u\in\mathbf{R}^n}F_k(u)$, 所以$F_k(u_k^*)\leq F_k(v)$.又$Av=b$, 故

$ F_k(v)=\lambda f_{\tau_k}(u)+\frac{1}{2}\|Av-b\|_2^2=\lambda f_{\tau_k}(v). $

因为$F_k(u)=\lambda f_{\tau_k}(u)+\frac{1}{2}\|Au-b\|_2^2$, $\|Au-b\|_2^2\geq0$, 故

$ F_k(u)\geq\lambda f_{\tau_k}(u). $

(3.1)

而由引理3.3有$f_{\tau_k}(v)\leq\|v\|_1$, 则

$ F_k(u_k^*)\leq F_k(v)=\lambda f_{\tau_k}(v)\leq\lambda\|v\|_1. $

(3.2)

结合式(3.1)、(3.2) 及$f_{\tau_k}(u)$的光滑性有

$ \|u_k\|_1\leq_{\tau_k}(u_k)+\frac{n\tau_k}{2}\leq\frac{1}{\lambda}F_k(u_k)+\frac{n\tau_k}{2}\leq\frac{1}{\lambda}\lambda\|v\|_1+\frac{n\tau_k}{2}=\|v\|_1+\frac{n\tau_k}{2}, $

故当光滑参数$\tau_k\rightarrow0$时, $\{u_k\}$是有界序列.

另一方面, 参考文献[8], 可证序列$\{u_k\}$的极限点$u^*$是可行点, 且极限点$u^*$是$\ell_1$范数最小化问题(1.2) 的KKT点.因为$\|u\|_1$是凸的, 对于问题(1.2) 是凸规划问题, 故$u^*$是$\ell_1$范数最小化问题(1.2) 的最优解.

在压缩感知理论中, 当$\ell_1$范数最小化问题$\min\{\|u\|_1:Au=b\}$有唯一解时, 精确恢复才会实现.参考文献[8], 有下面的推论.

推论3.5  设$\ell_1$范数最小化问题$\min\{\|u\|_1:Au=b\}$有唯一最优解, 令$\{u_k\}$表示由BFGS校正拟牛顿法迭代产生的序列, 当$\|u\|_1$的光滑参数$\tau_k\rightarrow0$时, 则$\lim\limits_{k\rightarrow\infty}u_k=u^*, $这里$u^*=\arg\min\{\|u\|_1:Au=b\}$.

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

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

本文用BFGS校正拟牛顿法与对称秩1校正拟牛顿法解决大规模信号恢复问题作比较, 其中对称秩1校正公式如下:

$ B_{k+1}=B_k+\frac{(y_k-B_ks_k)(y_k-B_ks_k)^T}{(y_k-B_ks_k)^Ts_k}. $

本文提出用BFGS校正拟牛顿法解决大规模信号恢复问题, 最后做了两类数值实验, 第一类反映了随着正则化参数$\lambda$的变化, 相对残差、相对误差、运行时间的变化, 以及随着原信号的维数的变化, 相对残差、相对误差、运行时间的变化; 第二类实验用本文提出的BFGS校正拟牛顿法与对称秩1校正拟牛顿法解决大规模信号恢复问题作对比, 实验结果表明用BFGS校正拟牛顿法解决大规模信号恢复问题是可行的.该算法在图像处理中有直接的应用价值, 相关研究见文献[11, 12].