：測量矩陣、稀疏正則與共軛梯度實戰(zhàn))
簡介這份Matlab單像素成像算法比較代碼包面向已掌握基本成像原理、希望系統(tǒng)對比不同重建方法的讀者圍繞單像素成像中的稀疏采樣與圖像重建問題提供DGI、梯度下降、共軛梯度下降、Poisson最大似然、交替投影、稀疏表示壓縮感知及全變分壓縮感知等七類方法的可運行示例。包內(nèi)共10個文件以9個.m腳本和1個txt說明為主腳本分別對應各算法實現(xiàn)及誤差計算、演示入口說明文件用于輔助理解調(diào)用關系壓縮包約12KB體量輕便便于快速查閱與二次修改。目前已有2762人學習適合課程實驗、算法復現(xiàn)和不同重建策略的初步對照。讀者可借此理解各方法在采樣率、重建質(zhì)量與計算開銷上的差異并基于統(tǒng)一誤差函數(shù)開展結(jié)果比較為單像素成像實驗選型與后續(xù)改進提供參考。1. 一個探測器的悖論為什么幾百次測量就能還原一張圖一份典型的 Matlab 單像素成像腳本擺在面前圖像 64×64 共 4096 個像素測量次數(shù)只給 800 次重構(gòu)那一段調(diào)的是共軛梯度下降。按線性代數(shù)直覺這是嚴重欠定的方程組解不唯一按壓縮感知的說法這叫欠采樣重構(gòu)聽起來像玄學。兩個東西拼在一起居然能出圖這件事本身值得停一秒想清楚到底是哪一方在補償缺失的那些方程如果以為共軛梯度下降天生就能做稀疏重構(gòu)那后面對參數(shù)的所有直覺都會跑偏。把這條鏈路拆開看過一次之后mask 怎么選、正則項怎么調(diào)、為什么有人跑 20 次迭代就收斂而有人跑到 500 次還在震蕩就都有解釋。2. 從光強到方程單像素成像的測量矩陣是怎么長出來的2.1 掩模調(diào)制與標量測量一步物理怎么變成一行矩陣單像素成像的硬件鏈路其實短得可憐光源照到場景上反射光經(jīng)過一組掩模圖案調(diào)制再由一個聚光鏡把所有能量收進單像素探測器。探測器沒有空間分辨能力它讀到的只是一個數(shù)——這一幀的總光強。想要二維圖像就必須在時間維度上換掩模換一次圖案讀一個數(shù)換 M 次得到 M 個數(shù)。關鍵在于每次讀到的這個數(shù)本質(zhì)上是場景光強分布與當前掩模的內(nèi)積。把場景離散成向量 x長度 n N×N第 i 次測量用的掩模離散成向量 φ_i同樣長度 n那么探測器讀數(shù) y_i ≈ Σ_j φ_ij · x_j。M 次測量堆起來就是 y ΦxΦ 是 M×n 的測量矩陣每一行對應一次掩模曝光。一段最樸素的 Matlab 正向仿真長這樣N 64; % 圖像邊長 n N * N; % 未知量個數(shù) 4096 img double(imread(target.png)); x img(:); % n x 1 列向量 M round(0.25 * n); % 采樣率 25%即 1024 次測量 rng(0); % 固定隨機源否則每次跑結(jié)果都不同 Phi randn(M, n) / sqrt(M); % 高斯隨機測量矩陣按行歸一 y Phi * x; % 探測器讀到的 M 個標量這段代碼里有兩處容易被忽略但影響極大。第一處是/ sqrt(M)把每一行的能量歸一到 1是為了讓 y 的量級與 x 處在同一水平否則后面正則項系數(shù) λ 的量綱會跟著采樣率一起漂移換個 M 就得重新調(diào)參。第二處是rng(0)隨機掩模一旦不固定種子同一份代碼兩次運行的重建質(zhì)量能差出好幾個 dB調(diào)試時會誤以為是算法本身不穩(wěn)定。這條鏈路還有個物理上的隱含條件常被跳過掩模刷新與探測器采樣必須同步。真實系統(tǒng)里 DMD 的刷新有建立時間探測器積分窗如果沒對齊讀到的就是兩張掩模的混合結(jié)果等效于 Φ 的每一行都被悄悄污染了一點點。仿真里看不出來一上硬件就是「圖總是糊一層」的源頭。2.2 稀疏基與不相干性為什么隨機掩模比哈達瑪掩模更容易命中有了 y Φx 還不夠M n 時它有無數(shù)個解。壓縮感知的解法是把 x 換到另一個域里表示x Ψθ其中 Ψ 是稀疏基DCT、小波、或者圖像本身θ 只有 K 個顯著非零元K n。代進去得到 y ΦΨθ AθA 是 M×n 的傳感矩陣。這里有個非常容易被跳過的約束Φ 和 Ψ 必須「不相干」。相干性的定義是把兩個基里最相似的一對原子的內(nèi)積挑出來再乘上 √n 歸一。相干性越高一個基里的稀疏信號在另一個基里就越「鋪得開」壓縮感知的恢復保證就越弱。這條約束直接決定了掩模怎么選。如果稀疏基用的是哈達瑪變換而測量掩模也用哈達瑪圖案兩個基完全重合相干性等于最大值這時候壓縮感知的理論保證直接失效——你等于在稀疏域里逐點采樣只是自己沒意識到而已。換成伯努利二值掩?；蛘吒咚闺S機掩模與 DCT、小波這些常見基的相干性都很低恢復才穩(wěn)。血淚經(jīng)驗是很多人第一次做單像素成像會選哈達瑪掩模因為它二值、好實現(xiàn)、有快速變換、重建還能用逆變換一步到位。它在高采樣率下效果確實漂亮但一旦把采樣率壓到 10% 以下重構(gòu)質(zhì)量的下降速度明顯比隨機掩模陡。原因不神秘哈達瑪掩模壓到低采樣率時是「均勻地少采」而隨機掩模是「不相干地亂采」后者的信息冗余更低。2.3 從 L0 到 L1凸松弛到底換來了什么理論上最想要的是 L0 范數(shù)最小在所有滿足 y Aθ 的解里挑非零元最少的那個。問題是這屬于組合優(yōu)化n 上千時窮舉不現(xiàn)實。L1 范數(shù)是 L0 在單位球上的凸包把它替換進去問題變成可解的凸優(yōu)化min ||θ||? s.t. y Aθ幾何上的直覺是L1 球在坐標軸方向有尖角當殘差約束的可行域高維空間里的一個仿射子空間與這個球相切時切點大概率落在坐標軸上也就是稀疏散點上而 L2 球是光滑的切點會落在球面上解里幾乎每個分量都不為零圖像就糊了。這就是為什么純最小二乘解出來的圖總是「整體形狀對、細節(jié)全丟」。實際有噪聲時等式約束要放松成不等式或懲罰形式min ?||y ? Aθ||?2 λ||θ||?到這里才出現(xiàn)了這份代碼里真正要解的東西。λ 大稀疏性強、解偏平滑、細節(jié)丟λ 小擬合狠、容易把噪聲也當信號擬合進去。這個式子不是二次的共軛梯度下降沒法直接吃它后面第 3 章要拆的就是這層關系。2.4 噪聲進來之后理想 RIP 在實驗臺上永遠不會成立教科書上的恢復保證依賴 RIP有限等距性質(zhì)對任意 K 稀疏向量 θ要求 (1?δ)||θ||2 ≤ ||Aθ||2 ≤ (1δ)||θ||2δ 足夠小。這個條件在仿真里勉強成立在真實裝置上幾乎不可能滿足。干擾源至少有四類。探測器有暗電流和讀出噪聲低光強下還疊加泊松噪聲表現(xiàn)為 y 的方差隨信號強度變化不是簡單的高斯白噪聲。掩模與場景的相對位置有微小位移等效于 A 的每一行都在抖。光源強度隨時間漂移等效于 y 整體乘了個慢變系數(shù)。掩模本身的二值化不理想灰度掩模的透過率不是嚴格的 0 和 1。這四類里最坑的是光源漂移因為它不改變圖像形狀只是讓 y 的整體尺度慢慢變重構(gòu)出來的圖看起來「形狀對、對比度不對」很容易被誤判成 λ 沒調(diào)好。常見的處理思路是在每一批測量里插一幀全通掩模所有元素為 1用它的讀數(shù)做歸一化基準把慢變系數(shù)除掉。這招在論文里不顯眼在實驗臺上能省掉大量無效調(diào)參。3. 共軛梯度下降它解的到底是不是壓縮感知那個問題3.1 Krylov 子空間與 A-共軛CG 比最速下降快在哪共軛梯度下降解決的是這類問題A 對稱正定求 Ax b。它的迭代不是在原始空間里沿負梯度走而是在一個逐步擴張的子空間里找最優(yōu)解這個子空間叫 Krylov 子空間 K_k(A, r?) span{r?, Ar?, A2r?, …}。每一步都把當前殘差正交化到之前所有搜索方向上于是搜索方向兩兩 A-共軛p_i? A p_j 0i ≠ j。這個性質(zhì)帶來的效果非常直觀。最速下降每步都沿當前負梯度走相鄰兩步的方向幾乎正交路徑是鋸齒狀的條件數(shù)大時收斂極慢。CG 因為方向互相 A-共軛不會重復走已經(jīng)走過的方向理論上在精確算術(shù)下 n 步內(nèi)必收斂——盡管浮點運算下遠達不到但實際往往幾十步就夠了。收斂速度有個常用估計誤差的 A-范數(shù)每步至少乘上 (√κ ? 1)/(√κ 1)κ 是 A 的條件數(shù)。這句話信息量很大條件數(shù) 100 時這個因子約 0.82收斂慢條件數(shù) 10 時約 0.52幾十步就夠。所以對 CG 來說比「迭代多少次」更重要的是「矩陣條件數(shù)是多少」。這也解釋了后面為什么一定要加正則項。3.2 正規(guī)方程、Tikhonov 與 CGLS三種寫法的數(shù)值差異最小二乘問題 min ||Aθ ? y||2 的最優(yōu)性條件是 A?Aθ A?y這就是正規(guī)方程。A?A 對稱半正定可以用 CG 解。但這里埋了個雷A?A 的條件數(shù)是 A 的條件數(shù)的平方。A 本身條件數(shù) 50 還算溫和平方之后 2500CG 的收斂因子就變成 0.96 左右跑幾百步都未必收斂。第二種寫法是 Tikhonov 正則把 (A?A λI)θ A?y 拿來做。加 λI 相當于給所有特征值抬高一個地板條件數(shù)從 κ2 降到 (σ_max2 λ)/(σ_min2 λ)小奇異值被壓住了收斂立刻變快。代價是解被往零的方向拉了一點圖像會略偏平滑——但這恰好和壓縮感知要的「抑制噪聲」方向一致很多時候不是壞事。第三種是 CGLS它不顯式構(gòu)造 A?A而是直接在 A 上做 Krylov 迭代每一步只做一次 A 乘和一次 A? 乘。數(shù)值上比正規(guī)方程穩(wěn)得多條件數(shù)不會平方放大。代價是實現(xiàn)稍復雜需要維護兩套遞推量。寫法求解對象有效條件數(shù)每步矩陣乘法適用場景正規(guī)方程 CGA?Aθ A?yκ22 次小規(guī)模、A 良態(tài)Tikhonov CG(A?AλI)θ A?y(σ2λ) 比值2 次有噪聲、需要正則CGLS / LSQR直接在 A 上迭代κ2 次大規(guī)模、病態(tài)有一處容易翻車在 Matlab 里顯式寫A * A會真的構(gòu)造出一個 n×n 矩陣。n 4096 時是 128 MB勉強能跑n 256×256 65536 時是 34 GB直接內(nèi)存爆炸。這就是下一節(jié)要講的函數(shù)句柄寫法為什么必要。3.3 手寫 CG 與 pcg從顯式矩陣到函數(shù)句柄先看手寫版本它把每一步在干什么攤開調(diào)試時能直接看到殘差曲線function [x, relres, k] cg_solve(Afun, b, x, tol, maxit) % Afun : 對稱正定算子句柄輸入列向量返回 A*v % b : 右端項 (n x 1) % x : 初值同時作為輸出 % tol : 相對殘差停止閾值常用 1e-6 % maxit : 最大迭代步數(shù) 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 會炸成天文數(shù)字圖直接變成一片雪花。加一句 warning 至少能讓你知道是自己參數(shù)設錯了而不是算法爛。實際工程里我更傾向直接用 Matlab 自帶的pcg把算子寫成匿名函數(shù)永遠不顯式構(gòu)造 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 是預條件子病態(tài)3 是迭代停滯4 是某步出現(xiàn)非有限值。調(diào)試時如果看到 flag 1別急著加迭代次數(shù)上限先看條件數(shù)——十有八九是 λ 太小或者干脆沒加?？吹?flag 4多半是 Phi 里有 Inf 或者 NaN往上游查數(shù)據(jù)歸一化。3.4 把 CG 塞進 L1IRLS 外層循環(huán)怎么搭到這一步必須把話說透CG 解的是二次型問題而壓縮感知的目標函數(shù)帶 L1 項不是二次的。CG 不能直接解壓縮感知它只能當內(nèi)層的一個線性求解器。這一層關系如果搞混調(diào)參會調(diào)得非常迷茫。把 L1 和 CG 接起來的經(jīng)典做法是迭代重加權(quán)最小二乘IRLS。核心思想是|θ_i| ≈ θ_i2 / |θ_i|用當前估計的 |θ_i| 當權(quán)重把 L1 問題轉(zhuǎn)成一串加權(quán) L2 問題每個 L2 子問題用 CG 解。lambda 1e-3; % 正則強度越大越稀疏 eps_reg 1e-6; % 防止 1/0常用 1e-6 ~ 1e-4 theta ones(n, 1); % 初值不要全零否則第一步權(quán)重會爆 Aop (v) Phi * (Phi * v); rhs Phi * y; for outer 1:15 w 1 ./ (abs(theta) eps_reg); % 越接近 0 的解權(quán)重越大 Rop (v) lambda * (w .* v); % 加權(quán) 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; % 回到圖像域三處參數(shù)值得單獨說。eps_reg決定了「多小算零」設太大會把真實的小系數(shù)一起壓掉圖像細節(jié)丟失設太小會讓權(quán)重矩陣條件數(shù)爆表內(nèi)層 CG 收斂極慢。經(jīng)驗值 1e-6 起步如果內(nèi)層 flag 經(jīng)常是 1就往 1e-4 調(diào)。lambda控制整體稀疏強度量級通常要和norm(Phi*y, inf)對齊看單純照抄別人的 1e-3 在你自己數(shù)據(jù)上很可能失效。外層循環(huán)次數(shù) 10~20 次足夠再往上收益遞減可以直接看非零元個數(shù)有沒有穩(wěn)定。注意這套 IRLS 只是眾多 L1 求解路線里的一種。ADMM 把它拆成「數(shù)據(jù)保真 軟閾值」兩個子問題交替解FISTA 用梯度步加軟閾值再加動量加速TV 正則則要在 L1 之外再套一層梯度算子。選哪條路線取決于你的先驗是「稀疏」還是「分片光滑」下一章對比的就是這件事。4. 四種重構(gòu)器同場對比采樣率降到多少開始崩4.1 統(tǒng)一實驗骨架掩模、采樣率與評價指標對比算法的第一原則是所有變量只動一個。固定同一張圖、同一個隨機種子、同一個 Φ只換重構(gòu)器否則你根本分不清差異來自算法還是來自掩模。骨架大概是這樣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; % 在這里依次調(diào)用 CG-Tikhonov / OMP / FISTA / TV-ADMM % 每個算法跑完記錄 psnr(x, x_rec)、ssim 和 toc end指標選 PSNR 和 SSIM 兩個就夠。PSNR 對整體亮度誤差敏感SSIM 對結(jié)構(gòu)相似度敏感兩者經(jīng)常給出不同排序——一個算法 PSNR 高但 SSIM 低通常意味著它把噪聲抹平了但邊緣也磨掉了。只看一個指標下結(jié)論容易得出偏頗的結(jié)論。還有一條必須統(tǒng)一每個算法的 λ 都要在自己的參數(shù)網(wǎng)格上取最優(yōu)再拿來對比。拿一個沒調(diào)參的 FISTA 去比一個調(diào)好的 CG這種對比沒有意義。代價是實驗量翻好幾倍但這是唯一能讓結(jié)論站住的做法。4.2 重建質(zhì)量隨采樣率的變化誰先崩趨勢大致是這樣的。在 50% 采樣率附近四種方法的 PSNR 差距不明顯往往在 1~2 dB 以內(nèi)因為信息本來就夠多什么方法都能湊出來。差距是從 20% 往下拉開的。純 CG 加 Tikhonov 屬于 L2 正則它給出的解在所有方向上都平滑本質(zhì)是「最小能量解」稀疏性完全沒有被利用。采樣率降到 20% 以下重建圖會明顯偏糊邊緣發(fā)軟紋理區(qū)域糊成一片。它的優(yōu)勢是快——迭代次數(shù)少每步只做兩次矩陣乘法。OMP 走的是另一條路貪心地每次挑一個與當前殘差最相關的原子然后在這一小撮原子上做最小二乘。它在極度稀疏K 很小且 M 不太小時效果很好但有兩個軟肋。一是需要預先知道稀疏度 K猜大了會過擬合噪聲猜小了會丟結(jié)構(gòu)二是每一步都要算 A?r復雜度隨原子數(shù)線性漲K 大時非常慢。FISTA 是 L1 正則的主力選手梯度步加軟閾值再疊動量收斂速率比 ISTA 快一個量級。它在 10%~30% 采樣率區(qū)間表現(xiàn)穩(wěn)定是中低采樣率下的性價比之選。軟閾值的副作用是會把小系數(shù)整體往零拉圖像會有輕微的「塊狀感」這是 L1 的固有偏置不是實現(xiàn) bug。TV 正則ADMM 或增廣拉格朗日求解在低采樣率下通常畫質(zhì)最好因為圖像的分片光滑先驗比純稀疏先驗更貼合自然圖像。代價是慢而且慢得不只是常數(shù)倍——它每次迭代內(nèi)部還要再解一個子問題整體比 FISTA 慢一個數(shù)量級很常見。提示這里的排序是趨勢層面的具體到某一張圖、某一組掩模名次可能互換。尤其是紋理極豐富的圖像TV 的分片光滑假設會失效反而不如 FISTA。4.3 運行時間與內(nèi)存TV 類方法慢在哪時間差異的量級值得單獨拎出來看因為它直接決定了你走不走得通高分辨率。方法先驗類型每次迭代代價典型迭代數(shù)主要失敗模式CG TikhonovL2 平滑2 次矩陣乘20~80低采樣率下糊OMP嚴格稀疏1 次相關 1 次 LSK 步K 估計錯就崩FISTAL1 稀疏2 次矩陣乘 軟閾值100~300塊狀偽影TV ADMM分片光滑內(nèi)層子問題求解100~500極慢、參數(shù)敏感從表里能讀出一個關鍵取舍CG 快但先驗弱TV 先驗強但慢。FISTA 卡在中間這也是它在實際項目里出場率最高的原因。內(nèi)存方面只要堅持用函數(shù)句柄四種方法的峰值內(nèi)存基本都在「Phi 本身 幾個 n 維向量」這個量級。真正的內(nèi)存殺手是顯式構(gòu)造 A?A 或者把 TV 的梯度算子展開成稠密矩陣。n 65536 時一個稠密 n×n 就是 34 GB這個坑踩一次能記一輩子。ADMM 實現(xiàn)里有一步常見的錯誤是「為了寫起來方便」把差分算子寫成diff(eye(n))那一行代碼就是內(nèi)存爆炸的直接原因正確寫法是把它寫成只做鄰域差分的函數(shù)句柄。5. 排查圖糊、圖碎、不收斂時先看哪幾個量5.1 圖像糊、碎、全黑先查信號鏈再查算法現(xiàn)象一重建圖整體形狀對但像蒙了一層霧邊緣發(fā)軟。原因大概率不是算法而是測量矩陣沒有按行歸一化。如果 Φ 的每一行能量不同y 的各個分量量級就會參差不齊最小二乘會自動更重視量級大的那些測量等價于給不同掩模加了不同權(quán)重。解決很直接構(gòu)造完 Phi 之后強制Phi Phi ./ vecnorm(Phi, 2, 2)再檢查vecnorm(Phi,2,2)是不是全 1?，F(xiàn)象二圖碎成馬賽克狀或者有明顯的方塊拼接感。這通常是稀疏基和掩模不匹配。用哈達瑪掩模配小波基或者用 DCT 基去配本身就帶 DCT 結(jié)構(gòu)的掩模相干性過高恢復出來的解在基里分布得很散重建圖上就表現(xiàn)為均勻的塊狀噪聲。解決是換隨機二值或高斯掩模重跑一次看看塊狀感有沒有消失。現(xiàn)象三重建全黑或者全白。這種八成是 λ 和數(shù)據(jù)量級嚴重失配。檢查一下norm(Phi*y)和lambda的數(shù)量級如果差了三四個數(shù)量級λ 就相當于無窮大解直接被壓到零。另一個可能是 y 里混進了 NaN往上順一遍數(shù)據(jù)讀入有沒有除零?，F(xiàn)象四低采樣率下某個固定區(qū)域永遠重建不出來。這往往是掩模覆蓋不均。隨機掩模雖然整體不相干但 M 很小時局部可能恰好沒有充分覆蓋某些位置尤其圖像邊角。解決是換一組種子多跑幾次取平均或者用結(jié)構(gòu)化隨機把區(qū)域分塊后每塊內(nèi)獨立隨機保證空間覆蓋均勻。5.2 不收斂、內(nèi)存爆、每次結(jié)果不同查參數(shù)與隨機源現(xiàn)象五pcg 返回 flag 1迭代跑滿還沒收斂。先看兩個量condest估一下有效算子的條件數(shù)以及 λ 到底設了多少。不加正則或 λ 極小的時候A?A 的最小奇異值接近零CG 收斂因子接近 1當然跑不滿。解決辦法是加 λ從 1e-3 起試或者換 CGLS/LSQR 這類不平方條件數(shù)的方法?，F(xiàn)象六跑到一半內(nèi)存爆掉。基本就兩種可能顯式構(gòu)造了 A?A或者 TV 的差分算子被展開成了稠密矩陣。前者改成函數(shù)句柄后者改成鄰域差分。還有個隱蔽的來源是spdiags用的位置不對——IRLS 里的權(quán)重矩陣 W 一定要用spdiags生成稀疏對角陣寫成diag(w)就是又一個 n×n 稠密矩陣?，F(xiàn)象七同一份代碼兩次運行結(jié)果差異很大。壓縮感知流程里至少有三個隨機源掩模生成、IRLS 初值、某些實現(xiàn)里的隨機打亂。在腳本開頭統(tǒng)一rng(0)并且在 IRLS 里固定初值結(jié)果就能復現(xiàn)。這一條看著瑣碎但凡是做算法對比不復現(xiàn)的結(jié)果等于沒做。注意rng要放在每次生成隨機量之前不能只放在腳本開頭一次。Matlab 的隨機狀態(tài)是全局的中間任何一處調(diào)用了隨機函數(shù)都會把它推進后面生成的東西就變了。現(xiàn)象八換成更高分辨率后代碼能跑但圖明顯變差。除了內(nèi)存和迭代次數(shù)還有個常被忽略的點采樣率的分母是 n分辨率一升同樣的 M 對應的采樣率就掉下來了。64×64 用 1000 次測量是 24%256×256 用同樣的 1000 次只有 1.5%完全不在一個難度級別。這類問題的排查方式就是打印一下實際的 M/n別憑感覺。6. 設計取舍先驗、采樣預算與驗證成本的三方博弈到這里可以把整條鏈路的選擇收斂成一個三角形采樣預算、先驗強度、驗證成本三者互相拉扯沒有一種重構(gòu)器能同時拿下三邊。采樣預算寬裕時M/n 在 30% 以上L2 加 Tikhonov 就夠用實現(xiàn)簡單、跑得快、參數(shù)只有一個 λ調(diào)試成本極低。這個區(qū)間里上 TV 或 L1 屬于過度設計畫質(zhì)提升 1~2 dB代價是代碼復雜度和運行時間翻幾倍。很多人一上來就上 TV其實是被「TV 更強」這句話帶跑了沒算過性價比。采樣預算緊的時候M/n 在 10% 以下先驗就是命。這時候純 L2 基本沒法看必須用 L1 或 TV。選 L1 還是 TV取決于圖像內(nèi)容自然場景、醫(yī)學影像這類分片光滑的TV 占優(yōu)稀疏點目標、編碼孔徑這類本身就是稀疏的L1 更合適。選錯先驗的后果比不調(diào)參數(shù)還嚴重。6.1 一套流程是否可信的最小驗證每次改完重構(gòu)代碼我一般會跑一個三板斧驗證全程不到五分鐘但能擋掉絕大多數(shù)低級錯誤。驗證項做法期望結(jié)果正向一致性用真實 x 算 y再用無噪聲的 CG 反解殘差應接近機器精度圖與原圖幾乎一致采樣率下降曲線50% → 20% → 10% → 5%固定其他變量PSNR 單調(diào)下降不應出現(xiàn)反彈參數(shù)敏感性λ 取 1e-5 到 1e0 掃一遍應存在明顯的最優(yōu)區(qū)間曲線呈單峰第一項最容易被跳過但它能擋掉最致命的一類錯誤正向模型和反向重構(gòu)用的是同一套 Φ 嗎有沒有在中間不小心轉(zhuǎn)置了有一次我把Phi * Phi里的轉(zhuǎn)置漏掉一個殘差怎么都降不下去查了一下午才發(fā)現(xiàn)在那一行。第二項能擋掉過擬合式的誤判——如果 5% 采樣率下的 PSNR 反而比 10% 高那一定是某處有信息泄漏。第三項的參數(shù)掃描值得多說一句。λ 的最優(yōu)區(qū)間通常在兩個數(shù)量級以內(nèi)如果掃出來曲線是平的或者一直往上走說明你的 λ 量綱和數(shù)據(jù)量級沒對齊得先做歸一化再掃。這一步做完再去做跨算法對比結(jié)論才站得住。從那以后我每次接手一份單像素成像的重構(gòu)代碼第一件事不是看它用了什么高級算法而是先打印出 Φ 的行范數(shù)、y 的均值方差、λ 的實際數(shù)量級這三個數(shù)。這三個數(shù)對不上后面所有調(diào)參都是白費。希望幫到你。本文還有配套的精品資源點擊獲取 單像素成像壓縮感知重構(gòu)：測量矩陣、稀疏正則與共軛梯度實戰(zhàn)
簡介這份Matlab單像素成像算法比較代碼包面向已掌握基本成像原理、希望系統(tǒng)對比不同重建方法的讀者圍繞單像素成像中的稀疏采樣與圖像重建問題提供DGI、梯度下降、共軛梯度下降、Poisson最大似然、交替投影、稀疏表示壓縮感知及全變分壓縮感知等七類方法的可運行示例。包內(nèi)共10個文件以9個.m腳本和1個txt說明為主腳本分別對應各算法實現(xiàn)及誤差計算、演示入口說明文件用于輔助理解調(diào)用關系壓縮包約12KB體量輕便便于快速查閱與二次修改。目前已有2762人學習適合課程實驗、算法復現(xiàn)和不同重建策略的初步對照。讀者可借此理解各方法在采樣率、重建質(zhì)量與計算開銷上的差異并基于統(tǒng)一誤差函數(shù)開展結(jié)果比較為單像素成像實驗選型與后續(xù)改進提供參考。1. 一個探測器的悖論為什么幾百次測量就能還原一張圖一份典型的 Matlab 單像素成像腳本擺在面前圖像 64×64 共 4096 個像素測量次數(shù)只給 800 次重構(gòu)那一段調(diào)的是共軛梯度下降。按線性代數(shù)直覺這是嚴重欠定的方程組解不唯一按壓縮感知的說法這叫欠采樣重構(gòu)聽起來像玄學。兩個東西拼在一起居然能出圖這件事本身值得停一秒想清楚到底是哪一方在補償缺失的那些方程如果以為共軛梯度下降天生就能做稀疏重構(gòu)那后面對參數(shù)的所有直覺都會跑偏。把這條鏈路拆開看過一次之后mask 怎么選、正則項怎么調(diào)、為什么有人跑 20 次迭代就收斂而有人跑到 500 次還在震蕩就都有解釋。2. 從光強到方程單像素成像的測量矩陣是怎么長出來的2.1 掩模調(diào)制與標量測量一步物理怎么變成一行矩陣單像素成像的硬件鏈路其實短得可憐光源照到場景上反射光經(jīng)過一組掩模圖案調(diào)制再由一個聚光鏡把所有能量收進單像素探測器。探測器沒有空間分辨能力它讀到的只是一個數(shù)——這一幀的總光強。想要二維圖像就必須在時間維度上換掩模換一次圖案讀一個數(shù)換 M 次得到 M 個數(shù)。關鍵在于每次讀到的這個數(shù)本質(zhì)上是場景光強分布與當前掩模的內(nèi)積。把場景離散成向量 x長度 n N×N第 i 次測量用的掩模離散成向量 φ_i同樣長度 n那么探測器讀數(shù) y_i ≈ Σ_j φ_ij · x_j。M 次測量堆起來就是 y ΦxΦ 是 M×n 的測量矩陣每一行對應一次掩模曝光。一段最樸素的 Matlab 正向仿真長這樣N 64; % 圖像邊長 n N * N; % 未知量個數(shù) 4096 img double(imread(target.png)); x img(:); % n x 1 列向量 M round(0.25 * n); % 采樣率 25%即 1024 次測量 rng(0); % 固定隨機源否則每次跑結(jié)果都不同 Phi randn(M, n) / sqrt(M); % 高斯隨機測量矩陣按行歸一 y Phi * x; % 探測器讀到的 M 個標量這段代碼里有兩處容易被忽略但影響極大。第一處是/ sqrt(M)把每一行的能量歸一到 1是為了讓 y 的量級與 x 處在同一水平否則后面正則項系數(shù) λ 的量綱會跟著采樣率一起漂移換個 M 就得重新調(diào)參。第二處是rng(0)隨機掩模一旦不固定種子同一份代碼兩次運行的重建質(zhì)量能差出好幾個 dB調(diào)試時會誤以為是算法本身不穩(wěn)定。這條鏈路還有個物理上的隱含條件常被跳過掩模刷新與探測器采樣必須同步。真實系統(tǒng)里 DMD 的刷新有建立時間探測器積分窗如果沒對齊讀到的就是兩張掩模的混合結(jié)果等效于 Φ 的每一行都被悄悄污染了一點點。仿真里看不出來一上硬件就是「圖總是糊一層」的源頭。2.2 稀疏基與不相干性為什么隨機掩模比哈達瑪掩模更容易命中有了 y Φx 還不夠M n 時它有無數(shù)個解。壓縮感知的解法是把 x 換到另一個域里表示x Ψθ其中 Ψ 是稀疏基DCT、小波、或者圖像本身θ 只有 K 個顯著非零元K n。代進去得到 y ΦΨθ AθA 是 M×n 的傳感矩陣。這里有個非常容易被跳過的約束Φ 和 Ψ 必須「不相干」。相干性的定義是把兩個基里最相似的一對原子的內(nèi)積挑出來再乘上 √n 歸一。相干性越高一個基里的稀疏信號在另一個基里就越「鋪得開」壓縮感知的恢復保證就越弱。這條約束直接決定了掩模怎么選。如果稀疏基用的是哈達瑪變換而測量掩模也用哈達瑪圖案兩個基完全重合相干性等于最大值這時候壓縮感知的理論保證直接失效——你等于在稀疏域里逐點采樣只是自己沒意識到而已。換成伯努利二值掩?；蛘吒咚闺S機掩模與 DCT、小波這些常見基的相干性都很低恢復才穩(wěn)。血淚經(jīng)驗是很多人第一次做單像素成像會選哈達瑪掩模因為它二值、好實現(xiàn)、有快速變換、重建還能用逆變換一步到位。它在高采樣率下效果確實漂亮但一旦把采樣率壓到 10% 以下重構(gòu)質(zhì)量的下降速度明顯比隨機掩模陡。原因不神秘哈達瑪掩模壓到低采樣率時是「均勻地少采」而隨機掩模是「不相干地亂采」后者的信息冗余更低。2.3 從 L0 到 L1凸松弛到底換來了什么理論上最想要的是 L0 范數(shù)最小在所有滿足 y Aθ 的解里挑非零元最少的那個。問題是這屬于組合優(yōu)化n 上千時窮舉不現(xiàn)實。L1 范數(shù)是 L0 在單位球上的凸包把它替換進去問題變成可解的凸優(yōu)化min ||θ||? s.t. y Aθ幾何上的直覺是L1 球在坐標軸方向有尖角當殘差約束的可行域高維空間里的一個仿射子空間與這個球相切時切點大概率落在坐標軸上也就是稀疏散點上而 L2 球是光滑的切點會落在球面上解里幾乎每個分量都不為零圖像就糊了。這就是為什么純最小二乘解出來的圖總是「整體形狀對、細節(jié)全丟」。實際有噪聲時等式約束要放松成不等式或懲罰形式min ?||y ? Aθ||?2 λ||θ||?到這里才出現(xiàn)了這份代碼里真正要解的東西。λ 大稀疏性強、解偏平滑、細節(jié)丟λ 小擬合狠、容易把噪聲也當信號擬合進去。這個式子不是二次的共軛梯度下降沒法直接吃它后面第 3 章要拆的就是這層關系。2.4 噪聲進來之后理想 RIP 在實驗臺上永遠不會成立教科書上的恢復保證依賴 RIP有限等距性質(zhì)對任意 K 稀疏向量 θ要求 (1?δ)||θ||2 ≤ ||Aθ||2 ≤ (1δ)||θ||2δ 足夠小。這個條件在仿真里勉強成立在真實裝置上幾乎不可能滿足。干擾源至少有四類。探測器有暗電流和讀出噪聲低光強下還疊加泊松噪聲表現(xiàn)為 y 的方差隨信號強度變化不是簡單的高斯白噪聲。掩模與場景的相對位置有微小位移等效于 A 的每一行都在抖。光源強度隨時間漂移等效于 y 整體乘了個慢變系數(shù)。掩模本身的二值化不理想灰度掩模的透過率不是嚴格的 0 和 1。這四類里最坑的是光源漂移因為它不改變圖像形狀只是讓 y 的整體尺度慢慢變重構(gòu)出來的圖看起來「形狀對、對比度不對」很容易被誤判成 λ 沒調(diào)好。常見的處理思路是在每一批測量里插一幀全通掩模所有元素為 1用它的讀數(shù)做歸一化基準把慢變系數(shù)除掉。這招在論文里不顯眼在實驗臺上能省掉大量無效調(diào)參。3. 共軛梯度下降它解的到底是不是壓縮感知那個問題3.1 Krylov 子空間與 A-共軛CG 比最速下降快在哪共軛梯度下降解決的是這類問題A 對稱正定求 Ax b。它的迭代不是在原始空間里沿負梯度走而是在一個逐步擴張的子空間里找最優(yōu)解這個子空間叫 Krylov 子空間 K_k(A, r?) span{r?, Ar?, A2r?, …}。每一步都把當前殘差正交化到之前所有搜索方向上于是搜索方向兩兩 A-共軛p_i? A p_j 0i ≠ j。這個性質(zhì)帶來的效果非常直觀。最速下降每步都沿當前負梯度走相鄰兩步的方向幾乎正交路徑是鋸齒狀的條件數(shù)大時收斂極慢。CG 因為方向互相 A-共軛不會重復走已經(jīng)走過的方向理論上在精確算術(shù)下 n 步內(nèi)必收斂——盡管浮點運算下遠達不到但實際往往幾十步就夠了。收斂速度有個常用估計誤差的 A-范數(shù)每步至少乘上 (√κ ? 1)/(√κ 1)κ 是 A 的條件數(shù)。這句話信息量很大條件數(shù) 100 時這個因子約 0.82收斂慢條件數(shù) 10 時約 0.52幾十步就夠。所以對 CG 來說比「迭代多少次」更重要的是「矩陣條件數(shù)是多少」。這也解釋了后面為什么一定要加正則項。3.2 正規(guī)方程、Tikhonov 與 CGLS三種寫法的數(shù)值差異最小二乘問題 min ||Aθ ? y||2 的最優(yōu)性條件是 A?Aθ A?y這就是正規(guī)方程。A?A 對稱半正定可以用 CG 解。但這里埋了個雷A?A 的條件數(shù)是 A 的條件數(shù)的平方。A 本身條件數(shù) 50 還算溫和平方之后 2500CG 的收斂因子就變成 0.96 左右跑幾百步都未必收斂。第二種寫法是 Tikhonov 正則把 (A?A λI)θ A?y 拿來做。加 λI 相當于給所有特征值抬高一個地板條件數(shù)從 κ2 降到 (σ_max2 λ)/(σ_min2 λ)小奇異值被壓住了收斂立刻變快。代價是解被往零的方向拉了一點圖像會略偏平滑——但這恰好和壓縮感知要的「抑制噪聲」方向一致很多時候不是壞事。第三種是 CGLS它不顯式構(gòu)造 A?A而是直接在 A 上做 Krylov 迭代每一步只做一次 A 乘和一次 A? 乘。數(shù)值上比正規(guī)方程穩(wěn)得多條件數(shù)不會平方放大。代價是實現(xiàn)稍復雜需要維護兩套遞推量。寫法求解對象有效條件數(shù)每步矩陣乘法適用場景正規(guī)方程 CGA?Aθ A?yκ22 次小規(guī)模、A 良態(tài)Tikhonov CG(A?AλI)θ A?y(σ2λ) 比值2 次有噪聲、需要正則CGLS / LSQR直接在 A 上迭代κ2 次大規(guī)模、病態(tài)有一處容易翻車在 Matlab 里顯式寫A * A會真的構(gòu)造出一個 n×n 矩陣。n 4096 時是 128 MB勉強能跑n 256×256 65536 時是 34 GB直接內(nèi)存爆炸。這就是下一節(jié)要講的函數(shù)句柄寫法為什么必要。3.3 手寫 CG 與 pcg從顯式矩陣到函數(shù)句柄先看手寫版本它把每一步在干什么攤開調(diào)試時能直接看到殘差曲線function [x, relres, k] cg_solve(Afun, b, x, tol, maxit) % Afun : 對稱正定算子句柄輸入列向量返回 A*v % b : 右端項 (n x 1) % x : 初值同時作為輸出 % tol : 相對殘差停止閾值常用 1e-6 % maxit : 最大迭代步數(shù) 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 會炸成天文數(shù)字圖直接變成一片雪花。加一句 warning 至少能讓你知道是自己參數(shù)設錯了而不是算法爛。實際工程里我更傾向直接用 Matlab 自帶的pcg把算子寫成匿名函數(shù)永遠不顯式構(gòu)造 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 是預條件子病態(tài)3 是迭代停滯4 是某步出現(xiàn)非有限值。調(diào)試時如果看到 flag 1別急著加迭代次數(shù)上限先看條件數(shù)——十有八九是 λ 太小或者干脆沒加?？吹?flag 4多半是 Phi 里有 Inf 或者 NaN往上游查數(shù)據(jù)歸一化。3.4 把 CG 塞進 L1IRLS 外層循環(huán)怎么搭到這一步必須把話說透CG 解的是二次型問題而壓縮感知的目標函數(shù)帶 L1 項不是二次的。CG 不能直接解壓縮感知它只能當內(nèi)層的一個線性求解器。這一層關系如果搞混調(diào)參會調(diào)得非常迷茫。把 L1 和 CG 接起來的經(jīng)典做法是迭代重加權(quán)最小二乘IRLS。核心思想是|θ_i| ≈ θ_i2 / |θ_i|用當前估計的 |θ_i| 當權(quán)重把 L1 問題轉(zhuǎn)成一串加權(quán) L2 問題每個 L2 子問題用 CG 解。lambda 1e-3; % 正則強度越大越稀疏 eps_reg 1e-6; % 防止 1/0常用 1e-6 ~ 1e-4 theta ones(n, 1); % 初值不要全零否則第一步權(quán)重會爆 Aop (v) Phi * (Phi * v); rhs Phi * y; for outer 1:15 w 1 ./ (abs(theta) eps_reg); % 越接近 0 的解權(quán)重越大 Rop (v) lambda * (w .* v); % 加權(quán) 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; % 回到圖像域三處參數(shù)值得單獨說。eps_reg決定了「多小算零」設太大會把真實的小系數(shù)一起壓掉圖像細節(jié)丟失設太小會讓權(quán)重矩陣條件數(shù)爆表內(nèi)層 CG 收斂極慢。經(jīng)驗值 1e-6 起步如果內(nèi)層 flag 經(jīng)常是 1就往 1e-4 調(diào)。lambda控制整體稀疏強度量級通常要和norm(Phi*y, inf)對齊看單純照抄別人的 1e-3 在你自己數(shù)據(jù)上很可能失效。外層循環(huán)次數(shù) 10~20 次足夠再往上收益遞減可以直接看非零元個數(shù)有沒有穩(wěn)定。注意這套 IRLS 只是眾多 L1 求解路線里的一種。ADMM 把它拆成「數(shù)據(jù)保真 軟閾值」兩個子問題交替解FISTA 用梯度步加軟閾值再加動量加速TV 正則則要在 L1 之外再套一層梯度算子。選哪條路線取決于你的先驗是「稀疏」還是「分片光滑」下一章對比的就是這件事。4. 四種重構(gòu)器同場對比采樣率降到多少開始崩4.1 統(tǒng)一實驗骨架掩模、采樣率與評價指標對比算法的第一原則是所有變量只動一個。固定同一張圖、同一個隨機種子、同一個 Φ只換重構(gòu)器否則你根本分不清差異來自算法還是來自掩模。骨架大概是這樣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; % 在這里依次調(diào)用 CG-Tikhonov / OMP / FISTA / TV-ADMM % 每個算法跑完記錄 psnr(x, x_rec)、ssim 和 toc end指標選 PSNR 和 SSIM 兩個就夠。PSNR 對整體亮度誤差敏感SSIM 對結(jié)構(gòu)相似度敏感兩者經(jīng)常給出不同排序——一個算法 PSNR 高但 SSIM 低通常意味著它把噪聲抹平了但邊緣也磨掉了。只看一個指標下結(jié)論容易得出偏頗的結(jié)論。還有一條必須統(tǒng)一每個算法的 λ 都要在自己的參數(shù)網(wǎng)格上取最優(yōu)再拿來對比。拿一個沒調(diào)參的 FISTA 去比一個調(diào)好的 CG這種對比沒有意義。代價是實驗量翻好幾倍但這是唯一能讓結(jié)論站住的做法。4.2 重建質(zhì)量隨采樣率的變化誰先崩趨勢大致是這樣的。在 50% 采樣率附近四種方法的 PSNR 差距不明顯往往在 1~2 dB 以內(nèi)因為信息本來就夠多什么方法都能湊出來。差距是從 20% 往下拉開的。純 CG 加 Tikhonov 屬于 L2 正則它給出的解在所有方向上都平滑本質(zhì)是「最小能量解」稀疏性完全沒有被利用。采樣率降到 20% 以下重建圖會明顯偏糊邊緣發(fā)軟紋理區(qū)域糊成一片。它的優(yōu)勢是快——迭代次數(shù)少每步只做兩次矩陣乘法。OMP 走的是另一條路貪心地每次挑一個與當前殘差最相關的原子然后在這一小撮原子上做最小二乘。它在極度稀疏K 很小且 M 不太小時效果很好但有兩個軟肋。一是需要預先知道稀疏度 K猜大了會過擬合噪聲猜小了會丟結(jié)構(gòu)二是每一步都要算 A?r復雜度隨原子數(shù)線性漲K 大時非常慢。FISTA 是 L1 正則的主力選手梯度步加軟閾值再疊動量收斂速率比 ISTA 快一個量級。它在 10%~30% 采樣率區(qū)間表現(xiàn)穩(wěn)定是中低采樣率下的性價比之選。軟閾值的副作用是會把小系數(shù)整體往零拉圖像會有輕微的「塊狀感」這是 L1 的固有偏置不是實現(xiàn) bug。TV 正則ADMM 或增廣拉格朗日求解在低采樣率下通常畫質(zhì)最好因為圖像的分片光滑先驗比純稀疏先驗更貼合自然圖像。代價是慢而且慢得不只是常數(shù)倍——它每次迭代內(nèi)部還要再解一個子問題整體比 FISTA 慢一個數(shù)量級很常見。提示這里的排序是趨勢層面的具體到某一張圖、某一組掩模名次可能互換。尤其是紋理極豐富的圖像TV 的分片光滑假設會失效反而不如 FISTA。4.3 運行時間與內(nèi)存TV 類方法慢在哪時間差異的量級值得單獨拎出來看因為它直接決定了你走不走得通高分辨率。方法先驗類型每次迭代代價典型迭代數(shù)主要失敗模式CG TikhonovL2 平滑2 次矩陣乘20~80低采樣率下糊OMP嚴格稀疏1 次相關 1 次 LSK 步K 估計錯就崩FISTAL1 稀疏2 次矩陣乘 軟閾值100~300塊狀偽影TV ADMM分片光滑內(nèi)層子問題求解100~500極慢、參數(shù)敏感從表里能讀出一個關鍵取舍CG 快但先驗弱TV 先驗強但慢。FISTA 卡在中間這也是它在實際項目里出場率最高的原因。內(nèi)存方面只要堅持用函數(shù)句柄四種方法的峰值內(nèi)存基本都在「Phi 本身 幾個 n 維向量」這個量級。真正的內(nèi)存殺手是顯式構(gòu)造 A?A 或者把 TV 的梯度算子展開成稠密矩陣。n 65536 時一個稠密 n×n 就是 34 GB這個坑踩一次能記一輩子。ADMM 實現(xiàn)里有一步常見的錯誤是「為了寫起來方便」把差分算子寫成diff(eye(n))那一行代碼就是內(nèi)存爆炸的直接原因正確寫法是把它寫成只做鄰域差分的函數(shù)句柄。5. 排查圖糊、圖碎、不收斂時先看哪幾個量5.1 圖像糊、碎、全黑先查信號鏈再查算法現(xiàn)象一重建圖整體形狀對但像蒙了一層霧邊緣發(fā)軟。原因大概率不是算法而是測量矩陣沒有按行歸一化。如果 Φ 的每一行能量不同y 的各個分量量級就會參差不齊最小二乘會自動更重視量級大的那些測量等價于給不同掩模加了不同權(quán)重。解決很直接構(gòu)造完 Phi 之后強制Phi Phi ./ vecnorm(Phi, 2, 2)再檢查vecnorm(Phi,2,2)是不是全 1?，F(xiàn)象二圖碎成馬賽克狀或者有明顯的方塊拼接感。這通常是稀疏基和掩模不匹配。用哈達瑪掩模配小波基或者用 DCT 基去配本身就帶 DCT 結(jié)構(gòu)的掩模相干性過高恢復出來的解在基里分布得很散重建圖上就表現(xiàn)為均勻的塊狀噪聲。解決是換隨機二值或高斯掩模重跑一次看看塊狀感有沒有消失。現(xiàn)象三重建全黑或者全白。這種八成是 λ 和數(shù)據(jù)量級嚴重失配。檢查一下norm(Phi*y)和lambda的數(shù)量級如果差了三四個數(shù)量級λ 就相當于無窮大解直接被壓到零。另一個可能是 y 里混進了 NaN往上順一遍數(shù)據(jù)讀入有沒有除零?，F(xiàn)象四低采樣率下某個固定區(qū)域永遠重建不出來。這往往是掩模覆蓋不均。隨機掩模雖然整體不相干但 M 很小時局部可能恰好沒有充分覆蓋某些位置尤其圖像邊角。解決是換一組種子多跑幾次取平均或者用結(jié)構(gòu)化隨機把區(qū)域分塊后每塊內(nèi)獨立隨機保證空間覆蓋均勻。5.2 不收斂、內(nèi)存爆、每次結(jié)果不同查參數(shù)與隨機源現(xiàn)象五pcg 返回 flag 1迭代跑滿還沒收斂。先看兩個量condest估一下有效算子的條件數(shù)以及 λ 到底設了多少。不加正則或 λ 極小的時候A?A 的最小奇異值接近零CG 收斂因子接近 1當然跑不滿。解決辦法是加 λ從 1e-3 起試或者換 CGLS/LSQR 這類不平方條件數(shù)的方法?，F(xiàn)象六跑到一半內(nèi)存爆掉。基本就兩種可能顯式構(gòu)造了 A?A或者 TV 的差分算子被展開成了稠密矩陣。前者改成函數(shù)句柄后者改成鄰域差分。還有個隱蔽的來源是spdiags用的位置不對——IRLS 里的權(quán)重矩陣 W 一定要用spdiags生成稀疏對角陣寫成diag(w)就是又一個 n×n 稠密矩陣?，F(xiàn)象七同一份代碼兩次運行結(jié)果差異很大。壓縮感知流程里至少有三個隨機源掩模生成、IRLS 初值、某些實現(xiàn)里的隨機打亂。在腳本開頭統(tǒng)一rng(0)并且在 IRLS 里固定初值結(jié)果就能復現(xiàn)。這一條看著瑣碎但凡是做算法對比不復現(xiàn)的結(jié)果等于沒做。注意rng要放在每次生成隨機量之前不能只放在腳本開頭一次。Matlab 的隨機狀態(tài)是全局的中間任何一處調(diào)用了隨機函數(shù)都會把它推進后面生成的東西就變了。現(xiàn)象八換成更高分辨率后代碼能跑但圖明顯變差。除了內(nèi)存和迭代次數(shù)還有個常被忽略的點采樣率的分母是 n分辨率一升同樣的 M 對應的采樣率就掉下來了。64×64 用 1000 次測量是 24%256×256 用同樣的 1000 次只有 1.5%完全不在一個難度級別。這類問題的排查方式就是打印一下實際的 M/n別憑感覺。6. 設計取舍先驗、采樣預算與驗證成本的三方博弈到這里可以把整條鏈路的選擇收斂成一個三角形采樣預算、先驗強度、驗證成本三者互相拉扯沒有一種重構(gòu)器能同時拿下三邊。采樣預算寬裕時M/n 在 30% 以上L2 加 Tikhonov 就夠用實現(xiàn)簡單、跑得快、參數(shù)只有一個 λ調(diào)試成本極低。這個區(qū)間里上 TV 或 L1 屬于過度設計畫質(zhì)提升 1~2 dB代價是代碼復雜度和運行時間翻幾倍。很多人一上來就上 TV其實是被「TV 更強」這句話帶跑了沒算過性價比。采樣預算緊的時候M/n 在 10% 以下先驗就是命。這時候純 L2 基本沒法看必須用 L1 或 TV。選 L1 還是 TV取決于圖像內(nèi)容自然場景、醫(yī)學影像這類分片光滑的TV 占優(yōu)稀疏點目標、編碼孔徑這類本身就是稀疏的L1 更合適。選錯先驗的后果比不調(diào)參數(shù)還嚴重。6.1 一套流程是否可信的最小驗證每次改完重構(gòu)代碼我一般會跑一個三板斧驗證全程不到五分鐘但能擋掉絕大多數(shù)低級錯誤。驗證項做法期望結(jié)果正向一致性用真實 x 算 y再用無噪聲的 CG 反解殘差應接近機器精度圖與原圖幾乎一致采樣率下降曲線50% → 20% → 10% → 5%固定其他變量PSNR 單調(diào)下降不應出現(xiàn)反彈參數(shù)敏感性λ 取 1e-5 到 1e0 掃一遍應存在明顯的最優(yōu)區(qū)間曲線呈單峰第一項最容易被跳過但它能擋掉最致命的一類錯誤正向模型和反向重構(gòu)用的是同一套 Φ 嗎有沒有在中間不小心轉(zhuǎn)置了有一次我把Phi * Phi里的轉(zhuǎn)置漏掉一個殘差怎么都降不下去查了一下午才發(fā)現(xiàn)在那一行。第二項能擋掉過擬合式的誤判——如果 5% 采樣率下的 PSNR 反而比 10% 高那一定是某處有信息泄漏。第三項的參數(shù)掃描值得多說一句。λ 的最優(yōu)區(qū)間通常在兩個數(shù)量級以內(nèi)如果掃出來曲線是平的或者一直往上走說明你的 λ 量綱和數(shù)據(jù)量級沒對齊得先做歸一化再掃。這一步做完再去做跨算法對比結(jié)論才站得住。從那以后我每次接手一份單像素成像的重構(gòu)代碼第一件事不是看它用了什么高級算法而是先打印出 Φ 的行范數(shù)、y 的均值方差、λ 的實際數(shù)量級這三個數(shù)。這三個數(shù)對不上后面所有調(diào)參都是白費。希望幫到你。本文還有配套的精品資源點擊獲取