近似消息传递(Approximate Message Passing)算法简介

近似消息传递(Approximate Message Passing)算法简介

1 前言

近似���息传递(Approximate Message Passing, AMP)算法是基于消息传递算法，也叫和-积算法(Sum-Product Algorithm, SPA)，还被称为置信传播(Belief Propagation, BP)算法，经过一系列假设与简化得来，这其中包括了中心极限定理(Central Limit Theory, CLT)和泰勒级数(Taylor Series)展开等[@zou_concise_2022]。

2 基础知识

这一节简单介绍基础知识，包括SPA、CLT、高斯分布乘以高斯分布、泰勒级数以及后验概率密度函数的偏导数求法的结果。

2.1 和-积算法(Sum-Product Algorithm, SPA)

消息传递算法是基于因子图(factor graph)的节点(node)之间的消息沿边(edge)相互传递而得名。

我们假设一个线性估计模型： y = H x + w ( 1 ) \mathbf{y = Hx + w} \quad(1) y=Hx+w(1)

其中， y ∈ C M \mathbf{y} \in \mathbb{C}^M y∈CM是已知的观测数据(observed data)， H ∈ C M × N \mathbf{H} \in \mathbb{C}^{M\times N} H∈CM×N是已知的观测矩阵， x ∈ C N \mathbf{x}\in \mathbb{C}^N x∈CN是待估计的信号，并且知道先验分布为 p ( x ) p(\mathbf{x}) p(x)， w ∈ C M \mathbf{w} \in \mathbb{C}^M w∈CM是复高斯白噪声，即 w ∼ C N ( w ; 0 , σ 2 I ) \mathbf{w} \sim \mathcal{CN}(\mathbf{w};\mathbf{0},\sigma^2 \mathbf{I}) w∼CN(w;0,σ2I)。

根据公式(1)存在一个全局概率密度分布(global probability density distribution,PDF) p ( y , x ∣ H ) p(\mathbf{y,x|H}) p(y,x∣H)，为了简单起见，我们省略 H \mathbf{H} H，写为 p ( y , x ) p(\mathbf{y,x}) p(y,x)。我们假设 y \mathbf{y} y和 x \mathbf{x} x的各自的元素都是分别各自独立同分布的。那么我们可以将全局PDF进行分解：

p ( y , x ) = p ( y ∣ x ) p ( x ) ( 2 ) = ∏ a M p ( y a ∣ x ) ∏ i N p ( x i ) ( 3 ) = ∏ a f a ( x ) ( 4 ) \begin{aligned} p(\mathbf{y,x}) & = p(\mathbf{y|x})p(\mathbf{x}) \quad(2) \\ & = \prod_a^M p(y_a |\mathbf{x}) \prod_i^N p(x_i) \quad(3)\\ & = \prod_a f_a(\mathbf{x}) \quad(4) \end{aligned} p(y,x)​=p(y∣x)p(x)(2)=a∏M​p(ya​∣x)i∏N​p(xi​)(3)=a∏​fa​(x)(4)​ 公式(2)是条件概率公式，公式(3)是根据独立假设而得来，公式(4)是因子形式分解，中的 x \mathbf{x} x是与该因子相连的全部变量，不是待估计的信号 x ∈ C N \mathbf{x}\in \mathbb{C}^N x∈CN。

从变量节点 x i x_i xi​传递到因子节点 f a ( x ) f_a(\mathbf{x}) fa​(x)的消息为： μ i → a ( x i ) ∝ ∏ b ≠ a μ i ← b ( x i ) ( 5 ) \mu_{i \rightarrow a}(x_i) \propto \prod_{b \neq a} \mu_{i \leftarrow b}(x_i) \quad(5) μi→a​(xi​)∝b=a∏​μi←b​(xi​)(5)

从因子节点 f a ( x ) f_a(\mathbf{x}) fa​(x)传递到变量节点 x i x_i xi​的消息为： μ i ← a ( x i ) ∝ ∫ x \ i f a ( x ) ∏ j ≠ i μ j → a ( x j ) d x \ i ( 6 ) \mu_{i \leftarrow a}(x_i) \propto \int_{\mathbf{x}_{\backslash i}} f_a(\mathbf{x}) \prod_{j \neq i}\mu_{j \rightarrow a}(x_j) \mathbf{d} \mathbf{x}_{\backslash i} \quad(6) μi←a​(xi​)∝∫x\i​​fa​(x)j=i∏​μj→a​(xj​)dx\i​(6)

j ≠ i j \neq i j=i和 b ≠ a b \neq a b=a的原因是避免传递自身的消息，即传递外信息，从而不会出现自己不断置信自己的情况。消息的本质是PDF，由于这里未进行归一化，所以使用正比于 ∝ \propto ∝。

对于我们处理的模型， M × N M \times N M×N的因子图为：

这张图省略了观测变量 y \mathbf{y} y，可以画出来，再用斜线涂黑，代表观测。整张图只有两种不同的节点，并且相同类型的节点不会连接，分别是变量节点(variable node)和因子节点(factor node)，又叫校验节点。消息在两类不同的节点之间相互传递。

那么对于上图来说，公式(5)和(6)可以分别写为： μ i → a ( x i ) ∝ p ( x i ) ∏ b ≠ a μ i ← b ( x i ) ( 7 ) μ i ← a ( x i ) ∝ ∫ x \ i p ( y a ∣ x ) ∏ j ≠ i N μ j → a ( x j ) d x \ i ( 8 ) \begin{aligned} \mu_{i \rightarrow a}(x_i) & \propto p(x_i) \prod_{b \neq a} \mu_{i \leftarrow b}(x_i) \quad(7) \\ \mu_{i \leftarrow a}(x_i) & \propto \int_{\mathbf{x}_{\backslash i}} p(y_a|\mathbf{x}) \prod_{j \neq i}^N \mu_{j \rightarrow a}(x_j) \mathbf{d} \mathbf{x}_{\backslash i} \quad(8) \end{aligned} μi→a​(xi​)μi←a​(xi​)​∝p(xi​)b=a∏​μi←b​(xi​)(7)∝∫x\i​​p(ya​∣x)j=i∏N​μj→a​(xj​)dx\i​(8)​

注意公式(5)和(6)是最一般的形式，即公式(4)，而公式(7)和(8)是针对公式(3)得来。可以看出来公式(7)和(8)只是变量节点 x i x_i xi​和因子节点 p ( y a ∣ x ) p(y_a|\mathbf{x}) p(ya​∣x)之间在进行消息传递，而没有 x i x_i xi​与 p ( x i ) p(x_i) p(xi​)之间的传递。这是因为 p ( x i ) p(x_i) p(xi​)是边缘节点，只与 x i x_i xi​相连，只能传递给 x i x_i xi​消息 p ( x i ) p(x_i) p(xi​)，所以不需要更新。或者说，因为传递的是外信息， p ( x i ) p(x_i) p(xi​)就算被更新了，也是传递最原始的先验知识 p ( x i ) p(x_i) p(xi​)。

由于这里的因子图不是可以得到精确解的树状结构，而是环状结构，消息会在两种节点间进行迭代，得到近似解。所以，我们把公式(7)和(8)重写为： μ i → a t + 1 ( x i ) ∝ p ( x i ) ∏ b ≠ a μ i ← b ( x i ) ( 9 ) μ i ← a t ( x i ) ∝ ∫ x \ i p ( y a ∣ x ) ∏ j ≠ i N μ j → a ( x j ) d x \ i ( 10 ) \begin{aligned} \mu_{i \rightarrow a}^{t+1}(x_i) & \propto p(x_i) \prod_{b \neq a} \mu_{i \leftarrow b}(x_i) \quad(9) \\ \mu_{i \leftarrow a}^{t}(x_i) & \propto \int_{\mathbf{x}_{\backslash i}} p(y_a|\mathbf{x}) \prod_{j \neq i}^N \mu_{j \rightarrow a}(x_j) \mathbf{d} \mathbf{x}_{\backslash i} \quad(10) \end{aligned} μi→at+1​(xi​)μi←at​(xi​)​∝p(xi​)b=a∏​μi←b​(xi​)(9)∝∫x\i​​p(ya​∣x)j=i∏N​μj→a​(xj​)dx\i​(10)​ 其中 t t t是第 t t t次迭代。

对于我们的需要的后验知识，可以利用下式求出后验分布： p ( x i t ) ∝ p ( x i ) ∏ a = 1 M μ i ← a t ( x i ) ( 11 ) p(x_i^t) \propto p(x_i) \prod_{a=1}^{M} \mu_{i \leftarrow a }^t (x_i) \quad(11) p(xit​)∝p(xi​)a=1∏M​μi←at​(xi​)(11)

即把所有与 x i x_i xi​相连的因子节点传递给它的消息进行合并，这被称为因子图的边缘分布，但并不是一般的概率分布的边缘分布(marginal

distribution)，它实际上是后验分布 p ( x i ∣ y ) p(x_i|\mathbf{y}) p(xi​∣y)。

至此，我们完成了简要讨论SPA。

2.2 中心极限定理(Central Limit Theory, CLT)

我们这里考虑最简单的独立同分布(independent identical distribution, i.i.d.)的形式，也叫Lindeberg-Feller中心极限定理。

定理 1 (独立同分布中心极限定理).

假设 { X n } n = 1 N \{X_n\}^N_{n=1} {Xn​}n=1N​的 N N N个样本是i.i.d.的，有 E ( X n ) = m \mathbb{E}(X_n) = m E(Xn​)=m且 V a r ( X n ) = σ 2 > 0 {\rm Var}(X_n) = \sigma^2 > 0 Var(Xn​)=σ2>0。当 N N N足够大时， X ˉ n = 1 N ∑ n = 1 N X n \bar{X}_n = \frac{1}{N}\sum_{n=1}^{N}X_n Xˉn​=N1​∑n=1N​Xn​近似服从正态分布 N ( m , σ 2 N ) \mathcal{N}(m, \frac{\sigma^2}{N}) N(m,Nσ2​)。

我们这里将正态分布重写为标准正态分布，即 X ˉ n − m σ / N ∼ N ( 0 , 1 ) \frac{\bar{X}_n - m}{\sigma / \sqrt{N} } \sim \mathcal{N}(0, 1) σ/N ​Xˉn​−m​∼N(0,1)。同样的，在复数域下，可以写为 X ˉ n − m σ / N ∼ C N ( 0 , 1 ) \frac{\bar{X}_n - m}{\sigma / \sqrt{N} } \sim \mathcal{CN}(0, 1) σ/N ​Xˉn​−m​∼CN(0,1)

2.3 高斯分布乘高斯分布

接下来介绍高斯分布乘高斯分布的引理：

定理 2 (Gaussian reproduction lemma).

对于两个高斯分布， N ( x ; a , A ) \mathcal{N}(x;a,A) N(x;a,A)和 N ( x ; b , B ) \mathcal{N}(x;b,B) N(x;b,B)，它们的乘积为：

N ( x ; a , A ) N ( x ; b , B ) = N ( x ; c , C ) N ( 0 ; a − b , A + B ) ∝ N ( x ; c , C ) ( 12 ) \mathcal{N}(x;a,A)\mathcal{N}(x;b,B) = \mathcal{N}(x;c,C)\mathcal{N}(0;a-b,A+B) \propto \mathcal{N}(x;c,C) \quad(12) N(x;a,A)N(x;b,B)=N(x;c,C)N(0;a−b,A+B)∝N(x;c,C)(12)其中 1 C = 1 A + 1 B \frac{1}{C} = \frac{1}{A} + \frac{1}{B} C1​=A1​+B1​, c C = a A + b B \frac{c}{C} = \frac{a}{A} + \frac{b}{B} Cc​=Aa​+Bb​。

这里的证明很简单，简单的想办法把两个高斯分布的乘积拼凑高斯分布的形式即可，可以很显式地得出均值和方差，可以将等式两边除以 N ( 0 ; a − b , A + B ) \mathcal{N}(0;a-b,A+B) N(0;a−b,A+B)以归一化，得到真正的PDF。同样的，给出复数域的情况：

C N ( x ; a , A ) C N ( x ; b , B ) = C N ( x ; c , C ) C N ( 0 ; a − b , A + B ) ∝ C N ( x ; c , C ) ( 13 ) \mathcal{CN}(x;a,A)\mathcal{CN}(x;b,B) = \mathcal{CN}(x;c,C)\mathcal{CN}(0;a-b,A+B) \propto \mathcal{CN}(x;c,C)\quad(13) CN(x;a,A)CN(x;b,B)=CN(x;c,C)CN(0;a−b,A+B)∝CN(x;c,C)(13)至此，我们完成了基础知识的简单介绍。

2.4 一阶泰勒级数展开

对于一个两变量的函数 f ( x , y ) f(x,y) f(x,y)，并假设其Lipschitz连续，我们一阶泰勒级数展开有：

f ( x + Δ x , y + Δ y ) ≈ f ( x , y ) + Δ x f x ′ ( x , y ) + Δ y f y ′ ( x , y ) ( 14 ) f(x+\Delta x ,y+\Delta y ) \approx f(x,y) + \Delta x f'_x(x,y)+ \Delta y f'_y(x,y)\quad(14) f(x+Δx,y+Δy)≈f(x,y)+Δxfx′​(x,y)+Δyfy′​(x,y)(14)其中， f x ′ f'_x fx′​和 f y ′ f'_y fy′​分别代表了对 x x x和 y y y分别求偏导。

2.5 后验概率密度函数的偏导数求法

对于一个任意有界非负的函数 f ( x ) f(x) f(x)（其实是加权），并定义一个分布 P ( x ) = f ( x ) C N ( x ; m , v ) ∫ f ( x ) C N ( x ; m , v ) d x \mathcal{P}(x) = \frac{f(x)\mathcal{CN}(x;m,v)}{\int f(x)\mathcal{CN}(x;m,v) {\rm d}x} P(x)=∫f(x)CN(x;m,v)dxf(x)CN(x;m,v)​。其均值和方差表示为 E ( x ) = ∫ x P ( x ) d x \mathbb{E}(x) = \int x \mathcal{P}(x) {\rm d}x E(x)=∫xP(x)dx和 V a r ( x ) = ∫ [ x − E ( x ) ] 2 P ( x ) d x {\rm Var}(x) = \int [x - \mathbb{E}(x)]^2 \mathcal{P}(x) {\rm d}x Var(x)=∫[x−E(x)]2P(x)dx。我们有：

∂ ∫ x P ( x ) d x ∂ m = ∫ x x − m v f ( x ) C N ( x ; m , v ) d x ⋅ ∫ f ( x ) C N ( x ; m , v ) d x [ ∫ f ( x ) C N ( x ; m , v ) d x ] 2 − ∫ x f ( x ) C N ( x ; m , v ) d x ⋅ ∫ x − m v f ( x ) C N ( x ; m , v ) d x [ ∫ f ( x ) C N ( x ; m , v ) d x ] 2 = V a r ( x ) v ( 15 ) \begin{split} \frac{\partial \int x \mathcal{P}(x) {\rm d}x }{\partial m} & = \frac{\int x \frac{x-m}{v} f(x) \mathcal{CN}(x;m,v){\rm d}x \cdot \int f(x) \mathcal{CN}(x;m,v){\rm d}x }{[\int f(x)\mathcal{CN}(x;m,v){\rm d}x ]^2} \\ & \quad - \frac{\int x f(x) \mathcal{CN}(x;m,v){\rm d}x \cdot \int \frac{x-m}{v} f(x) \mathcal{CN}(x;m,v){\rm d}x }{[\int f(x)\mathcal{CN}(x;m,v){\rm d}x ]^2} \\ & = \frac{{\rm Var}(x)}{v} \end{split}\quad(15) ∂m∂∫xP(x)dx​​=[∫f(x)CN(x;m,v)dx]2∫xvx−m​f(x)CN(x;m,v)dx⋅∫f(x)CN(x;m,v)dx​−[∫f(x)CN(x;m,v)dx]2∫xf(x)CN(x;m,v)dx⋅∫vx−m​f(x)CN(x;m,v)dx​=vVar(x)​​(15)

3 AMP推导

我们回到SPA的模型以及公式(9)和(10)，开始对它们进行简化。

3.1 从因子节点传递到变量节点的消息

对于公式(10)中的 p ( y a ∣ x ) p(y_a|\mathbf{x}) p(ya​∣x)而言，我们这里不再省略 H \mathbf{H} H，并用 H a ∼ \mathbf{H}_{a\sim} Ha∼​表示其第 a a a行：

p ( y a ∣ H , x ) = p w a ( y a − H a ∼ x ) = C N ( y a ; H a ∼ x , σ 2 ) = C N ( w a = y a − H a ∼ x ; 0 , σ 2 ) ( 16 ) = 1 π σ 2 exp ⁡ ( − ∣ y a − H a ∼ x ∣ 2 σ 2 ) \begin{split} p(y_a|\mathbf{H,x}) & = p_{w_a}(y_a - \mathbf{H}_{a\sim}\mathbf{x} ) \\ & = \mathcal{CN}(y_a;\mathbf{H}_{a\sim}\mathbf{x} , \sigma^2 ) \\ & = \mathcal{CN}(w_a = y_a - \mathbf{H}_{a\sim}\mathbf{x}; 0, \sigma^2 )\quad(16)\\ & = \frac{1}{\pi \sigma^2}\exp( \frac{-|y_a - \mathbf{H}_{a\sim}\mathbf{x}|^2}{\sigma^2} ) \end{split} p(ya​∣H,x)​=pwa​​(ya​−Ha∼​x)=CN(ya​;Ha∼​x,σ2)=CN(wa​=ya​−Ha∼​x;0,σ2)(16)=πσ21​exp(σ2−∣ya​−Ha∼​x∣2​)​

其中 p w a ( ⋅ ) p_{w_a}(\cdot) pwa​​(⋅)代表是随机变量 w a w_a wa​的PDF， H a ∼ x \mathbf{H}_{a\sim}\mathbf{x} Ha∼​x也可以被写作 ∑ k = 1 N h a k x k \sum_{k=1}^{N}h_{ak}x_k ∑k=1N​hak​xk​。

为了简化公式(10)，我们将其重写为： μ i ← a t ( x i ) ∝ ∫ x \ i p ( y a ∣ H a ∼ x ) ∏ j ≠ i N μ j → a t ( x j ) d x \ i ∝ ∫ x \ i ∫ z a p ( y a ∣ z a ) δ ( z a − ∑ k = 1 N h a k x k ) d z a ∏ j ≠ i N μ j → a t ( x j ) d x \ i ∝ ∫ z a p ( y a ∣ z a ) E { δ ( z a − ∑ j ≠ i N h a j x j − h a i x i ) } d z a ( 17 ) \begin{split} \mu_{i \leftarrow a}^{t}(x_i) & \propto \int_{\mathbf{x}_{\backslash i}} p(y_a|\mathbf{H}_{a\sim}\mathbf{x}) \prod_{j \neq i}^N \mu_{j \rightarrow a}^t(x_j) \mathbf{d} \mathbf{x}_{\backslash i} \\ & \propto \int_{\mathbf{x}_{\backslash i}} \int_{z_a} p(y_a|z_a) \delta(z_a - \sum_{k=1}^{N}h_{ak}x_k) {\rm d}z_a \prod_{j \neq i}^N \mu_{j \rightarrow a}^t(x_j) \mathbf{d} \mathbf{x}_{\backslash i}\\ & \propto \int_{z_a} p(y_a|z_a) \mathbb{E}\left\{\delta \left( z_a - \sum_{j \neq i}^{N}h_{aj}x_j - h_{ai}x_i \right) \right\} {\rm d}z_a \quad(17) \end{split} μi←at​(xi​)​∝∫x\i​​p(ya​∣Ha∼​x)j=i∏N​μj→at​(xj​)dx\i​∝∫x\i​​∫za​​p(ya​∣za​)δ(za​−k=1∑N​hak​xk​)dza​j=i∏N​μj→at​(xj​)dx\i​∝∫za​​p(ya​∣za​)E⎩ ⎨ ⎧​δ ​za​−j=i∑N​haj​xj​−hai​xi​ ​⎭ ⎬ ⎫​dza​(17)​

其中，期望是基于 ∏ j ≠ i N μ j → a t ( x j ) \prod_{j \neq i}^N \mu_{j \rightarrow a}^t(x_j) ∏j=iN​μj→at​(xj​)求的。我们定义了新的随机变量 z a = H a ∼ x z_a = \mathbf{H}_{a\sim}\mathbf{x} za​=Ha∼​x，以及与它相关的随机变量 ζ i ← a t \zeta_{i \leftarrow a}^{t} ζi←at​。再定义根据 μ j → a t ( x j ) \mu_{j \rightarrow a}^{t}(x_j) μj→at​(xj​)与 x j x_j xj​相关的 ξ j → a t \xi_{j \rightarrow a}^{t} ξj→at​。并且给出 ξ j → a t \xi_{j \rightarrow a}^{t} ξj→at​的均值和方差分别为 x ^ j → a t \hat{x}_{j \rightarrow a}^{t} x^j→at​和 v ^ j → a t \hat{v}_{j \rightarrow a}^{t} v^j→at​。由公式(17)，当 N N N趋近于无穷时，根据CLT，可以得到随机变量 ζ i ← a t \zeta_{i \leftarrow a}^{t} ζi←at​收敛到高斯随机变量，其均值和方差为：

E ( ζ i ← a t ) = Z i ← a t + h a i x i , V a r ( ζ i ← a t ) = V i ← a t ( 18 ) \mathbb{E}(\zeta_{i \leftarrow a}^{t}) = Z^t_{i \leftarrow a}+h_{ai}x_i , \quad {\rm Var}(\zeta_{i \leftarrow a}^{t}) = V_{i \leftarrow a}^t\quad(18) E(ζi←at​)=Zi←at​+hai​xi​,Var(ζi←at​)=Vi←at​(18)

其中 Z i ← a t = ∑ j ≠ i h a j x ^ j → a t , V i ← a t = ∑ j ≠ i ∣ h a j ∣ 2 v ^ j → a t ( 19 ) Z^t_{i \leftarrow a} = \sum_{j \neq i} h_{aj}\hat{x}^t_{j \rightarrow a}, \quad V_{i \leftarrow a}^t = \sum_{j \neq i} |h_{aj}|^2 \hat{v}^t_{j \rightarrow a } \quad(19) Zi←at​=j=i∑​haj​x^j→at​,Vi←at​=j=i∑​∣haj​∣2v^j→at​(19)

实际上公式(19)可以视为在 x ∼ i \mathbf{x}_{\sim i} x∼i​下对 y a y_a ya​的估计。

基于这里的高斯假设，我们可以将公式(17)中的 E { δ ( z a − ∑ j ≠ i N h a j x j − h a i x i ) } \mathbb{E}\{\delta ( z_a - \sum_{j \neq i}^{N}h_{aj}x_j - h_{ai}x_i )\} E{δ(za​−∑j=iN​haj​xj​−hai​xi​)}替换为， C N ( z a ; h a i x i + Z i ← a t , V i ← a t ) \mathcal{CN}(z_a ; h_{ai}x_i + Z^t_{i \leftarrow a}, V_{i \leftarrow a}^t) CN(za​;hai​xi​+Zi←at​,Vi←at​)，再根据高斯分布乘高斯分布的引理， μ i ← a t ( x i ) \mu_{i \leftarrow a}^{t}(x_i) μi←at​(xi​)被近似为：

μ i ← a t ( x i ) ∝ ∫ z a p ( y a ∣ z a ) C N ( z a ; h a i x i + Z i ← a t , V i ← a t ) d z a ∝ ∫ z a C N ( y a ; H a ∼ x , σ 2 ) C N ( z a ; h a i x i + Z i ← a t , V i ← a t ) d z a ∝ ∫ z a C N ( z a ; y a , σ 2 ) C N ( z a ; h a i x i + Z i ← a t , V i ← a t ) d z a ∝ C N ( 0 ; y a − h a i x i − Z i ← a t , σ 2 + V i ← a t ) ∝ C N ( x i ; y a − Z i ← a t h a i , σ 2 + V i ← a t ∣ h a i ∣ 2 ) ( 20 ) \begin{split} \mu_{i \leftarrow a}^{t}(x_i) & \propto \int_{z_a} p(y_a|z_a)\mathcal{CN}(z_a ; h_{ai}x_i + Z^t_{i \leftarrow a}, V_{i \leftarrow a}^t) {\rm d}z_a \\ & \propto \int_{z_a} \mathcal{CN}(y_a;\mathbf{H}_{a\sim}\mathbf{x} , \sigma^2 ) \mathcal{CN}(z_a ; h_{ai}x_i + Z^t_{i \leftarrow a}, V_{i \leftarrow a}^t) {\rm d}z_a \\ & \propto \int_{z_a} \mathcal{CN}(z_a; y_a , \sigma^2 ) \mathcal{CN}(z_a ; h_{ai}x_i + Z^t_{i \leftarrow a}, V_{i \leftarrow a}^t) {\rm d}z_a \\ & \propto \mathcal{CN}(0; y_a - h_{ai}x_i - Z^t_{i \leftarrow a} , \sigma^2 + V_{i \leftarrow a}^t) \\ & \propto \mathcal{CN}(x_i;\frac{y_a - Z_{i \leftarrow a}^t}{h_{ai}},\frac{\sigma^2 + V_{i \leftarrow a}^t}{|h_{ai}|^2})\quad(20) \end{split} μi←at​(xi​)​∝∫za​​p(ya​∣za​)CN(za​;hai​xi​+Zi←at​,Vi←at​)dza​∝∫za​​CN(ya​;Ha∼​x,σ2)CN(za​;hai​xi​+Zi←at​,Vi←at​)dza​∝∫za​​CN(za​;ya​,σ2)CN(za​;hai​xi​+Zi←at​,Vi←at​)dza​∝CN(0;ya​−hai​xi​−Zi←at​,σ2+Vi←at​)∝CN(xi​;hai​ya​−Zi←at​​,∣hai​∣2σ2+Vi←at​​)(20)​

这其中第二个 ∝ \propto ∝到第三个 ∝ \propto ∝之间，利用了高斯分布的性质，更换了PDF的变量，以保证两个复高斯分布都是基于随机变量 z a z_a za​，从而可以使用高斯分布乘高斯分布的引理。

注意到此时消息变成了复数高斯的，且计算时也需要传入消息的均值与方差。高斯消息仅存在两个充分统计量，会让参数传递变得方便。在此，定义此条消息的均值和方差为：

x ^ i ← a t = y a − Z i ← a t h a i ; v ^ i ← a t = σ 2 + V i ← a t ∣ h a i ∣ 2 ( 21 ) \hat{x}_{i \leftarrow a}^t = \frac{y_a - Z_{i \leftarrow a}^t}{h_{ai}}; \quad \hat{v}_{i \leftarrow a}^t = \frac{\sigma^2 + V_{i \leftarrow a}^t}{|h_{ai}|^2}\quad(21) x^i←at​=hai​ya​−Zi←at​​;v^i←at​=∣hai​∣2σ2+Vi←at​​(21)

值得注意的是，这里对与PDF中的方差和均值没有区分";“和”|“，实际上这里更严谨需要使用”|“，但为了表明均值和方差，我们采用”;"。并且，此处如果 H \mathbf{H} H中某个元素为0，会导致除以0的情况发生导致错误。事实上，可以在后面的推导中避免这一问题。

3.2 从变量节点传递到因子节点的消息

我们回到公式(9)，我们先对 ∏ b ≠ a M μ i ← b t ( x i ) \prod_{b \neq a}^{M} \mu_{i \leftarrow b}^t (x_i) ∏b=aM​μi←bt​(xi​)进行分析：

∏ b ≠ a M μ i ← b t ( x i ) = ∏ b ≠ a M C N ( x i ; x ^ i ← b t , v ^ i ← b t ) ∝ C N ( x i ; r ^ i → a t , Σ ^ i → a t ) ( 22 ) \begin{split} \prod_{b \neq a}^{M} \mu_{i \leftarrow b}^t (x_i) & = \prod_{b \neq a}^{M} \mathcal{CN}(x_i;\hat{x}_{i \leftarrow b}^t , \hat{v}_{i \leftarrow b}^t ) \\ & \propto \mathcal{CN}(x_i;\hat{r}_{i \rightarrow a}^t , \hat{\Sigma}_{i \rightarrow a}^t )\quad(22) \end{split} b=a∏M​μi←bt​(xi​)​=b=a∏M​CN(xi​;x^i←bt​,v^i←bt​)∝CN(xi​;r^i→at​,Σ^i→at​)(22)​ 其中 1 Σ ^ i → a t = ∑ b ≠ a M 1 v ^ i ← b t = ∑ b ≠ a M ∣ h b i ∣ 2 σ 2 + V i ← b t ( 23 ) r ^ i → a t Σ ^ i → a t = ∑ b ≠ a M x ^ i ← b t v ^ i ← b t = ∑ b ≠ a M y b − Z i ← b t h b i ∣ h b i ∣ 2 σ 2 + V i ← b t = ∑ b ≠ a M h b i ∗ ( y b − Z i ← b t ) σ 2 + V i ← b t ( 24 ) \begin{aligned} \frac{1}{\hat{\Sigma}_{i \rightarrow a}^t} & = \sum_{b \neq a}^{M} \frac{1}{\hat{v}_{i \leftarrow b}^t} \\ & = \sum_{b \neq a}^{M} \frac{|h_{bi}|^2 }{\sigma^2 + V_{i \leftarrow b}^t} \quad(23)\\ \frac{\hat{r}_{i \rightarrow a}^t}{\hat{\Sigma}_{i \rightarrow a}^t} & = \sum_{b \neq a}^{M} \frac{\hat{x}_{i \leftarrow b}^t}{ \hat{v}_{i \leftarrow b}^t } \\ & = \sum_{b \neq a}^{M} \frac{y_b - Z_{i \leftarrow b}^t}{h_{bi}} \frac{|h_{bi}|^2 }{\sigma^2 + V_{i \leftarrow b}^t}\\ & = \sum_{b \neq a}^{M} \frac{h_{bi}^* (y_b - Z_{i \leftarrow b}^t) }{\sigma^2 + V_{i \leftarrow b}^t} \quad(24) \end{aligned} Σ^i→at​1​Σ^i→at​r^i→at​​​=b=a∑M​v^i←bt​1​=b=a∑M​σ2+Vi←bt​∣hbi​∣2​(23)=b=a∑M​v^i←bt​x^i←bt​​=b=a∑M​hbi​yb​−Zi←bt​​σ2+Vi←bt​∣hbi​∣2​=b=a∑M​σ2+Vi←bt​hbi∗​(yb​−Zi←bt​)​(24)​

这里完全是利用的公式(13)。之后，我们再结合先验信息 p ( x i ) p(x_i) p(xi​)： μ i → a t + 1 ( x i ) ∝ p ( x i ) ∏ b ≠ a M μ i ← b ( x i ) ∝ p ( x i ) C N ( x i ; r ^ i → a t , Σ ^ i → a t ) ( 25 ) \begin{split} \mu_{i \rightarrow a}^{t+1}(x_i) & \propto p(x_i) \prod_{b \neq a}^M \mu_{i \leftarrow b}(x_i) \\ & \propto p(x_i)\mathcal{CN}(x_i ;\hat{r}_{i \rightarrow a}^t , \hat{\Sigma}_{i \rightarrow a}^t )\quad(25) \end{split} μi→at+1​(xi​)​∝p(xi​)b=a∏M​μi←b​(xi​)∝p(xi​)CN(xi​;r^i→at​,Σ^i→at​)(25)​

我们不需要知道 μ i → a t + 1 ( x i ) \mu_{i \rightarrow a}^{t+1}(x_i) μi→at+1​(xi​)的确切分布，只需要知道其均值和方差，并近似为一个新的复高斯分布 μ i → a t + 1 ( x i ) ∝ C N ( x i ; x ^ i → a t + 1 , v ^ i → a t + 1 ) \mu_{i \rightarrow a}^{t+1}(x_i) \propto \mathcal{CN}(x_i;\hat{x}^{t+1}_{i \rightarrow a}, \hat{v}^{t+1}_{i \rightarrow a}) μi→at+1​(xi​)∝CN(xi​;x^i→at+1​,v^i→at+1​)，其中

x ^ i → a t + 1 = E ( x i ) = ∫ x i x i 1 C p ( x i ) C N ( x i ; r ^ i → a t , Σ ^ i → a t ) d x i = d e f F ( x i ; r ^ i → a t , Σ ^ i → a t ) ( 26 ) v ^ i → a t + 1 = V a r ( x i ) = ∫ x i ( x i − x ^ i → a t + 1 ) 2 1 C p ( x i ) C N ( x i ; r ^ i → a t , Σ ^ i → a t ) d x i = d e f G ( x i ; r ^ i → a t , Σ ^ i → a t ) ( 27 ) \begin{aligned} \hat{x}^{t+1}_{i \rightarrow a} & = \mathbb{E}(x_i) \notag \\ & = \int_{x_i} x_i \frac{1}{C}p(x_i) \mathcal{CN}(x_i ;\hat{r}_{i \rightarrow a}^t , \hat{\Sigma}_{i \rightarrow a}^t) {\rm d}x_i \notag \\ & \overset{def}{=} F(x_i;\hat{r}_{i \rightarrow a}^t , \hat{\Sigma}_{i \rightarrow a}^t) \quad(26) \\ \hat{v}^{t+1}_{i \rightarrow a} & = {\rm Var}(x_i) \notag \\ & = \int_{x_i} (x_i - \hat{x}^{t+1}_{i \rightarrow a} )^2 \frac{1}{C}p(x_i) \mathcal{CN}(x_i ;\hat{r}_{i \rightarrow a}^t , \hat{\Sigma}_{i \rightarrow a}^t) {\rm d}x_i \notag \\ & \overset{def}{=} G(x_i;\hat{r}_{i \rightarrow a}^t , \hat{\Sigma}_{i \rightarrow a}^t) \quad(27) \end{aligned} x^i→at+1​v^i→at+1​​=E(xi​)=∫xi​​xi​C1​p(xi​)CN(xi​;r^i→at​,Σ^i→at​)dxi​=defF(xi​;r^i→at​,Σ^i→at​)(26)=Var(xi​)=∫xi​​(xi​−x^i→at+1​)2C1​p(xi​)CN(xi​;r^i→at​,Σ^i→at​)dxi​=defG(xi​;r^i→at​,Σ^i→at​)(27)​其中， C = ∫ x i p ( x i ) C N ( x i ; r ^ i → a t , Σ ^ i → a t ) d x i C = \int_{x_i} p(x_i) \mathcal{CN} (x_i;\hat{r}_{i \rightarrow a}^t , \hat{\Sigma}_{i \rightarrow a}^t){\rm d}x_i C=∫xi​​p(xi​)CN(xi​;r^i→at​,Σ^i→at​)dxi​是归一化常数。当 p ( x i ) p(x_i) p(xi​)为确定分布时，这两个积分可以直接计算或近似，其复杂度为 O ( 1 ) \mathcal{O}(1) O(1)。

至此，两种节点之间的消息均推导完毕。到这里，实际上是对循环置信传播(Loopy Belief Propagation, LBP)进行了近似，减少了传递的参数量，将消息传递变成了统计量传递，但到目前为止，我们已经完成高斯化所有信息的处理。

3.3 进一步的简化

上述消息的计算均涉及到 O ( N ) \mathcal{O}(N) O(N)或 O ( M ) \mathcal{O}(M) O(M)（两者同等级，以下统一为 O ( N ) \mathcal{O}(N) O(N)）复杂度的计算过程，同时，在系统中，每次迭代存在 2 M N 2MN 2MN条消息，故实际上以上算法的复杂度约为 O ( N 3 ) \mathcal{O}(N^3) O(N3).考虑到三次复杂度依然难以接受，AMP对上述两类消息进行了进一步的简化。

我们约定， ∣ h i j ∣ |h_{ij}| ∣hij​∣的数量级为 O ( 1 / N ) O(1/\sqrt{N}) O(1/N ​)，因此 ∣ h i j ∣ 2 |h_{ij}|^2 ∣hij​∣2的数量级为 O ( 1 / N ) O(1/N) O(1/N)；同时约定 x j x_j xj​的数量级为 O ( 1 ) O(1) O(1)，因此其估计量 x ^ i \hat{x}_i x^i​和 v ^ i \hat{v}_i v^i​数量级同样为 O ( 1 ) O(1) O(1)。在以上约定下， y i y_i yi​的数量级同样的，也是 O ( 1 ) O(1) O(1)， Z a , V a Z_a, V_a Za​,Va​也同样应为 O ( 1 ) O(1) O(1)。因此，我们认为在在大系统极限下，单独存在的 O ( 1 / N ) O(1/\sqrt{N}) O(1/N ​)数量级或更低的变量以及求和后数量级为 O ( 1 / N ) O(1/\sqrt{N}) O(1/N ​)或更低可作为无穷小量而被忽略。

回忆公式(11)，我们提到了我们想要知道的后验分布 p t + 1 ( x i ∣ y ) p^{t+1}(x_i|\mathbf{y}) pt+1(xi​∣y)，我们定义：

Σ ^ i t = ( ∑ a = 1 M 1 v ^ i ← a t ) − 1 = ( ∑ a = 1 M ∣ h a i ∣ 2 σ 2 + V i ← a t ) − 1 ( 28 ) r ^ i t = Σ ^ i t ∑ a = 1 M h a i ∗ ( y a − Z i ← a t ) σ 2 + V i ← a t ( 29 ) \begin{aligned} \hat{\Sigma}_{i}^t & = \left(\sum_{a=1}^{M} \frac{1}{\hat{v}_{i \leftarrow a}^t}\right)^{-1} \notag \\ & = \left(\sum_{a=1}^{M} \frac{|h_{ai}|^2 }{\sigma^2 + V_{i \leftarrow a}^t}\right)^{-1} \quad(28) \\ \hat{r}_{i}^t & = \hat{\Sigma}_{i}^t \sum_{a = 1}^{M} \frac{h_{ai}^* (y_a - Z_{i \leftarrow a}^t) }{\sigma^2 + V_{i \leftarrow a}^t} \quad(29) \end{aligned} Σ^it​r^it​​=(a=1∑M​v^i←at​1​)−1=(a=1∑M​σ2+Vi←at​∣hai​∣2​)−1(28)=Σ^it​a=1∑M​σ2+Vi←at​hai∗​(ya​−Zi←at​)​(29)​ 其中公式(11)的 ∏ a = 1 M μ i ← a t ( x i ) \prod_{a=1}^{M} \mu^t_{i \leftarrow a}(x_i) ∏a=1M​μi←at​(xi​)正比于 C N ( x i ; r ^ i t , Σ ^ i t ) \mathcal{CN}(x_i; \hat{r}_i^t,\hat{\Sigma}_i^t ) CN(xi​;r^it​,Σ^it​)。相应的， p t + 1 ( x i ∣ y ) p^{t+1}(x_i|\mathbf{y}) pt+1(xi​∣y)的近似后验均值和方差可以表示为：

x ^ i t + 1 = F ( x i ; r ^ i t , Σ ^ i t ) ( 30 ) v ^ i t + 1 = G ( x i ; r ^ i t , Σ ^ i t ) ( 31 ) \begin{aligned} \hat{x}^{t+1}_{i} & = F(x_i;\hat{r}_{i}^t , \hat{\Sigma}_{i}^t) \quad(30) \\ \hat{v}^{t+1}_{i} & = G(x_i;\hat{r}_{i}^t , \hat{\Sigma}_{i}^t) \quad(31) \end{aligned} x^it+1​v^it+1​​=F(xi​;r^it​,Σ^it​)(30)=G(xi​;r^it​,Σ^it​)(31)​

同样地，我们定义 Z a t = ∑ i = 1 N h a i x ^ i → a t ( 32 ) V a t = ∑ i = 1 N ∣ h a i ∣ 2 v ^ i → a t ≈ V i ← a t ( 33 ) \begin{aligned} Z_a^t & = \sum_{i = 1}^{N} h_{ai} \hat{x}^t_{i \rightarrow a} \quad(32) \\ V_a^t & = \sum_{i = 1}^{N} |h_{ai}|^2 \hat{v}^t_{i \rightarrow a} \approx V_{i \leftarrow a}^t \quad(33) \end{aligned} Zat​Vat​​=i=1∑N​hai​x^i→at​(32)=i=1∑N​∣hai​∣2v^i→at​≈Vi←at​(33)​ 其中 ≈ \approx ≈是忽略了无穷小量。

对公式(26)中的 x ^ i → a t + 1 \hat{x}^{t+1}_{i \rightarrow a} x^i→at+1​应用一阶泰勒级数展开，有：

x ^ i → a t + 1 ≈ x ^ i t + 1 + Δ r ^ ∂ ∂ r F ( x i ; r ^ i t , Σ ^ i t ) + Δ Σ ^ ∂ ∂ Σ F ( x i ; r ^ i t , Σ ^ i t ) ( 34 ) \hat{x}^{t+1}_{i \rightarrow a} \approx \hat{x}^{t+1}_{i} + \Delta \hat{r} \frac{\partial}{\partial r} F(x_i;\hat{r}_{i}^t , \hat{\Sigma}_{i}^t) + \Delta \hat{\Sigma} \frac{\partial}{\partial \Sigma} F(x_i;\hat{r}_{i}^t , \hat{\Sigma}_{i}^t)\quad(34) x^i→at+1​≈x^it+1​+Δr^∂r∂​F(xi​;r^it​,Σ^it​)+ΔΣ^∂Σ∂​F(xi​;r^it​,Σ^it​)(34)

其中 Δ Σ = Σ ^ i → a t − Σ i t ^ = ( ∑ b ≠ a M ∣ h b i ∣ 2 σ 2 + V i ← b t ) − 1 − ( ∑ b = 1 M ∣ h b i ∣ 2 σ 2 + V i ← b t ) − 1 = ∣ h a i ∣ 2 σ 2 + V i ← a t ( ∑ b ≠ a M ∣ h b i ∣ 2 σ 2 + V i ← b t ) ( ∑ b = 1 M ∣ h b i ∣ 2 σ 2 + V i ← b t ) ≈ ∣ h a i ∣ 2 σ 2 + V a t ( ∑ b ≠ a M ∣ h b i ∣ 2 σ 2 + V i ← b t ) ( ∑ b = 1 M ∣ h b i ∣ 2 σ 2 + V i ← b t ) ≈ 0 ( 35 ) Δ r = r ^ i → a t − r ^ i t = Σ ^ i → a t ∑ b ≠ a M h b i ∗ ( y b − Z i ← b t ) σ 2 + V i ← b t − Σ ^ i t ∑ b = 1 M h b i ∗ ( y b − Z i ← b t ) σ 2 + V i ← b t ≈ − Σ ^ i t h a i ∗ ( y a − Z i ← a t ) σ 2 + V i ← a t ≈ − Σ ^ i t h a i ∗ ( y a − Z i ← a t ) σ 2 + V a t ( 36 ) \begin{aligned} \Delta \Sigma & = \hat{\Sigma}^t_{i \rightarrow a } - \hat{\Sigma^t_{i}} \notag \\ & = \left(\sum_{b \neq a}^{M} \frac{|h_{bi}|^2 }{\sigma^2 + V_{i \leftarrow b}^t}\right)^{-1} - \left(\sum_{b=1}^{M} \frac{|h_{bi}|^2 }{\sigma^2 + V_{i \leftarrow b}^t}\right)^{-1} \notag \\ & = \frac{ \frac{|h_{ai}|^2 }{\sigma^2 + V_{i \leftarrow a}^t} }{\left(\sum_{b \neq a}^{M} \frac{|h_{bi}|^2 }{\sigma^2 + V_{i \leftarrow b}^t}\right)\left(\sum_{b=1}^{M} \frac{|h_{bi}|^2 }{\sigma^2 + V_{i \leftarrow b}^t}\right)} \notag \\ & \approx \frac{ \frac{|h_{ai}|^2 }{\sigma^2 + V_{a}^t} }{\left(\sum_{b \neq a}^{M} \frac{|h_{bi}|^2 }{\sigma^2 + V_{i \leftarrow b}^t}\right)\left(\sum_{b=1}^{M} \frac{|h_{bi}|^2 }{\sigma^2 + V_{i \leftarrow b}^t}\right)} \notag \\ & \approx 0 \quad(35) \\ \Delta r & = \hat{r}^t_{i \rightarrow a } - \hat{r}^t_{i} \notag \\ & = \hat{\Sigma}^t_{i \rightarrow a }\sum_{b \neq a}^{M} \frac{h_{bi}^* (y_b - Z_{i \leftarrow b}^t) }{\sigma^2 + V_{i \leftarrow b}^t} - \hat{\Sigma}_{i}^t \sum_{b = 1}^{M} \frac{h_{bi}^* (y_b - Z_{i \leftarrow b}^t) }{\sigma^2 + V_{i \leftarrow b}^t} \notag\\ & \approx -\hat{\Sigma}_{i}^t \frac{h_{ai}^* (y_a - Z_{i \leftarrow a}^t) }{\sigma^2 + V_{i \leftarrow a}^t} \notag\\ & \approx -\hat{\Sigma}_{i}^t \frac{h_{ai}^* (y_a - Z_{i \leftarrow a}^t) }{\sigma^2 + V_{a}^t} \quad(36) \end{aligned} ΔΣΔr​=Σ^i→at​−Σit​^​= ​b=a∑M​σ2+Vi←bt​∣hbi​∣2​ ​−1−(b=1∑M​σ2+Vi←bt​∣hbi​∣2​)−1=(∑b=aM​σ2+Vi←bt​∣hbi​∣2​)(∑b=1M​σ2+Vi←bt​∣hbi​∣2​)σ2+Vi←at​∣hai​∣2​​≈(∑b=aM​σ2+Vi←bt​∣hbi​∣2​)(∑b=1M​σ2+Vi←bt​∣hbi​∣2​)σ2+Vat​∣hai​∣2​​≈0(35)=r^i→at​−r^it​=Σ^i→at​b=a∑M​σ2+Vi←bt​hbi∗​(yb​−Zi←bt​)​−Σ^it​b=1∑M​σ2+Vi←bt​hbi∗​(yb​−Zi←bt​)​≈−Σ^it​σ2+Vi←at​hai∗​(ya​−Zi←at​)​≈−Σ^it​σ2+Vat​hai∗​(ya​−Zi←at​)​(36)​

其中利用到了 V a t = V i ← a t + O ( 1 / N ) V_a^t = V_{i \leftarrow a}^t + O(1/N) Vat​=Vi←at​+O(1/N)和 Σ ^ i t = Σ ^ i → a t + O ( 1 / N ) \hat{\Sigma}_i^t =\hat{\Sigma}_{i \rightarrow a}^t + O(1/N) Σ^it​=Σ^i→at​+O(1/N)的近似。利用后验概率密度函数的偏导数求法，可以得到 ∂ ∂ r F ( x i ; r , Σ i t ) ∣ r = r ^ i t = v ^ i t + 1 Σ ^ i t \frac{\partial}{\partial r} F(x_i;r,\Sigma^t_i)|_{r = \hat{r}^t_i} = \frac{\hat{v}_i^{t+1}}{\hat{\Sigma}_i^t} ∂r∂​F(xi​;r,Σit​)∣r=r^it​​=Σ^it​v^it+1​​，从而公式(34)可以化简为： x ^ i → a t + 1 ≈ x ^ i t + 1 − h a i ∗ ( y a − Z i ← a t ) σ 2 + V a t v ^ i t + 1 ( 37 ) \hat{x}^{t+1}_{i \rightarrow a} \approx \hat{x}_i^{t+1} - \frac{h_{ai}^*(y_a - Z^t_{i \leftarrow a})}{\sigma^2+V_a^t}\hat{v}_i^{t+1} \quad(37) x^i→at+1​≈x^it+1​−σ2+Vat​hai∗​(ya​−Zi←at​)​v^it+1​(37)

同样地，我们对公式(27)进行一阶泰勒展开，有： v ^ i → a t + 1 ≈ v ^ i t + 1 + Δ r ∂ ∂ r G ( x i ; r ^ i t , Σ ^ i t ) ( 38 ) \hat{v}^{t+1}_{i \rightarrow a} \approx \hat{v}_i^{t+1} + \Delta r \frac{\partial}{\partial r}G(x_i;\hat{r}_i^t,\hat{\Sigma}^t_i)\quad(38) v^i→at+1​≈v^it+1​+Δr∂r∂​G(xi​;r^it​,Σ^it​)(38)

将公式(36)和(38)带入(33)得到： V a t ≈ ∑ i = 1 N ∣ h a i ∣ 2 ( v ^ i t − Σ ^ i t h a i ∗ ( y a − Z i ← a t ) σ 2 + V a t × ∂ ∂ r G ( x i ; r ^ i t , Σ ^ i t ) ) ≈ ∑ i = 1 N ∣ h a i ∣ 2 v ^ i t − ∑ i = 1 N ∣ h a i ∣ 2 ( ∑ b = 1 M ∣ h b i ∣ 2 σ 2 + V i ← b t ) − 1 h a i ∗ ( y a − Z i ← a t ) σ 2 + V a t × ∂ ∂ r G ( x i ; r ^ i t , Σ ^ i t ) = ∑ i = 1 N ∣ h a i ∣ 2 v ^ i t + O ( 1 / N ) ≈ ∑ i = 1 N ∣ h a i ∣ 2 v ^ i t ( 39 ) \begin{split} V_a^t & \approx \sum_{i=1}^{N} |h_{ai}|^2 \left( \hat{v}_i^t -\hat{\Sigma}_{i}^t \frac{h_{ai}^* (y_a - Z_{i \leftarrow a}^t) }{\sigma^2 + V_{a}^t} \times \frac{\partial}{\partial r}G(x_i;\hat{r}_i^t,\hat{\Sigma}^t_i) \right) \\ & \approx \sum_{i=1}^{N} |h_{ai}|^2 \hat{v}_i^t - \sum_{i=1}^{N} |h_{ai}|^2 \left( \sum_{b=1}^{M}\frac{|h_{bi}|^2}{\sigma^2 + V^t_{i \leftarrow b}} \right)^{-1} \frac{h_{ai}^* (y_a - Z_{i \leftarrow a}^t)}{\sigma^2 + V_{a}^t} \times \frac{\partial}{\partial r}G(x_i;\hat{r}_i^t,\hat{\Sigma}^t_i) \\ & = \sum_{i=1}^{N} |h_{ai}|^2 \hat{v}_i^t + O(1/\sqrt{N}) \\ & \approx \sum_{i=1}^{N} |h_{ai}|^2 \hat{v}_i^t \quad(39) \end{split} Vat​​≈i=1∑N​∣hai​∣2(v^it​−Σ^it​σ2+Vat​hai∗​(ya​−Zi←at​)​×∂r∂​G(xi​;r^it​,Σ^it​))≈i=1∑N​∣hai​∣2v^it​−i=1∑N​∣hai​∣2(b=1∑M​σ2+Vi←bt​∣hbi​∣2​)−1σ2+Vat​hai∗​(ya​−Zi←at​)​×∂r∂​G(xi​;r^it​,Σ^it​)=i=1∑N​∣hai​∣2v^it​+O(1/N ​)≈i=1∑N​∣hai​∣2v^it​(39)​

其中 O ( 1 / N ) O(1/\sqrt{N}) O(1/N ​)的来源是 h a i ∗ h_{ai}^* hai∗​。再将公式(37)代入(32)中，可得： Z a t ≈ ∑ i = 1 N h a i x ^ i t − ∑ i = 1 N ∣ h a i ∣ 2 ( y a − Z i ← a t − 1 ) σ 2 + V a t − 1 v ^ i t = ∑ i = 1 N h a i x ^ i t − ∑ i = 1 N ∣ h a i ∣ 2 v ^ i t ( y a − Z a t − 1 + h a i x ^ i t − 1 ) σ 2 + V a t − 1 ≈ ∑ i = 1 N h a i x ^ i t − V a t ( y a − Z a t − 1 ) σ 2 + V a t − 1 ( 40 ) \begin{split} Z_a^t & \approx \sum_{i=1}^{N} h_{ai} \hat{x}_i^t - \sum_{i=1}^{N} \frac{|h_{ai}|^2(y_a - Z_{i \leftarrow a}^{t-1})}{\sigma^2 + V_{a}^{t-1}} \hat{v}_{i}^{t} \\ & = \sum_{i=1}^{N} h_{ai} \hat{x}_i^t - \sum_{i=1}^{N} \frac{|h_{ai}|^2\hat{v}_{i}^{t}(y_a - Z_{a}^{t-1} + h_{ai}\hat{x}_i^{t-1})}{\sigma^2 + V_{a}^{t-1}} \\ & \approx \sum_{i=1}^{N} h_{ai} \hat{x}_i^t - \frac{V_a^t(y_a - Z_{a}^{t-1})}{\sigma^2 + V_{a}^{t-1}} \quad(40) \end{split} Zat​​≈i=1∑N​hai​x^it​−i=1∑N​σ2+Vat−1​∣hai​∣2(ya​−Zi←at−1​)​v^it​=i=1∑N​hai​x^it​−i=1∑N​σ2+Vat−1​∣hai​∣2v^it​(ya​−Zat−1​+hai​x^it−1​)​≈i=1∑N​hai​x^it​−σ2+Vat−1​Vat​(ya​−Zat−1​)​(40)​

其中最后一步的化简是因为 h a i x ^ i t − 1 = O ( 1 / N ) h_{ai}\hat{x}_i^{t-1}=O(1/\sqrt{N}) hai​x^it−1​=O(1/N ​)。再将公式(37)代入(29)中： r ^ i t ≈ Σ ^ i t ∑ a = 1 M h a i ∗ ( y a − Z a t + h a i x ^ i t ) σ 2 + V a t = x ^ i t + Σ ^ i t ∑ a = 1 M h a i ∗ ( y a − Z a t ) σ 2 + V a t ( 41 ) \begin{split} \hat{r}_{i}^t & \approx \hat{\Sigma}_{i}^t \sum_{a = 1}^{M} \frac{h_{ai}^* (y_a - Z_{a}^{t} + h_{ai}\hat{x}_i^{t}) }{\sigma^2 + V_{a}^t} \\ & = \hat{x}^t_i + \hat{\Sigma}_i^t \sum_{a=1}^{M} \frac{h_{ai}^* (y_a - Z_{a}^{t} )}{\sigma^2 + V_{a}^t} \quad(41) \end{split} r^it​​≈Σ^it​a=1∑M​σ2+Vat​hai∗​(ya​−Zat​+hai​x^it​)​=x^it​+Σ^it​a=1∑M​σ2+Vat​hai∗​(ya​−Zat​)​(41)​ 到此为止，我们完成AMP的推导。

4 AMP伪算法及总结

最终，我们把算法表示在这里：

需要注意的是，此处的上标 T T T不是转置(transpose)，而是第 T T T迭代之后的结果。还需要强调的是(30)和(31)其实就是后验分布的均值和方差，在伪算法中给出了贝叶斯的形式，注意这里是对后验分布求的均值和方差。 r ^ i t \hat{r}_{i}^t r^it​和 Σ ^ i t \hat{\Sigma}_{i}^t Σ^it​是似然函数（复高斯分布）的均值和方差。还有一些文献会引入残差项（一般为 s ^ a t \hat{s}_a^t s^at​）来表示 ( y a − Z a t ) / ( σ 2 + V a t ) (y_a - Z_a^t)/(\sigma^2+V_a^t) (ya​−Zat​)/(σ2+Vat​)等。

本文大量参考了

https://etc.lmdewz.xyz/?p=94

https://arxiv.org/abs/2201.07487