0 引言

压缩感知理论近年来在各研究领域都得到了广泛的应用，例如医学成像、CT断层扫描、机器学习等。2006年Candès等^[1]指出：当信号具有稀疏特性时，可以通过远小于信号长度的少量观测值来精确地重构稀疏信号。但是一般的自然信号s本身并不是稀疏的，需要在某种稀疏基上进行稀疏表示。不妨给定有限长离散信号s，信号稀疏度为k(即含有k个非零值)，则s可以表示成一组正交基的线性组合：

$ \mathit{\boldsymbol{s}} = \sum\limits_{i = 1}^n {{\psi _i}{x_i}} = {\mathit{\boldsymbol{\psi }}^\rm{T}}\mathit{\boldsymbol{x}} $

(1)

其中ψ=[ψ₁, ψ₂, …, ψ_n]是一组正交基。

Candès等^[2]指出若压缩感知矩阵A满足有限等距性(Restricted Isometry Property, RIP)条件，并且x为稀疏度为k的信号，那么可以通过求解下面的l₀范数最小化问题便可以精确地恢复x:

$ \left\{ \begin{array}{l} \min {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\left\| \mathit{\boldsymbol{x}} \right\|_0}\\ {\rm{s}}{\rm{.t}}{\rm{.}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \mathit{\boldsymbol{Ax}} = \mathit{\boldsymbol{b}} \end{array} \right. $

(2)

其中：‖x‖₀指的是非零元素的个数，A是传感矩阵。

压缩感知主要包括信号的稀疏表示、观测矩阵的设计以及稀疏信号恢复三个方面^[3]，主要是通过非线性的重构算法(最优化方法)来恢复信号。

由于Ax=b是欠定的，‖·‖₀是非凸不可微函数，并且问题(2) 是NP难的，因此通常采用l₁范数来替代l₀范数，那么重构稀疏信号的问题即变成了求解优化问题(3)：

$ \left\{ \begin{array}{l} \min {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\left\| \mathit{\boldsymbol{x}} \right\|_1}\\ {\rm{s}}{\rm{.t}}{\rm{.}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \mathit{\boldsymbol{Ax}} = \mathit{\boldsymbol{b}} \end{array} \right. $

(3)

定义1  定义测量矩阵A的RIP参数δ_k满足下式的最小值δ：

$ {\kern 1pt} {\kern 1pt} (1 - {\delta _k}{\kern 1pt} )\left\| \mathit{\boldsymbol{x}} \right\|_2^2 \leqslant \left\| {\mathit{\boldsymbol{Ax}}} \right\|_2^2 \leqslant (1 + {\delta _k})\left\| \mathit{\boldsymbol{x}} \right\|_2^2 $

其中x为k阶稀疏信号。

若δ_k < 1，则称测量矩阵A满足k阶RIP。特别地，当δ_2k < $\sqrt 2 $-1 < 1时, l₀范数与l₁范数是等价的^[4]。

一般有两大类算法对问题(3) 进行近似求解，即贪婪追踪法和凸松弛法，文献[5-9]提出利用匹配追踪算法(Matching Pursuit, MP)和正交追踪算法(Orthogonal Matching Pursuit, OMP)来求解l₀最小范数问题，大大地提高了计算速度，且易于实现，但是恢复能力不强。由于传统贪婪算法在抗噪方面不是很强，文献[10]提出了具有较强鲁棒性的压缩采样匹配追踪(Compressive Sampling Matching Pursuit, CoSaMP)，但由于CoSaMP算法需要已知信号的稀疏度k，而实际应用中信号的稀疏度k往往是未知的。因此CoSaMP算法在解决稀疏信号重构问题上存在一些问题。MP算法来源于贪婪追踪算法，特点是计算复杂度低，但需要较多的观测值，重构精度低；而OMP算法和CoSaMP算法具有良好的稳定性和理论保证，但仅针对大规模的问题恢复率较高。凸松弛算法，这类方法通过将非凸问题转化为凸问题求解找到信号的逼近，其中最常用的方法就是基追踪(Basis Pursuit, BP)^[11]，该方法提出利用l₁范数替代l₀范数来解决最优化问题。这两大类算法在较高的稀疏度下或在较低的观测度下，都很难对高斯信号进行高效重构。

1 反馈神经网络算法 1.1 反馈神经网络

鉴于上述算法所存在的弊端，本文采用神经网络算法来解决压缩感知理论中的信号恢复问题。主要从以下几个方面展开工作：第一部分是建立非线性等式约束优化问题相应的反馈神经网络(Recurrent Neural Network, RNN)，为便于研究，接着选取合适的凸函数f(x)来逼近‖x‖₁, 随后根据梯度下降思想建立神经动力学方程；第二部分探讨所建立的反馈神经网络的稳定性、收敛性和收敛速度与步长的关系等因素；第三部分给出网络的停时条件以及算法具体步骤；第四部分是随机生成离散信号，利用神经网络重构离散信号，并且与已有的重构算法进行对比，从而说明反馈神经网络算法的有效性和准确性。

不论是l₀范数最小化还是l₁范数最小化问题都可以归结为带约束的优化问题。其中一类通过设计一系列反馈神经网络求解带约束的优化问题的方法统称为神经动力学优化方法。以下是一个非线性等式约束优化问题：

$\left\{ \begin{array}{l} \min\ f(\mathit{\boldsymbol{x}})\\ {\rm{s}}{\rm{.t}}{\rm{.}} \ h(\mathit{\boldsymbol{x}}) = 0 \end{array} \right. $

(4)

其中：x∈Rⁿ，f:Rⁿ → R是目标函数，h:Rⁿ → R^m(m < n)是等式约束条件，假设f是二次可微凸函数，h是Rⁿ上的线性函数。假设问题(4) 至少有一个可行解，由Weierstrass定理知，(4) 至少有一个最小值点。通过Matlab平台，可以利用反馈神经网络求解此类问题。

假设x^*是问题(4) 的一个最优解，M₁是目标函数f(x)的一个下界，即M₁≤f(x^*)。函数^[5]：

$ \begin{array}{l} F(\mathit{\boldsymbol{x}}, {M_1}) = {[{(f(\mathit{\boldsymbol{x}})-{M_1})^ + }]^2} = \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left\{ \begin{array}{l} \frac{1}{2}{(f(\mathit{\boldsymbol{x}}) - {M_1})^2},{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} f(\mathit{\boldsymbol{x}}) \geqslant {M_1}\\ 0,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} f(\mathit{\boldsymbol{x}}) < {M_1} \end{array} \right. \end{array} $

由于F(x, M₁)是连续可微的非减凸函数，并且:

$ \nabla F(\mathit{\boldsymbol{x}}, {M_1}) = \left\{ \begin{array}{l} (f(\mathit{\boldsymbol{x}}) - {M_1})\nabla f(\mathit{\boldsymbol{x}}),{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} f(\mathit{\boldsymbol{x}}) \geqslant {M_1}\\ 0,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} f(\mathit{\boldsymbol{x}}) < {M_1} \end{array} \right. $

作能量函数：$E(\mathit{\boldsymbol{x}})=F(\mathit{\boldsymbol{x}}, {{M}_{1}})+\frac{1}{2}h{{(\mathit{\boldsymbol{x}})}^{\rm{T}}}h(\mathit{\boldsymbol{x}})$，其中x=(x₁, x₂, …, x_n)，显然E(x, M₁)是Rⁿ上的凸函数，且∇E(x, M₁)=∇F(x, M₁)+h(x)^T∇h(x)。

考虑凸优化问题：

$ \min \ E(\mathit{\boldsymbol{x}}, {{M}_{1}})=F(\mathit{\boldsymbol{x}}, {{M}_{1}})+\frac{1}{2}h{{(\mathit{\boldsymbol{x}})}^{\rm{T}}}h(\mathit{\boldsymbol{x}}) $

(5)

根据梯度下降法，可以建立求解问题(5) 的神经网络：

$ \left\{ \begin{align} & \frac{\rm{d}\mathit{\boldsymbol{x}}}{\rm{d}\mathit{t}}=-\nabla E(\mathit{\boldsymbol{x}}, {{M}_{1}}) \\ & \mathit{\boldsymbol{x}}({{t}_{0}})={{\mathit{\boldsymbol{x}}}_{0}} \\ \end{align} \right. $

(6)

如果神经网络系统(6) 是稳定的，那么其平衡点就是优化问题(5) 的最优解。设系统(6) 的存在平衡点x₁^*，则E(x₁^*, M₁)≤E(x^*, M₁)，因为x^*是问题(4) 的一个最优解，所以有E(x^*, M₁)=[(f(x^*)-M₁)⁺]²。因为M₁是目标函数f(x)的一个下界，又不等式：

$ \sqrt{E(\mathit{\boldsymbol{x}}_{1}^{*}, {{M}_{1}})}\leqslant \sqrt{E({{\mathit{\boldsymbol{x}}}^{*}}, {{M}_{1}})}=\frac{1}{2}(f({{\mathit{\boldsymbol{x}}}^{*}})-{{M}_{1}}) $

令${{M}_{2}}={{M}_{1}}+2\sqrt{E(\mathit{\boldsymbol{x}}_{1}^{*}, {{M}_{1}})}\leqslant f({{\mathit{\boldsymbol{x}}}^{*}})$，则M₂是最优目标函数f(x^*)的一个下界。以M₂作为优化问题(4) 的一个新下界，即M₂≤f(x^*)。

类似地，再次作能量函数：

$ E(\mathit{\boldsymbol{x}}, {{M}_{2}})={{[{{(f(\mathit{\boldsymbol{x}})-{{M}_{2}})}^{+}}]}^{2}}+\frac{1}{2}\mathit{h}{{(\mathit{\boldsymbol{x}})}^{\rm{T}}}\mathit{h}(\mathit{\boldsymbol{x}}) $

(7)

得到相应问题(7) 的子凸优化问题：

$ \min \ E(\mathit{\boldsymbol{x}}, {{M}_{2}}) $

(8)

此时以x₁^*作为初始点，建立优化问题(8) 的神经网络：

$ \left\{ \begin{align} & \frac{\rm{d}\mathit{\boldsymbol{x}}}{\rm{d}\mathit{t}}=-\nabla E(\mathit{\boldsymbol{x}}, {{M}_{2}}) \\ & \mathit{\boldsymbol{x}}({{t}_{0}})=\mathit{\boldsymbol{x}}_{1}^{\rm{*}} \\ \end{align} \right.;k=0, 1, \cdots $

(9)

基于逐步迭代连续近似的思想，令：

$ {{M}_{k+1}}={{M}_{k}}+2\sqrt{E(\mathit{\boldsymbol{x}}_{k}^{*}, {{M}_{k}})} $

其中x_k^*是凸优化问题(10) 的最优解：

$ \min \ E(\mathit{\boldsymbol{x}}, {{M}_{k}}) $

(10)

再用x_k^*作为初始点，构造一系列子凸优化问题：

$ \min \ E(\mathit{\boldsymbol{x}}, {{M}_{k+1}}) $

(11)

相应的神经子网络模型：

$ \left\{ \begin{align} & \frac{\rm{d}\mathit{\boldsymbol{x}}}{\rm{d}\mathit{t}}=-\nabla E(\mathit{\boldsymbol{x}}, {{M}_{k+1}}) \\ & \mathit{\boldsymbol{x}}({{t}_{0}})=\mathit{\boldsymbol{x}}_\mathit{k}^{\rm{*}} \\ \end{align} \right.;k=0, 1, \cdots $

(12)

此神经网络的特点是用前一个子优化问题对应的神经网络的平衡点x_k^*构造新的能量函数E(x, M_k+1)，并建立初值为x_k^*子神经网络(12)。这一系列神经网络称为反馈神经网络(RNN)。

1.2 稀疏信号的恢复

结合式(3) 和式(4)，考虑建立反馈神经网络来解决问题(3)，由于算子‖·‖₁在x_i=0处不可微，因此不妨考虑凸函数f(x_i)=ln(cosh(ax_i))/a(其中a > 1) 来近似逼近‖x‖₁，从图 2可以看出a越大，该近似效果越精确：

$ {{\left\| \mathit{\boldsymbol{x}} \right\|}_{1}}\approx f(\mathit{\boldsymbol{x}})=\sum\limits_{i=1}^{n}{\frac{\rm{ln}(\rm{cosh}(\mathit{a}{{\mathit{x}}_{\mathit{i}}}))}{a}} $

对于函数y=ln(cosh(ax))/a，有y′=tanh(ax)，并且双曲正切函数是严格单调增函数。

值得注意的是，双曲正切函数在神经网络中经常被当作激活函数使用。此时问题(3) 变成了：

$ \left\{ \begin{align} & \min \frac{1}{a}\sum\limits_{i=1}^{n}{\ln (\cosh (\mathit{a}{{\mathit{x}}_{\mathit{i}}}))} \\ & \rm{s}\rm{.t}\rm{.}\ \ \ \ \mathit{\boldsymbol{b}}\rm{-}\mathit{\boldsymbol{Ax}}=0 \\ \end{align} \right. $

(13)

其中：目标函数是凸函数，约束函数b-Ax是连续函数，因此由凸分析理论知，优化问题(13) 的局部最优解即是全局最优解。相应可建立问题(13) 的Hopfield能量函数：$E(\mathit{\boldsymbol{x}}, {{M}_{1}})=F(\mathit{\boldsymbol{x}}, {{M}_{1}})+\frac{1}{2}\left\| \mathit{\boldsymbol{Ax}}-\mathit{\boldsymbol{b}} \right\|$，由于其中(ln(cosh(ax_i)))/a≥0(i=1, 2, …, n)，取M₁=0，则:

$ F(\mathit{\boldsymbol{x}}, {{M}_{1}})=\frac{1}{2}f{{(\mathit{\boldsymbol{x}})}^{2}} $

所以$E(\mathit{\boldsymbol{x}}, {{M}_{1}})=\frac{1}{2}f{{(\mathit{\boldsymbol{x}})}^{2}}+\frac{1}{2}\left\| \mathit{\boldsymbol{Ax}}-\mathit{\boldsymbol{b}} \right\|$

易知E(x, M₁)≥0是凸函数，则：

$ \nabla E(\mathit{\boldsymbol{x}}, {{M}_{1}})=f(\mathit{\boldsymbol{x}})\nabla f(\mathit{\boldsymbol{x}})+{{(\mathit{\boldsymbol{Ax}}-\mathit{\boldsymbol{b}})}^{\rm{T}}}\mathit{\boldsymbol{A}} $

那么对于凸优化问题：min E(x, M₁)

根据梯度下降法，建立神经动力学方程：

$ \left\{ \begin{align} & \frac{\rm{d}\mathit{\boldsymbol{x}}}{\rm{d}\mathit{t}}=-f(\mathit{\boldsymbol{x}})\nabla f(\mathit{\boldsymbol{x}})-{{(\mathit{\boldsymbol{Ax}}-\mathit{\boldsymbol{b}})}^{\rm{T}}}\mathit{\boldsymbol{A}} \\ & {\mathit{\boldsymbol{x}}}({{t}_{0}})={{\mathit{\boldsymbol{x}}}_{0}} \\ \end{align} \right. $

(14)

2 稳定性与收敛性分析

为了探讨RNN的稳定性与收敛性能，下面给出两个定理并证明：

定理1  对任何初值x₀∈Rⁿ，神经网络(14) 存在唯一最优解x^*的充分必要条件是^[9]：$\left\{ \begin{align} & \nabla E({{\mathit{\boldsymbol{x}}}^{*}}, M)=0 \\ & \mathit{\boldsymbol{b}}-\mathit{\boldsymbol{A}}{{\mathit{\boldsymbol{x}}}^{*}}=0 \\ \end{align} \right.$。

定理2  设x^*是问题(13) 的一个最优解，M是问题(13) 的一个下界，即M≤f(x^*)，定义一个能量函数序列$E(\mathit{\boldsymbol{x}}, {{M}_{k}})={{[{{(f(\mathit{\boldsymbol{x}})-{{M}_{k}})}^{+}}]}^{2}}+\frac{1}{2}{{\left\| \mathit{\boldsymbol{Ax}}-\mathit{\boldsymbol{b}} \right\|}^{2}}$，x_k^*是系统(14) 的一个平衡点，即为问题(13) 的最优解，且序列f(x_k^*)收敛到问题(13) 的一个最优值。

证明  设x^*是问题(13) 的一个最优解，$\mathit{\boldsymbol{\hat{x}}}$是问题(13) 的可行解，若目标函数f(x)≤f(x^*), 则$\mathit{\boldsymbol{\hat{x}}}$也是问题(13) 的最优解且f(x)=f(x^*)，下证x^*收敛到(13) 的一个最优解。

1) 先证f(x_k^*)≤f(x^*)。

因为${{M}_{k+1}}={{M}_{k}}+2\sqrt{E(\mathit{\boldsymbol{x}}_{k}^{*}, {{M}_{k}})}\geqslant {{M}_{k}}$，所以序列M_k是单调非减函数列，又有假设知M₁ < f(x^*)。由数学归纳法可证M_k+1≤f(x^*)，即序列M_k是单调非减有上界序列。即：

$ \begin{align} & f({{\mathit{\boldsymbol{x}}}^{*}})\geqslant {{M}_{k+1}}={{M}_{k}}+2\sqrt{E(\mathit{\boldsymbol{x}}_{k}^{*}, {{M}_{k}})}\geqslant \\ & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {{M}_{k}}+2{{[f(x_{k}^{*})-{{M}_{k}}]}^{+}} \\ \end{align} $

① 当f(x_k^*)≥M_k时，f(x^*)≥M_k+f(x_k^*)-M_k≥f(x_k^*)。

② 当f(x_k^*) < M_k时，f(x^*)≥M_k > f(x_k^*)，所以f(x_k^*)≤f(x^*)。

2) 下证序列{x_k^*}收敛到问题(13) 的一个可行解。

由1) 知，序列{M_k}是单调非减有上界序列，由极限存在定理知，序列{M_k}必存在极限，即存在M^*，使得$\begin{eqnarray*}\underset{k\to \infty }{\mathop{\lim }}\, {{M}_{k}}={{M}^{*}}\end{eqnarray*}$，且M^* < f(x^*)。

再对等式${{M}_{k+1}}={{M}_{k}}+2\sqrt{E(\mathit{\boldsymbol{x}}_{k}^{*}, {{M}_{k}})}$两边同时取极限，得$\begin{eqnarray*}\underset{k\to \infty }{\mathop{\lim }}\, E(\mathit{\boldsymbol{x}}_{k}^{*}, {{M}_{k}})=0\end{eqnarray*}$。

又因为$E(\mathit{\boldsymbol{x}}_{k}^{*}, {{M}_{k}})\rm{=}{{\left[{{\rm{(}\mathit{\boldsymbol{f}}(\mathit{\boldsymbol{x}}_{\mathit{k}}^{*})-{{M}_{\mathit{k}}}\rm{)}}^{+}} \right]}^{2}}+\frac{1}{2}||\mathit{\boldsymbol{Ax}}_{\mathit{k}}^{*}-\mathit{\boldsymbol{b}}|{{|}^{2}}$，所以$\mathop {\lim }\limits_{k \to \infty } {{\rm{[(}}f(\mathit{\boldsymbol{x}}_k^*)-{M_k}{{\rm{)}}^ + }]^2} = 0$且$\mathop {\lim }\limits_{k \to \infty } \frac{1}{2}||\mathit{\boldsymbol{Ax}}_\mathit{k}^* - \mathit{\boldsymbol{b}}|{|^2}{\rm{ = }}0$。即E(x_k^*, M_k) → 0(k → ∞)，即{x_k^*}收敛到问题(13) 的一个可行解。

所以序列f(x_k^*)收敛到问题(13) 的一个最优值。证毕。

3 停时条件

神经网络(14) 的平衡点序列x_k^*收敛到问题(13) 的一个最优解, 从定理2的证明可看出该系统的停时条件是$\mathop {\lim }\limits_{k \to \infty } E(\mathit{\boldsymbol{x}}_k^*, {M_k}) = 0$，即给定ε，E(x_k^*, M_k) < ε时，可得到问题(13) 的最优解。式(14) 中的-∇E(x, M₁)是一个非线性的函数，因此考虑使用数值解法，传统的定步长方法为了达到精度要求，其步长会设置得较小。而变步长的方法随着对解的搜索自适应地调整求解步长，每步迭代使用外推方法得到估计误差，当估计误差较大时减小步长，估计误差较小时则适当增加步长。考虑到速度和求解精度，本文选用BS23方法^[10]，这是一种带步长控制的龙格-库塔方法。在网格0=t₀ < t₁ < t₂… < t_n上逐步迭代计算x(t_k)，直到得到收敛值x(t_n)，每一个x(t_k+1)以及步长Δt_k+1的值都是基于上一个值x(t_k)和步长Δt_k计算得到。

重构稀疏信号x的主要思想是通过建立反馈神经网络模型，再由龙格-库塔方法，利用Matlab平台逐步迭代寻找问题(3) 的最优解和最优值，直到满足停时条件。下面给出RNN算法的具体步骤。

步骤1  初始化。设t=0，取初始点x(t₀)=x₀∈Rⁿ，给定步长Δ > 0和ε∈[10^-15, 10^-5]，k:=1, M₁≤f(x^*)。

步骤2  计算梯度：u(t)=(f(x)-M_k)∇f(x)+(Ax-b)^TA。

步骤3  状态更新：x(t+Δt)=x(t)-Δt·u(t)。

步骤4  计算：$l = {\left\| {\mathit{\boldsymbol{u}}(t)} \right\|_2} = \sum\limits_{i = 1}^n {u_i^2(t)} $。

步骤5  递归与停止：若l < ε，则令x_k=x(t)并输出x_k，停止计算，并输出最优解x(t)和最优值E(x_k)；否则令k=k+1，t=t+Δt，${M_{k + 1}} = {M_k} + 2\sqrt {E({\mathit{\boldsymbol{x}}_k}, {M_k})} $并转到步骤2继续迭代，直到满足停时条件。

4 仿真实验

通过生成仿真信号，采用神经网络模型(14) 求解状态向量$\mathit{\boldsymbol{\hat x}}$，来检验该模型的收敛性质。对于算法求解精度，本文采用总体相对误差R_total，定义如下：

$ {\kern 1pt} {R_{{\rm{total}}}} = {\left\| {\mathit{\boldsymbol{x}} - \mathit{\boldsymbol{\hat x}}} \right\|_2}/{\left\| \mathit{\boldsymbol{x}} \right\|_2} $

其中：x表示原始稀疏信号向量，$\mathit{\boldsymbol{\hat x}}$表示通过RNN求得的解。若R_total < 10^-3，则认为信号恢复正确。对于实验数据的生成方式如下：

1) 产生k稀疏信号x∈Rⁿ, k个非零元素的位置是随机产生的，满足[1, n]的均匀随机分布。相应的非零元素的大小也是随机产生的。

3) 计算观测向量：b=Ax，其中b∈R^m。

由文献[2]知，当矩阵A满足RIP条件时，精确恢复信号x所需的观测量的数量m满足：m≥Ck log(n/k)即可，其中C是常数，即m只与问题的规模参数组合(n, k)有关。取实验参数为：测量矩阵大小m=64，n=256，原始信号稀疏度k=10。根据上述数据的生成方法，可以产生如图 4和图 5所示的原始稀疏信号和观测向量。

2) 产生观测矩阵A∈R^m×n，矩阵的所有元素是随机生成的并且服从高斯分布N(0, 1)，rank(A)=m。

由图 6可知，在b与A给定的情况下，x在神经网络(14) 的不断反馈之下，最终重构并且非零元素均收敛到既定位置上的元素。为了检验不同稀疏度下和观测次数m与正确重构概率的关系，不妨取稀疏度k=[4,12,20,28,36]。

每一组(k, m, n)，执行200次随机实验，由图 7可知，对于不同稀疏度k，当观测点数m大于等于110时，RNN方法能正确恢复稀疏信号的概率极高。而由图 8可知，对于不同观测点数的m，当稀疏度k小于等于20时，RNN方法能正确恢复稀疏信号的概率极高。

由表 1知，模型(8) 只需要较少的观测次数就可以正确地恢复稀疏信号，也可以看出，当m大于100时，其他算法也获得较低的恢复误差，说明当观察次数m足够大时，通常适用的算法同样也能获得较精确的恢复效果。

数据产生之后根据观测向量b通过RNN(14) 进行稀疏信号恢复，选择初始误差容限ε=10^-3，为一般性考虑，取随机向量作为初始x₀，选取实验参数令k=10，m=64，n=256，经神经网络(14) 不断反馈求解状态向量${\mathit{\boldsymbol{\hat x}}}$，状态变化曲线如图 6所示。当相对误差R_total < 10^-3时, 认为这次实验成功, 由图 6知, RNN模型经过迭代之后都能得到收敛解, 证明了该系统有着良好的恢复性能。(其中${\mathit{\boldsymbol{\hat x}}}$是x的重构结果)。

5 结语

本文基于目标函数的性质，设计快速重构稀疏信号的反馈神经网络算法，通过讨论能量函数的梯度，证明该RNN模型最优解的存在性与原问题的近似等价性，讨论了网络的稳定性，并且通过实验验证了算法的有效性；但是该算法也有不足之处，反复实验之下，发现最佳的迭代次数会导致运算量增加，从而稀疏信号重构时间加长，因而该网络需要进一步改进。