在本文中，我将研究Python中的压缩感知（也称为压缩感知、压缩采样和稀疏采样）。由于压缩感知的概念可以应用于各种各样的主题，我将主要关注如何将其应用于声音和图像等一维和二维领域。具体来说，我将展示如何获取高度不完整的信号样本数据集，并重建底层的声音或图像。这是一种非常强大的技术。
你可能知道，有许多不同类型的规范。也许最常见和最广泛认可的是$L^2$norm：
$L^2$norm很好，因为它易于计算、易于区分，并且具有直观的吸引力（例如，向量的norm是其长度）。许多非常重要的算法和方法都依赖于$L^2$，包括最小二乘拟合。
也就是说，L^2$norm并不是所有东西的goto解决方案。其他规范也有许多有趣和有用的特性。考虑$^ ^ 1美元的标准：
它不求每个元素的平方，只取其绝对值。尽管绝对值令人恼火，因为它经常在其导数中引入不连续性，但与在$L^2$norm中发生的平方相比，它确实有一些独特的性质。压缩感知就是利用这些特性。
让我们用Python可视化一些数据，看看我在说什么。
#确保你'；我们安装了以下软件包import numpy as np import matplotlib as mpl import matplotlib。pyplot作为plt导入scipy。优化为spopt import scipy。FFT包装为spfft导入scipy。ndimage为spimg导入cvxpy为cvx
首先我们要做的是创建一些任意的线性数据，包括一些噪声。让我们使用合成方程式：
#生成一些噪声x=np的数据。排序（np.random.uniform（0,10,15））y=3+0.2*x+0.1*np。随机的兰登（兰（x））
现在让我们在数据样本上拟合两行。对于第一行，我们将使用$L^1$norm作为良好匹配的标准；对于第二行，我们将使用$L^2$norm。
#求L1直线拟合L1_拟合=λ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=lambda 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来调整数据。复制（）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=5000 t=np。linspace（0，1/8，n）y=np。sin（1394*np.pi*t）+np。sin（3266*np.pi*t）yt=spfft。dct（y，norm='；正交'；）
在上面的图中，我们看到信号有一个清晰的模式，但不是微不足道的。最上面一行的曲线图是不同尺度下时域中的信号。最下面一行中的图是频谱域中的信号（即信号的频率内容）。特别是考虑到频域，我们注意到除了代表两个正弦频率的两个尖峰之外，频谱大部分为零。
现在想象一下对10%的时间信号进行采样（见下文）。你会有一个数据集，在肉眼看来，看起来像胡说八道。基本信号仍然是相同的，其频率内容也是相同的（除两个尖峰外，大部分为零）。有人可能会问，是否有可能从不完整的数据中提取这两个主频，以便我们重建信号？答案是肯定的！
#提取信号m=500的小样本#10%样本ri=np。随机的选择（n，m，replace=False）#指数ri的随机样本。排序（）#排序不是严格必要的，但便于绘制t2=t[ri]y2=y[ri]
在这种情况下，由于信号的频率内容非常稀疏，压缩感知成为可能。这就是$L^1$norm发挥作用的地方。我们要做的是，在所有可能的信号中，找到与已知数据匹配的最简单的信号。换句话说，我们希望使用一个最小化例程来找到一组满足两个条件的频率：（a）基本信号与我们的数据完全匹配（或尽可能接近）；（b）频率的$L^1$norm最小化。这样的例行程序将产生一个稀疏的解决方案——这正是我们想要的。
在Python中，有几种方法可以实现这一点。也许最简单的方法是利用凸优化库CVXPY。使用下面的代码最小化信号频率的范数，并限制候选信号应与我们的不完整样本完全匹配。
#创建idct矩阵运算符A=spfft。idct（np.同一性（n），常模='；正交'，轴=0）A=A[ri]#进行L1优化vx=cvx。变量（n）目标=cvx。最小化（cvx.norm（vx，1））约束=[A*vx==y2]prob=cvx。问题（目标、约束）结果=问题。求解（verbose=True）
你可能会问：那该死的$A$matrix是什么？这是整个派对的关键。让我解释一下。
为了实现极小化，我们必须以某种方式将问题转化为线性方程组：
具体来说，我们想要导出一个矩阵$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的图像瀑布进行随机采样和重建（约1200×1600像素）。由于$A$矩阵带来的内存限制，我们将首先考虑缩小版本的图像（约50×65像素）。在接下来的部分中，我们将扩展例程以处理大型图像。
请注意，SciPy不提供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，范数='；正交'；，轴=0）。T，norm='；正交'，轴=0）def idct2（x）：返回spfft。idct（spfft.idct（x.T，norm='；正交'；，轴=0）。T，norm='；正交'，axis=0）#读取原始图像并缩小尺寸以获得速度xOrg=spimg。imread（'；escher#u瀑布.jpeg'；，flatte=True，mode='；L'；）#读取灰度X=spimg。缩放（Xorig，0.04）ny，nx=X。形状
与上一节一样，我们将随机抽取图像索引样本，形成我们的$b$矩阵。然后，我们将生成$A$矩阵。
为2D图像数据创建$A$矩阵需要比在1D情况下更巧妙的技巧。在接下来的推导中，我们将使用Kronecker乘积$\otimes$和2D离散余弦变换是可分离的这一事实来生成运算符$A$。
假设$X$是光谱域中的图像，$D_i=\textit{idct}（i_i）$，其中$i_i$是大小为$i$的单位矩阵。然后：
如果$\textit{vec}（X）$是将$X$的列堆叠在彼此之上的向量运算符，则：
显然，Kronecker产品是我们想要的转换矩阵$\psi$。因此，我们的矩阵$A$变成了$A=\phi（D_n\otimes D_m）$，其中$\phi$是采样矩阵。你可以用Numpy计算Kronecker积。克朗。这种方法的主要问题是Kronecker产品可以很快变得非常巨大。如果你的目标图像是$m$乘$n$，并且你要采集$k$样本，那么$A$矩阵的大小为$（mnk）^2$。这就是说，对于小图像来说，这是很好的。
#提取信号k=round（nx*ny*0.5）#50%样本ri=np的小样本。随机的选择（nx*ny，k，replace=False）#指数b=X的随机样本。T平坦[ri]b=np。展开_dims（b，轴=1）#使用kron（大ny*nx的内存错误）A=np创建dct矩阵运算符。kron（spfft.idct）（np.identity（nx），norm='；正交'，轴=0），spfft。idct（np.身份（ny），常模='；正交'，轴=0）A=A[ri，：]#与φ乘以kron相同#进行L1优化vx=cvx。变量（nx*ny）目标=cvx。最小化（cvx.norm（vx，1））约束=[A*vx==b]prob=cvx。问题（目标、约束）结果=问题。求解（verbose=True）Xat2=np。数组（vx.value）。挤压（）
#重建信号Xat=Xat2。重塑（nx，ny）。T#堆栈列Xa=idct2（Xat）#如果不是np，则确认解决方案。allclose（X.T.flat[ri]，Xa.T.flat[ri]）：打印（'；警告：样本索引处的值与原始值不匹配。'；）创建蒙版图像（用于可视化）蒙版=np。零（X形状）遮罩。T平坦[ri]=255 Xm=255*np。一个（X.形状）Xm。T平坦[ri]=X。T平[里]
好吧，结果不太好。最左边的原始图像几乎无法理解。分辨率很低，所以我们必须取50%的大ish样本（布尔掩码显示在左中；掩码图像显示在右中）。不管怎样，这一过程显然奏效了：最右边的重建图像肯定与原始图像很接近，尽管效果很差。
考虑到我们正在进行的概念验证，有很多方法可以改进。基于Kronecker的方法虽然易于实现，但对于大型图像来说是不可用的。还有什么其他方法？
使用CVXPY的凸优化不一定是找到$L^1$最小值的唯一方法。通过一点在线研究，我找到了L-BFGS算法6及其变体OWL-QN3。我们对OWL-QN算法特别感兴趣，因为它允许通过最小化以下形式的函数来拟合$L^1$正则化模型：
其中，$f$是可微凸损失函数，$C$是常数。在我们的例子中，我们可以将$f$定义为最小二乘目标函数，它只是剩余平方的$L^2$范数：
现在剩下的就是编码了！在尝试了几个不同的选项后，我最终决定使用libLBFGS（用C编写）来实现OWL-QN。为了从Python中访问它，我使用Python和Numpy的C API对其进行了包装。你可以在PyLBFGS找到我的实现。有关安装说明和基本用法，请参阅项目自述。如果你遇到虫子，请告诉我。
libLBFGS（扩展为PyLBFGS）的好处在于，您可以随意定义目标函数。换句话说，我们不必盲目遵循$Ax=b$模型。重要的是我们能够计算剩余平方的范数及其梯度。我们根本不需要生成$A$！
下面的代码比用文字更好地解释了我的意思。请特别注意传递给OWL-QN算法的evaluate回调。
从pylbfgs导入owlqn def评估（x，g，步骤）：&34""内存中的求值回调""" # 我们想要返回两件事：#（1）残差的范数平方，和（（Ax-b）^2） 和#（2）梯度2*A和#39；（Ax-b）#首先展开x列x2=x。重塑（（nx，ny））。T#Ax只是x2 Ax2=idct2（x2）#堆栈列和提取样本Ax=Ax2的逆2D dct。T平坦[ri]。重塑（b.形状）#计算剩余的Ax-b及其2-范数平方Axb=Ax-bfx=np。求和（np.power（Axb，2））#将剩余向量（kx1）投影到空白图像（ny x nx）上Axb2=np。零（x2.形状）Axb2。T扁平[ri]=Axb#先填充列#A和#39；（Ax-b）只是Axb2 AtAxb2=2*dct2（Axb2）AtAxb=AtAxb2的2D dct。T重塑（x.shape）#堆叠列#复制渐变向量np。copyto（g，AtAxb）返回fx#缩放图像的分数，以随机采样，采样大小=（0.1,0.01）#读取原始图像Xorig=spimg。imread（'；escher#u瀑布.jpeg'；）ny，nx，nchan=xOrg。形状#对于每个样本大小Z=[np.zero（Xorig.shape，dtype='；uint8'；）对于样本大小为[u]的s，掩码=[np.zero（Xorig.shape，dtype='；uint8'；）对于样本大小中的s]对于i，枚举中的s（样本大小）：#创建随机抽样索引向量k=round（nx*ny*s）ri=np。随机的选择（nx*ny，k，replace=False）#索引的随机样本#范围内j的每个颜色通道（nchan）：#提取通道X=Xorig[：，：，j]。挤压（）#创建遮罩图像（用于可视化）Xm=255*np。一个（X.形状）Xm。T平坦[ri]=X。T平面[ri]掩模[i][：，：，j]=Xm#对图像进行随机采样，将其存储在向量b=X中。T平坦[ri]。astype（float）#在内存Xat2=owlqn（nx*ny，evaluate，None，5）中执行L1最小化#将输出转换回空间域Xat=Xat2。重塑（nx，ny）。T#stack columns Xa=idct2（Xat）Z[i][：，：，j]=Xa。aType（'；uint8'；）
OWL-QN的快速C实现允许我们处理整个瀑布图像的样本，而无需任何缩放。我们现在可以处理图像的三个颜色通道中的每一个，而不是像以前那样在灰度中进行所有操作。上面显示的解决方案真正展示了压缩感知的威力。原始的全彩图像显示在左侧。中间的图像是随机的10%样本。解决方案图像位于右侧。虽然该解决方案包含一些明显的瑕疵，由于不良的颜色通道混合，整体精度是不可思议的。此外，只要格外小心和注意，这些瑕疵可能会被去除——要么通过后处理（例如高斯滤波器），要么通过考虑颜色通道的改进压缩感知实现。一种可能性是尝试预先确定可能的调色板，然后将其纳入压缩感知例程。
我还包括了一张根据1%的可用数据重建的图像，这只是为了好玩和咯咯笑。它肯定是模糊的，但仍然可以辨认！
数据驱动的建模和科学计算：整合复杂系统动力学和大数据的方法(2013). ↩
坎德、伊曼纽尔·J.和迈克尔·B·沃金。“压缩采样简介。”信号处理杂志，IEEE 25.2（2008）：21-30。 ↩
安德鲁、盖伦和高剑锋。“L1正则对数线性模型的可扩展训练。”第24届国际机器学习会议记录。ACM，2007年。 ↩
维基百科撰稿人。“压缩感知。”维基百科，免费的百科全书。维基百科，免费百科全书，2016年3月26日。网状物5月26日。2016. ↩
维基百科撰稿人。“克罗内克产品。”维基百科，免费的百科全书。维基百科，免费百科全书，1
......