压缩传感引论(shawei)

“压缩传感”引论香港大学电机电子工程学系高效计算方法研究小组沙威2008年11月20日Email:[e-mail]:~wsha到了香港做博士后，我的研究领域继续锁定计算电磁学。我远离小波和计算时谐分析有很长的时间了。尽管如此，我仍然陆续收到一些年轻学者和工程人员的来信，和我探讨有关的问题。对这个领域，我在理论上几乎有很少的贡献，唯一可以拿出手的就是几篇网上发布的帖子和一些简单的入门级的代码。但是，从小波和其相关的理论学习中，我真正懂得了一些有趣的知识，并获益良多。我深切地感到越来越多的学者和工程师开始使用这个工具解决实际的问题，我也发现互联网上关于这方面的话题多了起来。但是，我们永远不能只停留在某个阶段，因为当今学术界的知识更新实在太快。就像我们学习了一代小波，就要学习二代小波；学习了二代小波，就要继续学习方向性小波（X-let）。我也是在某个特殊的巧合下不断地学习某方面的知识。就像最近，我的一个友人让我帮她看看“压缩传感”（CompressiveSensing）这个话题的时候，我的兴趣又一次来了。我花了一个星期，阅读文献、思考问题、编程序、直到写出今天的帖子。我希望这篇帖子，能对那些没进入且迫切想进入这个领域的学者和工程师有所帮助。并且，我也希望和我一个星期前一样，对这个信号处理学界的“一个大想法”（ABigIdea）丝毫不了解的人，可以尝试去了解它。我更希望，大家可以和我探讨这个问题，因为我到现在甚至不完全确定我对压缩传感的某些观点是否正确，尽管我的简单的不到50行的代码工作良好。在这个领域中，华裔数学家陶哲轩和斯坦福大学的统计学家DavidDonoho教授做出了重要的贡献。在这个引言中，我用简单的关键字，给出我们为什么需要了解甚至是研究这个领域的原因。那是因为，我们从中可以学习到，下面的这些：矩阵分析、统计概率论、拓扑几何、优化与运筹学、泛函分析、时谐分析（傅里叶变换、小波变换、方向性小波变换、框架小波）、信号处理、图像处理等等。所以，我们有什么理由，拒绝这个有意思的东西呢？让我们开始吧。传统思路——正交变换对于一维的信号1×∈NRx，大多数情况下，信息是冗余的。我们可以通过正交变换的方法来压缩它。正变换：xyΨ=，反变换yxHΨ=。这里，IHH=ΨΨ=ΨΨ，NNC×∈Ψ，I是单位矩阵。对于1×∈NCy，能量较x集中，本质上去除了x中的相关性。因此，我们只保留K个较大的分量，而把其它KN−个置为零。通过反变换，我们能够近乎完美的重建原始信号。因为，那KN−个变换域系数的贡献，实在微乎其微。具有这样性质的信号被称为K“稀疏”（Sparsity）的。于是，我们有了如下编码解码的策略：编码：构造Ψ，做正变换xyΨ=，保留y中最重要的K个分量，和其对应的位置。解码：把K个分量放回到对应的位置，其它位置填0，构造HΨ，反变换yxHˆˆΨ=。而解码能否近乎得到原始信号呢？显然，我们希望δ≤−=−22||ˆ||||ˆ||yyxx，δ是一个小的常数。但更有效的是用相对误差δ≤−22||||/||ˆ||yyy。但这种编码解码方法有些缺点：1、考虑到香农（Shannon）采样定理，为了获得很好的信号分辨率，采样间隔会很小，造成了原始信号长度会很长，因此变换过程会消耗很长的时间。2、K个需要保留的重要分量的位置，是随着信号的不同而不同的。因此，这种策略是“自适应”（Adaptive）的，且需要分配多余的空间存储这些位置。3、一旦在传输过程中K个分量中的某几个丢失了，后果可想而知。如果我们制作一个音频设备，1将带来电力的消耗和用户的不满，2将带来存储空间的增加，3将带来较差的抗干扰能力。新的思路——压缩传感压缩传感（CompressiveSensing）是一个很有意思的新的方向。它也正成为信号处理领域的“ABigIdea”。对于信号1×∈NRx，我们可以找到它的M个线性测量（LinerMeasurement），xsΦ=。NMR×∈Φ。这里，Φ的每一行可以看作是一个传感器（Sensor）,它与信号相乘，拾取（Acquisition）了信号的一部分信息。拥有了这M个测量和Φ，我们就可以近乎完美的重构原始信号了。听起来“相当”传奇，事实上，它基于如下严格的数学最优化（Optimization）问题：目标函数0||ˆ||miny，且满足等式约束syH=ΦΨˆ或者，可以写成02||ˆ||||ˆ||minyysHλ+ΦΨ−求解该最优化问题，得到变换域的yˆ，然后反变换，便可以得到时域的xˆ。公式中的2是我们熟悉的2-范数，而0是什么呢？是0-范数，也就是向量yˆ中非零元素的个数。看起来很有道理，因为yˆ是待求的变换域向量，它是K稀疏的。使yˆ非零元素的个数尽量小，也就是保留了尽量少的重要的K个分量，显然这几个分量可以近乎完美重构x。我们回到传统的思路，这K个分量是我们在变换域“自适应”找的，而该优化算法也可以使我们找到这K个分量。这就足够了吗？显然不行，我们仍然没有探讨测量矩阵Φ需要满足的性质。我们用极限分析法。如果我们把Φ构造成和HΨ极端相似（Coherence）的矩阵，也就是拿出Ψ的前M行。用这个算法求yˆ，我们将得到⎟⎟⎠⎞⎜⎜⎝⎛=×−×1)(10ˆMNMsy，这显然是错误的。也就是说，你强迫的认为前M个变换域分量是重要的。而事实是，重要的K个分量的位置我们事先是不知道的，是随着信号的不同而不同的。当然，你可以将Φ恰好构造成对应最重要分量的K行，得到正确的结果。而这种的做法要付出的概率代价KNC1。也就是说，你必须穷举KNC次，才能得到你想要的结果。但是，即使你有幸碰到了它，也并不能肯定这个结果就是对的。因此，我们选择Φ和HΨ极端不相似（ExtremelyIncoherence）。于是，Φ很大程度上和随机（Randomness）这个词相联系，它可以是满足高斯分布的白噪声矩阵，或贝努里分布的1±矩阵（也称作Noiselet）等等。除此之外，我们希望线性测量有稳定的能量性质：δδ+≤ΦΨ≤−1||ˆ||||ˆ||122yyH，也就是它要保持K个重要分量的长度。综合上面的，我们有了如下编码解码的策略：编码：构造Φ，生成测量xsΦ=，保留s。解码：构造同样的Φ，构造任一种正交变换HΨ，根据s重构x。压缩传感的优势：1、非自适应（Non-Adaptive）的，一开始就可以传输长度较短的信号，甚至突破采样定理的极限。2、抗干扰，s中任何一项都是重要的，或者说不重要的。丢失了某几项，仍然可以完美重构。它的缺点：1、实际中，s的长度一般是重要分量长度的4倍，才能近乎完美重构。数学上更严格的，KM4≈或者)(log2KNKM≥。2、重构（恢复）算法是NP问题。即使将0-范数转化为1-范数，由于其不可微性（Indifferentiable），算法的计算复杂度仍然很高。它的应用前景广泛：低成本数码相机和音频采集设备；节电型音频和图像采集设备；天文（图像本身就稀疏，例如天空的星星）；网络；军事（用很简易的摄像机随机记录场景，可以完全重构军事地图）；超宽带（雷达信号处理）。这里，值得指出的是，美国的工程学家已经设计出了实际的产品。快速算法——正交匹配追踪对于0-范数的优化问题，实际上是NP问题，就是在多项式时间内难以求解，甚至无法验证解的可靠性。于是，我们必须将0-范数换一下，变成1-范数。为什么不是2-范数呢，那样就会简单多了。毕竟2-范数的优化问题可以转化成2次型问题，而1-范数，∑=iiyy|ˆ|||ˆ||1，在0点处不可导，因此无论是梯度算法、矩阵求导等等手段都变得相形见绌。因此，基于1-范数的优化算法需要特殊处理，且复杂度很高。下面我们来解释下为什么要用1-范数，而不是2-范数。我们令恢复（Recovery）矩阵HTΦΨ=，则等式约束可重写为：syT=ˆ。yˆ中未知数有N个，方程只有M个，且NM。因此，方程有无穷多解。从几何上说，0ˆ=−syT是一个超平面，为了简化，在2-D问题中（1=K，yˆ只有两个元素待求）可认为它就是一条直线。而范数约束呢？0-范数是一个十字架，因此它的最外侧（范数的最小值）是4个点。所以其和直线的交点，必然在坐标轴上。也就是说，能使yˆ产生更多的0，这正是我们想要的“稀疏”的结果。2范数是一个圆，因此它的最外侧边界和直线的交点（就是切线的概念），以压倒性的概率不在坐标轴上，除了直线的斜率恰好为0或者无穷大。其实直线的斜率恰好为0或者无穷大，是不可能的，因为Φ和HΨ极端不相似。只有Φ取Ψ的某一行时，两者相似，才会发生斜率恰好为0或者无穷大的情况（因此，你的胜算只有12/1C，但你不知道哪个是对的。）。依上所述，用2-范数优化的结果，使yˆ几乎没有0，这是我们不期望的。而1-范数是一个菱形，四个角都在坐标轴上，因此它和直线的交点以压倒性的概率落在坐标轴上。这就是我们要用1-范数的原因。事实上，p范数满足，10≤p，它的外边界都向中心凹（Concave）。而1p，外边界向外凸（Convex）。所以前者的外边界和直线的交点无疑落在坐标轴上。根据这个几何解释，我们可以将问题转化成：12||ˆ||||ˆ||minyyTsλ+−这显然是一个非线性（Non-Linear）的凸（Convex）优化问题。众所周知，对于优化问题，我们一般用梯度的方法来求解。而对1||ˆ||y，在0点导数不存在，因为这个点正好位于两条直线的交点上，左右导数不相等。这也正是很多数学方法的考虑。像子梯度（Subgradient）法、平滑近似法（SmoothApproximation）等等。我们暂不谈这些方法，因为它需要特殊的数学背景。我们谈一谈工程领域最常用的正交匹配追踪法（OrthogonalMatchingPursuit）。它的思想本质上还是来自于这个K“稀疏”。我们绕了一圈，还是为了找这K个关键的分量。既然是关键，显然它的系数的绝对值应该比其它KN−个分量大得多。为了简单起见，我们先假设1=K，唯一非零元素qyˆ在yˆ中对应的位置在q。于是yTˆ就是恢复矩阵T的第q列qT与yˆ中的非零元素qyˆ的乘积，即qqqsyT=ˆ。且δ−22||||/||||sssq。换句话说，T的第q列与s的相似程度最高，即qrsTsTsTsTrHrHqq≠==|,,||||||,|。所以，我们只要计算恢复矩阵T的所有列与s的内积，找到内积绝对值最大的那列就行了，该列对应的位置就是q。根据最小二乘法，sTTTyHqqHqq1)(ˆ−=，就是使2||ˆ||qqyTs−最小的那个qyˆ。呵呵，感觉到了吗，实际上这有点或者非常像施密特（Schimidt）正交化方法。余量qqqqnTTTsTsr−=,,，始终同qT正交。这也是为什么这个方法叫“正交”匹配追踪的意思了。而匹配，显然是你找到了最大的|,|sTq。对于1K呢，其实是差不多的。我们找到余量nr同T中所有列向量最大的那个即可（但第一次找到的那列要排除，因为它已经保留了下来。）。于是，找到使()21212||ˆˆ,||⎟⎟⎠⎞⎜⎜⎝⎛−qqqqyyTTs最小的那个⎟⎟⎠⎞⎜⎜⎝⎛12ˆˆqqyy。这里，1qT是我们第一次找到的那一列，2qT是我们新找到的那一列（也要记住它的列号2q）。可以看出，⎟⎟⎠⎞⎜⎜⎝⎛=12ˆˆˆqqqyyy被更新了，由原来的一个变成两个了，也就是我们找到两个在变换域最关键的元素和其在yˆ中对应的位置了。令()12,qqqTTT=，余量nr又一次被写为：qqqqnTTTsTsr−=,,。我们再一次看到了施密特正交化的影子。继续上面的步骤，直至找到变换域所有K个最重要的分量。也就是说，正交匹配追踪的迭代次数Km≥。实际操作上只要满足δ22||||/||||srn，迭代就可以中止了。最后，假设正交变换Ψ是傅里叶变换，我们来分析下计算复杂度。生成恢复矩阵T，采用FFT快速算法，需要计算复杂度))(log(2NMNO；找到T中匹配的列向量需要)(mNMO次操作，注意m是算法需要的迭代次