问题背景
图像反卷积的目标是从一张模糊图像中恢复出清晰的原始图像。在实际测量系统中，由于设备精度有限、响应时间非零等原因，获取的图像总是真实场景的变形版本。以对焦不准的相机为例，一个理想的点在成像后变成了一个扩散的斑点，整张图像就是所有点对应斑点的叠加结果，即发生了卷积模糊。
描述这一过程最简单的数学模型是线性时不变系统。在二维情况下，观测方程为：
其中 x_{n,m} 是原始清晰图像，h 是系统的脉冲响应（即点扩散函数），b_{n,m} 是加性噪声（假设为空间白噪声），y_{n,m} 是观测到的模糊图像。写成矩阵形式就是：
这里 H 是卷积矩阵，具有 Toeplitz 结构。反卷积问题的任务就是根据已知的 \boldsymbol{y} 和 H，估计出 \boldsymbol{x}。
卷积本质上是一个低通滤波过程——测量系统无法跟踪高频信息，高频分量在数据中已经被大幅衰减甚至消失。反卷积需要恢复这些丢失的高频成分，这使得问题非常困难
正则化最小二乘方法
从最小二乘到正则化
直接使用最小二乘估计（即最小化 \|\boldsymbol{y} - H\boldsymbol{x}\|^2）会导致噪声被急剧放大，得到的结果极不稳定。原因在于卷积矩阵 H 在高频方向上接近奇异，求逆时高频噪声会被放大到不可接受的程度。
为了克服这个问题，引入正则化项来约束解的光滑性。最自然的想法是惩罚相邻像素之间的大差异，由此得到正则化最小二乘准则：
其中 D 是一阶差分矩阵，\mu 是正则化参数，控制数据拟合与光滑性之间的折中。该准则关于 \boldsymbol{x} 是二次的，其最小值有闭式解：
\mu 太小则噪声抑制不足，\mu 太大则图像过度平滑、细节丢失。
循环近似与傅里叶域求解
直接计算上述闭式解需要对大矩阵求逆，计算代价极高。为了加速，利用循环矩阵在傅里叶基下可对角化这一性质。具体做法是将 Toeplitz 矩阵 H 近似为循环矩阵 \tilde{H}，当信号长度 N 远大于脉冲响应长度 P 时，这个近似非常精确。
循环矩阵 \tilde{H} 可以分解为：
其中 W 是傅里叶矩阵，\Lambda_h 是对角矩阵，其对角元素就是脉冲响应的 FFT。对差分矩阵 D 做同样的循环近似后，正则化最小二乘解在傅里叶域变为逐元素运算：
其中 \hat{X}_n、Y_n 分别是 \hat{\boldsymbol{x}} 和 \boldsymbol{y} 的傅里叶变换的第 n 个分量，\lambda_h^n 和 \lambda_d^n 分别是 \tilde{H} 和 \tilde{D} 的第 n 个特征值。
贝叶斯框架下的概率建模
正则化最小二乘方法的核心缺陷是需要手动调节 \mu。变分贝叶斯方法通过引入概率框架，可以实现对超参数的自动估计，从而实现非监督反卷积。
似然函数
由于噪声 \boldsymbol{b} 假设为高斯白噪声，给定 \boldsymbol{x} 后，观测数据 \boldsymbol{y} 的条件分布为：
用精度参数 \gamma_b = 1/\sigma_b^2（方差的倒数）改写，协方差矩阵变为 \frac{1}{\gamma_b}I：
先验分布
正则化项 \mu\,\boldsymbol{x}^t D^t D\,\boldsymbol{x} 对应的概率解释是对 \boldsymbol{x} 施加一个零均值的高斯先验分布，其协方差矩阵与 (D^t D)^{-1} 成正比：
同样改写为精度形式，令 \gamma_x = 1/\sigma_x^2：
D\boldsymbol{x} 衡量的是相邻像素的差异，先验分布倾向于让这些差异较小，从而鼓励光滑的解。
超参数的先验与联合后验
\gamma_b 和 \gamma_x 是超参数，为了实现完全贝叶斯估计，需要对它们也施加先验分布。由于图像像素数量很大，数据对超参数有充分的信息量，同时我们不想引入主观偏好，因此选择无信息的 Jeffreys 先验：
综合以上所有部分，图像和超参数的联合后验分布为：
正则化参数 \mu 在贝叶斯框架中隐含在超参数的比值中：
通过估计 \gamma_b 和 \gamma_x，\mu 被自动确定。
变分贝叶斯近似
可分离近似的基本思想
联合后验 p(\mathring{\boldsymbol{x}}, \gamma_b, \gamma_x|\mathring{\boldsymbol{y}}) 的精确计算是不可行的，因为 \boldsymbol{x}、\gamma_b、\gamma_x 之间存在耦合。变分贝叶斯方法的核心思想是用一个可分离的分布族来近似真实后验：
这意味着我们假设 \boldsymbol{x} 与超参数在近似分布下是独立的。通过最小化近似分布与真实后验之间的 KL 散度，可以得到每个因子的最优更新公式。由于使用了共轭先验，各因子的函数形式是已知的：
q(\mathring{\boldsymbol{x}}) 的推导
根据变分贝叶斯的一般公式，\boldsymbol{x} 的最优近似分布通过对联合对数后验关于其他变量的期望得到：
这里省略了 \int q(\gamma_x)\, d\gamma_x = 1，同时由于我们知道 \int \gamma_b\, q(\gamma_b)\, d\gamma_b = \widetilde{\gamma_b}，\gamma_x 同理，因此上式可以简化为
这是一个高斯分布的核，由此可以直接读出协方差和均值。在傅里叶域中，由于 \Lambda_h 和 \Lambda_d 都是对角矩阵，协方差矩阵也是对角的：
将 \widetilde{\gamma_b} 提出后，均值可以化简为：
这个形式和正则化最小二乘解在傅里叶域的表达完全一致，只不过正则化参数 \mu 被替换为 \widetilde{\gamma_x}/\widetilde{\gamma_b}，即超参数期望的比值。
其中超参数的期望由 Gamma 分布的均值给出：
q(\gamma_b) 的推导
对联合对数后验关于 \boldsymbol{x} 和 \gamma_x 求期望：
对 \gamma_x 的积分由于归一化条件直接消去相关项。对 \mathring{\boldsymbol{x}} 的积分需要计算 \|\mathring{\boldsymbol{y}} - \Lambda_h\mathring{\boldsymbol{x}}\|^2 在 q(\mathring{\boldsymbol{x}}) 下的期望。展开这个二次型：
回顾二阶矩公式，对于均值为 m_x、方差为 \sigma_x^2 的随机变量，有
利用这个关系，期望展开后可以识别出 Gamma 分布的标准形式 \gamma^{\alpha-1}e^{-\beta\gamma}，从而得到更新后的参数。经过整理，中间步骤为：
由此读出 Gamma 分布的参数：
第一项 \|\mathring{\boldsymbol{y}} - \Lambda_h m_x\|^2 是当前估计下的残差能量，第二项 \Sigma_x^2\,\Lambda_h^*\Lambda_h 是由 \boldsymbol{x} 的不确定性（方差）传播到观测空间后的贡献。两项之和衡量的是数据与模型之间的总体不匹配程度。
q(\gamma_x) 的推导
类似地，对联合对数后验关于 \boldsymbol{x} 和 \gamma_b 求期望，得到：
更新参数为：
第一项是当前均值估计的差分能量（衡量光滑程度），第二项是方差在差分空间中的贡献。
算法实现
整体流程
变分贝叶斯算法是一个交替迭代过程，每次迭代依次更新三组参数直到收敛。由于所有运算在傅里叶域中都是逐元素操作（对角矩阵），计算效率很高。
初始化设定 \widetilde{\gamma_b} 和 \widetilde{\gamma_x} 的初值为1，以及 Gamma 分布的形状参数 \alpha_b = \alpha_x = N/2（在整个迭代过程中保持不变）。
每次迭代的更新顺序为：首先根据当前的 \widetilde{\gamma_b} 和 \widetilde{\gamma_x} 计算 \boldsymbol{x} 的后验协方差和均值，然后利用更新后的 m_x 和 \Sigma_x^2 计算 \gamma_b 的 Gamma 参数并更新 \widetilde{\gamma_b}，最后用同样的方式更新 \gamma_x 的参数和 \widetilde{\gamma_x}。
核心更新公式与代码实现
协方差矩阵（傅里叶域逐元素求逆）：
SigmaX2 = 1 ./ (gbMean * FourRI2 + g1Mean * FourD);
均值（傅里叶域逐元素运算）：
mX_fourier = gbMean * SigmaX2 .* conj(FourRI) .* data;
\gamma_b 的 Gamma 分布参数 \beta_b，包含残差能量和方差传播两项：
residual = data - FourRI .* mX_fourier;
betaBruit = 0.5 * (sum2(abs(residual).^2) + sum2(FourRI2 .* SigmaX2));
\gamma_x 的 Gamma 分布参数 \beta_x，差分能量在空间域通过二维卷积计算
方差传播项在傅里叶域计算：
DmX_h = conv2(mX, [1, -1], 'same');
DmX_v = conv2(mX, [1; -1], 'same');
normDmX2 = sum2(DmX_h.^2) + sum2(DmX_v.^2);
betaX1 = 0.5 * (normDmX2 + sum2(FourD .* SigmaX2));
超参数期望的更新：
gbMean = alphaBruit / betaBruit;
g1Mean = alphaX1 / betaX1;
实验结果分析
上图展示了对 Donnees1 数据集运行变分贝叶斯算法的结果。左上为输入的模糊图像，右上为经过 100 次迭代后恢复的图像，下方两张图分别是超参数 \gamma_b 和 \gamma_1 随迭代次数的变化曲线。
可见重建图像很好的恢复了人物、摄像机、以及背景的细节，成功恢复了卷积过程中被低通滤波消除的高频信息
\gamma_b 和\gamma_1 的收敛曲线显示超参数快速收敛到稳定值，整个过程无需手动调试超参数，验证了变分贝叶斯方法的非监督特性。
总结
变分贝叶斯方法相比手动调参的正则化最小二乘方法，最大的优势在于超参数的自动估计。算法通过迭代，让噪声精度 \gamma_b 和先验精度 \gamma_x 自适应地从数据中学习，以确定正则化参数 \mu = \widetilde{\gamma_x}/\widetilde{\gamma_b} 。从计算效率角度看，由于全部在傅里叶域操作，每次迭代的复杂度主要由 FFT 决定，与基于 MCMC 采样的贝叶斯方法相比，收敛速度快、计算时间短，重建质量相同。