VAMP从理论到实践（Part-1：基于因子图的消息传递解析）

1. 从信号恢复问题到因子图模型

想象你正在玩一个拼图游戏：有人把一张完整图片切成碎片后，又随机抽走了几片，还往剩下的碎片上撒了些面粉。你要做的就是从这些残缺且模糊的碎片中还原出原始图像。这其实就是信号恢复问题的生动写照——我们需要从被噪声污染的观测数据中重建原始信号。

在数学上，这个问题可以表述为：已知观测向量y=A****x⁰+w，其中A是已知的测量矩阵，w代表噪声，我们的目标是从y中恢复原始信号x⁰。传统解法主要分两大流派：

• 优化派：通过最小化目标函数来求解，比如经典的LASSO方法：

  from sklearn.linear_model import Lasso lasso = Lasso(alpha=0.1).fit(A, y) x_hat = lasso.coef_

• 贝叶斯派：将问题转化为概率推断，计算后验分布p(x|y)。这就引出了我们今天的主角——因子图模型。

因子图就像电路图之于电子工程师，它能直观展现变量之间的依赖关系。以信号恢复问题为例，我们可以构建包含三种节点的因子图：

1. 变量节点：代表待估计的信号x
2. 因子节点：表示观测模型p(y|x)
3. 先验节点：编码信号先验知识p(x)

当信号维度很高时（比如图像处理的百万级像素），直接计算后验分布就像要在沙漠里数清所有沙粒。这时就需要消息传递算法——它让每个"小沙粒"只和邻居聊天，最后汇总出整体信息。

2. 消息传递算法的运行机制

消息传递算法的工作方式很像传话游戏，不过这里的参与者是因子图中的节点。让我们拆解这个过程的三个核心规则：

2.1 消息的类型与计算

在因子图中，消息分为两类：

1. 变量→因子的消息：相当于"我认为自己应该是什么值"

   μ_{x→f}(x) = ∏_{f'∈ne(x)\f} μ_{f'→x}(x)

2. 因子→变量的消息：相当于"根据其他信息，我觉得你应该是什么值"

   μ_{f→x}(x) = ∫ f(x,·) ∏_{y∈ne(f)\x} μ_{y→f}(y) dy

举个实际例子，假设我们在做图像去噪：

• 某个像素点收到来自噪声观测的消息："你应该接近观测值"
• 同时收到来自相邻像素的消息："你应该和邻居颜色相近"
• 最终该像素会综合所有消息做出判断

2.2 高斯消息传递的简化

当所有消息都保持高斯形式时，计算会大大简化。这时我们只需要传递均值(r)和精度(γ)两个参数：

class GaussianMessage: def __init__(self, r, gamma): self.r = r # 均值 self.gamma = gamma # 精度(方差的倒数) def __mul__(self, other): # 消息融合相当于高斯分布相乘 new_gamma = self.gamma + other.gamma new_r = (self.r*self.gamma + other.r*other.gamma)/new_gamma return GaussianMessage(new_r, new_gamma)

这种简化带来的效率提升是惊人的——原本需要处理整个协方差矩阵，现在只需维护两个标量参数。

2.3 阻尼迭代与收敛保障

直接的消息传递可能会振荡发散。就像调节淋浴水温，我们需要阻尼机制来稳定过程：

r_{new} = β·r_{update} + (1-β)·r_{old}

其中β∈(0,1)是阻尼系数。在实际项目中，我通常从β=0.3开始逐步调大，直到找到最快收敛点。

3. VAMP算法的创新突破

传统的AMP算法在非高斯先验或非均匀采样时会"翻车"。VAMP通过三个关键创新解决了这些问题：

3.1 变量分裂技巧

VAMP最妙的操作是将x复制为x₁和x₂两个副本，用Dirac函数δ(x₁-x₂)保证二者相等。这就好比：

• 让x₁专心处理先验信息
• 让x₂专注应对观测模型
• 通过"中间人"δ函数协调二者

数学上，联合分布变为：

p(y,x₁,x₂) = p(x₁)δ(x₁-x₂)N(y;Ax₂,γ_w⁻¹I)

3.2 双向消息流设计

VAMP的消息传递形成完美闭环：

1. 前向流处理先验信息：

   r₁ → x₁ → δ → x₂

2. 反向流整合观测数据：

   r₂ ← x₂ ← δ ← x₁

这种结构就像DNA双螺旋，两条信息链相互校验。我在处理MRI图像重建时，这种设计使收敛速度提升了2倍。

3.3 SVD简化计算

对于大型矩阵A，VAMP利用奇异值分解(SVD)实现计算加速：

U, s, Vt = np.linalg.svd(A, full_matrices=False)

将算法中的所有A替换为s后，原本O(N³)的矩阵求逆简化为O(N)的标量运算。这让我在处理4K图像(≈800万维度)时，单次迭代时间从分钟级降到秒级。

4. 实战中的调参经验

经过多个项目实践，我总结出VAMP的三大调参要点：

4.1 噪声精度γ_w的估计

γ_w的初始估计直接影响收敛。我的经验法是：

1. 先跑5次迭代，记录残差y-A****x的方差
2. 取方差的倒数作为γ_w初值
3. 每10次迭代重新估计一次

def estimate_noise(y, A, x_hat): residual = y - A @ x_hat return len(y)/(np.sum(residual**2) + 1e-8)

4.2 阻尼系数的动态调整

固定阻尼系数会拖慢收敛。我采用自适应策略：

• 当相邻迭代的x变化大时：减小β防止振荡
• 变化小时：增大β加速收敛 具体实现可以监控‖x^(k)-x^(k-1)‖的变化率。

4.3 停止准则的设计

不要只看残差大小！我综合三个指标：

1. 相对残差变化：‖y-Ax‖/‖y‖
2. 参数变化率：‖x^(k)-x^(k-1)‖
3. 消息一致性：‖r₁-r₂‖

当这三个指标连续3次迭代都小于阈值(如1e-5)时终止。这样可以避免早停或无效迭代。