一类信号恢复问题的投影迭代硬阈值算法❋

李庆龙， 曾雪迎

(中国海洋大学数学科学学院， 山东 青岛 266100)

稀疏恢复在信号和图像处理、机器学习等许多领域中都得到了广泛的关注。基于部分傅里叶系数恢复原始信号是稀疏恢复中的一类重要问题，在压缩感知、磁共振成像等领域中具有重要的应用价值[1]。该问题需要寻找一个可以在傅里叶变换域内拟合采样不足的傅里叶系数的信号。稀疏恢复是典型的不适定反问题。数据可稀疏逼近使得从部分傅里叶系数中恢复原始信号成为可能[2]。该问题可以建模为L0模最优化问题，但它是一个非凸不连续的NP难问题，模型的理论分析和数值求解均非常困难。在不增加测量数据的情形下，已有文献的解决方案通常是将L0模松弛为L1模[2-4]或其它非凸但连续的稀疏性度量函数[5]。

虽然L1模对应的凸优化问题具有快速算法[6-7],但其会在恢复信号中引入偏差，因此通常不能取得令人满意的恢复效果。对于低维问题，增加观测傅里叶系数的数量可以带来更准确的解决方案，但它耗时且耗力，而且高维问题将会挑战设备硬件采样和存储的时间和物理成本，特别是在诸如磁共振成像等具体应用背景中代价较大。将L0模松弛为非凸连续函数的方法虽然取得了很好的效果，但其适用范围相对有限，目前主要用于求解稀疏信号的变换域系数，一般不能在信号域直接解决问题。

已有文献在构造数据稀疏变换方面做了大量的工作。根据文献[8-10]的结果, 一般认为冗余表示系统对数据进行稀疏表示有很好的效果。这些冗余表示系统能够更好地捕捉不同的图像特征，并且比正交基具有更佳的稀疏逼近能力。冗余性也使得表示系统的设计和训练更加灵活，另外可以带来更为鲁棒的图像表示。

本文提出了一种简单有效的基于冗余系统的稀疏信号恢复模型，求解从部分傅里叶系数中恢复原始信号的问题。该方法将L0模与正交投影技术相结合，约束优化模型的目标函数可以有效利用原始数据在给定冗余表示系统中的稀疏性。我们选择紧小波框架来对所恢复信号进行稀疏逼近，并称该方法为投影迭代硬阈值算法(pIHTA)。L0模的邻近算子和正交投影技术被用来求解算法所产生的子问题。由于每个子问题都有闭式解，使得整个算法简洁高效。我们进一步证明了pIHTA算法的全局收敛性，并通过数值实验验证了其有效性。

1 相关数学背景

设D∈Rm×n(m≥n)，如果∃a,b>0，满足

则称D的行向量是Rn空间的一组框架，当a=b时，称其为紧框架。进一步，如果D的列向量规范正交，即DTD=I，则称D的行向量是Rn空间的一组紧小波框架。框架具有完备性，我们可以用其对Rn空间中的任意信号进行表示。由于冗余性的存在，信号表示方法并不唯一。稀疏表示旨在从所有表示方法中寻找所用非零元系数最少的一种，即最稀疏的一种表示方法。

本文主要考虑紧小波框架，设x∈Rn，则α=Dx称为x的小波框架系数，用其可精确的重建信号x=DTα。

设ψ：Rn→(-∞,+∞]是一个正常的下半连续函数，α>0，给定x∈Rn,则由ψ定义的临近算子proxαψ:Rn→Rn定义为

本文所考虑的不完整傅里叶数据的观测模型可描述为

y=PFx+η。

(1)

式中：F∈Cn×n表示离散傅里叶变换矩阵；P∈Rd×n(d

2 投影迭代硬阈值算法

由于d

(2)

由于L0模的非凸不连续性，加之与紧小波框架变换D的复合，模型(2)的数值求解异常困难。已有文献通常将L0模松弛为L1模，进而采用相对成熟的凸优化方法求解[1-3]。但L1模仅是L0模的凸近似，当二者之间的等价性条件不满足时，由此导致的模型误差不可避免。事实上，采用L1模度量稀疏性通常会给恢复信号带来偏差[11]。

本文首先基于约束条件，将L0模和D的复合运算进行分离，在紧小波框架变换域进行求解，给出模型(2)的如下的等价模型

(3)

式中：α是x的小波框架系数；R(D)表示D的值域。显然模型(2)和(3)的解具有如下关系：x=DTα。

为了处理模型(3)中的约束条件，我们引入指示函数lR(α)，

则模型(3)可等价变形为:

(4)

下面将着重讨论等价模型(4)的数值求解。为此我们令

显然g(α)非凸不连续，f(α)连续可微且其梯度具有常数为1的Lipschitz连续性，即：

模型(4)可由临近梯度算法求解[6]，其迭代格式为

αk+1=proxγg(αk-γ∇f(αk))=

+

(5)

式中γ为下降步长。虽然公式(5)是一种相对简单的迭代格式，但约束条件α∈R(D)使得αk+1并无闭式解，导致算法不高效。为了解决该问题，类似于文献[12]中的策略，我们采用交替优化技术，首先处理L0模对应项，然后再考虑约束条件，对应的迭代格式为

(6)

式中PR(D)是空间R(D)上的正交投影算子。

下面证明公式(6)中两个子问题均具有闭式解。对于第一个子问题由邻近算子的定义可知

(7)

对于第二个子问题，由于DTD=I，显然

注意到

∇f(α)=DFTPT(PFDTα-y)，

(8)

该迭代格式具有闭式解，其计算仅涉及傅里叶变换，紧小波框架变换和硬阈值算子。每步迭代过程中的信号可通过

给出。我们称该格式对应的算法为投影迭代硬阈值算法。

3 收敛性分析

对于向量α∈Rd,我们用N(α)表示α中零元素的索引集，即

N(α)∶={i|(α)i=0,1≤i≤d}。

我们引入函数

并且令

易证u(α)连续可微，并且

∇u(α)=DFTPT(PFDTα-y)+

为了估计∇u(α)的Lipschitz常数，引入常数

我们有下面的引理

式中λi(A)表示A的第i个特征值，由于F的正交性和D的列规范正交性，因此

故

故

综上可得

(9)

由公式(9)可知，对∀n成立

上式结合结论(1)可知结论(2)成立。

由结论(2)，∃K，对∀k>K，有

则Λ是闭凸集，记PΛ为Λ上的投影算子，当k>K时，我们有

4 数值实验

本节我们用不同的数据集，欠采样模式和紧小波框架来进行数值实验，展示pIHTA方法比之前较先进的pFISTA方法[12]所恢复的图像有更好的品质，更快的收敛速度。pFISTA方法在模型(4)中采用L1模度量所恢复信号在紧小波框架域的稀疏性，因而是凸优化问题，便于进行理论分析和数值求解。pFISTA方法是求解相应凸优化模型的一种快速算法，与本文方法的主要区别在于迭代过程中采用软阈值算子更新系数α。此外，我们采用相对误差(RLNE)：

我们需要设置几个算法参数。首先根据多次实验以及之前成果的经验，我们令γ=1，令最大权重λmax为1，λmin为0.000 1。运行算法包括内循环和外循环,内循环更新正则化参数λ，其收缩率为0.2。我们定义

并基于该相对误差设置内、外循环的停机准则。当Diff≤10-3时，λ就会根据收缩率变化，当Diff≤10-5，算法迭代终止。

平移不变离散小波变换(SIDWT)[13-14]，又称非抽取小波，因为可以较好地抑制重建伪影，一般作为典型的紧小波框架用在图像信号的恢复中，我们用具有四个分解层的Daubechies小波作为SIDWT。

本实验所用大脑图像切片(见图1(a))是3T西门子扫描仪用加权涡轮自旋回声序列扫描一个健康志愿者得到的。为了增加不同方法对比结果的说服力，我们同时选取图像处理领域中的标准图像之一“Cameraman”(见图1(b))进行数值实验。我们将根据图2中的掩模进行采样(白黑像素暗示着采样和未采样)。

图1 测试图像

图2 30%采样的二维高斯掩模

图3是在30%采样率的二维高斯掩模下分别使用pFISTA方法[12]和本文方法恢复大脑图片和Cameran图像的结果。如果仅从图像恢复视觉效果来比较，两种方法差别不大，我们重点关注两种算法所恢复图像的量化误差比较以及算法的运行时间。如图4和表1所示，在相同参数设置下，两组图像的数据实验均显示pIHTA方法的重构误差和运行时间(单位：s)优于pFISTA方法。

图3 不同方法的重构图像

图4 不同方法的相对误差和时间对比曲线结果(平移不变离散小波)

表1 两种方法在不同小波框架下的误差及运行时间比较

为进一步验证本文方法的效果，我们选择轮廓波(Contourlet)紧小波框架进行同样的数值实验[15-16]。在相同的参数设置下，运行时间和RLNE的实验结果见表1，同样表明了本文算法的优越性。

5 结语

本文提出了一种投影迭代硬阈值算法来求解基于不完全傅里叶观测数据的信号恢复问题。相应的数学模型采用L0模度量所恢复信号在紧小波框架域的稀疏性。投影迭代硬阈值方法可以有效克服L0稀疏优化模型的NP难问题，同时迭代过程的计算简洁高效。我们证明了该数值算法所产生序列的全局收敛性。数值结果表明与pFISTA方法相比，本文算法具有更好的信号恢复效果，收敛速度也相对更快。本文方法目前尚局限于从不完整的正交变换系数中恢复原始信号，我们将把该方法推广到一般的线性变换情形，用来求解诸如压缩感知、图像复原等更广泛的信号恢复问题。