3.2.1 DPS 或 PnPDP 這類算法是否基於原有 SDE 進行修改？#

DPS (Diffusion Posterior Sampling)、DiffPIR (Denoising Diffusion Models for Plug-and-Play Image Restoration)以及其他許多 PnPDP 算法，是基於擴散模型的逆向 SDE/ODE 進行修改的。

根據InverseBench 論文，擴散模型通常被表述為正向擴散和反向去噪過程的隨機微分方程 (SDEs)。擴散模型的訓練目標是學習得分函數 $(\nabla_{\mathbf{x}_t} \log p_t\mathbf{x}_t)$。一旦訓練完成，模型可以通過求解反向 SDE 來生成新樣本。

PnPDP 方法的基本思想就是修改或利用這個反向 SDE，來從後驗分佈 $P(\mathbf{x}|\mathbf{y})$ 而不是先驗分佈 $P(\mathbf{x})$ 中生成樣本。後驗分佈根據貝葉斯定理可以表示為 $P(\mathbf{x}|\mathbf{y}) \propto P(\mathbf{x})P(\mathbf{y}|\mathbf{x})$。其中 $P(\mathbf{x})$ 由擴散先驗建模，而 $P(\mathbf{y}|\mathbf{x})$ 則由前向模型 $G$ 計算。

具體來說，指導式方法 (Guidance-based methods)（如 DPS）通過在反向 SDE 中添加一個似然得分項 $\nabla_{\mathbf{x}_t} \log p_t(\mathbf{y}|\mathbf{x}_t)$ 來實現這一點。這個似然得分項與前向模型 $G$ 相關，通過查詢 $G$ 來估計 $p(\mathbf{y}|\mathbf{x}_0)$。

3.2.2 通用後驗算法是什麼？#

「後驗」本身是一個分佈，而「算法」是從中推斷或採樣的方法。不過，如果我們將其理解為從後驗分佈採樣的通用框架，那麼最常見的是基於馬爾可夫鏈蒙特卡洛 (MCMC) 方法的算法，例如：

1. Langevin Monte Carlo (LMC)：這是一種常見的 MCMC 方法，用於從目標分佈中採樣。它利用目標分佈的梯度（得分函數）來引導採樣過程。在擴散模型的背景下，它與得分匹配 (score matching) 方法緊密相關。
2. Hamiltonian Monte Carlo (HMC)：另一種更複雜的 MCMC 方法，利用哈密頓動力學來減少隨機遊走行為，提高採樣效率。

PnPDP（Plug-and-Play Diffusion Priors）方法的通用框架通常涉及在擴散模型的逆向採樣過程中引入一個數據一致性（或似然）項。核心思想是在每一步去噪時，除了考慮擴散先驗的去噪方向，還要引入一個梯度（或修正項）來引導樣本符合觀測數據。

以下是一個通用PnPDP框架的偽代碼，它基於逆向SDE/ODE的離散化，並融入了數據一致性項。這個框架是大多數PnPDP算法的基礎，包括DPS和DiffPIR。

1
Algorithm: Generic Plug-and-Play Diffusion Prior Inference
2

3
// y:觀測, G:前向模型, D_theta:去噪器, sigma_t:噪聲參數, T:總步數, lambda_t:似然權重
4
Require: y, G, D_theta, sigma_t, T, lambda_t
5

6
// 從標準高斯噪聲初始化 x_T
7
1: x_T = sample_from_prior_noise()
8
// 從總步數 T 倒序迭代到 1 (逆向採樣)
9
2: for t from T down to 1 do
10
// 擴散模型預測當前噪聲 epsilon_theta_pred
11
3:    epsilon_theta_pred = D_theta(x_t, t)
12
// 從噪聲預測中估計乾淨圖像 x_0_pred_from_prior
13
4:    x_0_pred_from_prior = (x_t - (1 - alpha_bar_t)**0.5 * epsilon_theta_pred) / (alpha_bar_t**0.5)
14
// 計算數據一致性損失
15
5:    loss_data_consistency = L_func(G_func(x_0_pred_from_prior), y)
16
// 計算似然項梯度
17
6:    gradient_likelihood = Compute_Gradient_of_Likelihood(x_0_pred_from_prior, y, G_func, L_func)
18
// 將似然梯度引入去噪器的噪聲預測進行修正
19
7:    epsilon_guided = epsilon_theta_pred - lambda_t * gradient_likelihood
20
// 使用修正後的噪聲預測執行一步去噪
21
8:    x_t_minus_1 = D_theta.denoise_step_with_guided_epsilon(x_t, epsilon_guided, t)
22
// 更新到下一個時間步
23
9:    x_t = x_t_minus_1
24
// 結束循環
25
10: end for
26
// 返回最終估計的乾淨圖像
27
11: return x_0_pred_from_prior

算法描述：

此算法描述了 PnPDP 的通用推斷框架，它在預訓練擴散模型的逆向採樣過程中，引入了數據一致性（似然）項來引導圖像重建。

輸入 (Input):

• y: 觀測數據 (例如，黑洞觀測值)
• G: 前向模型，定義了從原始圖像 $x_0$ 到觀測數據 $y$ 的映射，即 $y=G(x_0)$
• D_theta: 預訓練的擴散模型去噪器 (Diffusion model denoiser)，通常是一個像 U-Net 的神經網絡，用於估計噪聲 $\epsilon$ 或直接預測乾淨圖像 $x_0$
• sigma_t: 擴散過程中的噪聲標準差序列 (或其他擴散參數，例如 $\alpha_t$, $\bar{\alpha}_t$)
• T: 總的擴散步數
• lambda_t: 數據一致性項的權重或步長，用於調整先驗與似然之間的平衡，可以隨時間步 $t$ 變化

輸出 (Output):

• x_0_hat: 重建後的圖像 (例如，黑洞圖像)

算法步驟 (Algorithm Steps):

1. 初始化 (Initialization):

   • 從標準高斯噪聲分佈中採樣一個初始噪聲圖像：x_T = sample_from_prior_noise() (通常是 $N(0,I)$)

2. 迭代推斷 (Iterative Inference) - 從 T 倒數到 1:

   • 對於每個時間步 t (從 T 到 1):

     • 步驟 1: 先驗 (Prior) 部分 - 利用擴散模型去噪

       • 使用預訓練的擴散模型去噪器 D_theta 預測當前時間步的噪聲 epsilon_theta_pred:epsilon_theta_pred = D_theta(x_t, t)

       • 從預測的噪聲 epsilon_theta_pred 中，根據擴散模型的前向過程公式 $q(x_t|x_0) = \mathcal{N}(x_t; \sqrt{\bar{\alpha}_t}x_0, (1-\bar{\alpha}_t)I)$，逆推得到一個初步的乾淨圖像估計 x_0_pred_from_prior。其中 $\alpha_t = 1 - \beta_t$ 且 $\bar{\alpha}_t = \prod_{s=1}^t \alpha_s$。

         $x_{0,pred\_from\_prior} = \frac{x_t - \sqrt{1-\bar{\alpha}_t}\epsilon_{\theta,pred}}{\sqrt{\bar{\alpha}_t}}$

         (在程式碼實現中，這通常表示為 x_0_pred_from_prior = (x_t - (1 - alpha_bar_t)**0.5 * epsilon_theta_pred) / (alpha_bar_t**0.5))

     • 步驟 2: 似然 (Likelihood) / 數據一致性 (Data Consistency) 部分

       • 計算當前初步估計的乾淨圖像 x_0_pred_from_prior 經過前向模型 $G$ 轉換後，與實際觀測數據 $y$ 之間的損失。此處 $L$ 代表數據一致性損失函數 (例如，對於黑洞問題，可以是 $\chi^2$ 統計量)。

         $loss\_data\_consistency=L(G(x_{0,pred\_from\_prior}),y)$

         (在程式碼實現中，這通常表示為 loss_data_consistency = L_func(G_func(x_0_pred_from_prior), y))

       • 計算似然項關於 x_0_pred_from_prior 的梯度。這個梯度引導圖像向符合觀測數據的方向修正。

         $gradient\_likelihood=Compute\_Gradient\_of\_Likelihood(x_{0,pred\_from\_prior},y,G,L)$

         (請注意：Compute_Gradient_of_Likelihood 函數的實現取決於前向模型 $G$ 的性質，可能需要自動微分、數值近似或特定求解器。在程式碼中，這可能表示為 gradient_likelihood = Compute_Gradient_of_Likelihood(x_0_pred_from_prior, y, G_func, L_func))

     • 步驟 3: 結合先驗和似然進行更新

       • 將擴散模型預測的噪聲 epsilon_theta_pred 與來自似然項的梯度修正結合起來，得到引導後的噪聲 $\epsilon_{guided}$。lambda_t 作為指導強度，權衡兩者的影響。

         $\epsilon_{guided}=\epsilon_{\theta,pred}−\lambda_t\cdot gradient\_likelihood$

         (在程式碼實現中，這通常表示為 epsilon_guided = epsilon_theta_pred - lambda_t * gradient_likelihood)

       • 使用修正後的引導噪聲 $\epsilon_{guided}$ 執行擴散模型的一步去噪採樣，從 $x_t$ 迭代到 $x_{t-1}$。此步驟的具體公式依賴於所使用的擴散採樣器（例如 DDPM 或 DDIM）。

         $x_{t−1}=D_\theta.denoise\_step\_with\_guided\_epsilon(x_t,\epsilon_{guided},t)$

         (在程式碼實現中，這通常表示為 x_t_minus_1 = D_theta.denoise_step_with_guided_epsilon(x_t, epsilon_guided, t))

       • 更新當前時間步的圖像狀態：x_t = x_t_minus_1。

3. 返回結果 (Return Result):

   • 返回最後一步得到的 x_0_pred_from_prior，作為重建後的最終圖像 x_0_hat。

关键修正點 Compute_Gradient_of_Likelihood：

這是需要重點關注和修改的部分。對於黑洞成像這樣的前向模型 $G$ 是非線性且可能難以直接求導的情況，如何計算 $L(G(x), y)$ 對 $x$ 的梯度是挑戰。

3.2.3 DPS 與 DiffPir 在通用算法上的改動僞代碼#

TIP

DPS 和 DiffPIR 都屬於 PnPDP 框架，但它們在處理數據一致性項和整合方式上有所不同。

1
Algorithm: Diffusion Posterior Sampling (DPS)
2
// y:觀測, G:前向模型, D_theta:去噪器, alpha_bar_t_sequence:擴散參數, T:總步數, s:指導強度
3
Require: y, G, D_theta, alpha_bar_t_sequence, T, s
4
// 初始化 x_T 為標準高斯噪聲
5
1: x_T = sample_from_prior_noise()
6
// 從總步數 T 倒序迭代到 1
7
2: for t from T down to 1 do
8
// 獲取當前時間步的擴散參數 alpha_bar_t
9
3:    alpha_bar_t = alpha_bar_t_sequence[t]
10
// 擴散模型預測當前噪聲 epsilon_theta_pred
11
4:    epsilon_theta_pred = D_theta(x_t, t)
12
// 從噪聲預測中估計乾淨圖像 x_0_pred
13
5:    x_0_pred = (x_t - (1 - alpha_bar_t)**0.5 * epsilon_theta_pred) / (alpha_bar_t**0.5)
14
// 計算似然損失
15
6:    loss_likelihood = L_func(G_func(x_0_pred), y)
16
// 計算似然損失關於 x_0_pred 的梯度
17
7:    gradient_likelihood_x0 = Compute_Gradient_of_Loss_wrt_x0(loss_likelihood, x_0_pred)
18
// 將似然梯度投影到 x_t 空間
19
8:    gradient_likelihood_xt_for_epsilon = (1 - alpha_bar_t) * gradient_likelihood_x0
20
// 將似然梯度引入去噪器的噪聲預測進行修正
21
9:    epsilon_guided = epsilon_theta_pred - s * gradient_likelihood_xt_for_epsilon
22
// 使用修正後的噪聲預測執行一步去噪
23
10:   x_t_minus_1 = Standard_Diffusion_Denoise_Step(x_t, epsilon_guided, t, alpha_bar_t_sequence)
24
// 更新到下一個時間步
25
11:   x_t = x_t_minus_1
26
// 結束循環
27
12: end for
28
// 返回最終估計的乾淨圖像
29
13: return x_0_pred

1
Algorithm: DiffPIR
2

3
// y:觀測, G:前向模型, D_theta:去噪器, alpha_t_sequence:擴散參數, T:總步數, lambda_dc:一致性強度
4
Require: y, G, D_theta, alpha_t_sequence, alpha_bar_t_sequence, T, lambda_dc
5
// 初始化 x_T 為標準高斯噪聲
6
1: x_T = sample_from_prior_noise()
7
// 初始化當前處理的噪聲圖像
8
2: current_x = x_T
9
// 從總步數 T 倒序迭代到 1
10
3: for t from T down to 1 do
11
// 獲取當前時間步的擴散參數 alpha_t
12
4:    alpha_t = alpha_t_sequence[t]
13
// 獲取當前時間步的擴散參數 alpha_bar_t
14
5:    alpha_bar_t = alpha_bar_t_sequence[t]
15
// 擴散模型預測乾淨圖像 (去噪子問題)
16
6:    x_0_denoised_by_prior = D_theta.predict_x0(current_x, t)
17
// 求解數據一致性子問題
18
7:    x_0_data_consistent = Solve_Data_Consistency_Problem_with_Prior(y, G_func, x_0_denoised_by_prior, lambda_dc)
19
// 計算下一步需要添加的噪聲量
20
8:    noise_for_next_step = (1 - alpha_bar_t_sequence[t-1])**0.5 * torch.randn_like(x_0_data_consistent)
21
// 將處理後的 x_0 重新加噪，回到當前擴散步的噪聲水平
22
9:    current_x = (alpha_bar_t_sequence[t-1])**0.5 * x_0_data_consistent + noise_for_next_step
23
// 結束循環
24
10: end for
25
// 返回最終的數據一致性重建結果
26
11: return x_0_data_consistent

3.3.1 針對黑洞前向模型的 Compute_Gradient_of_Likelihood 或 Solve_Data_Consistency_Problem 進行特化和改進：#

• 可微分物理模擬 (Differentiable Physics Simulation):
  • 思想: 將黑洞成像的前向模型 $G$（包括從圖像到可見度再到閉合量的過程）實現為一個完全可微分的模塊。這可能需要用深度學習框架（如PyTorch, TensorFlow）重寫物理模型，或者利用像 JAX 這樣的庫進行自動微分。
  • 優勢: 獲得精確的梯度，避免了近似誤差，可能提高收斂性和圖像質量。
  • 挑戰: 實現複雜物理模型的可微分性非常困難，特別是 VLBI 數據中的非線性閉合量和非加性噪聲。
• 基於優化或數值梯度的似然修正：
  • 思想: 如果 $G$ 無法完全可微分，或者計算成本過高，可以使用數值方法（如有限差分）來近似似然梯度。
  • 改進點:
    • 自適應步長差分：根據當前 $x_t$ 或 $x_0$ 預測值的特性，動態調整差分計算中的步長，以平衡精度和穩定性。
    • 隨機近似梯度 (Stochastic Gradient Approximation)：特別對於非常高的維度，每次迭代都計算完整梯度可能太慢。可以考慮使用蒙特卡洛方法或類似於零階優化 (zeroth-order optimization) 的技術來隨機近似梯度。
    • 迭代優化器集成：在 DiffPIR 的數據一致性子問題中，您可以使用專為黑洞問題優化的迭代算法（如傳統的 CLEAN 或基於稀疏正則化的優化器），作為求解 Solve_Data_Consistency_Problem 的內部循環。
• 非加性噪聲的處理：
  • 思想: 黑洞觀測值中的非加性噪聲（來自增益和相位誤差）使得傳統的 $L2$ 損失不再是最佳選擇。
  • 改進點:
    • 貝葉斯似然建模：更精確地建模 $P(\mathbf{y}∣\mathbf{x})$，考慮到這些噪聲特性。例如，不是簡單的 $||\mathbf{G}(\mathbf{x})-\mathbf{y}||^2$，而是基於泊松分佈、伽馬分佈或其他更適合非加性噪聲的概率模型來構建似然函數。
    • 變分推斷似然：將噪聲參數也視為潛在變量，並在似然項中進行變分推斷。

3.3.2 結合多個先驗的策略：#

• 動態先驗融合：
  • 思想: 不僅是簡單的加權，而是在 PnP 迭代的不同階段或根據當前重建圖像的質量，動態調整通用先驗和領域特定先驗的影響力。
  • 改進點:
    • 信心分數 (Confidence Score)：訓練一個小型網絡來預測當前圖像與哪個先驗的「風格」更匹配，或者哪個先驗在此時更「可靠」，然後動態調整它們的權重。
    • 級聯/多階段PnP：設計一個分階段的生成過程。例如，在粗略階段使用通用先驗進行快速收斂，在精細階段切換到 GRMHD 先驗以捕捉黑洞的細節。

3.3.3 處理多模態後驗：#

• 多軌跡採樣 (Multi-trajectory Sampling)：
  • 思想: 運行多個並行的 PnPDP 採樣軌跡，每個軌跡從不同的隨機初始化開始，或在中間步驟引入隨機分岔，以探索後驗分佈的不同模式。
  • 改進點:
    • 有效的粒子管理：參考 SMC (Sequential Monte Carlo) 方法中的粒子濾波技術，管理和重採樣這些軌跡，確保它們能有效探索多模態空間。
    • 模式聚類與選擇：在生成結束後，對多個結果進行聚類，識別不同的模式，並提供多個有物理意義的重建結果。

3.3.4 針對穩定性條件的集成：#

• 引入約束或正則化：
  • 思想: 針對黑洞成像中潛在的數值穩定性問題（例如，如果將 PDE 求解器作為前向模型的一部分），在損失函數中加入懲罰項，或在優化過程中強制執行某些約束（如 CFL 條件）。
  • 改進點:
    • 可微分約束層：設計一個可微分層，檢查中間圖像是否滿足物理約束，並提供反向梯度來引導採樣避免不穩定區域。
    • 基於預測的穩定性判斷：訓練一個小型預測網絡，判斷當前 $x_t$ 如果通過前向模型是否會導致不穩定，並在此時調整 PnPDP 的步長或策略。

3.4.1 訓練 Prior (先驗)#

• Prior 1: 通用圖像先驗 (CIFAR10)

  • 數據集: CIFAR10 (或其他大型通用圖像數據集，如 ImageNet, LSUN 2)。

  • 目標: 訓練一個 Denoising Diffusion Probabilistic Model (DDPM)來學習通用自然圖像的低級別特徵（如邊緣、紋理）和一般結構 4。這類模型已被證明能生成高質量樣本。

  • 模型架構: 類似於 Ho et al.中使用的 U-Net 6，並結合 Transformer 的位置嵌入和自注意力機制。

  • 訓練目標: 採用 DDPM 的簡化訓練目標$L_{simple}(\theta) := \mathbb{E}_{t, x_0, \epsilon} \left[ \left\| \epsilon - \epsilon_\theta(\sqrt{\bar{\alpha}_t} x_0 + \sqrt{1 - \bar{\alpha}_t} \epsilon, t) \right\|^2 \right]$。

• Prior 2: 領域相關先驗 (TCIR 颱風衛星圖片)

  • 數據集: TCIR (已擁有此颱風衛星圖片數據集)。
  • 目標: 訓練另一個 DDPM 來學習衛星圖片特有的結構和模式，例如雲層的形態、氣旋的螺旋結構等。儘管不是直接的黑洞圖像，但它比 CIFAR10 更接近高空或天文視角的圖像數據，有助於過渡。
  • 模型架構: 同 Prior 1，或根據 TCIR 數據特性進行調整。
  • 訓練目標: 同 Prior 1。

3.4.2 PnP 後驗設計：結合多先驗與數值穩定性修正#

我們將在 PnP 框架的迭代採樣過程中，引入兩種主要修改：

• 多先驗融合策略 (Multi-Prior Fusion Strategy):

  • 思想: 在去噪步驟中，動態地結合兩個預訓練先驗（CIFAR10 先驗的去噪器 $D\_theta,CIFAR$ 和 TCIR 先亞的去噪器 $D\_theta,TCIR$）的輸出。
  • 方法: 可以採用加權平均的方式融合兩個去噪器的噪聲預測 $epsilon\_textpred$。權重可以固定，也可以根據時間步 t 或當前圖像的特徵動態調整。例如，在去噪早期（大 $t$），通用先驗可能更有用；在去噪後期（小 t），領域相關先驗可能更重要。

• 黑洞數值穩定性修正 (Numerical Stability Correction for Black Hole):

  • 問題: 黑洞成像的前向模型 $G$（特別是當它涉及複雜物理模擬或 PDE 求解時，如在InverseBench: Benchmarking Plug-and-Play Diffusion Priors for Inverse Problems in Physical Sciences 論文中提及的「Full Waveform Inversion」和「Navier-Stokes equation」中所述）可能對輸入有嚴格的數值穩定性條件（例如 CFL 條件）。如果傳遞給 $G$ 的圖像 $x_0$ 不符合這些條件，前向模擬可能會失敗或產生無意義的結果。傳統 PnPDP 方法可能未考慮這些穩定性條件。

  • 解決方法:

    1. 穩定性預測模塊 (Stability Prediction Module): 引入一個輕量級的模塊，該模塊能夠在將 $x_{0,pred}$ 傳遞給前向模型 $G$ 之前，預測其是否會導致數值不穩定。這個模塊可以是一個簡單的檢查器（如果穩定性條件可明確表達，如 CFL 條件），或者是一個在穩定/不穩定數據上訓練的分類器。

    2. 梯度修正與懲罰 (Gradient Correction and Penalty):

       • 懲罰項: 在似然損失 $L(G(x_{0,pred}),y)$ 中，添加一個穩定性懲罰項 $L\_stability(x_{0,pred})$。當 $x_{0,pred}$ 不滿足穩定性條件時，這個懲罰項會急劇增大，從而在梯度下降時引導 $x_{0,pred}$ 遠離不穩定區域。

       • 引導策略調整: 如果預測到不穩定，可以調整似然梯度 $gradient\_likelihood$ 的方向或強度，使其優先將 $x_0$ 推向穩定區域，或者降低當前步的學習率/指導強度，避免模型崩潰。

       • 回退機制: 如果在某一步驟中，即使經過修正仍然產生了不穩定的 $x_0$，可以考慮回退到上一步的 $x_t$ 或 $x_0$，並以更小的學習率或不同的隨機噪聲重新嘗試。

    3. 可微分物理約化 (Differentiable Physics Regularization): 若可能，將黑洞物理特性（如亮度分佈、光子環尺寸比例）納入可微分正則化項，直接約束 $x_0$ 的物理合理性，這有助於隱式地提升數值穩定性。

3.4.3 修改後的 PnP 後驗偽代碼#

我們將基於 DPS 的框架進行修改，因為它直接操作得分函數，方便引入額外指導。

新算法: Multi-Prior Guided DPS with Stability Correction (MP-DPS-SC)

1
Algorithm: Multi-Prior Guided DPS with Stability Correction (MP-DPS-SC)
2

3
 // y:觀測數據, G:前向模型, D_theta_CIFAR:CIFAR去噪器, D_theta_TCIR:TCIR去噪器, alpha_bar_t_sequence:擴散參數序列, T:總步數, s:指導強度, w_cifar_func:CIFAR權重函數, w_tcir_func:TCIR權重函數, lambda_stability:穩定性懲罰強度, Is_Stable_func:穩定性判斷函數
4
Require: y, G, D_theta_CIFAR, D_theta_TCIR, alpha_bar_t_sequence, T, s, w_cifar_func, w_tcir_func, lambda_stability, Is_Stable_func
5

6
// 初始化 x_T 為標準高斯噪聲分布 N(0, I)
7
1: x_T = sample_from_prior_noise()
8
// 從總步數 T 倒序迭代到 1 (逆向採樣過程)
9
2: for t = T down to 1 do
10
// 獲取當前時間步的累積擴散參數 alpha_bar_t
11
3:    alpha_bar_t = alpha_bar_t_sequence[t]
12
// CIFAR10 先驗擴散模型預測噪聲 epsilon_theta
13
4:    epsilon_theta_cifar = D_theta_CIFAR(x_t, t)
14
// TCIR 先驗擴散模型預測噪聲 epsilon_theta
15
5:    epsilon_theta_tcir = D_theta_TCIR(x_t, t)
16
// 根據權重函數融合兩個先驗的噪聲預測
17
6:    epsilon_theta_fused = w_cifar_func(t) * epsilon_theta_cifar + w_tcir_func(t) * epsilon_theta_tcir
18
// 從融合後的噪聲預測中估計乾淨圖像 x_0_pred
19
7:    x_0_pred = (x_t - (1 - alpha_bar_t)**0.5 * epsilon_theta_fused) / (alpha_bar_t**0.5)
20
// 初始化穩定性懲罰項為0
21
8:    stability_penalty = 0
22
// 初始化當前數據一致性指導強度
23
9:    current_guidance_strength = s
24
// 檢查當前估計的 x_0_pred 是否滿足前向模型 G 的穩定性條件
25
10:   if not Is_Stable_func(x_0_pred) then
26
// 如果不穩定，計算基於不穩定程度的懲罰值
27
11:      stability_penalty = lambda_stability * Compute_Instability_Measure(x_0_pred)
28
 // 可選步驟：在不穩定情況下降低數據一致性指導強度，以促進穩定性或避免發散
29
12:      current_guidance_strength = s * 0.5
30
13:   end if
31
// 計算總似然損失，包含原始數據一致性項和穩定性懲罰項
32
14:   loss_likelihood = L_func(G_func(x_0_pred), y) + stability_penalty
33
// 計算總似然損失關於 x_0_pred 的梯度
34
15:   gradient_likelihood_x0 = Compute_Gradient_of_Loss_wrt_x0(loss_likelihood, x_0_pred)
35
// 將 x_0 空間的似然梯度投影到 x_t 空間，以便修正 epsilon
36
16:   gradient_likelihood_xt_for_epsilon = (1 - alpha_bar_t) * gradient_likelihood_x0
37
// 將投影後的似然梯度引入融合的噪聲預測，得到引導後的噪聲 epsilon_guided
38
17:   epsilon_guided = epsilon_theta_fused - current_guidance_strength * gradient_likelihood_xt_for_epsilon
39
// 使用引導後的噪聲 epsilon_guided 執行標準擴散模型的一步去噪，得到 x_(t-1)
40
18:   x_t_minus_1 = Standard_Diffusion_Denoise_Step(x_t, epsilon_guided, t, alpha_bar_t_sequence)
41
// 將當前時間步的圖像狀態更新為 x_(t-1)
42
19:   x_t = x_t_minus_1
43
// 結束逆向採樣循環
44
20: end for
45
// 返回最終估計的乾淨圖像 x_0_pred
46
21: return x_0_pred

算法描述：

此算法融合了來自 CIFAR10 和 TCIR 數據集的兩個擴散模型先驗，並在逆向採樣過程中加入了針對黑洞前向模型數值穩定性的修正，以從物理觀測值生成黑洞圖像。

輸入 (Input):

• y: 黑洞的物理觀測值 (例如 VLBI 閉合量)
• G: 黑洞前向模型，定義了從圖像 $x_0$ 到觀測數據 $y$ 的映射 $(y=G(x_0))$
• D_theta_CIFAR: 預訓練的 CIFAR10 擴散模型去噪器
• D_theta_TCIR: 預訓練的 TCIR 擴散模型去噪器
• alpha_bar_t_sequence: 預計算的擴散參數序列 $\bar{\alpha}_t$
• T: 總的擴散步數
• s: 似然指導強度 (超參數)
• w_cifar(t): CIFAR10 先驗的權重函數 (可隨 $t$ 變化，例如 $w\_cifar(t)=t/T$)
• w_tcir(t): TCIR 先驗的權重函數 (可隨 $t$ 變化，例如 $w\_tcir(t)=1−t/T$)
• lambda_stability: 穩定性懲罰項的強度 (超參數)
• Is_Stable(x_0): 判斷 $x_0$ 是否滿足前向模型穩定性條件的函數 (返回布林值)

輸出 (Output):

• x_0_hat: 重建後的黑洞圖像

算法步驟 (Algorithm Steps):

1. 初始化 (Initialization):

   • x_T = sample_from_prior_noise() // 從 $N(0,I)$ 初始化 $x_T$

2. 迭代推斷 (Iterative Inference) - 從 T 倒數到 1:

   • 對於每個時間步 t (從 T 到 1):

     • alpha_bar_t = alpha_bar_t_sequence[t]

     • 步驟 1: 多先驗去噪預測 (Multi-Prior Denoising Prediction)

       • epsilon_theta_cifar = D_theta_CIFAR(x_t, t) // CIFAR10 先驗預測噪聲

       • epsilon_theta_tcir = D_theta_TCIR(x_t, t) // TCIR 先驗預測噪聲

       • epsilon_theta_fused = w_cifar(t) * epsilon_theta_cifar + w_tcir(t) * epsilon_theta_tcir // 融合兩個先驗的噪聲預測

       • 從融合後的噪聲預測中估計乾淨圖像 $x_0$:

         $x_{0,pred} = \frac{x_t - \sqrt{1-\bar{\alpha}_t}\epsilon_\theta}{\sqrt{\bar{\alpha}_t}}$

     • 步驟 2: 似然損失與穩定性修正 (Likelihood Loss with Stability Correction)

       • stability_penalty = 0 // 初始化穩定性懲罰

       • current_guidance_strength = s // 初始化當前指導強度

       • 如果 not Is_Stable(x_0_pred): // 檢查 $x_0$ 預測是否穩定

         • stability_penalty = lambda_stability * Compute_Instability_Measure(x_0_pred) // 計算不穩定性懲罰，例如基於偏離穩定區域的程度
         • current_guidance_strength = s * 0.5 // 可選: 降低指導強度以避免不穩定，或回退重試

       • 計算似然損失 (可選地包含穩定性懲罰):

         $loss\_likelihood=L(G(x_{0,pred}),y)+stability\_penalty$

       • 計算似然損失關於 $x_{0,pred}$ 的梯度:

         $gradient\_likelihood\_x0=Compute\_Gradient\_of\_Loss\_wrt\_x0(loss\_likelihood,x_{0,pred})$

     • 步驟 3: 投影梯度與修正噪聲 (Project Gradient and Correct Epsilon)

       • 將似然梯度從 $x_0$ 空間投影到 $x_t$ 空間:

         $gradient\_likelihood\_xt\_for\_epsilon=(1−\bar{\alpha}_t) \cdot gradient\_likelihood\_x0$

       • 將似然指導項加到融合後的噪聲預測中:

         $\epsilon_{guided}=\epsilon_{\theta,fused}−current\_guidance\_strength \cdot gradient\_likelihood\_xt\_for\_epsilon$

     • 步驟 4: 執行一步去噪 (Perform One Denoising Step)

       • 使用修正後的 $\epsilon_{guided}$ 計算下一步的 $x_{t−1}$:

         x_t_minus_1 = Standard_Diffusion_Denoise_Step(x_t, epsilon_guided, t, alpha_bar_t_sequence)

       • 更新當前時間步的圖像狀態：x_t = x_t_minus_1。

3. 返回結果 (Return Result):

   • 返回最後一步得到的 x_0_pred (或 x_0_pred_from_prior)，作為重建後的最終圖像 x_0_hat。