压缩感知

基本概念

关于压缩感知的一些基础概念，可以看 Baraniuk的这篇Lecture Note（想不到都有4000+的引用量），还有Baraniuk在Microsoft Research的一个报告，以及一个基本和报告内容一致的PDF文件

Shannon/Nyquist采样定律

Shannon/Nyquist采样定律 告诉我们，如果要在采样以后能够完整保留保证原始信号中信息，采样频率必须大于信号中最高频率的2倍。一直以来，声音、图像等传感器采样的方式都服从了Shannon/Nyquist采样定律，在保证了信号完整性的同时带来的是大量的数据。比如图像传感器采集的原始图像数据，每张图都有好几十兆。为了减小数据量，才出现了诸如JPEG，JPEG2000等压缩算法。

奈奎斯特采样和压缩采样
图像来源
这样的过程可以用上图的上半部分来表示，x为原始信号，经过Uniform sample 采样得到了一个N维度的向量，然后在经过DCT或者小波变换，得到了M个主要的成分。这里M远远小于N的维度，所以实现了压缩。压缩感知要实现的是上图的下半部分，直接从信号中采样M个成分，然后回复出近似与原始信号x的N向量。如果能够实现这一过程的话，中间的压缩过程就没有了，将会极大的方便数据存储和传输。

稀疏性/可压缩性

很多的信号都是具有稀疏性或者可压缩性。Baraniuk的Lecture Note中对此做了定义。如果存在一个Nx1维度的信号x，那么在N维空间中它必然可以表示：

image.png
其中，ψ_i为该空间中N个正交基向量中的第i个，并且为N维，对应的s_i为对应正交基向量的系数。如果一个信号x能够仅用k个基向量的线性组合来表示，那么这个信号x被称为k-sparse，也就是剩余的N-K个系数都为零。可压缩性和k-sparse的定义稍有不同，如果N个系数中，仅仅只有几个数字非常大，而其他的数字则非常小，那么信号x就具有压缩性。 信号的稀疏性
图像来源
上图中一张图像具有N个像素，虽然在时域这N个像素基本上都不等于零，但是如果对其进行小波变化的话，基本上除了前面的系数以为，其他的系数都很小，对于声音信号等其他信号，也具有这一特征。也就是说，总可以把信号弄到某个域中，让它在这个域具有稀疏性。

• 对于一个稀疏信号所在空间，其存在N个正交基向量，由于其只有K个位置存在系数，而其他位置系数为零，如果把这些点画在N维空间中，这些点是多个K维子空间的集合。比如在三维空间中，有一个维度的系数为零，那么结果就是几个面的组合，如下图的第一部分。
• 对于一个可压缩信号而言，则将构建一个类似于充满’毛刺‘的子空间，如下图第二部分所示。
• 对于具有稀疏性和可压缩性的信号，他们在空间中的所在位置都是非常靠近空间坐标系的，这是他们共同的特点。

稀疏性和可压缩性

压缩感知

假如已知一个信号x具有稀疏性，如下图所示，只有三个点非零（用绿、红、蓝三个点表示），如果使用一个MxN的测量矩阵（图中红色区域）对其进行测量，那么每次的测量就是x和矩阵某一行向量的内积，得到的测量结果就是具有M个维度的数据y。如果M很小的话，那么就实现了dimension reduction。同时，我们希望M跟K非常的接近，这样我们才能够尽量的保留信号中的信息。

压缩采样
图像来源

一个很明显的问题是，如果我们的MxN矩阵不是满秩（full rank）的，也就是说各个行向量不是独立的那么，我们必然会丢失一些信息。举个极端的例子，如果每个行的行向量是一样的，那么得到的y将是M个相同的值。这显然是无法让我们复原原始信号的。但是，如果我们的信号x是k-sparse的,那么其实我们只关心那个k个值，对应到矩阵乘法中，就是我们MxN这个测量矩阵中的某K列。只要这K列的数据是满秩的就OK了。也就是如果我们能够找到这样一个矩阵，其中任意MxK的子矩阵是full rank，我们就可以用这样的观测矩阵进行测量，并且又可能可以恢复数据。非常幸运的是，一个随机的矩阵就有很大的概率具有这样的性质。所以，我们可以使用下面的方法来对数据进行随机采样。其中采样矩阵为MxN维度随机矩阵。采样得到的每一个值都是信号x和随机矩阵一个行向量的内积。

随机采样

信号恢复

通过上面的操作，我们获得了y，并且也知道自己产生的MxN随机矩阵是啥样，剩下的就是重建我们的信号x。

信号恢复

因为可以将x看做在N维空间中，而这个MxN矩阵则可以看做一个在这个N维空间中的N-M维度的超平面，这个超平面穿过了真正的x，当然也穿越了一些别的x，我们要做的就是利用某种标准来找到真正的x。大牛们证明L2作为标准是不行滴，L0虽然是正确的标准，但是很难求解。如果将L1作为标准也能找到正确的解，但是求解就简单很多。

L1信号重建

很多的算法能够支持L1的求解，但是我们不用关心这些算法的处理过程，只需要对其进行使用即可。

处理在不稀疏的信号

前面我们提及的几个例子中，信号x只有几个地方有值，其余位置的值皆为零。但是如果我们的x是一张图像的话，大部分的位置都将不是零，这种情况我们要如何处理呢？图像虽然在时空域不是稀疏的，但是它在DCT、小波等变换后的域是稀疏的，所以我可以在这个域里面处理它。如下图所示，虽然x信号不稀疏，但是它可以表示为一个稀疏信号a和变换矩阵ψ的结果。如果将变换矩阵ψ和原来的矩阵结合为一个新的矩阵，那么就相当于在稀疏信号a上采样了。

处理任何域的稀疏信号

Python实现压缩感知

1D信号

在有了以上的基础知识以后，我们可以来做一些实验了。主要参考资料来源Imaging via Compressive Sampling、Compressive Sensing in Python 。

首先，我们产生两个sin信号组合而成的信号

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import scipy.optimize as spopt
import scipy.fftpack as spfft
import scipy.ndimage as spimg
import cvxpy as cvx


n = 5000
t = np.linspace(0, 0.125, n)
y = np.sin(1394 * np.pi * t) + np.sin(3266 * np.pi * t)
yt = spfft.dct(y, norm='ortho')

plt.subplot(221)
plt.plot(t,y,linewidth=0.35)
plt.title('Original signal')

plt.subplot(222)
plt.plot(yt,linewidth=0.5)
plt.title('Frequency Component')

plt.subplot(223)
plt.plot(t[:1000],y[:1000],linewidth=0.45)
plt.title('Original signal, first 1000 points ')

plt.subplot(224)
plt.plot(yt[:500],linewidth=0.5)
plt.title('Original signal, first 1000 points ')

plt.tight_layout()
plt.show()

上面代码生成的图像如下：

两个正弦信号组合

上面代码中，出现的下面这一行的目的是对一维数据进行DCT变换：

yt = spfft.dct(y, norm='ortho')

可以看到，在DCT域中，我们的信号其实是非常的稀疏的。根据压缩感知理论，对这样一个稀疏的信号，其实我们采用随机采样的方式就可以恢复出原始的信号。但是，很显然的一点是，我们无法直接在DCT域进行采样，只能在时空域对信号进行采样。所以， 我们对时空域的5000个信号点进行随机采样：

m = 500 # 10% sample
ri = np.random.choice(n, m, replace=False) # random sample of indices
ri.sort() # sorting not strictly necessary, but convenient for plotting
t2 = t[ri]
y2 = y[ri]
plt.plot(t,y,linewidth=0.35)
plt.scatter(t2,y2,marker='o',color='r',s=5)
plt.show()

从5000个数据点中随机采样500个点
从上图中可以看到，500个点随机分布，杂乱无章，如何从中提取出原始的信号呢？从处理在不稀疏的信号这一节可以看出，我们可以将这5000个时空域的数据看做经过反离散余弦变换（iDCT）变换得到的，而"真正的原始信号"则是‘余弦域’的。现在我们要找到一个变换矩阵，这个矩阵将余弦域信号转化为时空域。对一个对角矩阵进行反余弦变换，就能得到这样一个矩阵，由于我们对信号进行了随机采样，所以将随机采样的位置也考虑其中。最后得到的A矩阵就是我们整体的采样矩阵，包含了变换和随机采样的信息。

A = spfft.idct(np.identity(n), norm='ortho', axis=0)
A = A[ri]

上面的第一行，首先生成了这样一个矩阵，然后第二行按照随机采样的方式取出相应的行。

到这里，我们做好了所有的准备工作：

• 得到了采样后的数据点
• 得到了采样矩阵

剩下的就是在L1条件下，求解出我们的原始信号。然后在对得到的信号进行反余弦变换，就得到了时空域的信号了。

vx = cvx.Variable(n)
objective = cvx.Minimize(cvx.norm(vx, 1))
constraints = [A*vx == y2]
prob = cvx.Problem(objective, constraints)
result = prob.solve(verbose=True)

plt.plot(sig,linewidth=0.35)
plt.show()

重建后的信号

2D信号

在重建了一个1D信号以后，我们可以对2D信号进行压缩采样并重构。很显然，一般的图像信号在时空域也是非稀疏的，我们仍然需要设计一个转换矩阵，将其转化到一个稀疏域下，这里，我们继续使用DCT变换。根据2D-Signal Transforms and The Separability Property中的第7页的内容，可以看到，对于一个可分的图像变换，可以将其看做两个一维变换分别作用于图像矩阵的行和列。这一点类似高斯滤波的可分性。而行变换和列变换的克罗内克积则等效于这个2D变换。也就是说，如果将图像展开为一维数据以后，乘以这个克罗内克积得到的矩阵，就相当于对图像做了2D iDCT变换：

将2D DCT分离

下面我们来进行对图像随机采样，并且复原的工作，首先，我们读取lena图，然后随机取出里面50%的像素

Xorig = spimg.imread('lena_gray.gif',flatten=True, mode='L')
X = spimg.zoom(Xorig,0.25)
ny,nx=X.shape

k = round(nx*ny*0.5)
ri = np.random.choice(nx*ny,np.int(k), replace=False)


plt.subplot(121)
plt.imshow(Xorig,cmap='gray')
plt.axis('off')
plt.title('Original')

plt.subplot(122)
X1d = X.flatten()
b = X1d[ri]
Xsam = np.zeros(np.shape(X.flatten()))
Xsam[ri] = b;
Xsam = np.reshape(Xsam,[ny,nx])
plt.imshow(Xsam,cmap='gray')
plt.title('50% Sample')
plt.axis('off')
plt.show()

Lena原图和随机采样后的图

从上面的结果可以看到，通过对Lena图像进行50%随机的采样，最后得到的图像已经模糊不清了，但是我们人眼仍然是可以看出内容是Lena图像。

A = np.kron(
    spfft.idct(np.identity(nx), norm='ortho', axis=0),
    spfft.idct(np.identity(ny), norm='ortho', axis=0)
    )