简介这份Matlab单像素成像算法比较代码包面向已掌握基本成像原理、希望系统对比不同重建方法的读者围绕单像素成像中的稀疏采样与图像重建问题提供DGI、梯度下降、共轭梯度下降、Poisson最大似然、交替投影、稀疏表示压缩感知及全变分压缩感知等七类方法的可运行示例。包内共10个文件以9个.m脚本和1个txt说明为主脚本分别对应各算法实现及误差计算、演示入口说明文件用于辅助理解调用关系压缩包约12KB体量轻便便于快速查阅与二次修改。目前已有2762人学习适合课程实验、算法复现和不同重建策略的初步对照。读者可借此理解各方法在采样率、重建质量与计算开销上的差异并基于统一误差函数开展结果比较为单像素成像实验选型与后续改进提供参考。1. 一个探测器的悖论为什么几百次测量就能还原一张图一份典型的 Matlab 单像素成像脚本摆在面前图像 64×64 共 4096 个像素测量次数只给 800 次重构那一段调的是共轭梯度下降。按线性代数直觉这是严重欠定的方程组解不唯一按压缩感知的说法这叫欠采样重构听起来像玄学。两个东西拼在一起居然能出图这件事本身值得停一秒想清楚到底是哪一方在补偿缺失的那些方程如果以为共轭梯度下降天生就能做稀疏重构那后面对参数的所有直觉都会跑偏。把这条链路拆开看过一次之后mask 怎么选、正则项怎么调、为什么有人跑 20 次迭代就收敛而有人跑到 500 次还在震荡就都有解释。2. 从光强到方程单像素成像的测量矩阵是怎么长出来的2.1 掩模调制与标量测量一步物理怎么变成一行矩阵单像素成像的硬件链路其实短得可怜光源照到场景上反射光经过一组掩模图案调制再由一个聚光镜把所有能量收进单像素探测器。探测器没有空间分辨能力它读到的只是一个数——这一帧的总光强。想要二维图像就必须在时间维度上换掩模换一次图案读一个数换 M 次得到 M 个数。关键在于每次读到的这个数本质上是场景光强分布与当前掩模的内积。把场景离散成向量 x长度 n N×N第 i 次测量用的掩模离散成向量 φ_i同样长度 n那么探测器读数 y_i ≈ Σ_j φ_ij · x_j。M 次测量堆起来就是 y ΦxΦ 是 M×n 的测量矩阵每一行对应一次掩模曝光。一段最朴素的 Matlab 正向仿真长这样N 64; % 图像边长 n N * N; % 未知量个数 4096 img double(imread(target.png)); x img(:); % n x 1 列向量 M round(0.25 * n); % 采样率 25%即 1024 次测量 rng(0); % 固定随机源否则每次跑结果都不同 Phi randn(M, n) / sqrt(M); % 高斯随机测量矩阵按行归一 y Phi * x; % 探测器读到的 M 个标量这段代码里有两处容易被忽略但影响极大。第一处是/ sqrt(M)把每一行的能量归一到 1是为了让 y 的量级与 x 处在同一水平否则后面正则项系数 λ 的量纲会跟着采样率一起漂移换个 M 就得重新调参。第二处是rng(0)随机掩模一旦不固定种子同一份代码两次运行的重建质量能差出好几个 dB调试时会误以为是算法本身不稳定。这条链路还有个物理上的隐含条件常被跳过掩模刷新与探测器采样必须同步。真实系统里 DMD 的刷新有建立时间探测器积分窗如果没对齐读到的就是两张掩模的混合结果等效于 Φ 的每一行都被悄悄污染了一点点。仿真里看不出来一上硬件就是「图总是糊一层」的源头。2.2 稀疏基与不相干性为什么随机掩模比哈达玛掩模更容易命中有了 y Φx 还不够M n 时它有无数个解。压缩感知的解法是把 x 换到另一个域里表示x Ψθ其中 Ψ 是稀疏基DCT、小波、或者图像本身θ 只有 K 个显著非零元K n。代进去得到 y ΦΨθ AθA 是 M×n 的传感矩阵。这里有个非常容易被跳过的约束Φ 和 Ψ 必须「不相干」。相干性的定义是把两个基里最相似的一对原子的内积挑出来再乘上 √n 归一。相干性越高一个基里的稀疏信号在另一个基里就越「铺得开」压缩感知的恢复保证就越弱。这条约束直接决定了掩模怎么选。如果稀疏基用的是哈达玛变换而测量掩模也用哈达玛图案两个基完全重合相干性等于最大值这时候压缩感知的理论保证直接失效——你等于在稀疏域里逐点采样只是自己没意识到而已。换成伯努利二值掩模或者高斯随机掩模与 DCT、小波这些常见基的相干性都很低恢复才稳。血泪经验是很多人第一次做单像素成像会选哈达玛掩模因为它二值、好实现、有快速变换、重建还能用逆变换一步到位。它在高采样率下效果确实漂亮但一旦把采样率压到 10% 以下重构质量的下降速度明显比随机掩模陡。原因不神秘哈达玛掩模压到低采样率时是「均匀地少采」而随机掩模是「不相干地乱采」后者的信息冗余更低。2.3 从 L0 到 L1凸松弛到底换来了什么理论上最想要的是 L0 范数最小在所有满足 y Aθ 的解里挑非零元最少的那个。问题是这属于组合优化n 上千时穷举不现实。L1 范数是 L0 在单位球上的凸包把它替换进去问题变成可解的凸优化min ||θ||₁ s.t. y Aθ几何上的直觉是L1 球在坐标轴方向有尖角当残差约束的可行域高维空间里的一个仿射子空间与这个球相切时切点大概率落在坐标轴上也就是稀疏散点上而 L2 球是光滑的切点会落在球面上解里几乎每个分量都不为零图像就糊了。这就是为什么纯最小二乘解出来的图总是「整体形状对、细节全丢」。实际有噪声时等式约束要放松成不等式或惩罚形式min ½||y − Aθ||₂² λ||θ||₁到这里才出现了这份代码里真正要解的东西。λ 大稀疏性强、解偏平滑、细节丢λ 小拟合狠、容易把噪声也当信号拟合进去。这个式子不是二次的共轭梯度下降没法直接吃它后面第 3 章要拆的就是这层关系。2.4 噪声进来之后理想 RIP 在实验台上永远不会成立教科书上的恢复保证依赖 RIP有限等距性质对任意 K 稀疏向量 θ要求 (1−δ)||θ||² ≤ ||Aθ||² ≤ (1δ)||θ||²δ 足够小。这个条件在仿真里勉强成立在真实装置上几乎不可能满足。干扰源至少有四类。探测器有暗电流和读出噪声低光强下还叠加泊松噪声表现为 y 的方差随信号强度变化不是简单的高斯白噪声。掩模与场景的相对位置有微小位移等效于 A 的每一行都在抖。光源强度随时间漂移等效于 y 整体乘了个慢变系数。掩模本身的二值化不理想灰度掩模的透过率不是严格的 0 和 1。这四类里最坑的是光源漂移因为它不改变图像形状只是让 y 的整体尺度慢慢变重构出来的图看起来「形状对、对比度不对」很容易被误判成 λ 没调好。常见的处理思路是在每一批测量里插一帧全通掩模所有元素为 1用它的读数做归一化基准把慢变系数除掉。这招在论文里不显眼在实验台上能省掉大量无效调参。3. 共轭梯度下降它解的到底是不是压缩感知那个问题3.1 Krylov 子空间与 A-共轭CG 比最速下降快在哪共轭梯度下降解决的是这类问题A 对称正定求 Ax b。它的迭代不是在原始空间里沿负梯度走而是在一个逐步扩张的子空间里找最优解这个子空间叫 Krylov 子空间 K_k(A, r₀) span{r₀, Ar₀, A²r₀, …}。每一步都把当前残差正交化到之前所有搜索方向上于是搜索方向两两 A-共轭p_iᵀ A p_j 0i ≠ j。这个性质带来的效果非常直观。最速下降每步都沿当前负梯度走相邻两步的方向几乎正交路径是锯齿状的条件数大时收敛极慢。CG 因为方向互相 A-共轭不会重复走已经走过的方向理论上在精确算术下 n 步内必收敛——尽管浮点运算下远达不到但实际往往几十步就够了。收敛速度有个常用估计误差的 A-范数每步至少乘上 (√κ − 1)/(√κ 1)κ 是 A 的条件数。这句话信息量很大条件数 100 时这个因子约 0.82收敛慢条件数 10 时约 0.52几十步就够。所以对 CG 来说比「迭代多少次」更重要的是「矩阵条件数是多少」。这也解释了后面为什么一定要加正则项。3.2 正规方程、Tikhonov 与 CGLS三种写法的数值差异最小二乘问题 min ||Aθ − y||² 的最优性条件是 AᵀAθ Aᵀy这就是正规方程。AᵀA 对称半正定可以用 CG 解。但这里埋了个雷AᵀA 的条件数是 A 的条件数的平方。A 本身条件数 50 还算温和平方之后 2500CG 的收敛因子就变成 0.96 左右跑几百步都未必收敛。第二种写法是 Tikhonov 正则把 (AᵀA λI)θ Aᵀy 拿来做。加 λI 相当于给所有特征值抬高一个地板条件数从 κ² 降到 (σ_max² λ)/(σ_min² λ)小奇异值被压住了收敛立刻变快。代价是解被往零的方向拉了一点图像会略偏平滑——但这恰好和压缩感知要的「抑制噪声」方向一致很多时候不是坏事。第三种是 CGLS它不显式构造 AᵀA而是直接在 A 上做 Krylov 迭代每一步只做一次 A 乘和一次 Aᵀ 乘。数值上比正规方程稳得多条件数不会平方放大。代价是实现稍复杂需要维护两套递推量。写法求解对象有效条件数每步矩阵乘法适用场景正规方程 CGAᵀAθ Aᵀyκ²2 次小规模、A 良态Tikhonov CG(AᵀAλI)θ Aᵀy(σ²λ) 比值2 次有噪声、需要正则CGLS / LSQR直接在 A 上迭代κ2 次大规模、病态有一处容易翻车在 Matlab 里显式写A * A会真的构造出一个 n×n 矩阵。n 4096 时是 128 MB勉强能跑n 256×256 65536 时是 34 GB直接内存爆炸。这就是下一节要讲的函数句柄写法为什么必要。3.3 手写 CG 与 pcg从显式矩阵到函数句柄先看手写版本它把每一步在干什么摊开调试时能直接看到残差曲线function [x, relres, k] cg_solve(Afun, b, x, tol, maxit) % Afun : 对称正定算子句柄输入列向量返回 A*v % b : 右端项 (n x 1) % x : 初值同时作为输出 % tol : 相对残差停止阈值常用 1e-6 % maxit : 最大迭代步数 r b - Afun(x); p r; rs r * r; bn norm(b); for k 1:maxit Ap Afun(p); pAp p * Ap; if pAp 0 warning(算子非正定pAp %g检查正则项, pAp); break; end alpha rs / pAp; % 沿 p 方向走多远 x x alpha * p; r r - alpha * Ap; rs_new r * r; relres sqrt(rs_new) / bn; if relres tol, break; end p r (rs_new / rs) * p; % 新方向与旧方向 A-共轭 rs rs_new; end end这段里pAp 0那个检查不是摆设。压缩感知里如果 λ 设成 0 而 A 严重亏秩AᵀA 是半正定的pAp 可能接近 0alpha 会炸成天文数字图直接变成一片雪花。加一句 warning 至少能让你知道是自己参数设错了而不是算法烂。实际工程里我更倾向直接用 Matlab 自带的pcg把算子写成匿名函数永远不显式构造 AᵀAAfun (v) Phi * (Phi * v) lambda * v; % Tikhonov 算子 rhs Phi * y; [theta, flag, relres, iter] pcg(Afun, rhs, 1e-6, 200);flag的取值含义要记住0 是收敛1 是达到最大迭代仍没收敛2 是预条件子病态3 是迭代停滞4 是某步出现非有限值。调试时如果看到 flag 1别急着加迭代次数上限先看条件数——十有八九是 λ 太小或者干脆没加。看到 flag 4多半是 Phi 里有 Inf 或者 NaN往上游查数据归一化。3.4 把 CG 塞进 L1IRLS 外层循环怎么搭到这一步必须把话说透CG 解的是二次型问题而压缩感知的目标函数带 L1 项不是二次的。CG 不能直接解压缩感知它只能当内层的一个线性求解器。这一层关系如果搞混调参会调得非常迷茫。把 L1 和 CG 接起来的经典做法是迭代重加权最小二乘IRLS。核心思想是|θ_i| ≈ θ_i² / |θ_i|用当前估计的 |θ_i| 当权重把 L1 问题转成一串加权 L2 问题每个 L2 子问题用 CG 解。lambda 1e-3; % 正则强度越大越稀疏 eps_reg 1e-6; % 防止 1/0常用 1e-6 ~ 1e-4 theta ones(n, 1); % 初值不要全零否则第一步权重会爆 Aop (v) Phi * (Phi * v); rhs Phi * y; for outer 1:15 w 1 ./ (abs(theta) eps_reg); % 越接近 0 的解权重越大 Rop (v) lambda * (w .* v); % 加权 L2 正则项 theta pcg((v) Aop(v) Rop(v), rhs, 1e-6, 200, [], [], theta); fprintf(IRLS %2d: 非零元 %d\n, outer, nnz(abs(theta) 1e-3)); end x_rec Psi * theta; % 回到图像域三处参数值得单独说。eps_reg决定了「多小算零」设太大会把真实的小系数一起压掉图像细节丢失设太小会让权重矩阵条件数爆表内层 CG 收敛极慢。经验值 1e-6 起步如果内层 flag 经常是 1就往 1e-4 调。lambda控制整体稀疏强度量级通常要和norm(Phi*y, inf)对齐看单纯照抄别人的 1e-3 在你自己数据上很可能失效。外层循环次数 10~20 次足够再往上收益递减可以直接看非零元个数有没有稳定。注意这套 IRLS 只是众多 L1 求解路线里的一种。ADMM 把它拆成「数据保真 软阈值」两个子问题交替解FISTA 用梯度步加软阈值再加动量加速TV 正则则要在 L1 之外再套一层梯度算子。选哪条路线取决于你的先验是「稀疏」还是「分片光滑」下一章对比的就是这件事。4. 四种重构器同场对比采样率降到多少开始崩4.1 统一实验骨架掩模、采样率与评价指标对比算法的第一原则是所有变量只动一个。固定同一张图、同一个随机种子、同一个 Φ只换重构器否则你根本分不清差异来自算法还是来自掩模。骨架大概是这样rng(2024); img double(imread(target.png)); img imresize(img, [64 64]); x img(:); n numel(x); rates [0.05 0.10 0.20 0.35 0.50]; for r rates M round(r * n); Phi randn(M, n) / sqrt(M); y Phi * x; % 在这里依次调用 CG-Tikhonov / OMP / FISTA / TV-ADMM % 每个算法跑完记录 psnr(x, x_rec)、ssim 和 toc end指标选 PSNR 和 SSIM 两个就够。PSNR 对整体亮度误差敏感SSIM 对结构相似度敏感两者经常给出不同排序——一个算法 PSNR 高但 SSIM 低通常意味着它把噪声抹平了但边缘也磨掉了。只看一个指标下结论容易得出偏颇的结论。还有一条必须统一每个算法的 λ 都要在自己的参数网格上取最优再拿来对比。拿一个没调参的 FISTA 去比一个调好的 CG这种对比没有意义。代价是实验量翻好几倍但这是唯一能让结论站住的做法。4.2 重建质量随采样率的变化谁先崩趋势大致是这样的。在 50% 采样率附近四种方法的 PSNR 差距不明显往往在 1~2 dB 以内因为信息本来就够多什么方法都能凑出来。差距是从 20% 往下拉开的。纯 CG 加 Tikhonov 属于 L2 正则它给出的解在所有方向上都平滑本质是「最小能量解」稀疏性完全没有被利用。采样率降到 20% 以下重建图会明显偏糊边缘发软纹理区域糊成一片。它的优势是快——迭代次数少每步只做两次矩阵乘法。OMP 走的是另一条路贪心地每次挑一个与当前残差最相关的原子然后在这一小撮原子上做最小二乘。它在极度稀疏K 很小且 M 不太小时效果很好但有两个软肋。一是需要预先知道稀疏度 K猜大了会过拟合噪声猜小了会丢结构二是每一步都要算 Aᵀr复杂度随原子数线性涨K 大时非常慢。FISTA 是 L1 正则的主力选手梯度步加软阈值再叠动量收敛速率比 ISTA 快一个量级。它在 10%~30% 采样率区间表现稳定是中低采样率下的性价比之选。软阈值的副作用是会把小系数整体往零拉图像会有轻微的「块状感」这是 L1 的固有偏置不是实现 bug。TV 正则ADMM 或增广拉格朗日求解在低采样率下通常画质最好因为图像的分片光滑先验比纯稀疏先验更贴合自然图像。代价是慢而且慢得不只是常数倍——它每次迭代内部还要再解一个子问题整体比 FISTA 慢一个数量级很常见。提示这里的排序是趋势层面的具体到某一张图、某一组掩模名次可能互换。尤其是纹理极丰富的图像TV 的分片光滑假设会失效反而不如 FISTA。4.3 运行时间与内存TV 类方法慢在哪时间差异的量级值得单独拎出来看因为它直接决定了你走不走得通高分辨率。方法先验类型每次迭代代价典型迭代数主要失败模式CG TikhonovL2 平滑2 次矩阵乘20~80低采样率下糊OMP严格稀疏1 次相关 1 次 LSK 步K 估计错就崩FISTAL1 稀疏2 次矩阵乘 软阈值100~300块状伪影TV ADMM分片光滑内层子问题求解100~500极慢、参数敏感从表里能读出一个关键取舍CG 快但先验弱TV 先验强但慢。FISTA 卡在中间这也是它在实际项目里出场率最高的原因。内存方面只要坚持用函数句柄四种方法的峰值内存基本都在「Phi 本身 几个 n 维向量」这个量级。真正的内存杀手是显式构造 AᵀA 或者把 TV 的梯度算子展开成稠密矩阵。n 65536 时一个稠密 n×n 就是 34 GB这个坑踩一次能记一辈子。ADMM 实现里有一步常见的错误是「为了写起来方便」把差分算子写成diff(eye(n))那一行代码就是内存爆炸的直接原因正确写法是把它写成只做邻域差分的函数句柄。5. 排查图糊、图碎、不收敛时先看哪几个量5.1 图像糊、碎、全黑先查信号链再查算法现象一重建图整体形状对但像蒙了一层雾边缘发软。原因大概率不是算法而是测量矩阵没有按行归一化。如果 Φ 的每一行能量不同y 的各个分量量级就会参差不齐最小二乘会自动更重视量级大的那些测量等价于给不同掩模加了不同权重。解决很直接构造完 Phi 之后强制Phi Phi ./ vecnorm(Phi, 2, 2)再检查vecnorm(Phi,2,2)是不是全 1。现象二图碎成马赛克状或者有明显的方块拼接感。这通常是稀疏基和掩模不匹配。用哈达玛掩模配小波基或者用 DCT 基去配本身就带 DCT 结构的掩模相干性过高恢复出来的解在基里分布得很散重建图上就表现为均匀的块状噪声。解决是换随机二值或高斯掩模重跑一次看看块状感有没有消失。现象三重建全黑或者全白。这种八成是 λ 和数据量级严重失配。检查一下norm(Phi*y)和lambda的数量级如果差了三四个数量级λ 就相当于无穷大解直接被压到零。另一个可能是 y 里混进了 NaN往上顺一遍数据读入有没有除零。现象四低采样率下某个固定区域永远重建不出来。这往往是掩模覆盖不均。随机掩模虽然整体不相干但 M 很小时局部可能恰好没有充分覆盖某些位置尤其图像边角。解决是换一组种子多跑几次取平均或者用结构化随机把区域分块后每块内独立随机保证空间覆盖均匀。5.2 不收敛、内存爆、每次结果不同查参数与随机源现象五pcg 返回 flag 1迭代跑满还没收敛。先看两个量condest估一下有效算子的条件数以及 λ 到底设了多少。不加正则或 λ 极小的时候AᵀA 的最小奇异值接近零CG 收敛因子接近 1当然跑不满。解决办法是加 λ从 1e-3 起试或者换 CGLS/LSQR 这类不平方条件数的方法。现象六跑到一半内存爆掉。基本就两种可能显式构造了 AᵀA或者 TV 的差分算子被展开成了稠密矩阵。前者改成函数句柄后者改成邻域差分。还有个隐蔽的来源是spdiags用的位置不对——IRLS 里的权重矩阵 W 一定要用spdiags生成稀疏对角阵写成diag(w)就是又一个 n×n 稠密矩阵。现象七同一份代码两次运行结果差异很大。压缩感知流程里至少有三个随机源掩模生成、IRLS 初值、某些实现里的随机打乱。在脚本开头统一rng(0)并且在 IRLS 里固定初值结果就能复现。这一条看着琐碎但凡是做算法对比不复现的结果等于没做。注意rng要放在每次生成随机量之前不能只放在脚本开头一次。Matlab 的随机状态是全局的中间任何一处调用了随机函数都会把它推进后面生成的东西就变了。现象八换成更高分辨率后代码能跑但图明显变差。除了内存和迭代次数还有个常被忽略的点采样率的分母是 n分辨率一升同样的 M 对应的采样率就掉下来了。64×64 用 1000 次测量是 24%256×256 用同样的 1000 次只有 1.5%完全不在一个难度级别。这类问题的排查方式就是打印一下实际的 M/n别凭感觉。6. 设计取舍先验、采样预算与验证成本的三方博弈到这里可以把整条链路的选择收敛成一个三角形采样预算、先验强度、验证成本三者互相拉扯没有一种重构器能同时拿下三边。采样预算宽裕时M/n 在 30% 以上L2 加 Tikhonov 就够用实现简单、跑得快、参数只有一个 λ调试成本极低。这个区间里上 TV 或 L1 属于过度设计画质提升 1~2 dB代价是代码复杂度和运行时间翻几倍。很多人一上来就上 TV其实是被「TV 更强」这句话带跑了没算过性价比。采样预算紧的时候M/n 在 10% 以下先验就是命。这时候纯 L2 基本没法看必须用 L1 或 TV。选 L1 还是 TV取决于图像内容自然场景、医学影像这类分片光滑的TV 占优稀疏点目标、编码孔径这类本身就是稀疏的L1 更合适。选错先验的后果比不调参数还严重。6.1 一套流程是否可信的最小验证每次改完重构代码我一般会跑一个三板斧验证全程不到五分钟但能挡掉绝大多数低级错误。验证项做法期望结果正向一致性用真实 x 算 y再用无噪声的 CG 反解残差应接近机器精度图与原图几乎一致采样率下降曲线50% → 20% → 10% → 5%固定其他变量PSNR 单调下降不应出现反弹参数敏感性λ 取 1e-5 到 1e0 扫一遍应存在明显的最优区间曲线呈单峰第一项最容易被跳过但它能挡掉最致命的一类错误正向模型和反向重构用的是同一套 Φ 吗有没有在中间不小心转置了有一次我把Phi * Phi里的转置漏掉一个残差怎么都降不下去查了一下午才发现在那一行。第二项能挡掉过拟合式的误判——如果 5% 采样率下的 PSNR 反而比 10% 高那一定是某处有信息泄漏。第三项的参数扫描值得多说一句。λ 的最优区间通常在两个数量级以内如果扫出来曲线是平的或者一直往上走说明你的 λ 量纲和数据量级没对齐得先做归一化再扫。这一步做完再去做跨算法对比结论才站得住。从那以后我每次接手一份单像素成像的重构代码第一件事不是看它用了什么高级算法而是先打印出 Φ 的行范数、y 的均值方差、λ 的实际数量级这三个数。这三个数对不上后面所有调参都是白费。希望帮到你。本文还有配套的精品资源点击获取
简介这份Matlab单像素成像算法比较代码包面向已掌握基本成像原理、希望系统对比不同重建方法的读者围绕单像素成像中的稀疏采样与图像重建问题提供DGI、梯度下降、共轭梯度下降、Poisson最大似然、交替投影、稀疏表示压缩感知及全变分压缩感知等七类方法的可运行示例。包内共10个文件以9个.m脚本和1个txt说明为主脚本分别对应各算法实现及误差计算、演示入口说明文件用于辅助理解调用关系压缩包约12KB体量轻便便于快速查阅与二次修改。目前已有2762人学习适合课程实验、算法复现和不同重建策略的初步对照。读者可借此理解各方法在采样率、重建质量与计算开销上的差异并基于统一误差函数开展结果比较为单像素成像实验选型与后续改进提供参考。1. 一个探测器的悖论为什么几百次测量就能还原一张图一份典型的 Matlab 单像素成像脚本摆在面前图像 64×64 共 4096 个像素测量次数只给 800 次重构那一段调的是共轭梯度下降。按线性代数直觉这是严重欠定的方程组解不唯一按压缩感知的说法这叫欠采样重构听起来像玄学。两个东西拼在一起居然能出图这件事本身值得停一秒想清楚到底是哪一方在补偿缺失的那些方程如果以为共轭梯度下降天生就能做稀疏重构那后面对参数的所有直觉都会跑偏。把这条链路拆开看过一次之后mask 怎么选、正则项怎么调、为什么有人跑 20 次迭代就收敛而有人跑到 500 次还在震荡就都有解释。2. 从光强到方程单像素成像的测量矩阵是怎么长出来的2.1 掩模调制与标量测量一步物理怎么变成一行矩阵单像素成像的硬件链路其实短得可怜光源照到场景上反射光经过一组掩模图案调制再由一个聚光镜把所有能量收进单像素探测器。探测器没有空间分辨能力它读到的只是一个数——这一帧的总光强。想要二维图像就必须在时间维度上换掩模换一次图案读一个数换 M 次得到 M 个数。关键在于每次读到的这个数本质上是场景光强分布与当前掩模的内积。把场景离散成向量 x长度 n N×N第 i 次测量用的掩模离散成向量 φ_i同样长度 n那么探测器读数 y_i ≈ Σ_j φ_ij · x_j。M 次测量堆起来就是 y ΦxΦ 是 M×n 的测量矩阵每一行对应一次掩模曝光。一段最朴素的 Matlab 正向仿真长这样N 64; % 图像边长 n N * N; % 未知量个数 4096 img double(imread(target.png)); x img(:); % n x 1 列向量 M round(0.25 * n); % 采样率 25%即 1024 次测量 rng(0); % 固定随机源否则每次跑结果都不同 Phi randn(M, n) / sqrt(M); % 高斯随机测量矩阵按行归一 y Phi * x; % 探测器读到的 M 个标量这段代码里有两处容易被忽略但影响极大。第一处是/ sqrt(M)把每一行的能量归一到 1是为了让 y 的量级与 x 处在同一水平否则后面正则项系数 λ 的量纲会跟着采样率一起漂移换个 M 就得重新调参。第二处是rng(0)随机掩模一旦不固定种子同一份代码两次运行的重建质量能差出好几个 dB调试时会误以为是算法本身不稳定。这条链路还有个物理上的隐含条件常被跳过掩模刷新与探测器采样必须同步。真实系统里 DMD 的刷新有建立时间探测器积分窗如果没对齐读到的就是两张掩模的混合结果等效于 Φ 的每一行都被悄悄污染了一点点。仿真里看不出来一上硬件就是「图总是糊一层」的源头。2.2 稀疏基与不相干性为什么随机掩模比哈达玛掩模更容易命中有了 y Φx 还不够M n 时它有无数个解。压缩感知的解法是把 x 换到另一个域里表示x Ψθ其中 Ψ 是稀疏基DCT、小波、或者图像本身θ 只有 K 个显著非零元K n。代进去得到 y ΦΨθ AθA 是 M×n 的传感矩阵。这里有个非常容易被跳过的约束Φ 和 Ψ 必须「不相干」。相干性的定义是把两个基里最相似的一对原子的内积挑出来再乘上 √n 归一。相干性越高一个基里的稀疏信号在另一个基里就越「铺得开」压缩感知的恢复保证就越弱。这条约束直接决定了掩模怎么选。如果稀疏基用的是哈达玛变换而测量掩模也用哈达玛图案两个基完全重合相干性等于最大值这时候压缩感知的理论保证直接失效——你等于在稀疏域里逐点采样只是自己没意识到而已。换成伯努利二值掩模或者高斯随机掩模与 DCT、小波这些常见基的相干性都很低恢复才稳。血泪经验是很多人第一次做单像素成像会选哈达玛掩模因为它二值、好实现、有快速变换、重建还能用逆变换一步到位。它在高采样率下效果确实漂亮但一旦把采样率压到 10% 以下重构质量的下降速度明显比随机掩模陡。原因不神秘哈达玛掩模压到低采样率时是「均匀地少采」而随机掩模是「不相干地乱采」后者的信息冗余更低。2.3 从 L0 到 L1凸松弛到底换来了什么理论上最想要的是 L0 范数最小在所有满足 y Aθ 的解里挑非零元最少的那个。问题是这属于组合优化n 上千时穷举不现实。L1 范数是 L0 在单位球上的凸包把它替换进去问题变成可解的凸优化min ||θ||₁ s.t. y Aθ几何上的直觉是L1 球在坐标轴方向有尖角当残差约束的可行域高维空间里的一个仿射子空间与这个球相切时切点大概率落在坐标轴上也就是稀疏散点上而 L2 球是光滑的切点会落在球面上解里几乎每个分量都不为零图像就糊了。这就是为什么纯最小二乘解出来的图总是「整体形状对、细节全丢」。实际有噪声时等式约束要放松成不等式或惩罚形式min ½||y − Aθ||₂² λ||θ||₁到这里才出现了这份代码里真正要解的东西。λ 大稀疏性强、解偏平滑、细节丢λ 小拟合狠、容易把噪声也当信号拟合进去。这个式子不是二次的共轭梯度下降没法直接吃它后面第 3 章要拆的就是这层关系。2.4 噪声进来之后理想 RIP 在实验台上永远不会成立教科书上的恢复保证依赖 RIP有限等距性质对任意 K 稀疏向量 θ要求 (1−δ)||θ||² ≤ ||Aθ||² ≤ (1δ)||θ||²δ 足够小。这个条件在仿真里勉强成立在真实装置上几乎不可能满足。干扰源至少有四类。探测器有暗电流和读出噪声低光强下还叠加泊松噪声表现为 y 的方差随信号强度变化不是简单的高斯白噪声。掩模与场景的相对位置有微小位移等效于 A 的每一行都在抖。光源强度随时间漂移等效于 y 整体乘了个慢变系数。掩模本身的二值化不理想灰度掩模的透过率不是严格的 0 和 1。这四类里最坑的是光源漂移因为它不改变图像形状只是让 y 的整体尺度慢慢变重构出来的图看起来「形状对、对比度不对」很容易被误判成 λ 没调好。常见的处理思路是在每一批测量里插一帧全通掩模所有元素为 1用它的读数做归一化基准把慢变系数除掉。这招在论文里不显眼在实验台上能省掉大量无效调参。3. 共轭梯度下降它解的到底是不是压缩感知那个问题3.1 Krylov 子空间与 A-共轭CG 比最速下降快在哪共轭梯度下降解决的是这类问题A 对称正定求 Ax b。它的迭代不是在原始空间里沿负梯度走而是在一个逐步扩张的子空间里找最优解这个子空间叫 Krylov 子空间 K_k(A, r₀) span{r₀, Ar₀, A²r₀, …}。每一步都把当前残差正交化到之前所有搜索方向上于是搜索方向两两 A-共轭p_iᵀ A p_j 0i ≠ j。这个性质带来的效果非常直观。最速下降每步都沿当前负梯度走相邻两步的方向几乎正交路径是锯齿状的条件数大时收敛极慢。CG 因为方向互相 A-共轭不会重复走已经走过的方向理论上在精确算术下 n 步内必收敛——尽管浮点运算下远达不到但实际往往几十步就够了。收敛速度有个常用估计误差的 A-范数每步至少乘上 (√κ − 1)/(√κ 1)κ 是 A 的条件数。这句话信息量很大条件数 100 时这个因子约 0.82收敛慢条件数 10 时约 0.52几十步就够。所以对 CG 来说比「迭代多少次」更重要的是「矩阵条件数是多少」。这也解释了后面为什么一定要加正则项。3.2 正规方程、Tikhonov 与 CGLS三种写法的数值差异最小二乘问题 min ||Aθ − y||² 的最优性条件是 AᵀAθ Aᵀy这就是正规方程。AᵀA 对称半正定可以用 CG 解。但这里埋了个雷AᵀA 的条件数是 A 的条件数的平方。A 本身条件数 50 还算温和平方之后 2500CG 的收敛因子就变成 0.96 左右跑几百步都未必收敛。第二种写法是 Tikhonov 正则把 (AᵀA λI)θ Aᵀy 拿来做。加 λI 相当于给所有特征值抬高一个地板条件数从 κ² 降到 (σ_max² λ)/(σ_min² λ)小奇异值被压住了收敛立刻变快。代价是解被往零的方向拉了一点图像会略偏平滑——但这恰好和压缩感知要的「抑制噪声」方向一致很多时候不是坏事。第三种是 CGLS它不显式构造 AᵀA而是直接在 A 上做 Krylov 迭代每一步只做一次 A 乘和一次 Aᵀ 乘。数值上比正规方程稳得多条件数不会平方放大。代价是实现稍复杂需要维护两套递推量。写法求解对象有效条件数每步矩阵乘法适用场景正规方程 CGAᵀAθ Aᵀyκ²2 次小规模、A 良态Tikhonov CG(AᵀAλI)θ Aᵀy(σ²λ) 比值2 次有噪声、需要正则CGLS / LSQR直接在 A 上迭代κ2 次大规模、病态有一处容易翻车在 Matlab 里显式写A * A会真的构造出一个 n×n 矩阵。n 4096 时是 128 MB勉强能跑n 256×256 65536 时是 34 GB直接内存爆炸。这就是下一节要讲的函数句柄写法为什么必要。3.3 手写 CG 与 pcg从显式矩阵到函数句柄先看手写版本它把每一步在干什么摊开调试时能直接看到残差曲线function [x, relres, k] cg_solve(Afun, b, x, tol, maxit) % Afun : 对称正定算子句柄输入列向量返回 A*v % b : 右端项 (n x 1) % x : 初值同时作为输出 % tol : 相对残差停止阈值常用 1e-6 % maxit : 最大迭代步数 r b - Afun(x); p r; rs r * r; bn norm(b); for k 1:maxit Ap Afun(p); pAp p * Ap; if pAp 0 warning(算子非正定pAp %g检查正则项, pAp); break; end alpha rs / pAp; % 沿 p 方向走多远 x x alpha * p; r r - alpha * Ap; rs_new r * r; relres sqrt(rs_new) / bn; if relres tol, break; end p r (rs_new / rs) * p; % 新方向与旧方向 A-共轭 rs rs_new; end end这段里pAp 0那个检查不是摆设。压缩感知里如果 λ 设成 0 而 A 严重亏秩AᵀA 是半正定的pAp 可能接近 0alpha 会炸成天文数字图直接变成一片雪花。加一句 warning 至少能让你知道是自己参数设错了而不是算法烂。实际工程里我更倾向直接用 Matlab 自带的pcg把算子写成匿名函数永远不显式构造 AᵀAAfun (v) Phi * (Phi * v) lambda * v; % Tikhonov 算子 rhs Phi * y; [theta, flag, relres, iter] pcg(Afun, rhs, 1e-6, 200);flag的取值含义要记住0 是收敛1 是达到最大迭代仍没收敛2 是预条件子病态3 是迭代停滞4 是某步出现非有限值。调试时如果看到 flag 1别急着加迭代次数上限先看条件数——十有八九是 λ 太小或者干脆没加。看到 flag 4多半是 Phi 里有 Inf 或者 NaN往上游查数据归一化。3.4 把 CG 塞进 L1IRLS 外层循环怎么搭到这一步必须把话说透CG 解的是二次型问题而压缩感知的目标函数带 L1 项不是二次的。CG 不能直接解压缩感知它只能当内层的一个线性求解器。这一层关系如果搞混调参会调得非常迷茫。把 L1 和 CG 接起来的经典做法是迭代重加权最小二乘IRLS。核心思想是|θ_i| ≈ θ_i² / |θ_i|用当前估计的 |θ_i| 当权重把 L1 问题转成一串加权 L2 问题每个 L2 子问题用 CG 解。lambda 1e-3; % 正则强度越大越稀疏 eps_reg 1e-6; % 防止 1/0常用 1e-6 ~ 1e-4 theta ones(n, 1); % 初值不要全零否则第一步权重会爆 Aop (v) Phi * (Phi * v); rhs Phi * y; for outer 1:15 w 1 ./ (abs(theta) eps_reg); % 越接近 0 的解权重越大 Rop (v) lambda * (w .* v); % 加权 L2 正则项 theta pcg((v) Aop(v) Rop(v), rhs, 1e-6, 200, [], [], theta); fprintf(IRLS %2d: 非零元 %d\n, outer, nnz(abs(theta) 1e-3)); end x_rec Psi * theta; % 回到图像域三处参数值得单独说。eps_reg决定了「多小算零」设太大会把真实的小系数一起压掉图像细节丢失设太小会让权重矩阵条件数爆表内层 CG 收敛极慢。经验值 1e-6 起步如果内层 flag 经常是 1就往 1e-4 调。lambda控制整体稀疏强度量级通常要和norm(Phi*y, inf)对齐看单纯照抄别人的 1e-3 在你自己数据上很可能失效。外层循环次数 10~20 次足够再往上收益递减可以直接看非零元个数有没有稳定。注意这套 IRLS 只是众多 L1 求解路线里的一种。ADMM 把它拆成「数据保真 软阈值」两个子问题交替解FISTA 用梯度步加软阈值再加动量加速TV 正则则要在 L1 之外再套一层梯度算子。选哪条路线取决于你的先验是「稀疏」还是「分片光滑」下一章对比的就是这件事。4. 四种重构器同场对比采样率降到多少开始崩4.1 统一实验骨架掩模、采样率与评价指标对比算法的第一原则是所有变量只动一个。固定同一张图、同一个随机种子、同一个 Φ只换重构器否则你根本分不清差异来自算法还是来自掩模。骨架大概是这样rng(2024); img double(imread(target.png)); img imresize(img, [64 64]); x img(:); n numel(x); rates [0.05 0.10 0.20 0.35 0.50]; for r rates M round(r * n); Phi randn(M, n) / sqrt(M); y Phi * x; % 在这里依次调用 CG-Tikhonov / OMP / FISTA / TV-ADMM % 每个算法跑完记录 psnr(x, x_rec)、ssim 和 toc end指标选 PSNR 和 SSIM 两个就够。PSNR 对整体亮度误差敏感SSIM 对结构相似度敏感两者经常给出不同排序——一个算法 PSNR 高但 SSIM 低通常意味着它把噪声抹平了但边缘也磨掉了。只看一个指标下结论容易得出偏颇的结论。还有一条必须统一每个算法的 λ 都要在自己的参数网格上取最优再拿来对比。拿一个没调参的 FISTA 去比一个调好的 CG这种对比没有意义。代价是实验量翻好几倍但这是唯一能让结论站住的做法。4.2 重建质量随采样率的变化谁先崩趋势大致是这样的。在 50% 采样率附近四种方法的 PSNR 差距不明显往往在 1~2 dB 以内因为信息本来就够多什么方法都能凑出来。差距是从 20% 往下拉开的。纯 CG 加 Tikhonov 属于 L2 正则它给出的解在所有方向上都平滑本质是「最小能量解」稀疏性完全没有被利用。采样率降到 20% 以下重建图会明显偏糊边缘发软纹理区域糊成一片。它的优势是快——迭代次数少每步只做两次矩阵乘法。OMP 走的是另一条路贪心地每次挑一个与当前残差最相关的原子然后在这一小撮原子上做最小二乘。它在极度稀疏K 很小且 M 不太小时效果很好但有两个软肋。一是需要预先知道稀疏度 K猜大了会过拟合噪声猜小了会丢结构二是每一步都要算 Aᵀr复杂度随原子数线性涨K 大时非常慢。FISTA 是 L1 正则的主力选手梯度步加软阈值再叠动量收敛速率比 ISTA 快一个量级。它在 10%~30% 采样率区间表现稳定是中低采样率下的性价比之选。软阈值的副作用是会把小系数整体往零拉图像会有轻微的「块状感」这是 L1 的固有偏置不是实现 bug。TV 正则ADMM 或增广拉格朗日求解在低采样率下通常画质最好因为图像的分片光滑先验比纯稀疏先验更贴合自然图像。代价是慢而且慢得不只是常数倍——它每次迭代内部还要再解一个子问题整体比 FISTA 慢一个数量级很常见。提示这里的排序是趋势层面的具体到某一张图、某一组掩模名次可能互换。尤其是纹理极丰富的图像TV 的分片光滑假设会失效反而不如 FISTA。4.3 运行时间与内存TV 类方法慢在哪时间差异的量级值得单独拎出来看因为它直接决定了你走不走得通高分辨率。方法先验类型每次迭代代价典型迭代数主要失败模式CG TikhonovL2 平滑2 次矩阵乘20~80低采样率下糊OMP严格稀疏1 次相关 1 次 LSK 步K 估计错就崩FISTAL1 稀疏2 次矩阵乘 软阈值100~300块状伪影TV ADMM分片光滑内层子问题求解100~500极慢、参数敏感从表里能读出一个关键取舍CG 快但先验弱TV 先验强但慢。FISTA 卡在中间这也是它在实际项目里出场率最高的原因。内存方面只要坚持用函数句柄四种方法的峰值内存基本都在「Phi 本身 几个 n 维向量」这个量级。真正的内存杀手是显式构造 AᵀA 或者把 TV 的梯度算子展开成稠密矩阵。n 65536 时一个稠密 n×n 就是 34 GB这个坑踩一次能记一辈子。ADMM 实现里有一步常见的错误是「为了写起来方便」把差分算子写成diff(eye(n))那一行代码就是内存爆炸的直接原因正确写法是把它写成只做邻域差分的函数句柄。5. 排查图糊、图碎、不收敛时先看哪几个量5.1 图像糊、碎、全黑先查信号链再查算法现象一重建图整体形状对但像蒙了一层雾边缘发软。原因大概率不是算法而是测量矩阵没有按行归一化。如果 Φ 的每一行能量不同y 的各个分量量级就会参差不齐最小二乘会自动更重视量级大的那些测量等价于给不同掩模加了不同权重。解决很直接构造完 Phi 之后强制Phi Phi ./ vecnorm(Phi, 2, 2)再检查vecnorm(Phi,2,2)是不是全 1。现象二图碎成马赛克状或者有明显的方块拼接感。这通常是稀疏基和掩模不匹配。用哈达玛掩模配小波基或者用 DCT 基去配本身就带 DCT 结构的掩模相干性过高恢复出来的解在基里分布得很散重建图上就表现为均匀的块状噪声。解决是换随机二值或高斯掩模重跑一次看看块状感有没有消失。现象三重建全黑或者全白。这种八成是 λ 和数据量级严重失配。检查一下norm(Phi*y)和lambda的数量级如果差了三四个数量级λ 就相当于无穷大解直接被压到零。另一个可能是 y 里混进了 NaN往上顺一遍数据读入有没有除零。现象四低采样率下某个固定区域永远重建不出来。这往往是掩模覆盖不均。随机掩模虽然整体不相干但 M 很小时局部可能恰好没有充分覆盖某些位置尤其图像边角。解决是换一组种子多跑几次取平均或者用结构化随机把区域分块后每块内独立随机保证空间覆盖均匀。5.2 不收敛、内存爆、每次结果不同查参数与随机源现象五pcg 返回 flag 1迭代跑满还没收敛。先看两个量condest估一下有效算子的条件数以及 λ 到底设了多少。不加正则或 λ 极小的时候AᵀA 的最小奇异值接近零CG 收敛因子接近 1当然跑不满。解决办法是加 λ从 1e-3 起试或者换 CGLS/LSQR 这类不平方条件数的方法。现象六跑到一半内存爆掉。基本就两种可能显式构造了 AᵀA或者 TV 的差分算子被展开成了稠密矩阵。前者改成函数句柄后者改成邻域差分。还有个隐蔽的来源是spdiags用的位置不对——IRLS 里的权重矩阵 W 一定要用spdiags生成稀疏对角阵写成diag(w)就是又一个 n×n 稠密矩阵。现象七同一份代码两次运行结果差异很大。压缩感知流程里至少有三个随机源掩模生成、IRLS 初值、某些实现里的随机打乱。在脚本开头统一rng(0)并且在 IRLS 里固定初值结果就能复现。这一条看着琐碎但凡是做算法对比不复现的结果等于没做。注意rng要放在每次生成随机量之前不能只放在脚本开头一次。Matlab 的随机状态是全局的中间任何一处调用了随机函数都会把它推进后面生成的东西就变了。现象八换成更高分辨率后代码能跑但图明显变差。除了内存和迭代次数还有个常被忽略的点采样率的分母是 n分辨率一升同样的 M 对应的采样率就掉下来了。64×64 用 1000 次测量是 24%256×256 用同样的 1000 次只有 1.5%完全不在一个难度级别。这类问题的排查方式就是打印一下实际的 M/n别凭感觉。6. 设计取舍先验、采样预算与验证成本的三方博弈到这里可以把整条链路的选择收敛成一个三角形采样预算、先验强度、验证成本三者互相拉扯没有一种重构器能同时拿下三边。采样预算宽裕时M/n 在 30% 以上L2 加 Tikhonov 就够用实现简单、跑得快、参数只有一个 λ调试成本极低。这个区间里上 TV 或 L1 属于过度设计画质提升 1~2 dB代价是代码复杂度和运行时间翻几倍。很多人一上来就上 TV其实是被「TV 更强」这句话带跑了没算过性价比。采样预算紧的时候M/n 在 10% 以下先验就是命。这时候纯 L2 基本没法看必须用 L1 或 TV。选 L1 还是 TV取决于图像内容自然场景、医学影像这类分片光滑的TV 占优稀疏点目标、编码孔径这类本身就是稀疏的L1 更合适。选错先验的后果比不调参数还严重。6.1 一套流程是否可信的最小验证每次改完重构代码我一般会跑一个三板斧验证全程不到五分钟但能挡掉绝大多数低级错误。验证项做法期望结果正向一致性用真实 x 算 y再用无噪声的 CG 反解残差应接近机器精度图与原图几乎一致采样率下降曲线50% → 20% → 10% → 5%固定其他变量PSNR 单调下降不应出现反弹参数敏感性λ 取 1e-5 到 1e0 扫一遍应存在明显的最优区间曲线呈单峰第一项最容易被跳过但它能挡掉最致命的一类错误正向模型和反向重构用的是同一套 Φ 吗有没有在中间不小心转置了有一次我把Phi * Phi里的转置漏掉一个残差怎么都降不下去查了一下午才发现在那一行。第二项能挡掉过拟合式的误判——如果 5% 采样率下的 PSNR 反而比 10% 高那一定是某处有信息泄漏。第三项的参数扫描值得多说一句。λ 的最优区间通常在两个数量级以内如果扫出来曲线是平的或者一直往上走说明你的 λ 量纲和数据量级没对齐得先做归一化再扫。这一步做完再去做跨算法对比结论才站得住。从那以后我每次接手一份单像素成像的重构代码第一件事不是看它用了什么高级算法而是先打印出 Φ 的行范数、y 的均值方差、λ 的实际数量级这三个数。这三个数对不上后面所有调参都是白费。希望帮到你。本文还有配套的精品资源点击获取
