CS,全称一般被认为是 Compressive/Compressed Sensing/Sampling,中文叫做压缩感知,于 2004 至 2006 年之间由 David Donoho、Emamnuel Candes、Terry Tao 等人提出来之后,迅速发展壮大,虽然里组成 CS 理论的基本元素往前追溯都都有很多前人的研究,但是也是在把整个东西系统地提出来了之后,才引发了后续的广泛关注,从纯理论到纯应用甚至是新硬件的构造等等各方面。而 CS 本身也是认为是信号处理领域自 Shannon/Nyquist 采样理论以来的重大突破。

因为 Shannon/Nyquist 采样理论使得我们可以用离散信号去 perfectly 重构原来的(如果是 band-limited 并且采样率足够高的话)连续信号,而离散时间信号处理借助于计算机的力量使得人们可以实现现在的各种复杂的信号处理和滤波,如果我们还在用 Analog 设备直接对 Analog 信号进行各种处理的话,这些东西估计都难以想象吧。然而有些问题里的信号并不是完全 band-limited 的,或者要求的 Nyquist 频率过高使得硬件限制等各种原因无法完成有效地采样,然后压缩感知的理论跳出来说,如果信号在某个已知的基下面是稀疏的的话,我们可以使用远低于 Nyquist 频率的代价来完成采样并完美地恢复出原始信号,与此同时,在实际应用中的各种信号,由于具有其各自的结构性,所以通常在合适的 basic 或者 dictionary 下面都会表现出稀疏性(例如自然景观图像在小波域下是稀疏的),另外,Machine Learning 还专门有一个 topic 叫做 sparse learning,其重要的一个目标就是给定一个数据集,构造一个 dictionary 使得该数据集在这个 dictionary 下具有稀疏的线性表达。所以压缩感知的提出被认为是又一次重大革命。

当然把 CS 和 Shannon/Nyquist 采样理论直接并列起来比较的话,还是有一些设定不一样的地方,比如经典的采样理论一般讨论的需要采样的信号是连续时间的,无穷维的信号,而 CS 中虽然也有扩展理论尝试处理无穷维的信号,但是主要的着眼点还是在于长度为 NN 的有限维信号;此外,经典的采样理论中所谓的“采样”通常是指在某一固定的频率下沿着时间轴进行……“采样”,然而 CS 中的采样则是指去测量一个线性函数作用在原始信号上的值,换句话说是计算一个内积。

具体来说,CS 着眼于一个 NN 维的向量 xx,采样的过程是使用一个线性函数,也就是另一个 NN 维向量 aia_i 作用上去得到 yi=ai,xy_i=\langle a_i, x\rangle。将 nn 个采样线性函数按行排列成一个 n×Nn\times N 的矩阵 AA,于是问题就变成了已知 AAy=Axy=Ax,求 xx

线性代数的知识告诉我们,如果 n<Nn<N 的话,AA 的秩必定小于 NN,此时该线性系统有无穷多个解,如果 x0x_0 是其中任意一个解的话,那么所有

x=x0+xˉ,xˉN x' = x_0 + \bar{x},\quad \bar{x}\in\mathcal{N}

都是该问题的解,其中

N={xˉ:Axˉ=0} \mathcal{N}=\{\bar{x}: A\bar{x} = 0\}

图 1三维空间下的 2-sparse 子空间示意图。图片来自于《Compressed Sensing: Theory and Applications》。

是矩阵 AA 的 Null Space。这代表我们采样的数量还不够 unique 地决定原始的信号。但是如果我们已知 xx 是稀疏的,情况就不一样了。从上面的解的形式也让我们认识到 AA 的 Null Space 的结构将会是影响该问题的重要因素。不过,再进一步之前,我们需要先定义一些符号。首先,一个向量 xx,如果其中非零元素个数不超过 kk 的话,也就是说,如果 x0k\|x\|_0\leq k,那么我们称其为 kk-sparse 的,所有 kk-sparse 的向量构成一个集合,记为 Σk\Sigma_k,显然 ΣkΣk+1\Sigma_k\subset\Sigma_{k+1}

直观来说,Σk\Sigma_k 是由 (Nk)\binom{N}{k} 个线性子空间的并构成的一个空间,它本身并不是一个线性子空间,换句话说,两个 kk-sparse 的和通常并不再是 kk-sparse 的。R3\mathbb{R}^3 下的 Σ2\Sigma_2 的示意图如图 (img: 1) 所示。

如果我们已知了 xxkk-sparse 的,那么的可行解的范围就缩小了,如果刚才提到的 Null Space 性质良好,就能够保证我们有唯一解。具体来说,考虑一下如果 NΣ2k={0}\mathcal{N}\cap \Sigma_{2k}=\{0\} 的话,会怎样?假设 x1x_1x2x_2 同时是原问题的解,那么我们有

A(x1x2)=Ax1Ax2=yy=0 A(x_1-x_2) = Ax_1 - Ax_2 = y - y = 0

换句话说,x1x2Nx_1-x_2\in \mathcal{N},反过来,由于 x1x_1x2x_2 都是原问题的解,它们必定各自是 kk-sparse 的,所以 x1x2Σ2kx_1-x_2\in\Sigma_{2k},但由于 Σ2k\Sigma_{2k}N\mathcal{N} 的交集里只有零元素,因此 x1x2=0x_1-x_2=0,也就是说,原问题在这种情况下的解是唯一的。

为了更方便地描述这种情况,人们模仿矩阵的 rank 构造出一个叫做 spark 的量,一个矩阵 AA 的 spark 是最小的数 kk 使得 AA 存在 kk 列是线性相关的。对比一下,rank 则是最小的 kk 使得所有 k+1k+1 列都是线性相关的。根据定义容易看到,如果 spark(A)>2k\text{spark}(A)>2k 的话,那么 AA 的 Null Space 和 Σ2k\Sigma_{2k} 的交集必然只有零元素。

也就是说,如果我们已知原始的信号是 kk-sparse 的,那么用一个 spark 大于 2k2k 的感知矩阵 AA 来进行感知,就能保证所对应的原始信号是唯一的(并且这是充要条件)。不过光靠上面的分析只是证明了解唯一地存在,并没有告诉我们要如何去求得这个解,如果暴力枚举的话,将会变成一个 NP Hard 问题。所以在这里有必要简单地小结一下,CS 的基本理论中所研究的问题大概分为以下几块

  1. AA 需要满足什么样的性质可以保证压缩感知问题的解是唯一的;又有什么样的性质可以保证这个唯一的解是可以有效地(例如,多项式时间算法)解出来的。这里除了刚才提到的 spark 以外,还有许多诸如 Null Space Property (NSP), Restricted Isometry Property (RIP), Coherence 之类的。
  2. 构造具体的算法去进行 Decoding,这里通常分为三类:1\ell_1-norm 优化算法,将原始问题进行 convex relaxation;贪心算法,迭代求解;组合优化算法,通常属于理论计算机科学所研究的范围。
  3. 如何去验证一个给定的 AA 是否满足上面提到的要求,或者给定上面的要求如何去构造一个符合要求的 AA。前面提到的各种性质除了 Coherence 之外基本上都无法暴力直接验证,因为需要进行子集枚举;而满足性质的矩阵构造方面,现在已经普遍接受的是使用随机矩阵进行构造,并证明构造出来的矩阵以压倒性的概率 (Overwhelming Probability ;p) 满足给定的性质。

于是再回到刚才关于 Null Space 的讨论中,虽然 Spark 的条件给出了在 xΣkx\in\Sigma_k 时候解唯一的充要条件,但是仅考虑 kk-sparse 信号有时候还是过于局限,因为实际问题中有很多信号本身并不是严格的 kk-sparse,而仅仅是近似稀疏的,也就是是说,它们可以通过一个稀疏的向量来进行近似。具体地,我们可以定义如下的 kk-term approximation error

σk(x)X=minzΣkxzX \sigma_k(x)_X = \min_{z\in\Sigma_k}\|x-z\|_X

当然,如果 xx 本身就是 kk-sparse 的,那么 approximation error 就为零。此外,当 X\|\cdot\|_Xp\ell_p-norm 的时候,xx 的最佳 kk-term 近似其实就是保留绝对值最大的 kk 个分量,将剩余的分量全部置零。图 (fig: 1) 以最近的一幅合作水彩作为例子,可以看到在 Wavelet 系数域里,只要保留 5% 那么多非零系数已经可以达到相当令人眼满意的近似结果

(a)
(b)
(c)
图 1(a) Original Image. (b) Approximation with 5-percent wavelet coefficients. (c) Approximation with 1-percent wavelet coefficients.

回到 CS 的问题,在 xx 本身并不是 kk-sparse 的时候,我们一般不指望能够完美地恢复出 xx 来,但是通常可以希望做到和 best kk-term approximation 差不多好,具体来说,我们希望能够实现

xΔ(Ax)XC0σk(x)X(1) \|x-\Delta(A x)\|_X \leq C_0 \sigma_k(x)_X \label{25c95221d2ba2dfa343638b1b0b578404c184ee6}\tag{1}

其中 Δ\Delta 表示我们的 decoding/reconstruction 算法。直观来讲,等式的右边,就类似于图 (fig: 1) 中的近似结果,这是我们在看到完整的原图之后,把所有的系数从大到小排序然后只保留前 kk 个的结果;而左边则是只观察到了矩阵 AA 所对应的 nn 个 sampling 值之后进行重构的结果,我们希望的是在这样的情况下得到的结果和先观察到完整图之后再做近似的结果“差不多”。另外可以看到这个结果是把之前的情况包含进来作为特殊情况的:如果 xx 本身就是 kk-sparse 的,那么显然 σk(x)X=0\sigma_k(x)_X=0,所以我们可以保证无损地恢复出原来的信号来。

和之前一样,为了能够达到 (eq: 1) 中的目标,我们从 AA 的 Null Space 入手。之前的分析中我们要求 AA 的 Null Space 中不要有除了 0 以外的稀疏向量,现在我们考虑的对象变成了近似稀疏的向量,于是我们类似地要求 Null Space 中不要存在 0 以外的近似稀疏的向量。具体来说,我们将要求所有 hNh\in\mathcal{N} 满足

hXC0σ2k(h)X(2) \|h\|_X \leq C_0 \sigma_{2k}(h)_X \label{bd57b00b7d41c31dc22bec44c1fdc0a158e7a461}\tag{2}

这个不等式直观上来说,就是在说 AA 的 Null Space N\mathcal{N} 里的向量 hh 不应该将值“稀疏地”集中在某 2k2k 项以内。比如说,hh 如果是 2k2k-sparse 的话,那么式 (eq: 2) 右边将会等于零,于是左边也必须等于零,所以和之前一样,严格 2k2k-sparse 的向量只有零向量;而近似 2k2k-sparse 的向量也无法存在(注意这里的常数 C0C_0 是和 (eq: 1) 中对应起来的同一个常数)。

为了证明这一点,我们注意到对于任意的 hNh\in\mathcal{N},我们可以将它分解为三个部分: h=h1+h2+h3h=h_1+h_2+h_3,其中 h1h_1 由绝对值最大的 kk 项组成,h2h_2 由剩下的绝对值最大的 kk 项组成,而 h3=hh1h2h_3=h-h_1-h_2。首先,由于 h1Σk-h_1\in\Sigma_k,由 (eq: 1) 我们可以保证无损重建,也就是说 h1=Δ(A(h1))-h_1 = \Delta(A(-h_1))。另一方面,由于 hNh\in\mathcal{N},我们有

0=Ah=A(h1+h2+h3)A(h2+h3)=A(h1) 0 = Ah = A(h_1 + h_2 + h_3)\quad\Rightarrow\quad A(h_2+h_3) = A(-h_1)

因此同样地,我们有 h1=Δ(A(h2+h3))-h_1 = \Delta(A(h_2+h_3)),由此,注意到

hX=h2+h3Δ(A(h2+h3))XC0σk(h2+h3)C0h3X=C0σ2k(h) \begin{aligned} \|h\|_X &= \|h_2+h_3 - \Delta(A(h_2+h_3))\|_X \\ &\leq C_0\sigma_k(h_2+h_3) \\ &\leq C_0\|h_3\|_X \\ &= C_0\sigma_{2k}(h) \end{aligned}

即证。也就是说,为了实现 (eq: 1),必要条件是 AA 的 Null Space 里的向量满足 (eq: 2),该性质又被称为 Null Space Property (NSP)。实际上,通过改变一下常数,我们可以证明该条件同时是充分条件。具体来说,如果

hXC02σ2k(h)X,hN(3) \|h\|_X\leq \frac{C_0}{2}\sigma_{2k}(h)_X, \quad \forall h\in\mathcal{N} \label{6b34eab95966ec74518ca926454607979fdd229f}\tag{3}

那么我们可以令 decoder 为

Δ(y)=arg minAz=yσk(z)X \Delta(y) = \operatorname*{arg\,min}_{Az = y}\sigma_k(z)_X

则,由 decoder 的定义知道,xΔ(Ax)x-\Delta(Ax) 是属于 Null Space 的,于是根据 (eq: 3),我们有

xΔ(Ax)XC02σ2k(xΔ(Ax))XC02(σk(x)X+σk(Δ(Ax))X)C0σk(x)X \begin{aligned} \|x - \Delta(Ax)\|_X &\leq \frac{C_0}{2}\sigma_{2k}(x-\Delta(Ax))_X \\ &\leq \frac{C_0}{2}\left( \sigma_k(x)_X + \sigma_k(\Delta(Ax))_X \right) \\ &\leq C_0\sigma_k(x)_X \end{aligned}

即证。其中最后一个不等式是由于我们所定义的 decoder 是对其参数的 σk()X\sigma_k(\cdot)_X 进行最小化的缘故。当然和其他具体的 1\ell_1 最小化等 decoder 不一样,这个 decoder 也并不确定是否有有效地算法可以去进行求解的样子。

需要注意的是,我们上面的结论其实并没有明确地要求 xx 是怎样地“近似 kk-sparse”,实际上 xx 可以完全不 sparse,上面的结论仍然不会受到影响,但是结论本身可能就没有什么用处了,因为 (eq: 1) 右边本身就很大的话,这个 bound 就没有任何意义了。不过接下来我们还要再将我们的目标扩充一下:将测量误差考虑进来。换句话说,现在我们的 sample 结果将是

y=Ax+e y = Ax + e

其中 ee 代表测量误差。为了讨论简单起见,我们暂时回到 kk-sparse 的信号,此时我们希望我们的 sensing + decoding 过程是 stable 的,具体来说,我们希望对于 xΣkx\in\Sigma_k,有

Δ(Ax+e)x2Ce2 \|\Delta(Ax+e) - x\|_2 \leq C \|e\|_2

为了达到这个 stable 的要求,我们必须要有,对任意的 xΣ2kx\in\Sigma_{2k}

1Cx2Ax2(4) \frac{1}{C}\|x\|_2 \leq \|Ax\|_2 \label{c15050c56f4fcae5497334d62a94f987f382659e}\tag{4}

为了证明这一必要条件,我们将 xx 分解为 x=x1x2x=x_1-x_2,且 x1,x2Σkx_1,x_2\in\Sigma_k,并定义

e1=A(x2x1)2,e2=A(x1x2)2 e_1 = \frac{A(x_2-x_1)}{2},\quad e_2 = \frac{A(x_1-x_2)}{2}

Ax1+e1=Ax2+e2=A(x1+x2)2 Ax_1 + e_1 = Ax_2 + e_2 = \frac{A(x_1+x_2)}{2}

x^=Δ(Ax1+e1)=Δ(Ax2+e2)\hat{x}=\Delta(Ax_1+e_1)=\Delta(Ax_2+e_2),则

x2=x1x22=x1x^+x^x22x1x^2+x^x22Ce12+Ce22=CA(x1x2)2=CAx2 \begin{aligned} \|x\|_2 &= \|x_1-x_2\|_2 \\ &= \|x_1-\hat{x} + \hat{x}-x_2\|_2 \\ &\leq \|x_1-\hat{x}\|_2 + \|\hat{x}-x_2\|_2 \\ &\leq C\|e_1\|_2 + C\|e_2\|_2 \\ &= C\|A(x_1-x_2)\|_2 \\ &= C\|Ax\|_2 \end{aligned}

但是如果仅仅是 (eq: 4) 的话,我们可以仅仅通过对 AA 进行放大而达到任意想要的 CC stability。当然如果真的能够放大 sensing 矩阵而不同时增大测量误差的话,这确实是有效地消除测量误差所带来的影响的有效途径,但是实际中通常对 sensing 进行这种 naive 的放大之后相应的误差也会跟着放大,所以为了回避这个 trivial 的情况,我们再对 (eq: 4) 的右边也进行一下限制。于是有了下面这个性质。

1

定义 1(Restricted Isometry Property (RIP)). 对于矩阵 AA,如果存在常数 δk(0,1)\delta_k\in(0,1),使得对任意 xΣkx\in\Sigma_k,都有 (1δk)x22Ax22(1+δk)x22 (1-\delta_k)\|x\|_2^2 \leq \|Ax\|_2^2 \leq (1+\delta_k)\|x\|_2^2 那么我们称 AA 满足 kk 阶 RIP。

这是一个比之前更强的性质,由 δ(0,1)\delta \in (0,1) 很容易知道,如果 AA 满足 kk 阶 RIP,那么显然 AA 的 spark 是大于 kk 的,否则就存在非零 xΣkx\in\Sigma_k 使得

0=022=Ax22(1δk)x22>0 0 = \|0\|_2^2 = \|Ax\|_2^2 \geq (1-\delta_k)\|x\|_2^2 > 0

除此之外,RIP 也比 NSP 要强。具体来说,我们有如下的定理。简单起见,在接下来的讨论中,我们将 NSP 限制为 X\|\cdot\|_X1\ell_1-norm 的情况。

1

定理 1. 如果 AA 满足 2k2k 阶 RIP,并且 δ2k<1/3\delta_{2k}<1/3,那么 AA 也满足 2k2k 阶的 NSP,并且对应的常数为

C0=1δ2k13δ2k C_0 = \frac{1-\delta_{2k}}{1-3\delta_{2k}}

T0{1,,N}T_0\subset\{1,\ldots,N\}hh 中绝对值最大的 kk 项的下标集,T1T_1 为 除去 T0T_0 之后的绝对值最大的 kk 项的下标集,依此类推。记 T=T0T1T=T_0\cup T_1,显然,我们有

σ2k(h)1=hTc1 \sigma_{2k}(h)_1 = \|h_{T^c}\|_1

接下来我们先证明对于 hNh\in \mathcal{N},有

hT1C~hTc1 \|h_T\|_1 \leq \tilde{C}\|h_{T^c}\|_1

于是

h1=hT1+hTc1(1+C~)hTc1C0hTc1(5) \|h\|_1 = \|h_T\|_1 + \|h_{T^c}\|_1 \leq (1+\tilde{C})\|h_{T^c}\|_1\triangleq C_0\|h_{T^c}\|_1 \label{6e22ac41d7058bfa493e8f8fe3127633d4290918}\tag{5}

在证明中我们还需有用到如下的引理。

1

引理 1. 若 AA 满足 2k2k 阶 RIP,令 hhAA 的 Null Space 中的一个向量,集合 T0,T1,T2,T_0,T_1,T_2,\ldots 定义和刚才一样,并且 T=T0T1T=T_0\cup T_1,则

hT2αhT0c1k \|h_T\|_2 \leq \alpha \frac{\|h_{T_0^c}\|_1}{\sqrt{k}}

其中

α=2δ2k1δ2k \alpha = \frac{\sqrt{2}\delta_{2k}}{1-\delta_{2k}}

引理(的更 general 的情况,不要求 hh 属于 Null Space 时)的证明可以参见《Compressive Sensing: Theory and Applications》第一章中的引理 1.3。继续我们定理的证明,根据引理,我们有

hT12khT22αhT0c1=2α(hTc1+hT11)2α(hTc1+hT1) \begin{aligned} \|h_T\|_1 &\leq \sqrt{2k}\|h_T\|_2 \\ &\leq \sqrt{2}\alpha\|h_{T_0^c}\|_1 \\ &= \sqrt{2}\alpha \left( \|h_{T^c}\|_1 + \|h_{T_1}\|_1 \right) \\ &\leq \sqrt{2}\alpha \left( \|h_{T^c}\|_1 + \|h_T\|_1 \right) \end{aligned}

整理得

(12α)hT12αhTc1 (1-\sqrt{2}\alpha)\|h_T\|_1 \leq \sqrt{2}\alpha \|h_{T^c}\|_1

δ2k<1/3\delta_{2k}<1/3 时,我们有 12α>01-\sqrt{2}\alpha > 0,于是可以将系数除到右边而不改变不等号方向,从而得到:

hT12α12αhTc1C~hTc1 \|h_T\|_1 \leq \frac{\sqrt{2}\alpha}{1-\sqrt{2}\alpha}\|h_{T^c}\|_1 \triangleq \tilde{C} \|h_{T^c}\|_1

再带入相应的项即证。不过,RIP 既然是更强的条件,它自然也有自己的长处。我们刚才证明了为了能够让压缩感知在有感知误差的时候也表现的 stable,RIP 是必要条件。实际上,RIP 同时也是充分条件。

定理. 如果 AA 满足 2k2k 阶 RIP,并且 δ2k<21\delta_{2k}<\sqrt{2}-1,令 x^\hat{x} 是如下凸优化问题的最优解:

minx1,s.t.Axyϵ \min \|x'\|_1,\quad s.t. \|Ax'-y\|\leq \epsilon

其中 y=Ax+ey = Ax + e,而 ϵe2\epsilon \geq \|e\|_2 是感知误差的一个上界估计。则

x^x2C0σk(x)1k+C2ϵ \|\hat{x}-x\|_2 \leq C_0\frac{\sigma_k(x)_1}{\sqrt{k}} + C_2 \epsilon

其中

C0=21(12)δ2k1(1+2)δ2k,C2=41+δ2k1(1+2)δ2k C_0 = 2\frac{1-(1-\sqrt{2})\delta_{2k}}{1-(1+\sqrt{2})\delta_{2k}}, \quad C_2 = 4\frac{\sqrt{1+\delta_{2k}}}{1-(1+\sqrt{2})\delta_{2k}}

证明可以参见 (Candes, 2008) 或者《Compressive Sensing: Theory and Applications》第一章中整理过的定理 1.9 的证明。关于这个定理,有几点需要注意的:

小结一下,如果我们保证感知矩阵 AA 满足 2k2k 阶的 RIP,那么就能通过求解 1\ell_1-norm 优化问题来进行 decoding。不过,在满足 RIP 的情况下,除了 1\ell_1-norm 优化之外还有没有其他行之有效的 CS decoding 算法呢?要满足特定的 RIP 条件的时候,对于感知样本的数目 nn 有什么样的要求呢?具体的 AA 应该如何构造呢?虽然在本文开始的时候已经有一些剧透了,不过由于长度限制,具体的内容还是未完待续吧!