在这篇文章中，我将研究Python中的压缩感知(也称为压缩感知、压缩采样和稀疏采样)。由于压缩感知的概念可以广泛应用于各种学科，我将主要集中在如何将其在一维和二维上应用到声音和图像等方面。具体地说，我将展示如何获取高度不完整的信号样本数据集，并重建潜在的声音或图像。这是一项非常强大的技术。
正如你可能知道的，有许多不同类型的规范。也许最常见、最广为人知的是$L^2$规范：
$L^2$范数很好，因为它易于计算、易于区分，并且具有直观的吸引力(例如，向量的范数就是它的长度)。很多非常重要的算法和方法都依赖于$L^2$，包括最小二乘拟合。
也就是说，$L^2$规范并不是万能的Goto解决方案。其他规范也有许多有趣和有用的性质。考虑$L^1$规范：
它不是对每个元素平方，而是简单地取其绝对值。虽然绝对值是恼人的，因为它经常在其导数中引入不连续性，但与发生在$L^2$范数中的平方相比，它确实有一些独特的性质。压缩传感就是利用这些特性。
让我们用Python可视化一些数据，看看我在说什么。
#确保您已经安装了以下软件包import numpy as np import matplotlib as MPL import matplotlib.pylot as PLT import scipy.Optimize as spopt import scipy.fftpack as spfft import scipy.ndimage as spimg import cvxpy as CVX。
首先，我们要做的是创建一些任意的线性数据，包括一些噪声。让我们使用虚构的方程式：
#生成一些噪声x=np的数据。排序(NP。随机的。均匀(0，10，15)y=3+0.2*x+0.1*Np。随机的。Randn(len(X))。
现在，让我们为数据样本拟合两行。对于第一行，我们将使用$L^1$范数作为良好拟合的标准；对于第二行，我们将使用$L^2$范数。
#查找L1直线拟合L1_FIT=λx0，x，y：np。总和(np.。ABS(x0[0]*x+x0[1]-y))xopt1=spopt.。Fmin(func=L1_Fit，X0=[1，1]，args=(x，y))#查找L2直线拟合L2_Fit=λX0，x，y：np。总和(np.。幂(x0[0]*x+x0[1]-y，2)xopt2=spopt。Fmin(func=L2_Fit，x0=[1，1]，args=(x，y))。
请注意，这两个拟合似乎都很好地拟合了数据。当然，它们并不完全一致，但考虑到噪音，它们都是合理的近似值。
现在，让我们变得有点疯狂，并添加一些异常值。换句话说，让我们扰乱几个点，使它们远离线条。如果你仔细想想，这实际上并不是那么不寻常。在现实世界的数据中，异常值经常出现，引起各种令人头疼的问题。
#通过添加离线y2=y来调整数据。Copy()y2[3]+=4 y2[13]-=3#重新装配线xopt12=spopt。Fmin(func=L1_Fit，x0=[1，1]，args=(x，y2))xopt22=spopt。Fmin(func=L2_Fit，x0=[1，1]，args=(x，y2))。
当我们重新绘制$L^1$和$L^2$拟合时，我们会看到一些有趣的事情：$L^1$拟合仍然符合数据中的总体趋势，而$L^2$拟合似乎被异常值“破坏”了。这一切为什么要发生？归根结底是$L^2$ERROR平方，而$L^1$ERROR不平方。当您使用$L^2$误差解释将直线拟合到数据时，离群值的位移会产生不成比例的影响，因为它们已经很大的误差会被平方。只要看看我们示例中两个异常值的距离，想象一下将它们平方-当然，$L^2$线倾斜并不奇怪！
但是，当使用$L^1$解释错误时，离群值的贡献不会超过它们的位移。结果是一个更干净的合身，更符合我们对好合身应该是什么样子的直觉。正是这种有趣的特性为压缩传感打开了大门。
在此示例中(借用自Kutz 1)，我们将创建一个人造声波，对其进行10%的采样，并从10%的采样重建原始信号。这是一维压缩传感。
#两个正弦波之和n=5000t=Np。Linspace(0，1/8，n)y=Np。SIN(1394*np.)。Pi*t)+Np。SIN(3266*np.)。Pi*t)yt=spfft。DCT(y，Norm='；Ortho'；)。
在上面的曲线图中，我们看到信号有一个清晰的模式，但不是平凡的。顶行的曲线图是不同尺度的时域信号的曲线图。底行中的曲线图是谱域中的信号(即信号的频率内容)。特别考虑到频域，我们注意到，除了代表两个正弦频率的两个尖峰之外，频谱几乎为零。
现在想象一下对时间信号的10%进行采样(见下文)。你会有一个数据集，在肉眼看来，这是无稽之谈。基本信号IS仍然是相同的，其频率内容也是相同的(除了两个尖峰之外，大部分是零)。有人可能会问，是否有可能从不完整的数据中提取这两个主要频率，以便我们可以重建信号？答案是肯定的！
#提取信号m=500的小样本#10%样本ri=np。随机的。CHOICE(n，m，REPLACE=FALSE)#索引ri的随机样本。Sorte()#排序不是严格必需的，但便于绘制t2=t[ri]y2=y[ri]。
这种情况下的压缩感知是由于信号的频率内容高度稀疏这一事实而成为可能的。这就是$L^1$规范发挥作用的地方。我们要做的是，在所有可能的信号中，找到与已知数据匹配的最简单的信号。换句话说，我们希望使用最小化例程来找到一组满足两个条件的频率：(A)底层信号与我们数据的信号完全匹配(或尽可能接近)；以及(B)频率的$L^1$范数最小化。这样的例程将产生稀疏的解决方案-这正是我们想要的。
在Python中，有几种方法可以实现这一点。也许最简单的方法是利用凸优化库CVXPY。使用下面的代码最小化信号频率的范数，并约束候选信号应与我们的不完整样本完全匹配。
#创建idct矩阵运算符A=spfft。Idct(np.)。恒等式(N)，范数=#39；正交，轴=0)A=A[ri]#DO L1优化VX=CVX。变量(N)目标=CVX。最小化(CVX.。Norm(Vx，1))约束=[A*VX==y2]Prob=CVX。问题(目标，约束)结果=问题。求解(详细=True)
您可能会问：那个该死的$A$矩阵是什么？嗯，这是整个派对的关键。听我解释。
为了实现最小化，我们必须以某种方式将我们的问题归结为线性方程组：
具体地说，我们希望导出一个矩阵$A$，该矩阵可以与候选解$x$相乘，以产生包含数据样本的向量$b$。在我们当前问题的上下文中，候选解$x$存在于频域，而已知数据$b$存在于时域。显然，矩阵$A$执行从谱域到时域的采样和变换。
压缩传感实际上归结为能够正确地导出$A$运算符。幸运的是，有一种方法论。首先，让$f$作为矢量形式的目标信号(如果您的信号是二维或更高的，则将其展平)，并将$\phi$作为采样矩阵。然后：
现在，设$\psi$是将信号从谱域转换到时域的矩阵。给定频域中的解决方案$x$，它如下所示：
因此，$A$简单地由从域转换矩阵$\psi$采样的行组成。$\psi$矩阵很容易构造-它是作用于单位矩阵列的离散余弦逆变换。矩阵乘积$\psi x$等同于做idct(X)。
现在我们已经构造了$A$矩阵并运行最小化，我们可以通过将解从频域转换回时间域来重构信号。左下方是原始信号及其频率内容。右边是我们的$L^1$近似值。我想说，这对于只使用10%的数据来说是相当不错的！
一个突出的问题是，重建质量在$t=0$及其附近显著降低。这可能是因为我们的采样间隔违反了余弦变换的周期性边界条件要求。当然，给定一个事先不知道其性质的任意信号样本，很难不违反周期性边界条件。好消息是，现在我们对真实信号的频率有了一些非常明确的指示。如果需要，我们可以返回并在满足周期边界的间隔内重新采样信号。
现在让我们用我们从一维案例中学到的东西来做二维压缩传感。这才是真正有趣的开始，因为我们现在可以尝试重建图像了。
下面，我们将使用与前面完全相同的方法随机采样和重建M.C.Escher(约.。1200x1600像素)。由于$A$矩阵施加的内存限制，我们将首先考虑图像的缩小版本(约。50x65像素)。在接下来的部分中，我们将扩展该例程以处理大图像。
请注意，本网站不提供DCT或IDCT的2D版本。但是，通过认识到2D离散余弦变换只不过是作用于$x$行的DCT，然后是作用于其列的第二个DCT操作(反之亦然)，就可以很容易地构建它们：
出于个人喜好，我喜欢告诉SciPy的dct和idct方法对矩阵的列进行操作(而不是对行进行操作的默认行为)。首先，这使Python代码与MATLAB的代码保持一致。其次，它使构建矩阵运算符更直观(至少对我来说是这样)。例如，如果我们设$Y$是$m$×$n$矩阵，其中$i_m$和$i_n$分别是大小为$m$和$n$的单位矩阵，则。
这两个版本都可以工作，但我觉得第一个更干净，因为它自然地将矩阵运算符放在操作数前面。每当我提到dct或idct时，假设我的意思是轴=0的变化。
Def dct2(X)：返回spfft。DCT(spfft.。DCT(x.。T，Norm='；正交，轴=0)。T，Norm='；正交，轴=0)def idct2(X)：返回spfft。Idct(spfft.)。Idct(x.。T，Norm='；正交，轴=0)。T，Norm='；Ortho'；，Axis=0)#读取原始图像并缩小大小，速度Xorig=spimg。Imread('；esscher_water fall.jpeg'；，flatten=True，mode='；L'；)#以灰度X=spimg读取。Zoom(Xorig，0.04)NY，NX=X。形状。
与上一节一样，我们将随机抽取图像索引样本，形成$b$矩阵。然后，我们将生成我们的$A$矩阵。
与1D情况相比，为2D图像数据创建$A$矩阵需要更多的巧妙之处。在下面的推导中，我们将使用Kronecker乘积$\oTimes$和2D离散余弦变换是可分离的这一事实来产生我们的运算符$A$。
设$X$是谱域中的图像，$D_i=\textit{idct}(I_I)$，其中$i_i$是大小为$i$的单位矩阵。然后：
如果$\textit{vec}(X)$是将$X$列堆叠在一起的向量运算符，则：
显然，Kronecker产品是我们所需的转换矩阵$\psi$。因此，我们的矩阵$A$变成$A=\φ(D_n\o乘D_m)$，其中$\φ$是抽样矩阵。您可以用numpy.kron计算Numpy中的Kronecker乘积。这种方法的主要问题是，Kronecker产品可能很快就会变得非常庞大。如果您的目标图像是$m$x$n$，并且您要采集$k$样本，则$A$矩阵的大小为$(MNK)^2$。这就是说，对于小图像来说，这是可以接受的。
#提取信号k=round(nx*ny*0.5)的小样本#50%样本ri=np。随机的。CHOICE(NX*NY，k，REPLACE=FALSE)#指数b=X的随机样本。托特。平坦[ri]b=Np。EXPAND_DIMS(b，AXIS=1)#使用KRON(大NY*NX的内存错误)A=Np创建DCT矩阵运算符。克朗(Spfft.。Idct(np.)。标识(NX)，规范=正交，轴=0)，spfft。Idct(np.)。恒等式(NY)，范数=39；正交，轴=0))A=A[ri，：]#与φ乘以Kron#做L1优化Vx=cvx相同。变量(NX*NY)目标=CVX。最小化(CVX.。范数(vx，1))约束=[A*vx==b]prob=cvx。问题(目标，约束)结果=问题。解算(Verbose=True)Xat2=NP。阵列(VX.。值)。挤压()。
#重构信号Xat=Xat2。重塑(NX，NY)。T#堆栈列Xa=idct2(Xat)#如果不是NP，则确认解。全部关闭(X。托特。扁平[ri]，Xa.。托特。平面[ri])：print('；警告：样本索引处的值与原始不匹配。'；)#创建掩码的图像(用于可视化)掩码=np。零(X。形状)遮罩。托特。平坦[ri]=255 Xm=255*Np。一(X。Shape)XM。托特。平坦[ri]=X。托特。平坦[ri]。
好吧，结果不是很好。最左边的原始图像几乎看不懂。分辨率很低，所以我们必须取50%的大样本(布尔蒙版显示在左中；蒙版图像显示在右中)。无论如何，很明显这一过程是有效的：最右边的重建图像绝对接近于原始图像，尽管效果不佳。
考虑到我们的工作概念验证，有很多方法可以改进它。基于Kronecker的方法虽然易于实现，但被证明不适用于大图像。还有其他的方法吗？
使用CVXPY的凸优化不一定是求解$L^1$最小值的唯一方法。通过一些在线研究，我找到了L-BFGS算法6及其变体OWL-QN 3。我们对OWL-QN算法特别感兴趣，因为它允许通过最小化以下形式的函数来拟合$L^1$正则化模型：
其中$f$是可微凸损失函数，$C$是常数。在我们的例子中，我们可以将$f$定义为最小二乘目标函数，它只是残差平方的$L^2$范数：
现在剩下的就是编写代码了！在尝试了几个不同的选项之后，我最终决定使用libLBFGS(用C编写)来实现其OWL-QN。为了能够从Python访问它，我使用针对Python和Numpy的CAPI对其进行了包装。您可以在PyLBFGS找到我的实现。有关安装说明和基本用法，请参阅项目自述文件。如果你遇到虫子就告诉我。
关于libLBFGS(扩展为PyLBFGS)的好处是您可以随意定义目标函数。换句话说，我们不会被限制盲目地遵循$Ax=b$模型。重要的是我们能计算出残差平方的范数和它的梯度。我们根本不需要生成$A$！
下面的代码比文字更好地解释了我的意思。请特别注意传递给OWL-QN算法的EVALUE回调。
从pylbfgs导入owlqn def valuate(x，g，step)："；"；"；"；"；"；#我们希望返回两件事：#(1)残差的范数平方SUM((Ax-b).^2)和#(2)梯度2*A'；(Ax-b)#展开x列-first x2=x。重塑((NX，NY))。T#Ax正好是x2ax2=idct2(X2)#堆栈列和提取样本Ax=ax2的逆2DCT。托特。平坦[ri]。重塑(b.。形状)#计算剩余Ax-b及其2范数平方Axb=Ax-b Fx=Np。总和(np.。POWER(Axb，2))#将残差向量(K X 1)投影到空白图像(Ny X Nx)axb2=Np。零(x2.。形状)AXB2。托特。平面[ri]=Axb#Fill Columns-First#A'；(Ax-b)正好是axb2 AtAxb2=2*dct2(Axb2)AtAxb=AtAxb2的二维DCT。托特。重塑(x.。形状)#个堆栈列#复制梯度矢量Np。Copto(g，AtAxb)返回缩放图像的fx#个分数，以便在sample_size=(0.1，0.01)#read Original image Xorig=spimg处随机采样。Imread('；esch_water fall.jpeg'；)NY，NX，Nchan=Xorig。每个样本量Z=[Np的形状#。零(Xorig.。Shape，dtype='；uint8&39；)for s in Sample_Size]掩码=[np。零(Xorig.。Shape，dtype='；uint8&39；)for s in sample_size]for i，s in Enumerate(Sample_Size)：#创建随机采样索引向量k=round(nx*ny*s)ri=np。随机的。CHOICE(NX*NY，k，REPLACE=FALSE)#范围内j的每个颜色通道的索引随机样本#(Nchan)：#提取通道X=Xorig[：，：，j]。Squeeze()#创建遮罩的图像(用于可视化)xm=255*np。一(X。Shape)XM。托特。平坦[ri]=X。托特。平面[ri]掩码[i][：，：，j]=Xm#对图像进行随机采样，将其存储在向量b b=X中。托特。平坦[ri]。Astype(Float)#在内存中执行L1最小化Xat2=owlqn(NX*NY，Evaluate，None，5)#将输出转换回空间域Xat=Xat2。重塑(NX，NY)。T个堆栈列Xa=idct2(Xat)Z[i][：，：，j]=Xa。ASTYPE('；uint8'；)。
OWL-QN的快速C实现允许我们处理整个瀑布图像的样本，而不需要任何缩放。而且，我们现在可以处理图像的三个颜色通道中的每一个，而不是像以前那样一切都是灰度的。上面显示的解决方案真正展示了压缩感知的威力。左侧显示的是原始的全色图像。中间的图像是随机的10%样本。解决方案图像在右侧。尽管该解决方案由于颜色通道混合不良而包含一些明显的瑕疵，但总体精确度是不可思议的。此外，只要格外小心和注意，这些瑕疵可以通过后处理(例如，高斯滤波器)或通过考虑颜色通道的改进的压缩感测实现来去除。一种可能性是尝试并预先确定可能的调色板，然后将其合并到压缩感测例程中。
为了好玩，我还包括了一张从1%的可用数据重建的图像。它绝对是模糊的，但仍然可以辨认！
作者：Kutz，J.N.。“数据驱动的建模和科学计算：集成复杂系统和大数据的动力学的方法”(2013年)。↩
题名/责任者：Emmanuel J.。“关于压缩抽样的介绍。”信号处理杂志，IEEE25.2(2008年)：21-30。↩。
安德鲁、加伦和高剑锋。“L1正则化对数线性模型的可伸缩训练。”第24届机器学习国际会议论文集。ACM，2007。↩。
维基百科的贡献者。“压缩传感。”维基百科，免费的百科全书。维基百科，免费百科全书，2016年3月26日。万维网。5月26日。2016年，中国↩。
维基百科的贡献者。“Kronecker产品。”维基百科，免费的百科全书。维基百科，免费百科全书，1。
.