图像反卷积实验报告：正则化最小二乘、MCMC与变分贝叶斯方法

实验背景与问题描述

本实验的目标是消除图像的模糊，即从模糊图像中恢复清晰图像，这是一个反问题

观测方程的数学建模

描述这种变换的最简单模型是线性时不变滤波器，即卷积。系统可以用以下示意图表示：真实图像 x 通过系统 H 后，加上噪声 b，得到观测数据 y。

二维情况下的观测方程为：

其中 x_{n,m} 表示真实图像，y_{n,m} 表示可用的观测数据，h 是系统的脉冲响应，b_{n,m} 假设为白噪声。

这里涉及的是低通滤波器，这导致高频分量在数据中已经消失了，而反问题需要重建这些高频分量，在没有先验信息的情况下，导致反问题很困难。与此同时，由于噪声存在，当恢复被滤波的高频信号时，会导致噪声爆炸

矩阵形式与卷积矩阵

将上述方程矩阵形式如下

矩阵 \boldsymbol{H} 为具有Toeplitz结构的卷积矩阵：

监督方法

普通最小二乘

最简单的思路是普通最小二乘方法，即最小化数据拟合误差：

其解为

在频域中等价于

然而，由于反卷积问题是病态的（条件数不足），普通最小二乘法下的低通滤波器 H(\omega) 在高频处接近零，因此会出现噪音爆炸的问题

正则化准则

为了克服病态问题，使用正则化最小二乘准则：

• 第一项是数据拟合项，衡量重建图像与观测数据的相似性
• 第二项是正则化项，其中 \boldsymbol{D} 是差分算子，用于衡量信号相邻样本的差异，惩罚不平滑的解
• 参数 \mu 控制数据拟合与平滑先验之间的权衡

差分矩阵 \boldsymbol{D} 的形式为：

带有正则化的最小二乘法的闭式解：

循环近似与频域求解

直接计算上述闭式解需要对大矩阵求逆，这个计算量非常大，因此为了加速计算，我们采用循环矩阵在傅里叶基中可对角化的性质。

将Toeplitz矩阵 \boldsymbol{H} 近似为循环矩阵 \tilde{\boldsymbol{H}}：

循环矩阵假设信号是周期性的，首尾相连，当 N \gg P（信号长度远大于脉冲响应长度）时，这种近似对整体重建质量影响很小。

循环矩阵可以在傅里叶基中对角化：

其中 \boldsymbol{W} 是傅里叶矩阵，\boldsymbol{\Lambda}_h 是对角矩阵，其对角元素是 \tilde{\boldsymbol{H}} 的特征值（对循环矩阵第一行做FFT得到）

频域解的推导

将对角化代入正则化最小二乘解，利用傅里叶矩阵的酉性质 \boldsymbol{W}\boldsymbol{W}^\dagger = \boldsymbol{I}，

最终我们可以推导出频域解表达式为：

具体推导过程见附录1

下面是附录1 的内容

推导频域解表达式:

首先计算 \tilde{\boldsymbol{H}}^\dagger\tilde{\boldsymbol{H}}：

同理

因此：

对这个矩阵求逆，利用

和

可得:

将上式乘以

可得

即

附录1结束

记 \hat{\boldsymbol{X}} = \boldsymbol{W}\hat{\boldsymbol{x}} 和 \boldsymbol{Y} = \boldsymbol{W}\boldsymbol{y} 分别为重建图像和观测数据的傅里叶变换，对上式两边左乘 \boldsymbol{W}：

最后可得频域解的表达式

频率传递函数

由于 \boldsymbol{\Lambda}_h 和 \boldsymbol{\Lambda}_d 都是对角矩阵，设：

对角矩阵的共轭转置仍是对角矩阵，对角元素取共轭：

两个对角矩阵相乘，结果仍是对角矩阵，对角元素逐个相乘：

同理：

两者相加：

对角矩阵求逆就是每个对角元素取倒数：

再乘以 \boldsymbol{\Lambda}_h^\dagger：

因此，第 n 个对角元素就是频率传递函数：

由于整个矩阵是对角的，矩阵乘法变成逐元素乘法，重建图像的频域表示为：

写成向量形式：

计算复杂度

使用循环近似后，计算复杂度从 O(N^3) 降低到 O(N\log N)。

二维

扩展到二维时，需要同时考虑水平方向和垂直方向的差分，水平方向差分算子 \boldsymbol{D}_1 计算同一行相邻像素的差值：

垂直方向差分算子 \boldsymbol{D}_2 计算同一列相邻像素的差值

二维正则化项变为两个方向差分的平方和：

这等价于用二维拉普拉斯核进行卷积。拉普拉斯核检测图像中像素与其四邻域均值的偏差，中心像素权重为4，上下左右四个邻居权重各为-1：

二维实现步骤

具体实现步骤如下：

第一步，构造 \boldsymbol{\Lambda}_h，即脉冲响应在 N \times N 个点上的FFT-2D

第二步，构造 \boldsymbol{\Lambda}_d，即拉普拉斯核在 N \times N 个点上的FFT-2D，代码中使用的核为：

第三步，根据频率传递函数公式构造二维频率传递矩阵 \boldsymbol{G}_{\text{MCR}}：

• 其中 (u,v) 是二维频率坐标。

第四步，构造 \boldsymbol{Y}，即观测数据（模糊图像）的FFT-2D。

第五步，计算 \hat{\boldsymbol{X}}，即频率传递矩阵 \boldsymbol{G}_{\text{MCR}} 与 \boldsymbol{Y} 的逐元素乘积：

第六步，通过对 \hat{\boldsymbol{X}} 进行IFFT-2D回到空间域，得到重建图像 \hat{\boldsymbol{x}}。

Donnees1实验结果

当 \mu = 0.001 时，正则化强度不足，高频噪声被严重放大，图像噪声严重。

当 \mu = 1.0 时，正则化过强，图像过度平滑，边缘和细节严重丢失

\mu = 0.01 结果最好，MSE最小，既能够抑制噪声，又能够保证较好的图像清晰度

Donnees2实验结果

Donnees2的观测数据尺寸为252×252，与真实图像256×256不匹配，因此无法计算MSE。但是我们清楚的可见，\mu = 0.01 到 \mu = 0.1 之间的重建效果较好

监督方法的局限性

正则化参数 \mu 的选择很困难，需要长时间手动调试以找到最优值，因此下面引入非监督方法，即让算法自动从数据中学习最优的参数

非监督方法

贝叶斯框架

从优化到概率

正则化最小二乘方法需要手动设定参数 \mu，为了实现非监督估计，引入完全贝叶斯框架。正则化最小二乘可以从贝叶斯角度重新诠释：数据拟合项对应高斯似然函数，正则化项对应图像的高斯先验分布。

假设噪声服从高斯分布 \boldsymbol{b} \sim \mathcal{N}(0, \gamma_b^{-1}\boldsymbol{I})，其中 \gamma_b = 1/\sigma_b^2 是噪声精度（方差的倒数），那么观测数据 \boldsymbol{y} 的分布是以 \boldsymbol{H}\boldsymbol{x} 为均值的高斯分布，似然函数为：

先验分布

假设图像是平滑的，图像的先验分布设为：

• 其中 \gamma_1 = 1/\sigma_1^2 是先验精度，控制图像的平滑程度

超参数

在贝叶斯框架中，\gamma_b 和 \gamma_1 被称为超参数，其比值 \gamma_1/\gamma_b 对应正则化方法中的参数 \mu，非监督方法的核心思想是让算法自动从数据中估计这两个超参数

超参数的Jeffreys先验

这里我们对超参数的值选择为无信息先验分布，即没有特定取值偏好，即Jeffreys先验，其公式为：

换句话说，该先验在对数尺度上是均匀的，即 \log\gamma 服从均匀分布

联合后验分布

根据贝叶斯定理，后验分布正比于似然函数乘以先验分布：

将各项代入，得到图像和超参数的联合后验分布：

推导过程见附录2

这里 k 是归一化常数，吸收了所有与变量无关的项，包括分母 f(\boldsymbol{y})。这个联合后验分布描述了在给定观测数据 \boldsymbol{y} 后，图像 \boldsymbol{x} 和超参数 \gamma_b、\gamma_1 的联合概率分布。

附录2

联合后验分布的详细推导

根据贝叶斯定理，后验分布等于似然函数乘以先验分布再除以证据：

似然函数是给定图像 \boldsymbol{x} 和噪声精度 \gamma_b 时观测数据的条件分布：

图像先验是给定先验精度 \gamma_1 时图像的分布：

超参数的Jeffreys无信息先验：

将这四项相乘：

整理幂次项：

合并指数项（指数相加）：

将所有与变量无关的常数项（包括 (2\pi)^{-N} 和分母 f(\boldsymbol{y})）吸收到归一化常数 k 中，得到联合后验分布：

后验分布的复杂性

联合后验分布的形式虽然可以写出，但无法直接计算。主要困难在于分母 f(\boldsymbol{y}) 是一个高维积分:

涉及对所有可能的图像和超参数值积分，计算上不可行。此外，即使不考虑归一化常数，变量之间的耦合也使得直接从联合分布采样或求期望变得困难。

解决这个问题有两种主要方法：MCMC采样和变分贝叶斯近似。前者通过构造马尔可夫链来采样后验分布，后者通过可分离的近似分布来逼近真实后验。

MCMC方法：Gibbs采样

算法原理

Gibbs采样器是一种MCMC（马尔可夫链蒙特卡洛）方法，用于从复杂的联合后验分布中采样。其核心思想是：虽然无法直接从联合分布采样，但可以依次从每个变量的条件分布中采样。算法流程如下：

1. 初始化 \boldsymbol{x}^0，\gamma_b^0 和 \gamma_1^0
2. 从条件分布 f(\boldsymbol{x}|\gamma_b^k, \gamma_1^k, \boldsymbol{y}) 中抽取 \boldsymbol{x}^{k+1}
3. 从条件分布 f(\gamma_b|\boldsymbol{x}^{k+1}, \boldsymbol{y}) 中抽取 \gamma_b^{k+1}
4. 从条件分布 f(\gamma_1|\boldsymbol{x}^{k+1}) 中抽取 \gamma_1^{k+1}
5. 返回步骤2

每次迭代中，固定其他变量，更新一个变量。

条件分布的形式

从联合后验分布可以推导出各变量的条件分布。

图像 \boldsymbol{x} 的条件分布

固定 \gamma_b 和 \gamma_1 后，只看与 \boldsymbol{x} 相关的项：

指数项是 \boldsymbol{x} 的二次型，因此这是一个高斯分布。

噪声精度 \gamma_b 的条件分布

固定 \boldsymbol{x} ：

这是Gamma分布的形式，Gamma分布 \text{Gamma}(\alpha, \beta) 的密度函数正比于 \gamma^{\alpha-1}e^{-\beta\gamma}，对比可知

先验精度 \gamma_1 的条件分布，固定 \boldsymbol{x} 后：

同样是Gamma分布，\alpha_1 = N/2，\beta_1 = \|\boldsymbol{D}\boldsymbol{x}\|^2/2。

图像采样的计算困难

虽然 \boldsymbol{x} 的条件分布是高斯分布，但采样存在严重的计算困难。从高斯分布 \mathcal{N}(\boldsymbol{m}, \boldsymbol{\Sigma}) 采样的标准方法是 \boldsymbol{x} = \boldsymbol{m} + \boldsymbol{L}\boldsymbol{z}，其中 \boldsymbol{z} 是标准正态随机向量，\boldsymbol{L} 是协方差矩阵的Cholesky分解 \boldsymbol{\Sigma} = \boldsymbol{L}\boldsymbol{L}^t。

在我们的情况下，协方差矩阵为：

对于 256 \times 256 的图像，这是一个 65536 \times 65536 的矩阵，对其求逆和Cholesky分解在计算上完全不可行。

傅里叶域采样

解决方案是利用循环近似，在傅里叶域中进行采样。傅里叶变换后，卷积变成逐元素乘法，矩阵变成对角矩阵，高维耦合问题变成独立的一维问题。

傅里叶域中的模型为：

\mathring{\boldsymbol{x}} 的条件分布在傅里叶域中写为：

通过配方法整理指数项，可以得到：

Parseval关系

超参数的条件分布需要计算残差的范数 \|\boldsymbol{y} - \boldsymbol{H}\boldsymbol{x}\|^2 和 \|\boldsymbol{D}\boldsymbol{x}\|^2。由于Parseval关系，傅里叶变换保持范数：

因此可以直接在频域中计算这些范数，而无需变换回空域。这进一步提高了计算效率。

预热期（Burn-in）

Gibbs采样的初始值可能远离后验分布的高概率区域，前面若干次迭代的样本不具有代表性，需要丢弃。这段时间称为burn-in。当轨迹稳定在某个值附近波动时，说明马尔可夫链已经收敛到平稳分布。

实验结果

运行400次迭代，burn-in期设为100次。

从超参数链的轨迹可以观察到，\gamma_b 在约50次迭代后稳定在0.037附近，\gamma_1 在约50次迭代后稳定在0.02附近。burn-in = 100的设置是合适的，确保用于计算后验均值的样本来自平稳分布。

MCMC方法的MSE为339.09，与手动调参的最佳结果（335.99，\mu = 0.01）非常接近，证明了MCMC方法结果的正确性

变分贝叶斯方法

从采样到近似的思路转变

MCMC方法通过采样来探索后验分布，理论上可以达到任意精度，但需要大量迭代且每次运行结果不同。变分贝叶斯方法采用不同的思路：用一个简单的近似分布来逼近真实后验分布，将推断问题转化为优化问题。

可分离近似分布族

变分贝叶斯的核心假设是用可分离（factorizable）的近似分布 q 来代替真实后验分布 f：

可分离意味着假设变量之间相互独立，联合分布可以写成各自边缘分布的乘积。真实后验分布中变量之间有复杂的依赖关系，无法写成乘积形式，而近似分布通过独立性假设大大简化了问题

在所有可分离分布构成的集合中，我们寻找最接近真实后验的那一个，接近程度通常用KL散度衡量：

共轭先验与分布形式

由于先验分布与似然函数共轭，最优近似分布的形式可以确定。在我们的情况下，\boldsymbol{x} 的近似分布是高斯分布，\gamma_b 和 \gamma_1 的近似分布是Gamma分布：

变分更新的一般形式

变分贝叶斯的标准更新规则是：更新某个变量的分布时，对联合分布的对数关于其他变量取期望。对于 \mathring{\boldsymbol{x}} 的更新：

这里 \langle \cdot \rangle_{q_b^k q_1^k} 表示关于当前的 q_b^k(\gamma_b) 和 q_1^k(\gamma_1) 求期望，其核心思想是，更新 \boldsymbol{x} 时，我们不知道 \gamma_b 和 \gamma_1 的确切值，所以用它们当前的近似分布来平均化。

变分贝叶斯更新公式

变分贝叶斯的核心是迭代更新三个分布的参数。每次迭代中，固定其他变量的分布，更新一个变量的分布，具体是对联合分布的对数关于其他变量取期望。

对于图像 \mathring{\boldsymbol{x}} 的更新，取期望后随机变量 \gamma_b 和 \gamma_1 被替换为它们的期望值 \widetilde{\gamma_b} = \alpha_b/\beta_b 和 \widetilde{\gamma_1} = \alpha_1/\beta_1。结果是一个高斯分布，参数为：

对于噪声精度 \gamma_b 的更新，需要计算残差范数关于 \mathring{\boldsymbol{x}} 的期望。利用高斯分布的性质 \langle|\mathring{\boldsymbol{x}}|^2\rangle = |\boldsymbol{m}_x|^2 + \boldsymbol{\Sigma}_x^2，结果是一个Gamma分布，参数为：

\beta_b 中的第二项 |\boldsymbol{\Lambda}_h|^2(\boldsymbol{\Sigma}_x^2)^{k+1} 是图像估计不确定性带来的额外贡献，这是变分贝叶斯与点估计方法的关键区别。

对于先验精度 \gamma_1 的更新，类似地得到Gamma分布参数：

算法初始化 \widetilde{\gamma_b} = 1，\widetilde{\gamma_1} = 1，然后循环执行上述三步更新直到收敛，最终输出 \hat{\boldsymbol{x}} = \text{IFFT2}(\boldsymbol{m}_x)。

算法流程总结

变分贝叶斯算法的完整流程为：

初始化 \widetilde{\gamma_b} = 1，\widetilde{\gamma_1} = 1，预计算

循环直到收敛，每次迭代执行三步更新：

第一步更新图像分布参数：

第二步更新噪声精度分布参数：

第三步更新先验精度分布参数：

输出重建图像 \hat{\boldsymbol{x}} = \text{IFFT2}(\boldsymbol{m}_x)

与MCMC的区别

MCMC通过随机采样来探索后验分布，输出是样本序列，估计量通过样本均值计算，每次运行结果不同，理论上可以精确采样后验分布。变分贝叶斯通过优化来找到最佳近似分布，输出是确定性的分布参数，估计量直接是近似分布的均值，每次运行结果相同，但存在近似误差。

实验结果

运行100次迭代。

从收敛曲线可以观察到，\gamma_b 在约60次迭代后稳定在0.037附近，\gamma_1 在约60次迭代后稳定在0.02附近。与MCMC方法不同，变分贝叶斯的收敛曲线是平滑的确定性曲线，没有随机波动。

变分贝叶斯方法的MSE为338.96，超参数收敛值与MCMC方法高度一致。

三种方法对比

三种方法的重建质量非常接近，MSE差异不大。变分贝叶斯方法相比MCMC具有显著的计算效率优势：迭代次数更少（100次对比400次），不需要burn-in期，且每次运行结果确定可重复。MCMC方法的理论优势在于可以精确采样后验分布并提供不确定性量化，而变分贝叶斯只是近似。

结论

本实验系统研究了图像反卷积问题的三种方法。正则化最小二乘方法通过循环近似和FFT实现了高效计算，但需要手动选择正则化参数。贝叶斯框架将正则化参数的选择转化为超参数估计问题，MCMC方法通过Gibbs采样从后验分布中采样，变分贝叶斯方法通过可分离近似分布逼近后验分布。实验结果表明，三种方法的重建质量相近，非监督方法在无需手动调参的情况下达到了与最优监督方法相当的效果，变分贝叶斯在计算效率上优于MCMC。

​