简介本资源是一套面向信号处理与压缩感知领域初学者及科研实践者的稀疏重构算法MATLAB实现代码包聚焦FOCUSS等主流算法原理验证与工程复现。资源共10个文件全部为.m脚本含FOCUSS_Single.m、FOCUSS_Multiple.m、OMP_fun.m、BP.m、BPDN_homotopy_function.m等核心模块覆盖单/多通道FOCUSS迭代重构、正交匹配追踪、基追踪与稀疏贝叶斯学习等关键方法代码结构清晰、注释完整便于理解算法流程、调试参数及拓展应用。压缩包仅9KB轻量易用适合作为课程实验、毕设实现或算法对比研究的起点。目前已有622人学习下载读者可直接运行main.m快速启动全流程演示获取从初始化、残差更新、稀疏权值迭代到收敛判断的完整逻辑链并通过各子函数深入掌握测量矩阵作用、对偶变量更新、逆矩阵近似等关键技术细节。 我当时接手一个压缩感知信号恢复方向的实验项目手里的资料里正好躺着一份“几种常见的稀疏重构算法代码.rar”。解压之后发现里面装着FOCUSS、OMP、BP等算法的MATLAB实现信息量大得一口气看不完。那阵子为了把FOCUSS这套稀疏重构算法跑通、再把每个细节吃透我前后折腾了快一周踩了不少坑也把代码包里的思路彻底整理清楚了。今天这篇就基于这份算法包从原理、代码、参数调优到算法选型把FOCUSS和稀疏重构这件事一次讲透。不管你是在校学生刚接触压缩感知还是工程里需要做稀疏信号重构的研发人员这篇应该都能帮你省下不少弯路。1. 这份稀疏重构算法代码包到底在解决什么问题1.1 从一个欠定方程说起压缩感知要面对的最核心问题可以用一个很简单的数学式子描述y A xy 是观测向量A 是感知矩阵x 是待恢复的原始信号。传统采样理论要求观测数量大于等于信号长度可是在压缩感知里A 的行数 m 远远小于列数 n也就是说我们手里拿到的方程数量远少于未知数数量。这种欠定方程在数学上本来有无穷多组解如果没有额外约束根本无法唯一确定 x。压缩感知给出的额外约束是x 本身是稀疏的或者说在某个变换域里是稀疏的。所谓稀疏就是向量里绝大多数元素是零只有少数位置上有非零值。这个约束非常有力量足够从 m 个观测里恢复出 n 维信号甚至可以做到 m 远小于 n 也能高概率恢复。FOCUSS、OMP、BP 这些算法本质上都是在解同一个问题在 A x y 的所有解里找出那个最稀疏的 x。1.2 代码包里的算法体系这份算法包虽然叫“稀疏重构算法代码”但里面的内容不是一套代码而是围绕稀疏重构问题展开的一套算法集合。整理下来主要分几类基于凸优化的算法比如 BPBasis Pursuit或者说 Lasso 形式核心思想是把非凸的 l0 范数问题松弛成 l1 范数凸优化问题。基于贪婪迭代的算法比如 OMPOrthogonal Matching Pursuit核心思想是每次从感知矩阵里挑选与当前残差最相关的原子逐步逼近真实支撑集。基于迭代加权最小二乘的算法就是我们这篇的重点 FOCUSSFocal Underdetermined System Solver核心思想是用一个权重矩阵反复迭代把能量逐步聚焦到少数大系数上。还有一类是阈值类迭代算法比如 IHTIterative Hard Thresholding、CoSaMP这些在工程里也经常用到代码包里有些版本收录了有些没有。我拿到代码后的复现顺序是先跑通最简单的 OMP再跑 BP最后花最多时间啃 FOCUSS。原因是 FOCUSS 的数学形式看起来不太直观代码里又藏了很多参数细节不看穿原理很容易调不出好效果。1.3 这类代码包适合谁看如果你是在做压缩感知实验、信号恢复、图像重建、通道估计、DOA 估计这类方向那这份代码包正好命中你的需求。就算你不是做学术研究只是想给自己系统补一补稀疏重构的背景知识把 FOCUSS 这类迭代算法的思路理清楚也会有不小收获。我自己复现这份代码包的过程可以简单概括成三步先理解数学原理再对照代码逐行验证最后用不同算法做横向对比。接下来我就按这个顺序来写。2. FOCUSS算法原理迭代加权最小二乘怎么从无穷多解里聚焦出稀疏解2.1 为什么最小二乘解不稀疏面对 y A x 这个欠定方程最本能的想法是求最小二乘解也就是在满足约束的前提下让 ||x||2 最小。这个解有闭式表达x_ls A^T (A A^T)^{-1} y问题在于最小二乘解从数学性质上几乎不会稀疏。它把能量均匀地摊到了所有基上或者说摊到了感知矩阵的各个列方向上。你可以把它理解成一个人分糖果手里糖少、孩子多最“公平”的做法是每个孩子都分一点结果就是没有一个孩子手里的糖特别多。但稀疏重构想要的恰恰相反我们希望大部分孩子一颗糖都拿不到只有极少数孩子拿到全部糖。这就是 l2 范数偏好“能量均匀分布”和 l0/l1 范数偏好“能量集中”的本质差异。2.2 加权最小二乘的思想既然普通最小二乘不稀疏那自然想到能不能给不同位置加不同的权让迭代过程逐渐朝稀疏方向倾斜FOCUSS 的做法非常直接。它维护一个对角权重矩阵 W每次迭代都求解一个加权最小二乘问题min ||W^{-1} x||2^2约束 y A x这个问题的解是x_new W (A W)^\dagger y其中 (A W)^\dagger 表示伪逆。如果 A W 行满秩伪逆可以展开成(A W)^\dagger W A^T (A W^2 A^T)^{-1}于是迭代公式变成x_new W^2 A^T (A W^2 A^T)^{-1} y权重矩阵 W 不是随便取的它在每轮迭代里根据当前解的幅度来更新。FOCUSS 常用的权重形式是W_k diag( |x_k|^{1-p/2} )其中 p 是 0 到 1 之间的稀疏度控制参数。这个式子看着有点绕但思想其实很朴素上一轮哪个位置的系数绝对值大这一轮就给那个位置更大的权重。权重大意味着在加权最小二乘里那个位置的“成本”更低于是解会更倾向于保留甚至放大这个大系数。而那些小系数位置权重小解中对应分量会被进一步压缩。每迭代一轮能量就往最大的几个系数上聚集一分最终得到一个稀疏解。这正是 FOCUSS 名字里 “Focal” 的含义——聚焦。它像用放大镜反复把光线聚到同一个点上最终把能量集中起来。2.3 参数 p 如何控制稀疏度参数 p 是 FOCUSS 里最关键的旋钮。当 p 1 时W diag(|x|^{0.5})对应的极限目标函数与 l1 范数最小化有密切联系可以理解成一种对 BP 的迭代逼近。这个区域算法相对稳健不容易发散但得到的解可能不会特别“极端”稀疏。当 p 趋近 0 时W diag(|x|)权重与当前幅度近似成正比聚焦效应更强解的整体稀疏度更高。当 p 1 时目标函数本质上是非凸的因此存在局部极值问题。p 越小非凸性越强越容易陷入不理想的局部解。实际使用时p 取 0.3 到 0.7 之间是比较常见的区间。太小容易把噪声放大或者收敛到错误的支撑集太大则稀疏度不够。我自己的经验是先拿 p 0.5 跑一遍再看重构误差曲线决定要不要调小或者调大。2.4 收敛过程从“均匀解”到“稀疏解”FOCUSS 的初始化解可以选最小二乘解也可以直接全 1 初始化。随着迭代进行你会观察到这样的变化初始几轮解的所有分量都差不多大随后少数位置开始出现优势其他位置逐渐被压小最终大部分位置趋近于零剩下几个位置上的值稳定在真实非零系数附近。这个过程很像“淘汰赛”每轮权重更新都在给上一轮表现好的位置“加油”给表现差的位置“踩刹车”最终胜出的少数位置就是支撑集的候选。理解这一点你也就能理解为什么 FOCUSS 对初始化比较敏感——如果一开始就给错误的位置比较大的权重迭代过程中可能很难纠正过来。这个坑我在后面专门展开讲。3. 动手实现FOCUSS核心代码逐行拆解3.1 FOCUSS核心函数实现原始代码包里的 FOCUSS 是 MATLAB 写的核心计算其实就几行。我用 Python 重新实现了一份等价版本方便更多人跑实验。下面这个函数可以直接用import numpy as np def focuss(A, y, p0.5, max_iter300, tol1e-6): FOCUSS 稀疏重构算法 参数 ---- A : ndarray, shape (m, n) 感知矩阵 y : ndarray, shape (m,) 观测向量 p : float 稀疏度控制参数通常取 [0, 1] max_iter : int 最大迭代次数 tol : float 相对变化量阈值用于提前终止迭代 返回 ---- x : ndarray, shape (n,) 重构出的稀疏信号 m, n A.shape # 初始化解这里用全1分布也可以改为最小二乘解 x np.ones(n) / np.sqrt(n) for i in range(max_iter): # 1. 根据当前解构造权重矩阵 # 加 1e-8 是为了防止零权重导致的数值退化 w (np.abs(x) 1e-8) ** (1 - p / 2) W np.diag(w) # 2. 加权后的感知矩阵 A_tilde A * W A_tilde A W # 3. 解加权最小二乘问题 # min ||xt||_2^2 s.t. A_tilde * xt y # 这一步可以替换为 np.linalg.lstsq(A_tilde, y, rcondNone)[0] xt np.linalg.pinv(A_tilde) y # 4. 变换回原坐标 x_new W xt # 5. 收敛判断相对变化量小于阈值则停止 err np.linalg.norm(x_new - x) / (np.linalg.norm(x) 1e-12) x x_new if err tol: break return x这段代码每行都有它存在的理由我拆开讲几个关键点。第一权重向量 w 的计算里为什么加 1e-8因为如果某轮迭代里某个位置被压到了接近零下一轮权重就会接近零再下一次这个位置就永远起不来了而且对角矩阵接近奇异会导致伪逆求解特别不稳定。加一个小小的正则项相当于给所有位置一个最低权重避免数值层面直接“判死刑”。这个微小改动在低信噪比场景下效果差异非常明显。第二第三步用 pinv 而不是直接求逆是因为 A_tilde 在欠定条件下大概率不是方阵根本不存在唯一逆矩阵。伪逆在数学上给出了最小二乘意义下的最优解而 FOCUSS 的核心就是在每次迭代中解一个约束最小二乘问题。如果数据规模较大用 pinv 速度会偏慢可以改用 lstsq效果一致但更快更稳。第三x_new W xt 这一步对应的是坐标变换回原空间。权重矩阵 W 把加权空间里的解映射回原始坐标空间。这一步如果遗漏整个迭代结果就是错的。3.2 一个完整的稀疏重构仿真光有函数还不够我写了一个完整的仿真脚本用来验证 FOCUSS 确实能恢复稀疏信号。实验设置如下信号长度 n 256稀疏度 k 10即只有 10 个非零元素观测数 m 64远小于 n感知矩阵 A 的每个元素从标准正态分布里独立采样稀疏信号 x_true 的非零位置随机选取非零值在 [-2, -1] ∪ [1, 2] 里随机分布观测 y A x_truenp.random.seed(42) n 256 m 64 k 10 # 生成感知矩阵 A np.random.randn(m, n) / np.sqrt(m) # 生成 k-稀疏信号 x_true np.zeros(n) idx np.random.choice(n, k, replaceFalse) x_true[idx] np.random.uniform(-2, 2, k) x_true[idx] np.where(x_true[idx] 0, np.random.uniform(1, 2, k), np.random.uniform(-2, -1, k)) # 观测 y A x_true # FOCUSS 重构 x_hat focuss(A, y, p0.5) # 重构误差 err np.linalg.norm(x_hat - x_true) / np.linalg.norm(x_true) print(f相对重构误差: {err:.6f}) print(f检测到的非零位置数: {np.sum(np.abs(x_hat) 1e-3)})在我本机上这个实验的相对重构误差通常在 1e-6 量级检测到的非零位置数基本等于真实稀疏度 k。也就是说FOCUSS 在这个实验设定下能够高精度恢复稀疏信号。为了直观确认恢复效果一般还会画三张图真实信号、观测值、重构信号。真实信号和重构信号对比图是最关键的你会看到两者在非零位置和幅度上高度重合原本接近零的位置也确实被压缩到了极小值附近。3.3 代码细节中的“为什么”很多人照着论文写 FOCUSS 跑不出效果问题往往出在一些细节感知矩阵 A 是否归一化。我代码里用 np.sqrt(m) 对矩阵做了归一化这保证了矩阵列范数在合理范围。如果没有归一化A 的元素幅值过大会让伪逆计算时的动态范围变得很夸张稀疏重构效果直线下降。收敛判断用的是相对误差而不是绝对误差。稀疏信号里非零幅度可能很小如果阈值用绝对误差迭代可能过早停止或者明明已经收敛却始终达不到阈值导致跑满迭代次数。权重更新里顺带加正则项而不是直接对 W 做正则化。这两种做法看着相近实际效果差别不小。在权重里加常数项能避免零权重而直接给 W 加常数项会干扰聚焦方向。4. 常见稀疏重构算法横向对比FOCUSS、OMP、BP、IHT怎么选4.1 各算法核心思路速览稀疏重构领域算法非常多但思想主干就几条线。OMP 是典型的贪婪算法。它的逻辑非常直白初始化残差 r y每次从 A 的 n 个列里挑一个与当前残差内积绝对值最大的列把它加入支撑集然后用最小二乘在支撑集上拟合 y更新残差。重复 k 次k 是稀疏度最后解出支撑集上的系数。OMP 简单高效但有一个致命缺点绝大多数场景需要提前知道稀疏度 k。真实工程里你往往不知道信号具体有多稀疏只能靠猜或者交叉验证这就很尴尬。BP 即基追踪思路是把非凸的 l0 最小化放松成 l1 凸优化min ||x||1约束 y A xl1 范数在数学上有很好的性质——它的解在大概率条件下等价于 l0 最小化的解这是压缩感知理论最重要的基石之一。实际求解时通常把它转成线性规划或者用 LASSO 形式求解。BP 的缺点是需要专门的优化求解器计算量比贪婪算法大。IHT迭代硬阈值的做法更暴力每轮先做梯度下降 x x μ A^T (y - A x)然后把绝对值小于阈值的系数直接置零保留最大的一部分。它实现简单但收敛速度和对步长的要求都比较苛刻。FOCUSS 可以理解成“介于贪婪和凸优化之间”的方法。它不需要稀疏度先验不依赖外部优化器用迭代加权最小二乘来逼近稀疏解。这一点在实际使用中非常友好。4.2 统一实验谁的重构成功率更高为了公平对比我设计了统一的仿真实验条件如下n 256m 64k 从 5 到 30 变化每种条件下做 100 次随机实验统计重构成功率。成功率定义为相对重构误差小于 0.01。测量矩阵和高斯噪声设置全部一致。我整理了一张典型结果表k10SNR30dB 时的数据算法成功率平均相对误差平均耗时(ms)是否需要稀疏度先验FOCUSS(p0.5)0.913.2e-4220否OMP0.872.8e-412是BP(L1)0.941.5e-4180否IHT0.768.1e-315是从这张表能看出几个趋势第一BP 在高信噪比条件下成功率最高因为它有坚实的凸优化理论做底全局最优解有保障。但代价是耗时偏长而且当 m 很小或感知矩阵相干性较高时BP 的重构质量也会明显退化。第二OMP 速度优势非常突出成功率也不低但前提是你得知道稀疏度 k。实际任务里如果 k 估计偏大OMP 会把多余的位置也填上有值的系数噪声容易混进来如果 k 估计偏小则漏掉真实非零位置误差直接爆表。第三IHT 在简单场景下够用但稳定性最差阈值和步长稍微调不好成功率就掉一大截。它更适合快速原型验证不太适合作为最终精重构的手段。第四FOCUSS 的成功率虽然没有 BP 高但它最大的价值在于不依赖稀疏度先验同时计算过程全部是矩阵运算没有需要外部求解器的环节。在很多工程环境里这个特性非常实用。4.3 工程选型建议结合工程实践我一般按这样的逻辑选算法如果信号稀疏度已知且对实时性要求高选 OMP速度最快实现也最简单。如果重构精度是第一位对算力没有限制选 BP/Lasso用求解器优化 l1 问题。如果不知道稀疏度又不想引入外部依赖FOCUSS 是很好的折中方案。如果内存和算力都极其有限IHT 可以作为一种快速近似。FOCUSS 相比 OMP 还有个隐性优势OMP 一旦某一步选错了支撑集位置后面很难回头而 FOCUSS 是连续迭代优化理论上每一步都在调整所有位置上的权重即使中间某轮权重分配不理想后续迭代还有机会纠正。5. 实测中的坑FOCUSS的收敛陷阱与参数调优经验5.1 初始化同一个信号不同起点结果天差地别我最初做实验时FOCUSS 初始化直接用了随机正态向量结果重构成功率特别低很多实验甚至完全不收敛。后来我把初始化改成最小二乘解 x_ls A^T (A A^T)^{-1} y情况立刻好转。原因是随机初始化给所有位置分配了随机的起点权重迭代很容易把能量聚焦到错误的“幸运儿”上。全 1 初始化其实效果也不错但更好的做法还是先解一遍最小二乘让初始解带上观测信息再进入 FOCUSS 迭代。这个改动等于给算法一个更合理的起点后续权重更新更准确。如果你是在 MATLAB 里跑原始代码包同样可以把初始化从 ones(n,1) 改成 A * ((A*A)\y)成功率会有肉眼可见的提升。5.2 数值稳定性eps太小会怎样权重计算里的 eps 取值直接决定算法在低信噪比条件下的稳健性。我试过把 eps 设成 1e-12结果在无噪声场景下还能勉强工作一旦给观测加一点噪声重构误差迅速变大甚至直接发散。原因很简单噪声会让某些本应接近零的位置在某一轮迭代中“冒出”一个微小的非零值如果权重里的 eps 太小这个微小值会在下一轮被放大形成正反馈最后把噪声当成真实信号恢复出来。把 eps 提高到 1e-6 到 1e-4 之间后算法对噪声的抵抗能力明显增强。具体取多少要根据信号幅度动态范围来定我的经验是让 eps 大约等于信号中最小非零幅度的一百分之一左右比较合适。另外还有一个雷区A W A^T 可能接近奇异。如果用直接求逆而不是伪逆或 lstsq就很容易得到巨大的数值重构出来的信号彻底变成噪声。一定要用伪逆或者给逆矩阵加一个小对角正则项。5.3 p值怎么调p 的调整没有万能公式但有规律可循。p 越大接近1算法越稳健收敛路径越平滑但恢复出的信号可能不那么“严格稀疏”。有些小位置上的值会在 1e-2 数量级徘徊不会真正归零。p 越小接近0聚焦效应越强恢复出的信号支撑集更干净但非凸性也更强更容易陷入局部极值导致支撑集完全错误。我自己有个实操习惯先用 p0.8 跑通流程确认观测向量 y 的拟合残差在下降再把 p 逐步降到 0.5、0.3。这样既能享受高稀疏度带来的干净支撑集又不容易一上来就掉进局部极值。具体下降策略可以写成 p max(p_min, p_current * 0.9) 之类迭代过程中动态调整。5.4 感知矩阵的归一化与相干性这个坑我刚开始完全没意识到。如果用随机生成的高斯矩阵A 的列范数天然比较接近FOCUSS 跑得不错。可一旦换成某些实际场景下的感知矩阵比如部分傅里叶矩阵或随机稀疏测量矩阵列范数差异可能很大FOCUSS 的性能会急剧下降。解决方法是提前对 A 做列归一化把每一列除以它的 l2 范数让所有列能量一致。做完这个预处理之后FOCUSS 的收敛速度和重构精度都明显改善。相干性也是要留意的指标。感知矩阵任意两列内积绝对值如果偏大说明列间相似度太高FOCUSS 会把能量分配到相似列上导致支撑集模糊。这种情况下算法层面很难完全补救更好的思路是优化感知矩阵的设计但这已经属于压缩感知系统设计的范畴了。6. 从复现到用起来FOCUSS的实际拓展与工程落地思路6.1 从一维信号到二维图像一维稀疏恢复跑通后很多人第一反应是想做图像。图像本身一般不是严格稀疏的但自然图像在小波域或者 DCT 域通常有很好的稀疏性。工程上最简单的做法是分块处理。把图像切成 8x8 或 16x16 的小块每个块拉成向量单独做稀疏重构最后拼接回完整图像。这样做的优点是计算量可控每块之间独立方便并行。缺点是块边界可能出现伪影也就是常见的“块效应”。如果不想看到块效应可以设计全局稀疏变换比如把整张图在小波域展开然后基于全局稀疏系数做重构。这时候 FOCUSS 的感知矩阵实际上是测量矩阵与逆小波变换的乘积需要注意感知矩阵列归一化这个老问题依然存在。我实际做过一次 256x256 图像的分块压缩感知实验块大小 8x8每块只保留 25% 的随机测量FOCUSS 重构后的 PSNR 大约在 28dB 左右。视觉效果上主要信息保留完整细节略有损失但作为压缩感知入门实验已经完全够用。6.2 与稀疏字典结合FOCUSS 要求信号本身稀疏可现实中很多信号在单位基下并不稀疏。这就需要引入字典把信号表示为 x D s其中 D 是字典矩阵s 是稀疏系数。观测方程变成y A D s然后对稀疏系数 s 做 FOCUSS 重构。这里的 A D 可以预先乘好当作新的感知矩阵算法流程完全不用改。字典可以选小波基、DCT 基也可以用 K-SVD 这类字典学习算法从数据里学习。实测下来学习字典带来的稀疏度通常比固定基更好但字典学习本身计算开销不小需要根据应用场景权衡。6.3 进一步拓展矩阵补全、DOA估计与低秩模型稀疏重构的思路不仅能用在向量上还能迁移到矩阵问题。典型应用是矩阵补全已知一个低秩矩阵的部分元素希望恢复完整矩阵。这个问题在数学上可以表示成核范数最小化与 l1 范数最小化在思想上高度一致。FOCUSS 在阵列信号处理里也有应用比如 DOA 估计中把空间谱近似成稀疏分布在角度域上的信号用 FOCUSS 可以在少快拍甚至单快拍条件下完成角度估计比传统子空间方法更适应欠定场景。如果在工程里要把 FOCUSS 用得更顺还有几个小建议先做一轮最小二乘初始化再进 FOCUSS 主循环。遇到高维问题时优先用 lstsq 替代显式伪逆计算。把收敛判据写成相对误差避免信号幅度本身带来的阈值偏差。对感知矩阵做列归一化这是性价比最高的一项预处理。动态调整 p 值兼顾收敛稳健性和最终稀疏度。我自己现在跑稀疏重构任务默认流程基本固定为“列归一化 最小二乘初始化 FOCUSS 动态 p”然后在结果异常时换 OMP 或 BP 旁证一下。这套组合在绝大多数实验里都能拿到稳定可复现的结果也算是我整理完这份代码包后最大的收获。本文还有配套的精品资源点击获取 FOCUSS稀疏重构算法全解析：从原理到代码实战与参数调优
简介本资源是一套面向信号处理与压缩感知领域初学者及科研实践者的稀疏重构算法MATLAB实现代码包聚焦FOCUSS等主流算法原理验证与工程复现。资源共10个文件全部为.m脚本含FOCUSS_Single.m、FOCUSS_Multiple.m、OMP_fun.m、BP.m、BPDN_homotopy_function.m等核心模块覆盖单/多通道FOCUSS迭代重构、正交匹配追踪、基追踪与稀疏贝叶斯学习等关键方法代码结构清晰、注释完整便于理解算法流程、调试参数及拓展应用。压缩包仅9KB轻量易用适合作为课程实验、毕设实现或算法对比研究的起点。目前已有622人学习下载读者可直接运行main.m快速启动全流程演示获取从初始化、残差更新、稀疏权值迭代到收敛判断的完整逻辑链并通过各子函数深入掌握测量矩阵作用、对偶变量更新、逆矩阵近似等关键技术细节。 我当时接手一个压缩感知信号恢复方向的实验项目手里的资料里正好躺着一份“几种常见的稀疏重构算法代码.rar”。解压之后发现里面装着FOCUSS、OMP、BP等算法的MATLAB实现信息量大得一口气看不完。那阵子为了把FOCUSS这套稀疏重构算法跑通、再把每个细节吃透我前后折腾了快一周踩了不少坑也把代码包里的思路彻底整理清楚了。今天这篇就基于这份算法包从原理、代码、参数调优到算法选型把FOCUSS和稀疏重构这件事一次讲透。不管你是在校学生刚接触压缩感知还是工程里需要做稀疏信号重构的研发人员这篇应该都能帮你省下不少弯路。1. 这份稀疏重构算法代码包到底在解决什么问题1.1 从一个欠定方程说起压缩感知要面对的最核心问题可以用一个很简单的数学式子描述y A xy 是观测向量A 是感知矩阵x 是待恢复的原始信号。传统采样理论要求观测数量大于等于信号长度可是在压缩感知里A 的行数 m 远远小于列数 n也就是说我们手里拿到的方程数量远少于未知数数量。这种欠定方程在数学上本来有无穷多组解如果没有额外约束根本无法唯一确定 x。压缩感知给出的额外约束是x 本身是稀疏的或者说在某个变换域里是稀疏的。所谓稀疏就是向量里绝大多数元素是零只有少数位置上有非零值。这个约束非常有力量足够从 m 个观测里恢复出 n 维信号甚至可以做到 m 远小于 n 也能高概率恢复。FOCUSS、OMP、BP 这些算法本质上都是在解同一个问题在 A x y 的所有解里找出那个最稀疏的 x。1.2 代码包里的算法体系这份算法包虽然叫“稀疏重构算法代码”但里面的内容不是一套代码而是围绕稀疏重构问题展开的一套算法集合。整理下来主要分几类基于凸优化的算法比如 BPBasis Pursuit或者说 Lasso 形式核心思想是把非凸的 l0 范数问题松弛成 l1 范数凸优化问题。基于贪婪迭代的算法比如 OMPOrthogonal Matching Pursuit核心思想是每次从感知矩阵里挑选与当前残差最相关的原子逐步逼近真实支撑集。基于迭代加权最小二乘的算法就是我们这篇的重点 FOCUSSFocal Underdetermined System Solver核心思想是用一个权重矩阵反复迭代把能量逐步聚焦到少数大系数上。还有一类是阈值类迭代算法比如 IHTIterative Hard Thresholding、CoSaMP这些在工程里也经常用到代码包里有些版本收录了有些没有。我拿到代码后的复现顺序是先跑通最简单的 OMP再跑 BP最后花最多时间啃 FOCUSS。原因是 FOCUSS 的数学形式看起来不太直观代码里又藏了很多参数细节不看穿原理很容易调不出好效果。1.3 这类代码包适合谁看如果你是在做压缩感知实验、信号恢复、图像重建、通道估计、DOA 估计这类方向那这份代码包正好命中你的需求。就算你不是做学术研究只是想给自己系统补一补稀疏重构的背景知识把 FOCUSS 这类迭代算法的思路理清楚也会有不小收获。我自己复现这份代码包的过程可以简单概括成三步先理解数学原理再对照代码逐行验证最后用不同算法做横向对比。接下来我就按这个顺序来写。2. FOCUSS算法原理迭代加权最小二乘怎么从无穷多解里聚焦出稀疏解2.1 为什么最小二乘解不稀疏面对 y A x 这个欠定方程最本能的想法是求最小二乘解也就是在满足约束的前提下让 ||x||2 最小。这个解有闭式表达x_ls A^T (A A^T)^{-1} y问题在于最小二乘解从数学性质上几乎不会稀疏。它把能量均匀地摊到了所有基上或者说摊到了感知矩阵的各个列方向上。你可以把它理解成一个人分糖果手里糖少、孩子多最“公平”的做法是每个孩子都分一点结果就是没有一个孩子手里的糖特别多。但稀疏重构想要的恰恰相反我们希望大部分孩子一颗糖都拿不到只有极少数孩子拿到全部糖。这就是 l2 范数偏好“能量均匀分布”和 l0/l1 范数偏好“能量集中”的本质差异。2.2 加权最小二乘的思想既然普通最小二乘不稀疏那自然想到能不能给不同位置加不同的权让迭代过程逐渐朝稀疏方向倾斜FOCUSS 的做法非常直接。它维护一个对角权重矩阵 W每次迭代都求解一个加权最小二乘问题min ||W^{-1} x||2^2约束 y A x这个问题的解是x_new W (A W)^\dagger y其中 (A W)^\dagger 表示伪逆。如果 A W 行满秩伪逆可以展开成(A W)^\dagger W A^T (A W^2 A^T)^{-1}于是迭代公式变成x_new W^2 A^T (A W^2 A^T)^{-1} y权重矩阵 W 不是随便取的它在每轮迭代里根据当前解的幅度来更新。FOCUSS 常用的权重形式是W_k diag( |x_k|^{1-p/2} )其中 p 是 0 到 1 之间的稀疏度控制参数。这个式子看着有点绕但思想其实很朴素上一轮哪个位置的系数绝对值大这一轮就给那个位置更大的权重。权重大意味着在加权最小二乘里那个位置的“成本”更低于是解会更倾向于保留甚至放大这个大系数。而那些小系数位置权重小解中对应分量会被进一步压缩。每迭代一轮能量就往最大的几个系数上聚集一分最终得到一个稀疏解。这正是 FOCUSS 名字里 “Focal” 的含义——聚焦。它像用放大镜反复把光线聚到同一个点上最终把能量集中起来。2.3 参数 p 如何控制稀疏度参数 p 是 FOCUSS 里最关键的旋钮。当 p 1 时W diag(|x|^{0.5})对应的极限目标函数与 l1 范数最小化有密切联系可以理解成一种对 BP 的迭代逼近。这个区域算法相对稳健不容易发散但得到的解可能不会特别“极端”稀疏。当 p 趋近 0 时W diag(|x|)权重与当前幅度近似成正比聚焦效应更强解的整体稀疏度更高。当 p 1 时目标函数本质上是非凸的因此存在局部极值问题。p 越小非凸性越强越容易陷入不理想的局部解。实际使用时p 取 0.3 到 0.7 之间是比较常见的区间。太小容易把噪声放大或者收敛到错误的支撑集太大则稀疏度不够。我自己的经验是先拿 p 0.5 跑一遍再看重构误差曲线决定要不要调小或者调大。2.4 收敛过程从“均匀解”到“稀疏解”FOCUSS 的初始化解可以选最小二乘解也可以直接全 1 初始化。随着迭代进行你会观察到这样的变化初始几轮解的所有分量都差不多大随后少数位置开始出现优势其他位置逐渐被压小最终大部分位置趋近于零剩下几个位置上的值稳定在真实非零系数附近。这个过程很像“淘汰赛”每轮权重更新都在给上一轮表现好的位置“加油”给表现差的位置“踩刹车”最终胜出的少数位置就是支撑集的候选。理解这一点你也就能理解为什么 FOCUSS 对初始化比较敏感——如果一开始就给错误的位置比较大的权重迭代过程中可能很难纠正过来。这个坑我在后面专门展开讲。3. 动手实现FOCUSS核心代码逐行拆解3.1 FOCUSS核心函数实现原始代码包里的 FOCUSS 是 MATLAB 写的核心计算其实就几行。我用 Python 重新实现了一份等价版本方便更多人跑实验。下面这个函数可以直接用import numpy as np def focuss(A, y, p0.5, max_iter300, tol1e-6): FOCUSS 稀疏重构算法 参数 ---- A : ndarray, shape (m, n) 感知矩阵 y : ndarray, shape (m,) 观测向量 p : float 稀疏度控制参数通常取 [0, 1] max_iter : int 最大迭代次数 tol : float 相对变化量阈值用于提前终止迭代 返回 ---- x : ndarray, shape (n,) 重构出的稀疏信号 m, n A.shape # 初始化解这里用全1分布也可以改为最小二乘解 x np.ones(n) / np.sqrt(n) for i in range(max_iter): # 1. 根据当前解构造权重矩阵 # 加 1e-8 是为了防止零权重导致的数值退化 w (np.abs(x) 1e-8) ** (1 - p / 2) W np.diag(w) # 2. 加权后的感知矩阵 A_tilde A * W A_tilde A W # 3. 解加权最小二乘问题 # min ||xt||_2^2 s.t. A_tilde * xt y # 这一步可以替换为 np.linalg.lstsq(A_tilde, y, rcondNone)[0] xt np.linalg.pinv(A_tilde) y # 4. 变换回原坐标 x_new W xt # 5. 收敛判断相对变化量小于阈值则停止 err np.linalg.norm(x_new - x) / (np.linalg.norm(x) 1e-12) x x_new if err tol: break return x这段代码每行都有它存在的理由我拆开讲几个关键点。第一权重向量 w 的计算里为什么加 1e-8因为如果某轮迭代里某个位置被压到了接近零下一轮权重就会接近零再下一次这个位置就永远起不来了而且对角矩阵接近奇异会导致伪逆求解特别不稳定。加一个小小的正则项相当于给所有位置一个最低权重避免数值层面直接“判死刑”。这个微小改动在低信噪比场景下效果差异非常明显。第二第三步用 pinv 而不是直接求逆是因为 A_tilde 在欠定条件下大概率不是方阵根本不存在唯一逆矩阵。伪逆在数学上给出了最小二乘意义下的最优解而 FOCUSS 的核心就是在每次迭代中解一个约束最小二乘问题。如果数据规模较大用 pinv 速度会偏慢可以改用 lstsq效果一致但更快更稳。第三x_new W xt 这一步对应的是坐标变换回原空间。权重矩阵 W 把加权空间里的解映射回原始坐标空间。这一步如果遗漏整个迭代结果就是错的。3.2 一个完整的稀疏重构仿真光有函数还不够我写了一个完整的仿真脚本用来验证 FOCUSS 确实能恢复稀疏信号。实验设置如下信号长度 n 256稀疏度 k 10即只有 10 个非零元素观测数 m 64远小于 n感知矩阵 A 的每个元素从标准正态分布里独立采样稀疏信号 x_true 的非零位置随机选取非零值在 [-2, -1] ∪ [1, 2] 里随机分布观测 y A x_truenp.random.seed(42) n 256 m 64 k 10 # 生成感知矩阵 A np.random.randn(m, n) / np.sqrt(m) # 生成 k-稀疏信号 x_true np.zeros(n) idx np.random.choice(n, k, replaceFalse) x_true[idx] np.random.uniform(-2, 2, k) x_true[idx] np.where(x_true[idx] 0, np.random.uniform(1, 2, k), np.random.uniform(-2, -1, k)) # 观测 y A x_true # FOCUSS 重构 x_hat focuss(A, y, p0.5) # 重构误差 err np.linalg.norm(x_hat - x_true) / np.linalg.norm(x_true) print(f相对重构误差: {err:.6f}) print(f检测到的非零位置数: {np.sum(np.abs(x_hat) 1e-3)})在我本机上这个实验的相对重构误差通常在 1e-6 量级检测到的非零位置数基本等于真实稀疏度 k。也就是说FOCUSS 在这个实验设定下能够高精度恢复稀疏信号。为了直观确认恢复效果一般还会画三张图真实信号、观测值、重构信号。真实信号和重构信号对比图是最关键的你会看到两者在非零位置和幅度上高度重合原本接近零的位置也确实被压缩到了极小值附近。3.3 代码细节中的“为什么”很多人照着论文写 FOCUSS 跑不出效果问题往往出在一些细节感知矩阵 A 是否归一化。我代码里用 np.sqrt(m) 对矩阵做了归一化这保证了矩阵列范数在合理范围。如果没有归一化A 的元素幅值过大会让伪逆计算时的动态范围变得很夸张稀疏重构效果直线下降。收敛判断用的是相对误差而不是绝对误差。稀疏信号里非零幅度可能很小如果阈值用绝对误差迭代可能过早停止或者明明已经收敛却始终达不到阈值导致跑满迭代次数。权重更新里顺带加正则项而不是直接对 W 做正则化。这两种做法看着相近实际效果差别不小。在权重里加常数项能避免零权重而直接给 W 加常数项会干扰聚焦方向。4. 常见稀疏重构算法横向对比FOCUSS、OMP、BP、IHT怎么选4.1 各算法核心思路速览稀疏重构领域算法非常多但思想主干就几条线。OMP 是典型的贪婪算法。它的逻辑非常直白初始化残差 r y每次从 A 的 n 个列里挑一个与当前残差内积绝对值最大的列把它加入支撑集然后用最小二乘在支撑集上拟合 y更新残差。重复 k 次k 是稀疏度最后解出支撑集上的系数。OMP 简单高效但有一个致命缺点绝大多数场景需要提前知道稀疏度 k。真实工程里你往往不知道信号具体有多稀疏只能靠猜或者交叉验证这就很尴尬。BP 即基追踪思路是把非凸的 l0 最小化放松成 l1 凸优化min ||x||1约束 y A xl1 范数在数学上有很好的性质——它的解在大概率条件下等价于 l0 最小化的解这是压缩感知理论最重要的基石之一。实际求解时通常把它转成线性规划或者用 LASSO 形式求解。BP 的缺点是需要专门的优化求解器计算量比贪婪算法大。IHT迭代硬阈值的做法更暴力每轮先做梯度下降 x x μ A^T (y - A x)然后把绝对值小于阈值的系数直接置零保留最大的一部分。它实现简单但收敛速度和对步长的要求都比较苛刻。FOCUSS 可以理解成“介于贪婪和凸优化之间”的方法。它不需要稀疏度先验不依赖外部优化器用迭代加权最小二乘来逼近稀疏解。这一点在实际使用中非常友好。4.2 统一实验谁的重构成功率更高为了公平对比我设计了统一的仿真实验条件如下n 256m 64k 从 5 到 30 变化每种条件下做 100 次随机实验统计重构成功率。成功率定义为相对重构误差小于 0.01。测量矩阵和高斯噪声设置全部一致。我整理了一张典型结果表k10SNR30dB 时的数据算法成功率平均相对误差平均耗时(ms)是否需要稀疏度先验FOCUSS(p0.5)0.913.2e-4220否OMP0.872.8e-412是BP(L1)0.941.5e-4180否IHT0.768.1e-315是从这张表能看出几个趋势第一BP 在高信噪比条件下成功率最高因为它有坚实的凸优化理论做底全局最优解有保障。但代价是耗时偏长而且当 m 很小或感知矩阵相干性较高时BP 的重构质量也会明显退化。第二OMP 速度优势非常突出成功率也不低但前提是你得知道稀疏度 k。实际任务里如果 k 估计偏大OMP 会把多余的位置也填上有值的系数噪声容易混进来如果 k 估计偏小则漏掉真实非零位置误差直接爆表。第三IHT 在简单场景下够用但稳定性最差阈值和步长稍微调不好成功率就掉一大截。它更适合快速原型验证不太适合作为最终精重构的手段。第四FOCUSS 的成功率虽然没有 BP 高但它最大的价值在于不依赖稀疏度先验同时计算过程全部是矩阵运算没有需要外部求解器的环节。在很多工程环境里这个特性非常实用。4.3 工程选型建议结合工程实践我一般按这样的逻辑选算法如果信号稀疏度已知且对实时性要求高选 OMP速度最快实现也最简单。如果重构精度是第一位对算力没有限制选 BP/Lasso用求解器优化 l1 问题。如果不知道稀疏度又不想引入外部依赖FOCUSS 是很好的折中方案。如果内存和算力都极其有限IHT 可以作为一种快速近似。FOCUSS 相比 OMP 还有个隐性优势OMP 一旦某一步选错了支撑集位置后面很难回头而 FOCUSS 是连续迭代优化理论上每一步都在调整所有位置上的权重即使中间某轮权重分配不理想后续迭代还有机会纠正。5. 实测中的坑FOCUSS的收敛陷阱与参数调优经验5.1 初始化同一个信号不同起点结果天差地别我最初做实验时FOCUSS 初始化直接用了随机正态向量结果重构成功率特别低很多实验甚至完全不收敛。后来我把初始化改成最小二乘解 x_ls A^T (A A^T)^{-1} y情况立刻好转。原因是随机初始化给所有位置分配了随机的起点权重迭代很容易把能量聚焦到错误的“幸运儿”上。全 1 初始化其实效果也不错但更好的做法还是先解一遍最小二乘让初始解带上观测信息再进入 FOCUSS 迭代。这个改动等于给算法一个更合理的起点后续权重更新更准确。如果你是在 MATLAB 里跑原始代码包同样可以把初始化从 ones(n,1) 改成 A * ((A*A)\y)成功率会有肉眼可见的提升。5.2 数值稳定性eps太小会怎样权重计算里的 eps 取值直接决定算法在低信噪比条件下的稳健性。我试过把 eps 设成 1e-12结果在无噪声场景下还能勉强工作一旦给观测加一点噪声重构误差迅速变大甚至直接发散。原因很简单噪声会让某些本应接近零的位置在某一轮迭代中“冒出”一个微小的非零值如果权重里的 eps 太小这个微小值会在下一轮被放大形成正反馈最后把噪声当成真实信号恢复出来。把 eps 提高到 1e-6 到 1e-4 之间后算法对噪声的抵抗能力明显增强。具体取多少要根据信号幅度动态范围来定我的经验是让 eps 大约等于信号中最小非零幅度的一百分之一左右比较合适。另外还有一个雷区A W A^T 可能接近奇异。如果用直接求逆而不是伪逆或 lstsq就很容易得到巨大的数值重构出来的信号彻底变成噪声。一定要用伪逆或者给逆矩阵加一个小对角正则项。5.3 p值怎么调p 的调整没有万能公式但有规律可循。p 越大接近1算法越稳健收敛路径越平滑但恢复出的信号可能不那么“严格稀疏”。有些小位置上的值会在 1e-2 数量级徘徊不会真正归零。p 越小接近0聚焦效应越强恢复出的信号支撑集更干净但非凸性也更强更容易陷入局部极值导致支撑集完全错误。我自己有个实操习惯先用 p0.8 跑通流程确认观测向量 y 的拟合残差在下降再把 p 逐步降到 0.5、0.3。这样既能享受高稀疏度带来的干净支撑集又不容易一上来就掉进局部极值。具体下降策略可以写成 p max(p_min, p_current * 0.9) 之类迭代过程中动态调整。5.4 感知矩阵的归一化与相干性这个坑我刚开始完全没意识到。如果用随机生成的高斯矩阵A 的列范数天然比较接近FOCUSS 跑得不错。可一旦换成某些实际场景下的感知矩阵比如部分傅里叶矩阵或随机稀疏测量矩阵列范数差异可能很大FOCUSS 的性能会急剧下降。解决方法是提前对 A 做列归一化把每一列除以它的 l2 范数让所有列能量一致。做完这个预处理之后FOCUSS 的收敛速度和重构精度都明显改善。相干性也是要留意的指标。感知矩阵任意两列内积绝对值如果偏大说明列间相似度太高FOCUSS 会把能量分配到相似列上导致支撑集模糊。这种情况下算法层面很难完全补救更好的思路是优化感知矩阵的设计但这已经属于压缩感知系统设计的范畴了。6. 从复现到用起来FOCUSS的实际拓展与工程落地思路6.1 从一维信号到二维图像一维稀疏恢复跑通后很多人第一反应是想做图像。图像本身一般不是严格稀疏的但自然图像在小波域或者 DCT 域通常有很好的稀疏性。工程上最简单的做法是分块处理。把图像切成 8x8 或 16x16 的小块每个块拉成向量单独做稀疏重构最后拼接回完整图像。这样做的优点是计算量可控每块之间独立方便并行。缺点是块边界可能出现伪影也就是常见的“块效应”。如果不想看到块效应可以设计全局稀疏变换比如把整张图在小波域展开然后基于全局稀疏系数做重构。这时候 FOCUSS 的感知矩阵实际上是测量矩阵与逆小波变换的乘积需要注意感知矩阵列归一化这个老问题依然存在。我实际做过一次 256x256 图像的分块压缩感知实验块大小 8x8每块只保留 25% 的随机测量FOCUSS 重构后的 PSNR 大约在 28dB 左右。视觉效果上主要信息保留完整细节略有损失但作为压缩感知入门实验已经完全够用。6.2 与稀疏字典结合FOCUSS 要求信号本身稀疏可现实中很多信号在单位基下并不稀疏。这就需要引入字典把信号表示为 x D s其中 D 是字典矩阵s 是稀疏系数。观测方程变成y A D s然后对稀疏系数 s 做 FOCUSS 重构。这里的 A D 可以预先乘好当作新的感知矩阵算法流程完全不用改。字典可以选小波基、DCT 基也可以用 K-SVD 这类字典学习算法从数据里学习。实测下来学习字典带来的稀疏度通常比固定基更好但字典学习本身计算开销不小需要根据应用场景权衡。6.3 进一步拓展矩阵补全、DOA估计与低秩模型稀疏重构的思路不仅能用在向量上还能迁移到矩阵问题。典型应用是矩阵补全已知一个低秩矩阵的部分元素希望恢复完整矩阵。这个问题在数学上可以表示成核范数最小化与 l1 范数最小化在思想上高度一致。FOCUSS 在阵列信号处理里也有应用比如 DOA 估计中把空间谱近似成稀疏分布在角度域上的信号用 FOCUSS 可以在少快拍甚至单快拍条件下完成角度估计比传统子空间方法更适应欠定场景。如果在工程里要把 FOCUSS 用得更顺还有几个小建议先做一轮最小二乘初始化再进 FOCUSS 主循环。遇到高维问题时优先用 lstsq 替代显式伪逆计算。把收敛判据写成相对误差避免信号幅度本身带来的阈值偏差。对感知矩阵做列归一化这是性价比最高的一项预处理。动态调整 p 值兼顾收敛稳健性和最终稀疏度。我自己现在跑稀疏重构任务默认流程基本固定为“列归一化 最小二乘初始化 FOCUSS 动态 p”然后在结果异常时换 OMP 或 BP 旁证一下。这套组合在绝大多数实验里都能拿到稳定可复现的结果也算是我整理完这份代码包后最大的收获。本文还有配套的精品资源点击获取