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