简介：这份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 = 0，i ≠ 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 还算温和，平方之后 2500，CG 的收敛因子就变成 0.96 左右，跑几百步都未必收敛。

第二种写法是 Tikhonov 正则：把 (AᵀA + λI)θ = Aᵀy 拿来做。加 λI 相当于给所有特征值抬高一个地板，条件数从 κ² 降到 (σ_max² + λ)/(σ_min² + λ)，小奇异值被压住了，收敛立刻变快。代价是解被往零的方向拉了一点，图像会略偏平滑——但这恰好和压缩感知要的「抑制噪声」方向一致，很多时候不是坏事。

第三种是 CGLS，它不显式构造 AᵀA，而是直接在 A 上做 Krylov 迭代，每一步只做一次 A 乘和一次 Aᵀ 乘。数值上比正规方程稳得多，条件数不会平方放大。代价是实现稍复杂，需要维护两套递推量。

有一处容易翻车：在 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 可能接近 0，alpha 会炸成天文数字，图直接变成一片雪花。加一句 warning 至少能让你知道是自己参数设错了，而不是算法烂。

实际工程里我更倾向直接用 Matlab 自带的pcg，把算子写成匿名函数，永远不显式构造 AᵀA：

Afun = @(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 塞进 L1：IRLS 外层循环怎么搭

到这一步必须把话说透：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 快但先验弱，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 一套流程是否可信的最小验证

每次改完重构代码，我一般会跑一个三板斧验证，全程不到五分钟，但能挡掉绝大多数低级错误。

第一项最容易被跳过，但它能挡掉最致命的一类错误：正向模型和反向重构用的是同一套 Φ 吗？有没有在中间不小心转置了？有一次我把Phi' * Phi里的转置漏掉一个，残差怎么都降不下去，查了一下午才发现在那一行。第二项能挡掉过拟合式的误判——如果 5% 采样率下的 PSNR 反而比 10% 高，那一定是某处有信息泄漏。

第三项的参数扫描值得多说一句。λ 的最优区间通常在两个数量级以内，如果扫出来曲线是平的或者一直往上走，说明你的 λ 量纲和数据量级没对齐，得先做归一化再扫。这一步做完再去做跨算法对比，结论才站得住。

从那以后我每次接手一份单像素成像的重构代码，第一件事不是看它用了什么高级算法，而是先打印出 Φ 的行范数、y 的均值方差、λ 的实际数量级这三个数。这三个数对不上，后面所有调参都是白费。希望帮到你。

本文还有配套的精品资源，点击获取