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

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

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

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

1 相关数学背景

设D∈R^m×nm≥n，如果∃a, b>0，满足

$ a\|\boldsymbol{x}\|_{2}^{2} \leqslant\|\boldsymbol{D} \boldsymbol{x}\|_{2}^{2} \leqslant b\|\boldsymbol{x}\|_{2}^{2}, \forall \boldsymbol{x} \in {\bf{R}}^{n}, $

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

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

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

$ \operatorname{prox}_{\alpha \psi}(\boldsymbol{x}):=\underset{z \in {\bf{R}}^{n}}{\operatorname{argmin}} \frac{1}{2}\|\boldsymbol{z}-\boldsymbol{x}\|_{2}^{2}+\alpha \psi(\boldsymbol{z}) 。$

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

$ \boldsymbol{y}=\boldsymbol{P F} \boldsymbol{x}+\boldsymbol{\eta} \text { 。} $

(1)

式中：F∈C^n×n表示离散傅里叶变换矩阵；P∈R^d×n(d < n)是由n阶单位矩阵的其中d行组成的用于描述观测位置的选择矩阵；y∈C^d表示观测到的部分不准确的傅里叶系数；η∈R^d表示观测噪声，通常假设其服从零均值正态分布。

2 投影迭代硬阈值算法

由于d < n，问题(1)是典型的不适定反问题。因此我们需要在最大后验概率估计的框架下对恢复信号x做出合理的先验假设。信号在紧小波框架下具有稀疏性，被广泛地应用于求解信号和图像处理邻域内的诸多不适定问题^[8]。基于稀疏正则化求解公式(1)的数学模型可描述为

$ \min \limits_{\boldsymbol{x} \in {\bf{R}}^{n}} \frac{1}{2}\|\boldsymbol{P F} \boldsymbol{x}-\boldsymbol{y}\|_{2}^{2}+\lambda\|\boldsymbol{D} \boldsymbol{x}\|_{0} \text { 。} $

(2)

式中：λ为正则化参数，用来平衡数据拟合项和稀疏正则化项的权重；||·||₀为度量稀疏性的L₀模，用于统计向量中非零元的个数。

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

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

$ \min \limits_{\boldsymbol{\alpha} \in R(\boldsymbol{D})} \frac{1}{2}\left\|\boldsymbol{P F D}^{\mathrm{T}} \boldsymbol{\alpha}-\boldsymbol{y}\right\|_{2}^{2}+\lambda\|\boldsymbol{\alpha}\|_{0} 。$

(3)

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

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

$ l_{R}(\boldsymbol{\alpha}):=\left\{\begin{array}{l} 0, \boldsymbol{\alpha} \in R(\boldsymbol{D}) \\ +\infty, \boldsymbol{\alpha} \notin R(\boldsymbol{D}) \end{array}\right.。$

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

$ \min \limits_{\boldsymbol{\alpha} \in {\bf{R}}^{d}} \frac{1}{2}\left\|\boldsymbol{P F D}^{\mathrm{T}} \boldsymbol{\alpha}-\boldsymbol{y}\right\|_{2}^{2}+\lambda\|\boldsymbol{\alpha}\|_{0}+l_{R}(\boldsymbol{\alpha}) 。$

(4)

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

$ \left\{\begin{array}{l} f(\boldsymbol{\alpha})=\frac{1}{2}\left\|\boldsymbol{P F D}^{\mathrm{T}} \boldsymbol{\alpha}-\boldsymbol{y}\right\|_{2}^{2} \\ g(\boldsymbol{\alpha})=\lambda\|\boldsymbol{\alpha}\|_{0}+l_{R}(\boldsymbol{\alpha}) \end{array}\right.。$

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

$ \left\|\nabla f\left(\boldsymbol{\alpha}_{1}\right)-\nabla f\left(\boldsymbol{\alpha}_{2}\right)\right\|_{2} \leqslant\left\|\boldsymbol{\alpha}_{1}-\boldsymbol{\alpha}_{2}\right\|_{2}, \forall \boldsymbol{\alpha}_{1}, \boldsymbol{\alpha}_{2} \in {\bf{R}}^{d} 。$

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

$ \begin{array}{l} \ \ \ \ \ \ \ \ \boldsymbol{\alpha}_{k+1}=\operatorname{prox}_{\gamma{g}}\left(\boldsymbol{\alpha}_{k}-\gamma \nabla f\left(\boldsymbol{\alpha}_{k}\right)\right)= \\ \underset{\boldsymbol{\alpha}}{\operatorname{argmin}} \frac{1}{2}\left\|\boldsymbol{\alpha}-\left(\boldsymbol{\alpha}_{k}-\gamma \nabla f\left(\boldsymbol{\alpha}_{k}\right)\right)\right\|_{2}^{2}+\gamma \lambda\|\boldsymbol{\alpha}\|_{0}+ \\ l_{R}(\boldsymbol{\alpha})=\underset{\boldsymbol{\alpha} \in R(\boldsymbol{D})}{\operatorname{argmin}} \frac{1}{2}\left\|\boldsymbol{\alpha}-\left(\boldsymbol{\alpha}_{k}-\gamma \nabla f\left(\boldsymbol{\alpha}_{k}\right)\right)\right\|_{2}^{2}+ \\ \gamma \lambda\|\boldsymbol{\alpha}\|_{0} 。\end{array} $

(5)

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

$ \begin{aligned} \tilde{\boldsymbol{\alpha}}_{k+1} & \in \underset{\boldsymbol{\alpha}}{\operatorname{argmin}} \frac{1}{2}\left\|\boldsymbol{\alpha}-\left(\boldsymbol{\alpha}_{k}-\gamma \nabla f\left(\boldsymbol{\alpha}_{k}\right)\right)\right\|_{2}^{2}+\\ \gamma \lambda\|\boldsymbol{\alpha}\|_{0} &, \\ \boldsymbol{\alpha}_{k+1} &=P_{R(\boldsymbol{D})}\left(\tilde{\boldsymbol{\alpha}}_{k+1}\right) 。\end{aligned} $

(6)

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

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

$ \tilde{\boldsymbol{\alpha}}_{k+1} \in \operatorname{prox}_{\gamma \lambda\|\cdot\|_{0}}\left(\boldsymbol{\alpha}_{k}-\gamma \nabla f\left(\boldsymbol{\alpha}_{k}\right)\right) 。$

由于||·||₀和||·||₂²的可分性，ã_k+1可以逐元素由硬阈值算子给出，即

$ \begin{aligned} &\ \ \ \ \ \ \ \ \left(\tilde{\boldsymbol{\alpha}}_{k+1}\right)_{i} \in \\ &\left\{\begin{array}{l} \left(\boldsymbol{\alpha}_{k}-\gamma \nabla f\left(\boldsymbol{\alpha}_{k}\right)\right)_{i},\left|\left(\boldsymbol{\alpha}_{k}-\gamma \nabla f\left(\boldsymbol{\alpha}_{k}\right)\right)_{i}\right|>\sqrt{2 \gamma \lambda} \\ \left\{\left(\boldsymbol{\alpha}_{k}-\gamma \nabla f\left(\boldsymbol{\alpha}_{k}\right)\right)_{i}, 0\right\},\left|\left(\boldsymbol{\alpha}_{k}-\gamma \nabla f\left(\boldsymbol{\alpha}_{k}\right)\right)_{i}\right|=\sqrt{2 \gamma \lambda} \\ 0, \text { 其它 } \end{array}\right.。\end{aligned} $

(7)

对于第二个子问题，由于D^TD=I，显然

$ \boldsymbol{\alpha}_{k+1}=\boldsymbol{D} \boldsymbol{D}^{\mathrm{T}} \tilde{\boldsymbol{\alpha}}_{k+1}。$

注意到

$ \nabla f(\boldsymbol{\alpha})=\boldsymbol{D} \boldsymbol{F}^{\mathrm{T}} \boldsymbol{P}^{\mathrm{T}}\left(\boldsymbol{P} \boldsymbol{F} \boldsymbol{D}^{\mathrm{T}} \boldsymbol{\alpha}-\boldsymbol{y}\right), $

所以公式(6)中的迭代可以合并成关于ã_k+1的格式

$ \begin{aligned} &\tilde{\boldsymbol{\alpha}}_{k+1}=\operatorname{prox}_{\gamma \lambda\|\cdot\|_{0}}\left(\boldsymbol { D } \boldsymbol { D } ^ { \mathrm { T } } \tilde { \boldsymbol { \alpha } } _ { k } \gamma \boldsymbol { D } \boldsymbol { F } ^ { \mathrm { T } } \boldsymbol { P } ^ { \mathrm { T } } \left(\boldsymbol{P} \boldsymbol{F} \boldsymbol{D}^{\mathrm{T}} \tilde{\boldsymbol{\alpha}}_{k}-\right.\right.\\ \boldsymbol{y}) &。\end{aligned} $

(8)

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

$ \boldsymbol{x}_{k}=\boldsymbol{D}^{\mathrm{T}} \boldsymbol{\alpha}_{k}=\boldsymbol{D}^{\mathrm{T}} \tilde{\boldsymbol{\alpha}}_{k} $

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

3 收敛性分析

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

$ N(\boldsymbol{\alpha}):=\left\{i \mid(\boldsymbol{\alpha})_{i}=0,1 \leqslant i \leqslant d\right\} 。$

为了证明迭代序列{ã_k}的收敛性，我们首先通过下面的引理给出该序列的一个重要性质。

引理1   设{ã_k}是由公式(8)产生的迭代序列，则下列结论成立：

(1) 只要i∉N(ã_k)，必有|(ã_k)_i|≥ $\sqrt {{\rm{2}}\lambda \gamma } $；

(2)只要N(ã_k)≠N(ã_k+1)，必有||ã_k+1-ã_k||₂≥$\sqrt {{\rm{2}}\lambda \gamma } $。

证明  结论(1)可直接由ã_k的闭式解公式(7)得出。假设N(ã_k)≠N(ã_k+1)，则至少存在一个1≤j≤d，使得(ã_k)_j=0, |(ã_k+1)_j|≥ $\sqrt {{\rm{2}}\lambda \gamma } $和(ã_k+1)_j=0, |(ã_k)_j|≥ $\sqrt {{\rm{2}}\lambda \gamma } $两种情形有且必有其一成立，对这两种情形我们都有|(ã_k+1)_j-(ã_k)_j|≥ $\sqrt {{\rm{2}}\lambda \gamma } $, 因此||ã_k+1-ã_k||₂≥|(ã_k+1)_j-(ã_k)_j|≥ $\sqrt {{\rm{2}}\lambda \gamma } $。

我们引入函数

$ u(\boldsymbol{\alpha})=\frac{1}{2}\left\|\boldsymbol{P F D}^{\mathrm{T}} \boldsymbol{\alpha}-\boldsymbol{y}\right\|_{2}^{2}+\frac{1}{2 \gamma}\left\|\left(\boldsymbol{I}-\boldsymbol{D} \boldsymbol{D}^{\mathrm{T}}\right) \boldsymbol{\alpha}\right\|_{2}^{2}, $

并且令

$ E(\boldsymbol{\alpha})=u(\boldsymbol{\alpha})+\lambda\|\boldsymbol{\alpha}\|_{0} 。$

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

$ \begin{aligned} &\ \ \ \ \nabla u(\boldsymbol{\alpha})=\boldsymbol{D} \boldsymbol{F}^{\mathrm{T}} \boldsymbol{P}^{\mathrm{T}}\left(\boldsymbol{P} \boldsymbol{F} \boldsymbol{D}^{\mathrm{T}} \boldsymbol{\alpha}-\boldsymbol{y}\right)+\\ \frac{1}{\gamma}&\left(\boldsymbol{I}-\boldsymbol{D} \boldsymbol{D}^{\mathrm{T}}\right) \boldsymbol{\alpha} 。\end{aligned} $

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

$ L(\gamma)=\left\|\boldsymbol{D} \boldsymbol{F}^{\mathrm{T}} \boldsymbol{P}^{\mathrm{T}} \boldsymbol{P F} \boldsymbol{D}^{\mathrm{T}}+\frac{1}{\gamma}\left(\boldsymbol{I}-\boldsymbol{D} \boldsymbol{D}^{\mathrm{T}}\right)\right\|_{2}。$

我们有下面的引理

引理2   当0 < γ≤1时，L(γ)= $\frac{1}{\gamma }$；当γ>1时，L(γ)≤1。

证明  令$\boldsymbol{A} = \boldsymbol{D}{\boldsymbol{F}^{\rm{T}}}\left( {{\boldsymbol{P}^{\rm{T}}}\boldsymbol{P} - \frac{1}{\gamma }\boldsymbol{I}} \right)\boldsymbol{F}{\boldsymbol{D}^{\rm{T}}}$，则

$ L(\gamma)=\left\|\boldsymbol{A}+\frac{1}{\gamma} \boldsymbol{I}\right\|_{2}=\max \limits_{i}\left\{\left|\lambda_{i}(\boldsymbol{A})+\frac{1}{\gamma}\right|\right\} 。$

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

$ \lambda_{i}(\boldsymbol{A}) \in\left\{-\frac{1}{\gamma}, 1-\frac{1}{\gamma}, 0\right\} 。$

故

$ L(\gamma)=\max \left\{1, \frac{1}{\gamma}\right\} 。$

下面的定理给出了{ã_k}和{E(ã_k)}的收敛性。

定理1   设0 < γ≤1，{ã_k}由公式(8)产生，则下列结论成立：

(1) {E(ã_k)}严格单调递减且收敛；

(2) $\mathop {{\rm{lim}}}\limits_{k \to + \infty } ||{\boldsymbol{\tilde \alpha }_{k + 1}} - {\boldsymbol{\tilde \alpha }_k}|{|_2} = 0;$

(3) 存在正整数K>0，对∀k≥K，均有N(ã_k)=N(ã_K)。

证明   结合公式(8)可知，ã_k+1的迭代格式可等价变为

$ \begin{gathered} \tilde{\boldsymbol{\alpha}}_{k+1} \in \operatorname{prox}_{\gamma \lambda\|\cdot\|_{0}}\left(\tilde{\boldsymbol{\alpha}}_{k}-\gamma \nabla u\left(\tilde{\boldsymbol{\alpha}}_{k}\right)\right)= \\ \operatorname{argmin} \frac{1}{2}\left\|\boldsymbol{\alpha}-\left(\tilde{\boldsymbol{\alpha}}_{k}-\gamma \nabla u\left(\tilde{\boldsymbol{\alpha}}_{k}\right)\right)\right\|_{2}^{2}+\gamma \lambda\|\boldsymbol{\alpha}\|_{0} 。\end{gathered} $

将上式中的平方项打开，并丢掉常数项，我们有ã_k+1∈argmin $\frac{1}{2}||\boldsymbol{\alpha } - {\boldsymbol{\tilde \alpha }_k}||_{_2}^2\; + \;\gamma \left\langle {\boldsymbol{\alpha } - {{\boldsymbol{\tilde \alpha }}_k}, \;\;\nabla u\left( {{{\boldsymbol{\tilde \alpha }}_k}} \right)} \right\rangle + \gamma \lambda ||\boldsymbol{\alpha }|{|_0}$。

故

$ \begin{aligned} &\qquad \frac{1}{2}\left\|\tilde{\boldsymbol{\alpha}}_{k+1}-\tilde{\boldsymbol{\alpha}}_{k}\right\|_{2}^{2}+\gamma\left\langle\tilde{\boldsymbol{\alpha}}_{k+1}-\tilde{\boldsymbol{\alpha}}_{k}, \nabla u\left(\tilde{\boldsymbol{\alpha}}_{k}\right)\right\rangle+ \\ &\gamma \lambda\left\|\tilde{\boldsymbol{\alpha}}_{k+1}\right\|_{0} \leqslant \gamma \lambda\left\|\tilde{\boldsymbol{\alpha}}_{k}\right\|_{0} 。\end{aligned} $

另一方面，由引理2可知L(γ)=$\frac{1}{\gamma }$1γ，基于文献[6]中的引理5.7，我们有

$ \begin{aligned} &\qquad u\left(\tilde{\boldsymbol{\alpha}}_{k+1}\right) \leqslant u\left(\tilde{\boldsymbol{\alpha}}_{k}\right)+\left\langle\nabla u\left(\tilde{\boldsymbol{\alpha}}_{k}\right), \tilde{\boldsymbol{\alpha}}_{k+1}-\tilde{\boldsymbol{\alpha}}_{k}\right\rangle+\\ &\frac{1}{2 \gamma}\left\|\tilde{\boldsymbol{\alpha}}_{k+1}-\tilde{\boldsymbol{\alpha}}_{k}\right\|_{2}^{2} 。\end{aligned} $

综上可得

$ E\left(\tilde{\boldsymbol{\alpha}}_{k+1}\right) \leqslant E\left(\tilde{\boldsymbol{\alpha}}_{k}\right)-\frac{1}{2 \gamma}\left\|\tilde{\boldsymbol{\alpha}}_{k+1}-\tilde{\boldsymbol{\alpha}}_{k}\right\|_{2}^{2} 。$

(9)

因此，{E(ã_k)}严格单调递减且非负，故其收敛，结论(1)证毕。

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

$ \frac{1}{2 \gamma} \sum\limits_{k=1}^{n}\left\|\tilde{\boldsymbol{\alpha}}_{k+1}-\tilde{\boldsymbol{\alpha}}_{k}\right\|_{2}^{2} \leqslant E\left(\tilde{\boldsymbol{\alpha}}_{n+1}\right)-E\left(\widetilde{\boldsymbol{\alpha}}_{1}\right) 。$

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

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

$ \left\|\tilde{\boldsymbol{\alpha}}_{k}-\tilde{\boldsymbol{\alpha}}_{K}\right\|_{2}<\sqrt{2 \lambda \gamma}, $

由引理1可知N(ã_k)=N(ã_K)。

定理2   {ã_k}收敛到E(α)的局部最优解α^*，且

$ \lim \limits_{k \rightarrow+\infty} E\left(\tilde{\boldsymbol{\alpha}}_{k}\right)=E\left(\boldsymbol{\alpha}^{*}\right) 。$

证明  由定理 1，◫K，对∀k>K，N(ã_k)=N(ã_K+1)。记

$ \Lambda=\left\{\boldsymbol{\alpha} \in {\bf{R}}^{d} \mid N(\boldsymbol{\alpha}) \subset N\left(\tilde{\boldsymbol{\alpha}}_{K}\right)\right\}, $

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

$ \tilde{\boldsymbol{\alpha}}_{k+1}=P_{\Lambda}\left(\tilde{\boldsymbol{\alpha}}_{k}-\gamma \nabla u\left(\tilde{\boldsymbol{\alpha}}_{k}\right)\right) 。$

由文献[10]中定理1可知ã_k收敛于如下凸优化问题$\mathop {{\rm{min}}}\limits_{\boldsymbol{\alpha } \in \wedge } u\left( \boldsymbol{\alpha } \right)$的全局最优解α^*。由N(ã_K)的不变性可知$\mathop {{\rm{min}}}\limits_{k \to + \infty } E\left( {{{\boldsymbol{\tilde \alpha }}_k}} \right) = E\left( {\boldsymbol{\alpha } * } \right)$显然成立。

下证α^*是E(α)的局部极小值点。设Δα是α^*处的增量且||Δα||_∞≤||α^*||_∞，当Δα ∈  Λ时，显然E(α^*)≤E(α^*+Δα)；当Δα∉Λ时，||α^*+Δα||₀≥||α^*||₀+1，由u(α)的连续性可知，◫ε>0，当||Δα||_∞≤min{||α^*||_∞, ε}时，E(α^*)≤E(α^*+Δα)成立。因此，α^*是E(α)的局部极小值点。

4 数值实验

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

$ \text { RLNE: }=\frac{\|\check{\boldsymbol{x}}-\boldsymbol{x}\|_{2}}{\|\boldsymbol{x}\|_{2}} $

来进行定量比较，其中x表示真实图像，$\boldsymbol{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\smile$}} \over x} }$是重构图像。RLNE越小，图像恢复效果越好。本节所有的数值实验都是在装有win10操作系统的华硕笔记本电脑上进行的，CPU为Intel Core i5-8250U。

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

$ \text { Diff: }=\frac{\left\|\boldsymbol{x}_{n}-\boldsymbol{x}_{n+1}\right\|_{2}}{\left\|\boldsymbol{x}_{n}\right\|_{2}}, $

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

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

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

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

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

5 结语

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