压缩感知理论( CS)是近年提出的一种可将模拟信号以较低采样频率转化为数字形式压缩信号的有效途径。压缩感知理论指出, 若采样信号在某个变换域是稀疏的或者本身就是较稀疏的, 则可用一个不相干的观测矩阵, 直接采集少量的观测数据, 通过解一个优化问题便可从少量的采集数据中以较高的概率重构出原始信号。CS理论从另一个角度开辟了一种新型的高速采样理论, 它的出现改变了传统先高速采样后低码率压缩的信息采集模式, 可显著降低数据存储和传输的代价, 以及信号处理时间和计算成本。

从压缩感知理论可知, 要成功实现信号的低速采样及重构, 必须满足: 1)信号是稀疏的或者可稀疏的; 2)观测矩阵必须满足有限等距性质(RIP), 即与稀疏化表示矩阵之间的不相干性。从目前的研究结果来看, 利用Fourier变换、离散余弦变换(DCT)变换、小波变换等较好地实现了信号的稀疏化表示。此外, 利用正交基级联字典、字典学习等方法同样能够较为有效地实现复杂信号的稀疏化表示。而由于观测矩阵必须满足RIP条件, Candes曾在文献[1-2]中指出:在独立同分布的高斯随机变量形成的观测矩阵与任意稀疏化表示矩阵存在较强的不相干性, 压缩观测信号可较高概率进行精确恢复。由于随机矩阵具备优良的不相干特性, 目前的压缩感知研究过程中均被作为观测矩阵而得到广泛使用。作为一种普适的观测矩阵, 随机矩阵在实际使用过程中, 其计算复杂度较高, 特别是针对高维数据采集时, 存在占用内存较多, 无法适用于大规模应用等问题^[3]。

面对这一问题, 文献[4]提出了一种结构化随机矩阵, 这种矩阵所需的存储量小, 但构造结构随机矩阵过程涉及到随机排序、快速正交变换、随机下采样等操作, 存在不易实现等问题。文献[5-7]分别提出了一种基于分块压缩感知(block CS)的图像重构模型, 这种分块压缩感知提供了一种可以有效降低观测矩阵所需存储空间大的问题, 但对于稀疏度不均匀的原始图像, 需要改变块的大小和测量率, 这在某种程度上引入了一定的复杂度。同样为了降低观测矩阵的存储空间占用量, 文献[8]提出了一种将观测矩阵数值进行二值化处理的方法, 即将原先采用浮点数表示的观测矩阵元素采用二值化进行表示, 大大降低了观测矩阵所需的存储空间。

同样为了尝试降低观测矩阵所需存储空间, 本文提出一种使用半张量积(STP)的压缩感知方法。利用该方法, 成倍降低观测矩阵维数, 可有效降低观测矩阵所需的存储空间。

为验证基于半张量积压缩感知模型的有效性, 结合l^q(0 < q < 1)范数的迭代重加权最小二乘解的方法, 讨论了观测矩阵进行降维后的l^q(0 < q < 1)范数迭代重加权最小二乘重构方法, 并对算法的收敛性进行了讨论。最后针对1维稀疏信号和2维图像信号进行了关于重构概率、重构信号的峰值信噪比等指标的测试和比较。实验结果表明, 采用本文所述方法, 能够利用降维观测矩阵实现对信号的重构, 且重构信号的峰值信噪比并没有明显下降, 但观测矩阵所占的存储空间却可以成倍地降低。

半张量积是一种新型的矩阵乘法^[9]。对于传统矩阵乘法, 当A的列数(C(A))与B的行数(R(B))相等, 则矩阵A、B可乘, 而矩阵半张量积是将传统矩阵乘法推广到阶数不等情形, 即C(A)≠R(B), 从而半张量积运算可以实现前边矩阵的列数和后边矩阵的行数不等时的矩阵乘法运算。更为重要的是, 半张量积保持了传统矩阵乘法的全部主要性质, 所以说它是一种有力的数学工具, 对于实现不同阶的高维矩阵数字信号处理提供了一个非常好的途径。

定义1设T是一个np维行向量, X是一个p维列向量。将T分割成p个等长的块T¹, T², …, T^p, 它们都是n维行向量, 则定义左半张量积为

$ T \ltimes X = \sum\limits_{i = 1}^p {T^i}{x_i} \in {\mathbb{R}^n} $

(1)

式中，符号$\ltimes$表示向量的左半张量积。

定义2设X=[x₁, …, x_s]是一个行向量, Y=[y₁, …, y_t]^T是一个列向量, 根据两个向量的维度, 分别有以下两种情形:

1)如果t是s的因子, 即s=t×n, 则n维行向量

$ {\langle X, Y\rangle _L}: = \sum\limits_{k = 1}^t {{X^k}} {y_k} \in {\mathbb{R}^n} $

(2)

称为X和Y的左半张量内积, 式中$X = [{X^1}, \cdots, {X^t}], {X^i} \in {\mathbb{R}^n}, i = 1, 2, \cdots, t$。

2)如果s是t的因子, 即t=s×n, 则n维列向量

$ {\langle X, Y\rangle _L}: = {({\langle {Y^T}, {X^T}\rangle _L})^T} \in {\mathbb{R}^n} $

(3)

也称为X和Y的左半张量内积。

定义3设$M \in {\mathbb{R}^{m \times n}}, N \in {\mathbb{R}^{p \times q}}$。如果n是p的因子或者p是n的因子, 则有

$ C = M \ltimes N $

(4)

称C是M和N的左半张量积。式中$C = ({C^{ij}}) = {R_i}(M) \ltimes {C_j}(N), i = 1, 2, \cdots, m, j = 1, 2, \cdots, q$。

1)对于定义3所示的M和N, 若n=p, 则有$M \ltimes N = MN$, 即半张量积退化称为普通矩阵积。

2)两个矩阵的左半张量积的维数可以很容易地根据前一个矩阵列数和后一个矩阵行数的公因子消去来计算得到, 如

${F_{p \times qr}} \ltimes{G_{r \times s}} \ltimes{H_{qst \times l}} = {(F \ltimes G)_{p \times qs}} \ltimes{H_{qst \times l}} = {(F \ltimes G \ltimes H)_{pt \times l}} $

(5)

性质1只要$\ltimes$有定义, 即矩阵有合适的维数, 则$\ltimes$满足:

1)分配律

$\left\{ \begin{gathered} F \ltimes (aG \pm bH){\text{ = }}aF \ltimes G \pm bF \ltimes H \hfill \\ (aF \pm bG) \ltimes H = aF \ltimes H \pm bG \ltimes H \hfill \\ \end{gathered} \right. $

(6)

式中，$a, b \in \mathbb{R}$。

2)结合律

$(F \ltimes G) \ltimes H = F \ltimes (G \ltimes H) $

(7)

性质2

1)若F、G均可逆, 则

${(F \ltimes G)^{-1}} = {F^{-1}} \ltimes {G^{-1}} $

(8)

2)

${(F \ltimes G)^{\text{T}}} = {F^T} \ltimes {G^{\text{T}}} $

(9)

现有的压缩感知模型可以简单描述如下:设信号$x \in {\mathbb{R}^{N \times 1}}$在某个正交基矩阵Ψ_N×N下能够进行稀疏化表示, 即

${x_{N \times 1}} = {\psi _{N \times N}}{\theta _{N \times 1}} $

(10)

式中，正交基矩阵满足$\psi {\psi ^{\text{T}}} = {\psi ^{\text{T}}}\psi = I$, 信号经稀疏化表示后的系数${\theta _i} = \langle x{\text{, }}{\psi _i}\rangle = \psi _i^{\text{T}}x$。

若系数向量θ_N×1中仅有K个非零分量, 且K≤N, 则称信号x或系数向量θ是K-稀疏的。若存在一个与稀疏化表示矩阵Ψ_N×N不相干的矩阵Φ_M×N(其中$M \ll N$)作为压缩感知的观测矩阵, 则压缩感知模型可简单表示为

${y_{M \times 1}} = {\Phi _{M \times N}}{\psi _{N \times N}}{\theta _{N \times 1}} $

(11)

如式(11)所示, 对信号进行一次压缩观测, 便可以得到M个线性观测值$y \in {\mathbb{R}^{M \times 1}}$, 利用这些少量信号和系数向量θ是K-稀疏的约束条件, 通过凸优化或线性规划的方法可重构信号x。

在这个模型的基础上, 观测矩阵参与了CS处理方法中压缩观测与信号恢复这两个重要的数据处理过程, 且对整个方法的性能有着很大的影响, 因而观测矩阵的设计是压缩感知的一个核心问题。从式(11)可知, 观测矩阵的大小M×N仍与原始数据的维数N相关, 因此当面对高维数据时, 即N特别大, 观测矩阵Φ所占的存储空间较大, 信号恢复过程的计算量为维持在较高的水平。

本着为降低压缩感知过程中观测矩阵所占的存储空间以及降低计算量的角度, 提出一种应用半张量积的压缩感知模型, 使其有效降低观测矩阵在压缩感知以及重构过程中所占的存储空间容量。半张量积压缩感知模型定义为

${y_{M \times 1}} = {\Phi _{(M/t) \times (N/t)}} \ltimes {\psi _{N \times N}}{\theta _{N \times 1}} $

(12)

式中，${\Phi _{(M/t) \times (N/t)}}$是降维后的随机观测矩阵, 其维数为$(M/t) \times (N/t) \in {\mathbb{Z}^ + }, t \in {\mathbb{Z}^ + }$且t < M。此外, 从式(12)所示的压缩感知模型可知, 在相同数据类型的情形下, 其观测矩阵Φ占用的存储空间是式(11)所示观测矩阵所占存储空间的1/t²倍。

经过上述简单定义后, 利用式(5)(7), 可以验证式(12)是满足式(11)所示的压缩感知模型的, 同样用M个观测值可以对N维信号进行表示，即

${\Phi _{(M/t) \times (N/t)}} \ltimes ({\psi _{N \times N}} \ltimes {\theta _{N \times 1}}) = {\Phi _{(M/t) \times (N/t)}} \ltimes {(\psi \ltimes \theta)_{N \times 1}} = \\{(\Phi \ltimes \psi \ltimes \theta)_{M \times 1}} = {y_{M \times 1}} $

或

$({\Phi _{(M/t) \times (N/t)}} \ltimes {\psi _{N \times N}}) \ltimes {\theta _{N + 1}} = {(\Phi \ltimes \psi)_{M \times N}} \ltimes {\theta _{N + 1}} = \\ {\rm{ }}{(\Phi \ltimes \psi \ltimes \theta)_{M \times 1}} = {y_{M \times 1}} $

至此, 本文所述问题的关键在于, 如何采用式(12)所示的压缩感知模型来对N维信号进行恢复或重构。

若x_N×1本身是K-稀疏的, 则式(12)所示模型可用矩阵乘积形式表示为

$\left( \begin{matrix} {{y}_{1}} \\ {{y}_{2}} \\ \vdots \\ {{y}_{M}} \\ \end{matrix} \right)=\left( \begin{matrix} {{\phi }_{1,1}} & {{\phi }_{1,2}} & \cdots & {{\phi }_{1,N/t}} \\ {{\phi }_{\text{2},1}} & {{\phi }_{\text{2},\text{2}}} & \cdots & {{\phi }_{2,N/t}} \\ \vdots & \vdots & {} & \vdots \\ {{\phi }_{M/t,1}} & {{\phi }_{M/t,\text{2}}} & \cdots & {{\phi }_{M/t,N/t}} \\ \end{matrix} \right)\ltimes \left( \begin{matrix} {{x}_{1}} \\ {{x}_{2}} \\ \vdots \\ {{x}_{N}} \\ \end{matrix} \right)$

(13)

式中, ${\phi _{i, j}} \in {\Phi _{(M/t) \times (N/t)}}, i = 1, 2, \cdots, M/t, j = 1, 2, \cdots, N/t, $利用式(1)—(4)所示的半张量积定义, 式(13)表示为

$\left( {\begin{array}{*{20}{c}} {{y_1}}\\ {{y_2}}\\ \vdots \\ {{y_M}} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} {{\xi ^{1,1}} + \cdots + {\xi ^{1,j}} + \cdots + {\xi ^{1,\frac{N}{t}}}}\\ {{\xi ^{i,1}} + \cdots + {\xi ^{i,j}} + \cdots + {\xi ^{i,N/t}}}\\ \vdots \\ {{\xi ^{M/t,1}} + \cdots + {\xi ^{M/t,j}} + \cdots + {\xi ^{M/t,N/t}}} \end{array}} \right)$s

(14)

式中, ${\xi ^{i, j}} = {\phi _{i, j}}{[{X_{(j-1)t + 1}} \cdots {x_{jt}}]^{\rm{T}}}, i = 1, 2, \cdots, M/t, j = 1, 2, \cdots, N/t, $且ξ^{i, j}表示了一个t维的列向量。任意行元素ξ^{i, 1}+…+ξ^{i, j}+…+ξ^{i, N/t}的和也是一个t维的列向量。当观测矩阵维数为M×N时, 即t=1, 则式(14)可表示成${{\rm{y}}_{M \times 1}} = {\Phi _{M \times N}} \cdot {x_{N \times 1}}, $即满足式(11)所示压缩感知模型。

综上所述, 本文所提的基于半张量积的压缩感知模型是适用于式(10)所示的压缩感知模型的。为此, 借鉴目前已得到广泛研究的基于l^q(0 < q < 1)范数的迭代重加权最小二乘法^[10-11](IRLS)来讨论上述基于半张量积的压缩感知模型的重构方法。

相关研究已表明, 基于l^q(0 < q < 1)范数最小化求解, 在无噪声的情况下, 相比l¹范数最小化求解, 有更强的稀疏信号恢复能力, 且l^q范数最小化问题在有噪声干扰的情况下也是稳定的^[12], 在这些理论的基础上, 本文将探讨基于半张量积的l^q(0 < q < 1)范数的迭代重加权最小二乘解的方法。

令观测矩阵${\Phi _{(M/t) \times (N/t)}}$是行满秩矩阵, 即$rank({\Phi _{(M/t) \times (N/t)}}) = M/t, \parallel x{\parallel _0}$表示x的l⁰范数, 即表示x中非零元素的个数。则引入${\Phi _{(M/t) \times (N/t)}}$的所有K-个分量的值域可表示为

$ {{\mathfrak{R}}_{k}}=\left\{ y\in {{\mathbb{R}}^{M}}\left| y= \right.{{\Phi }_{\left( M/t \right)\times \left( N/t \right)}}\ltimes {{x}_{N\times 1}},x\in {{\mathbb{R}}^{N}},{{\left\| x \right\|}_{0}}\le K \right\} $

(15)

同样, 因为M < N, 所以对于任意的$y \in {\mathbb{R}^M}$, 方程${\Phi _{M/t \times N/t}}{x_{N \times 1}} = {y_{M \times 1}}$同样存在无穷多个解, 且令式(15)所有解的集合为$f(y) = {\Phi ^{-1}}(y)$。

若给定$y \in {\mathfrak{R} _k}, K < M/t, $考虑l^q(0 < q < 1)范数的最小化问题

$ \min \left\{ {{\left\| x \right\|}_{q}}\left| x\in {{\mathbb{R}}^{N}} \right.,\Phi \ltimes x=y \right\} $

(16)

式中, $\parallel x{\parallel _q} = {(\sum\limits_{i = 1}^N {|{x_i}{|^q}})^{\frac{1}{q}}}, \Phi \in {\mathbb{R}^{(M/t) \times (N/t)}}$。

在l^q(0 < q < 1)范数最小化的前提下, 考虑式(16)的近似

$ \min \left\{ \frac{1}{q}\sum\limits_{i=1}^{N}{{{\left( x_{i}^{2}+{{\varepsilon }^{1+q}} \right)}^{\frac{q}{2}}}\left| x\in {{\mathbb{R}}^{N}} \right.},\Phi \ltimes x=y \right\} $

(17)

当ε→0时, 式(17)所示的最小化问题即为式(16)的解，即

$ x=\min \left\{ \frac{1}{q}\sum\limits_{i=1}^{N}{{{\left( {{x}^{2}}+{{\varepsilon }^{1+q}} \right)}^{\frac{q}{2}}}\left| x\in {{\mathbb{R}}^{N}} \right.},\Phi \ltimes x=y \right\} $

(18)

对式(18), 引入Largrange乘子$\lambda \in {\mathbb{R}^M}$和Largrange函数, 有

$ F\left( x,\lambda \right)=\frac{1}{q}\sum\limits_{i=1}^{N}{{{\left( {{x}^{2}}+{{\varepsilon }^{1+q}} \right)}^{\frac{q}{2}}}}+{{\lambda }^{T}}\left( \Phi \ltimes x-y \right) $

(19)

则式(18)的解满足

$ \left\{ \begin{array}{l} \frac{{\partial F\left( {x,\lambda } \right)}}{{\partial {x_i}}} = \frac{{{x_i}}}{{{{\left( {x_i^2 + {\varepsilon ^{1 + q}}} \right)}^{\left( {2 - q} \right)/q}}}} + {\left( {{\Phi ^T}\lambda } \right)_i} = 0\\ \frac{{\partial F\left( {x,\lambda } \right)}}{{\partial {\lambda _i}}} = {\left( {\Phi \ltimes x - y} \right)_j} = 0 \end{array} \right. $

(20)

式中，i=1, 2, …, N; j=1, 2, …, M。对于式(20)采用迭代法进行求解，即

$ \left\{ \begin{array}{l} \frac{{x_i^{\left( {n + 1} \right)}}}{{{{\left( {{{\left( {x_i^{\left( n \right)}} \right)}^2} + \varepsilon _n^{1 + q}} \right)}^{\left( {2 - q} \right)/q}}}} + {\left( {{\Phi ^T}{\lambda ^{\left( {n + 1} \right)}}} \right)_i} = 0\\ \Phi \ltimes {x^{n + 1}} - y = 0 \end{array} \right. $

(21)

式中，i=1, 2, …, N。且在迭代过程中, 对于上述l^q(0 < q < 1)范数的最小化问题, 定义函数

$ \begin{array}{l} \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{T_q}\left( {x,w,\varepsilon } \right) = \\ \frac{q}{2}\left[ {\sum\limits_{i = 1}^N {x_i^2{w_i}} + \sum\limits_{i = 1}^N {\left( {{\varepsilon ^2}{w_i} + \frac{{2 - q}}{q}w_i^{\frac{q}{{\left( {q - 2} \right)}}}} \right)} } \right] \end{array} $

(22)

式中，$x \in {\mathbb{R}^N}, w \in \mathbb{R}_ + ^N, \varepsilon \in {\mathbb{R}_ + }, $且迭代权重w定义为

$w_i^{\left( n \right)} = {\left( {{{\left( {x_i^{\left( n \right)}} \right)}^2} + \varepsilon _n^{1 + q}} \right)^{\left( {2 - q} \right)/q}}$

(23)

迭代开始时, 令${w^{(0)}} = (1, \cdots, 1) \in {\mathbb{R}^N}, {\varepsilon _0} = 1$则x⁽ⁿ⁺¹⁾有

${x^{\left( {n + 1} \right)}} = \mathop {\arg \;\min }\limits_{x \in f\left( y \right)} {T_q}\left( {{x^{\left( n \right)}},{w^{\left( n \right)}},{\varepsilon _n}} \right)$

(24)

且迭代过程中, ε_n的更新采用

${\varepsilon _{n + 1}} = \min \left( {{\varepsilon _n},\frac{{r{{\left( {{x^{\left( {n + 1} \right)}}} \right)}_{K + 1}}}}{N}} \right)$

(25)

式中，K表示信号的稀疏度, r(x)表示将向量x中的每个分量取绝对值, 并进行单调递减排序后得到的向量, 即$r{(x)_1} \geqslant r{(x)_2} \geqslant \cdots r{(x)_N} \geqslant 0$。r(x⁽ⁿ⁺¹⁾)_K+1表示x⁽ⁿ⁺¹⁾经单调排序后的第K+1个分量。迭代过程中, 当ε_n=0时, 终止迭代, 并得到相应的稀疏解。

为了得到更明确的表达形式, 定义D_n为N×N的对角矩阵, 其矩阵中第i个对角分量表示为

${D_{i, i}} = \frac{1}{{w_i^{(n)}}}, i = 1, 2, \cdots, N $

则上述线性系统的解可表示为

${x^{\left( {n + 1} \right)}} = {D_n}\ltimes {\Phi ^T}\ltimes {\left( {{\Phi ^T}{D_n}\ltimes {\Phi ^T}} \right)^{ - 1}}\ltimes y$

(26)

式中，$\Phi \in {\mathbb{R}^{(M/t) \times (N/t)}}, {D_n} \in \mathbb{R}_ + ^{N \times N}$。

从而本文所涉及的半张量积l^q(0 < q < 1)范数的迭代重加权最小二乘解的算法如下:

1)初始化

${w^{(0)}} = (1, \cdots, 1), {\varepsilon _0} = 1, {x^{(0)}} = (1, \cdots, 1); $

2)根据w⁽ⁿ⁾得到N×N对角矩阵D_n;

3)根据式(26)计算并更新x⁽ⁿ⁺¹⁾;

4)根据式(25)计算并更新ε_n值;

5)根据式(23)计算并更新w⁽ⁿ⁾值;

6)更新迭代次数, 并判断ε_n是否满足条件, 若不满足, 返回第2)步, 继续迭代计算; 若满足则退出迭代, 并输出稀疏解x⁽ⁿ⁺¹⁾。

上述算法中, 对于第n次迭代(n≥0), T_q(x⁽ⁿ⁾, w⁽ⁿ⁾, ε_n)表示为

$\begin{array}{l} \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{T_q}\left( {{x^{\left( n \right)}},{w^{\left( n \right)}},{\varepsilon _n}} \right) = \\ \frac{q}{2}\left[ {\sum\limits_{i = 1}^N {{{\left( {x_i^{\left( n \right)}} \right)}^2}w_i^{\left( n \right)} + \sum\limits_{i = 1}^N {\left( {\varepsilon _n^2w_i^{\left( n \right)} + \frac{{2 - q}}{q}{{\left( {w_i^{\left( n \right)}} \right)}^{\frac{q}{{\left( {q - 2} \right)}}}}} \right)} } } \right] \end{array}$

(27)

式中，$w_i^{(n)} > 0, {(x_i^{(n)})^2} \geqslant 0, {\varepsilon _n} > 0$。由式(23)—(25)所示的更新方式, 则对于函数${T_q}({x^{(x)}}, {w^{(n)}}, {\varepsilon _n})$具有单调性质

$\begin{gathered} {T_q}({x^{(n + 1)}}, {w^{(n + 1)}}, {\varepsilon _{n + 1}}) \leqslant \hfill \\ {T_q}({x^{(n + 1)}}, {w^{(n)}}, {\varepsilon _{n + 1}}) \leqslant \hfill \\ {T_q}({x^{(n + 1)}}, {w^{(n)}}, {\varepsilon _n}) \leqslant \hfill \\ {T_q}({x^{(n)}}, {w^{(n)}}, {\varepsilon _n}) \hfill \\ \end{gathered} $

则对任意的n≥1, 可以得到

$\left\| {{x^{\left( n \right)}}} \right\|_q^q \le {T_q}\left( {{x^{\left( 1 \right)}},{w^{\left( 0 \right)}},{\varepsilon _0}} \right) = A$

(28)

且

$ w_i^{(n)} \geqslant {A^{-1}}, i = 1, 2, \cdots, N $

式中, A是正的常数。

此外, 对于n≥1, 同样可以得到

$\begin{gathered} {T_q}({x^{(n)}}, {w^{(n)}}, {\varepsilon _n})-{T_q}({x^{(n + 1)}}, {w^{(n + 1)}}, {\varepsilon _{n + 1}}) \geqslant \hfill \\ {T_q}({x^{(n)}}, {w^{(n)}}, {\varepsilon _n})-{T_q}({x^{(n + 1)}}, {w^{(n)}}, {\varepsilon _n}) = \hfill \\ \frac{q}{2}[\sum\limits_{i = 1}^N {{{(x_i^{(n)})}^2}w_i^{(n)}} ]-\frac{q}{2}[\sum\limits_{i = 1}^N {{{(x_i^{(n + 1)})}^2}w_i^{(n)}} ] = \hfill \\ \frac{q}{2}\sum\limits_{i = 1}^N {\{ [{{(x_i^{(n)})}^2}-} {(x_i^{(n + 1)})^2}]w_i^{(n)}\} \geqslant \hfill \\ \frac{q}{2}{A^{-1}}\sum\limits_{i = 1}^N {[{{(x_i^{(n)})}^2}-} {(x_i^{(n + 1)})^2}] \hfill \\ \end{gathered} $

由式(28), 便可得到

$\sum\limits_{i = 1}^N {[{{(x_i^{(n)})}^2}-} {(x_i^{(n + 1)})^2}] \leqslant \frac{{2{A^2}}}{q} $

式中，0 < q < 1, 从而可以得到迭代序列是收敛的, 即

$\mathop {\lim }\limits_{x \to \infty } ({x^{(n)}}-{x^{(n + 1)}}) = 0 $

为验证上述基于半张量积压缩感知算法的迭代收敛速度及精确恢复稀疏解的能力, 分别设计了两组实验进行验证。第1组选择1维稀疏信号进行重构, 并从1维信号在不同稀疏度情况下的恢复概率, 恢复信号的信噪比以及相同稀疏度情形下, 重构信号的收敛速度等角度进行比较。第2组实验选择了2维图像信号进行重构, 并从重构信号峰值信噪比的角度对重构性能与观测矩阵的降维比例进行比较。实验平台配备了Intel i7-4600 CPU, 2.1 GHz主频, 8 GB内存, 64位Windows 8操作系统, 仿真软件采用Matlab R2010b。

实验时, 构建1维稀疏信号

$\begin{gathered} x = zeros(N, 1) \hfill \\ p = randperm(N) \hfill \\ x(p(1:K)) = randn(K, 1) \hfill \\ \end{gathered} $

式中，K=1, 2, …, M/2, M为信号的观测次数, N=256, 则x便是表示一个长度为256, 稀疏度为K, 且非零位置为随机的1维稀疏信号。此处设定K的最大值为M/2, 原因是文献[13]中已指出, 一个向量要保证稀疏恢复的唯一性, 其非零元素的个数最多不超过一个极限值M/2。

为验证对上述1维稀疏信号重构性能, 设定测量次数M为128。仿真次数为500次。实验时构建了满足$N(0, \frac{1}{{(M/t)}})$高斯独立分布的$(M/t) \times (N/t)$维高斯随机观测矩阵${\Phi _{(M/t) \times (N/t)}}, $如表 1所示。

由表 1可知, 当t=1时, 其观测矩阵维数满足式(11)所示的压缩感知模型, 并没有实现降低观测矩阵所占存储空间的目的, 但当t=2和4时, 其观测矩阵大小分别是t=1时所对应观测矩阵的1/4和1/16。若随机观测矩阵中每一元素采用单精度浮点数进行表示, 即每个元素用4个字节表示, 则对于表 1中t=1所对应的观测矩阵需要128 KB字节的存储空间, 而t=2或4时, 对应观测矩阵仅需32 KB或8 KB字节, 从而大大降低了观测矩阵所需的存储空间以及重构过程中的内存占用率。

文献[11]中证明了基于l^q(0 < q < 1)范数的迭代重加权最小二乘解在q值较大时, 其恢复稀疏解的能力相对较强, 因此, 后续实验中均选择q=0.8来验证半张量积压缩感知重构方法对稀疏解的恢复情况。

图 1表示了观测次数M=128时, 不同稀疏度的原始信号x在不同观测矩阵的情况下其正确恢复的概率。实验中, 对每个稀疏度K(K=1, 2, …, M/2), 都随机生成了500次的稀疏向量x, 并计算其恢复的成功率。实验时, 选取算法的终止条件是ε_n < 10^－8, 即在每次重构过程中, 设ε₀=1, 若经过若干次迭代后, 有ε_n < 10^－8, 则退出重构。重构结束后, 计算重构信号 ${\hat x}$与原始信号x的相对误差E, 若相对误差E < 10^－5, 则认为重构是成功的, 反之则认为本次重构是错误的。1维信号的相对误差E以及重构概率为

$E = \frac{{\left\| {x - \hat x} \right\|}}{{{{\left\| x \right\|}_2}}}$

(29)

$pr = \frac{{{T_{correct}}}}{{{T_{total}}}}$

(30)

式中, T_correct表示正确重构的次数, T_total表示重构的总次数, x表示原始信号, ${\hat x}$表示重构信号。

图 1中, t=1表示采用的观测矩阵为Φ_128×256, t=2和t=4分别表示观测矩阵为Φ_64×128和Φ_32×64。从图 1可知, 当观测矩阵维数降低后, 即t=2和4时, 且当K值较大时, 对稀疏信号的重构概率相比t=1时有一定的下降, 但其整体上保持了较高的重构概率。此外从图 1中可知, 对观测矩阵进行降维后, 采用上述基于半张量积压缩感知模型对信号进行观测, 是可以实现对原始信号的精确重构。至于图 1中重构概率的降低, 其主要原因在于, 随着观测矩阵维数的降低, 对于每一个观测值y_i, 仅用了N/t个系数及未知量x进行表示, 且在相邻的t个观测值中(对于t个相邻观测值, 定义y₁, …, y_t表示为第1组相邻的t个观测值, y_(i－1)×t+1, …, y_i×t表示为第i组的t个相邻观测值, 则对于M个观测值y, 可分为M/t组), 所表示的系数都是相同的, 即

$\begin{array}{*{20}{c}} {\left( {\begin{array}{*{20}{c}} {{y_{\left( {i - 1} \right) \times t + 1}}}\\ {{y_{\left( {i - 1} \right) \times t + 2}}}\\ \vdots \\ {{y_{i \times t}}} \end{array}} \right)}\\ {\left( {\begin{array}{*{20}{c}} {{\phi _{i,1}}{x_1} + \cdots + {\phi _{i,j}}{x_{\left( {j - 1} \right) \times t + 1}} + \cdots + {\phi _{i,N/t}}{x_{\left( {N/t - 1} \right) \times t + 1}}}\\ {{\phi _{i,1}}{x_2} + \cdots + {\phi _{i,j}}{x_{\left( {j - 1} \right) \times t + 2}} + \cdots + {\phi _{i,N/t}}{x_{\left( {N/t - 1} \right) \times t + 2}}}\\ \vdots \\ {{\phi _{i,1}}{x_t} + \cdots + {\phi _{i,j}}{x_{\left( {j - 1} \right) \times t + t}} + \cdots + {\phi _{i,N/t}}{x_{\left( {N/t - 1} \right) \times t + t}}} \end{array}} \right)} \end{array}$

(31)

式中，y_(i－1)×t+1, …, y_i×t表示相邻的t个观测值, 每一个观测值表示未知量x的系数都是(ø _{i, 1}, …, ø _{i, N/t})。若N=10, M=6, t=2时, 构建的随机观测矩阵为Φ_3×5, 便有

$\left( \begin{matrix} {{y}_{1}} \\ {{y}_{2}} \\ {{y}_{3}} \\ {{y}_{4}} \\ {{y}_{5}} \\ {{y}_{6}} \\ \end{matrix} \right)\text{=}\left( \begin{matrix} {{\phi }_{1, 1}}{{x}_{1}}+{{\phi }_{1, 2}}{{x}_{3}}+{{\phi }_{1, 3}}{{x}_{5}}+{{\phi }_{1, 4}}{{x}_{7}}+{{\phi }_{1, 5}}{{x}_{9}} \\ {{\phi }_{1, 1}}{{x}_{2}}+{{\phi }_{1, 2}}{{x}_{4}}+{{\phi }_{1, 3}}{{x}_{6}}+{{\phi }_{1, 4}}{{x}_{8}}+{{\phi }_{1, 5}}{{x}_{10}} \\ {{\phi }_{2, 1}}{{x}_{1}}+{{\phi }_{2, 2}}{{x}_{3}}+{{\phi }_{2, 3}}{{x}_{5}}+{{\phi }_{2, 4}}{{x}_{7}}+{{\phi }_{2, 5}}{{x}_{9}} \\ {{\phi }_{2, 1}}{{x}_{2}}+{{\phi }_{2, 2}}{{x}_{4}}+{{\phi }_{2, 3}}{{x}_{6}}+{{\phi }_{2, 4}}{{x}_{8}}+{{\phi }_{2, 5}}{{x}_{9}} \\ {{\phi }_{3, 1}}{{x}_{1}}+{{\phi }_{3, 2}}{{x}_{3}}+{{\phi }_{3, 3}}{{x}_{5}}+{{\phi }_{3, 4}}{{x}_{7}}+{{\phi }_{3, 5}}{{x}_{9}} \\ {{\phi }_{3, 1}}{{x}_{2}}+{{\phi }_{3, 2}}{{x}_{4}}+{{\phi }_{3, 3}}{{x}_{6}}+{{\phi }_{3, 4}}{{x}_{8}}+{{\phi }_{3, 5}}{{\phi }_{10}} \\ \end{matrix} \right)$

式中，相邻两个观测值分别为[y₁, y₂]^T、[y₃, y₄]^T及[y₅, y₆]^T。从上式中可以明显看到, 相邻两个观测值在表示原始信号x时所采用的系数是一样的, 在用最小二乘法进行迭代时, 当原始信号x的稀疏度K较大时, 需要相对较多的迭代次数才能收敛, 即存在收敛速度相对较慢的问题, 从而造成了图 1中所示的重构概率的下降。但当原始信号x的稀疏度K较小时, 即原始信号x中存在大量的零元素, 观测矩阵维数下降对迭代收敛速度影响不大, 其结果如图 2所示。

图 2表示了在相同观测次数及稀疏度的情形下, 即M=128, N=256, K=30时, 采用不同观测矩阵, 经历相同次数迭代后, 其稀疏解的收敛速度。图 2中采用式(29)所示的相对误差E表示稀疏解的收敛速度。在迭代次数相同时, 其相对误差减小的越快表示收敛速度越快, 图 2中数据是经500次实验取相对误差平均值后得到的结果。

由图 2可知, 当稀疏度较低时, 将观测矩阵维数进行降维后, 即t=2或4时, 经过相同的迭代次数, 重构信号的相对误差均能达到相同的数量级, 即收敛速度大致相当。

为进一步验证观测矩阵降维后的重构性能, 进行针对固定原始信号长度N, 固定稀疏度K, 经不同的观测次数, 来比较其正确重构的概率。实验中, 设定N=256, K=30, 重构算法终止条件及判断重构是否成功的条件与图 1所设条件一致, 同样进行了500次实验, 其结果如图 3所示。

从图 3可知, 在稀疏度K相对较小的情形下, 经过相同次数的观测, 其正确重构概率大致相当。至此, 通过上述实验, 充分说明了利用本文所述的基于半张量积压缩感知模型, 对观测矩阵进行降维处理后, 虽然重构概率较t=1时有一定下降, 但对1维信号还是能够进行高概率重构。

在压缩感知与重构过程中, 其对象必须是稀疏信号才能进行精确的重构。2维图像信号在很多变换域下, 能够进行稀疏化表示, 如小波变换, DCT, 快速傅里叶变换(FFT)等等。实验中, 选取小波变换矩阵对原始2维图像进行稀疏化表示, 并将得到的稀疏信号投影到观测矩阵, 得到相应的观测信号y。通过基于l^q(0 < q < 1)范数的迭代重加权最小二乘法, 重构出原始图像的小波系数, 最后在通过小波系数反变换重构出原始图像的重构图像。

实验中采用4幅原始图像作为实验载体, 大小分别是为256×256像素的Lena、Peppers, 以及大小为512×512像素的Mandrill和DICOM医学灰度图像OT-colon, 其中OT-colon医学灰度图像可以从文献[14]得到。

为比较采用降维观测矩阵(t>1)对重构质量的影响, 进一步验证本文算法的可行性, 对4幅原始图像采用相同的采样率, 即M/N=0.5, 同时在设定观测矩阵维数时, 分别取t=1, 2和4。在得到各自的观测信号y后, 利用本文重构算法对原始图像进行50次重构, 分别观测其重构图像和重构图像峰值信噪比。图 4列出了Lena和Peppers的重构结果。

图 4(b)—(d)分别表示了对2维原始图像的重构结果。对应图 4(b)—(d)的高斯随机观测矩阵Φ的大小分别为128×256, 64×128和32×64。从图 4中可以较为直观地看到, 在采用降维后的高斯随机观测矩阵对信号进行采样并重构后, 其重构的2维图像质量并没有明显的下降, 这表明了本文基于半张量积的压缩感知算法对观测矩阵进行降维, 是同样可以对原始信号进行有效重构。且该方法非常直接有效地降低了观测矩阵的大小, 也降低了观测矩阵的存储空间占用量。

为进一步衡量2维图像的重构效果, 计算了各重构图像的峰值信噪比(PSNR), 结果如表 2所示。

表 2中列出了4幅原始图像采用3种不同大小的高斯随机观测矩阵分别经50次采样和重构后, 其重构图像峰值信噪比的最大值, 最小值和平均值。从表 2可知, 对高斯随机观测矩阵进行降维后(如t=2或4), 随着观测矩阵维数的减小, 其重构图像的峰值信噪比也逐渐降低, 但其总体上仍具备较高的重构质量, 且在主观视觉上很难觉察重构图像质量的下降。从而进一步验证了可以采用本文的方法, 利用降维后的观测矩阵对原始信号进行采样和重构。

实验中, 为进一步验证降维观测矩阵对重构性能的影响, 采用不同的采样率对原始图像进行了采样, 且在采样率相同时, 利用不同维数的高斯随机观测矩阵进行观测。实验同样进行了50次, 取重构图像峰值信噪比的平均值, 其结果如表 3所示。

如表 3所示, 在相同采样率的情形时, 随着高斯随机观测矩阵维数的降低, 经其采样重构后的峰值信噪比呈一定的降低趋势, 但下降趋势非常平缓。特别地, 当t=2或4时, 其重构信号的峰值信噪比与t=1时几乎相等, 即当观测矩阵大小降低至1/4时, 对2维图像信号的重构质量几乎没有影响。当观测矩阵维数继续降低时, 其重构图像的质量稍有下降, 且下降程度较t=2或4时更为明显。由此, 可以清楚地看到本文所提的基于半张量积的压缩感知模型, 有效降低了观测矩阵所占存储空间的。故在实际应用中, 可以根据实际情况, 找到重构质量与观测矩阵存储空间的折中点。

结合半张量积理论, 提出了一种新型的压缩感知模型。虽然所提出的压缩感知模型没有有效提升重构图像的质量, 但有效降低了观测矩阵在采样和重构过程中所占的存储空间。实验仿真结果表明, 在观测矩阵降维的过程中, 重构信号的质量并没有非常明显的下降, 且当选择适当大小的观测矩阵, 其重构图像的质量几乎没有下降。利用所述模型, 完全遵循了现有压缩感知的理论框架, 在采样和重构过程中不增加任何附加操作。

本文借鉴基于l^q(0 < q < 1)范数的迭代重加权最小二乘重构方法对采样信号进行了重构, 在重构过程中, 当原始信号的稀疏度K较大时, 存在重构时间增加, 实时性变差的现象。这一点是后续将要继续着重解决的一个重点。此外, 在后续工作中, 将继续尝试其他的重构方法, 以期在降低观测矩阵维数的情形下, 在提升或不降低重构信号质量的同时, 进一步加快重构速度, 以实现压缩感知的低存储, 快速化重构。