深入浅出GAMP算法（下）：MMSE估计和AWGN场景

前言

在前两篇博客中， 我们分别讲述了消息传递算法的来龙去脉 和 利用 高斯及泰勒展开近似得到的最大后验估计的GAMP版本。 这一篇博客，我们使用类似的推导，整理了在实际中可能更常用的， MMSE 最小均方误差估计版本的 GAMP 算法。

模型背景

我们旨在解决上图这样的问题， 已知输入 q \mathbf{q} q （先验信息）， 已知输出 y \mathbf{y} y（后验信息)， 已知变换矩阵 A \mathbf{A} A， 反推出变量 x \mathbf{x} x。 以AWGN信道举例： y = z + w = A x + w \mathbf{y}=\mathbf{z}+\mathbf{w}=\mathbf{A} \mathbf{x}+\mathbf{w} y=z+w=Ax+w y \mathbf{y} y已知而我们试图反推出 x \mathbf{x} x。 此时， 如果 x \mathbf{x} x有一些先验分布信息， 如稀疏分布， 即 x \mathbf{x} x中只有少量元素非零， 那么， 便可以通过概率的方式， 以GAMP算法进行求解。

MMSE

上一篇博客中我们写到的 MAP 版本的 GAMP 旨在给出以下的估计量： x ^ map  : = arg ⁡ max ⁡ x ∈ R n F ( x , z , q , y ) , z ^ = A x ^ \widehat{\mathbf{x}}^{\text {map }}:=\underset{\mathbf{x} \in \mathbb{R}^{n}}{\arg \max } F(\mathbf{x}, \mathbf{z}, \mathbf{q}, \mathbf{y}), \quad \widehat{\mathbf{z}}=\mathbf{A} \widehat{\mathbf{x}} x map :=x∈Rnargmax​F(x,z,q,y),z =Ax 即最大化后验概率。 我举一个例子， 比如经过计算得到， x = 0.5 x=0.5 x=0.5 的概率是 0.2 0.2 0.2， 而 x = 0.1 x=0.1 x=0.1的概率是 0.19 0.19 0.19， x = 0.11 x=0.11 x=0.11的概率是 0.18 0.18 0.18， x = 0.09 x=0.09 x=0.09的概率是 0.17 0.17 0.17。 此时， 如果使用最大后验 MAP 估计， 得到的 x x x估计结果就是 x = 0.5 x=0.5 x=0.5。 但在这种情况下这就有失偏颇， 因为他完全没有考虑其他情况的概率。 这就有点像通信中的硬判决，一刀切， 因此，综合了所有可能情况的软判决， 可能更显客观， 这也是我们要介绍的MMSE估计， 即： x ^ m m s e : = E [ x ∣ y , q ] \widehat{\mathrm{x}}^{\mathrm{mmse}}:=\mathbb{E}[\mathrm{x} \mid \mathrm{y}, \mathbf{q}] x mmse:=E[x∣y,q] 以期望作为估计值。下面， 我们就推导以MMSE为估计准则的GAMP版本。

回顾下MAP中的消息传递： 在每次迭代中， 实际要计算两个消息算子： Δ i → j ( t , x j ) = c o n s t + max ⁡ x \ x j f out  ( z i , y i ) + ∑ r ≠ j Δ i ← r ( t , x r ) (1) \begin{aligned} {\Delta}_{i \rightarrow j}\left(t, x_{j}\right)=\mathrm{const} +\max _{\mathbf{x}\backslash x_j} f_{\text {out }}\left(z_{i}, y_{i}\right)+\sum_{r \neq j} \Delta_{i \leftarrow r}\left(t, x_{r}\right) \end{aligned}\tag{1} Δi→j​(t,xj​)=const+x\xj​max​fout ​(zi​,yi​)+r​=j∑​Δi←r​(t,xr​)​(1) 和 Δ i ← j ( t + 1 , x j ) = c o n s t + f i n ( x j , q j ) + ∑ ℓ ≠ i Δ ℓ → j ( t , x j ) (2) \begin{aligned} {\Delta}_{i \leftarrow j}\left(t+1, x_{j}\right)=\mathrm{const} +f_{\mathrm{in}}\left(x_{j}, q_{j}\right)+\sum_{\ell \neq i} \Delta_{\ell \rightarrow j}\left(t, x_{j}\right) \end{aligned}\tag{2} Δi←j​(t+1,xj​)=const+fin​(xj​,qj​)+ℓ​=i∑​Δℓ→j​(t,xj​)​(2)

注意到(1)式中， 由于是MAP准则， 因此传递的消息里也是求了 max ⁡ \max max后的结果。 然而在MMSE准则下， 我们要传递的消息就变成了： Δ i → j ( t , x j ) = log ⁡ E ( p Y ∣ Z ( y i , z i ) ∣ x j ) +  const  (3) \Delta_{i \rightarrow j}\left(t, x_{j}\right)=\log \mathbb{E}\left(p_{Y \mid Z}\left(y_{i}, z_{i}\right) \mid x_{j}\right)+\text { const }\tag{3} Δi→j​(t,xj​)=logE(pY∣Z​(yi​,zi​)∣xj​)+ const (3) 和 Δ i ← j ( x j ) =  const  + log ⁡ p X ∣ Q ( x j ∣ q j ) + ∑ l ≠ i Δ l → j ( x j ) (4) \Delta_{i \leftarrow j}\left(x_{j}\right)=\text { const }+\log p_{X \mid Q}\left(x_{j} \mid q_{j}\right)+\sum_{l \neq i} \Delta_{l \rightarrow j}\left(x_{j}\right) \tag{4} Δi←j​(xj​)= const +logpX∣Q​(xj​∣qj​)+l​=i∑​Δl→j​(xj​)(4) 注意到， 在这个消息传递中， 传递的都是似然信息了， 也就是说有： p i ← r ( x r ) ∝ exp ⁡ Δ i ← r ( t , x r ) p_{i \leftarrow r}\left(x_{r}\right) \propto \exp \Delta_{i \leftarrow r}\left(t, x_{r}\right) pi←r​(xr​)∝expΔi←r​(t,xr​) 值得一提的是， 在(3)中的期望是对随机变量 z i = a i T x z_{i}=\mathbf{a}_{i}^{T} \mathbf{x} zi​=aiT​x 进行求取， 此时 x j x_j xj​ 是 固定的， 而 x \mathbf{x} x 中的其余项， 则服从独立分布 p i ← r ( x r ) ∝ exp ⁡ Δ i ← r ( t , x r ) p_{i \leftarrow r}\left(x_{r}\right) \propto \exp \Delta_{i \leftarrow r}\left(t, x_{r}\right) pi←r​(xr​)∝expΔi←r​(t,xr​)。

和 MAP估计一样，我们最后要得到的估计值也是由下式得出： p j ( x j ) ∝ exp ⁡ Δ j ( t , x j ) p_{j}\left(x_{j}\right) \propto \exp \Delta_{j}\left(t, x_{j}\right) pj​(xj​)∝expΔj​(t,xj​) 其中 Δ j ( t + 1 , x j ) = f i n ( x j , q j ) + ∑ i Δ i → j ( t , x j ) \Delta_{j}\left(t+1, x_{j}\right)=f_{\mathrm{in}}\left(x_{j}, q_{j}\right)+\sum_{i} \Delta_{i \rightarrow j}\left(t, x_{j}\right) Δj​(t+1,xj​)=fin​(xj​,qj​)+i∑​Δi→j​(t,xj​)

那么，后续我们就是通过合理的近似， 对消息传递算法进行简化。 类似地， 我们先定义如下变量： 我们需要定义的变量也变为： x ^ j ( t ) : = E [ x j ∣ Δ j ( t , ⋅ ) ] x ^ i ← j ( t ) : = E [ x j ∣ Δ i ← j ( t , ⋅ ) ] τ j x ( t ) : = var ⁡ [ x j ∣ Δ j ( t , ⋅ ) ] τ i ← j x ( t ) : = var ⁡ [ x j ∣ Δ i ← j ( t , ⋅ ) ] , \begin{aligned} \widehat{x}_{j}(t) &:=\mathbb{E}\left[x_{j} \mid \Delta_{j}(t, \cdot)\right] \\ \widehat{x}_{i \leftarrow j}(t) &:=\mathbb{E}\left[x_{j} \mid \Delta_{i \leftarrow j}(t, \cdot)\right] \\ \tau_{j}^{x}(t) &:=\operatorname{var}\left[x_{j} \mid \Delta_{j}(t, \cdot)\right] \\ \tau_{i \leftarrow j}^{x}(t) &:=\operatorname{var}\left[x_{j} \mid \Delta_{i \leftarrow j}(t, \cdot)\right], \end{aligned} x j​(t)x i←j​(t)τjx​(t)τi←jx​(t)​:=E[xj​∣Δj​(t,⋅)]:=E[xj​∣Δi←j​(t,⋅)]:=var[xj​∣Δj​(t,⋅)]:=var[xj​∣Δi←j​(t,⋅)],​ 其中 x ^ j \hat{x}_j x^j​ 就是我们最后要得到的估计量。

我们先对 (3) 进行近似， 这可以写为： Δ i → j ( t , x j ) =  const  + log ⁡ ∫ { x r } r ≠ j p Y ∣ Z ( y i ∣ a i j x j + ∑ r ≠ j a i r x r ⏟ ≜ z i ) ∏ r ≠ j e Δ i ← r ( t , x r ) . \begin{aligned} &\Delta_{i \rightarrow j}\left(t, x_{j}\right) \\ &=\text { const }+\log \int_{\left\{x_{r}\right\}_{r \neq j}} p_{Y \mid Z}(y_{i} \mid \underbrace{a_{i j} x_{j}+\sum_{r \neq j} a_{i r} x_{r}}_{\triangleq z_{i}}) \prod_{r \neq j} e^{\Delta_{i \leftarrow r}\left(t, x_{r}\right)} . \end{aligned} ​Δi→j​(t,xj​)= const +log∫{xr​}r​=j​​pY∣Z​(yi​∣≜zi​ aij​xj​+r​=j∑​air​xr​​​)r​=j∏​eΔi←r​(t,xr​).​ 而当矩阵维度较大时， 那么，根据中心极限定理， 相互独立的随机变量， 其和 服从 高斯分布， 且均值就是 各自均值之和， 而方差则是各自方差之和， 因此， 有： z i ∣ x j ∼ N ( a i j x j + p ^ i ← j ( t ) , τ i ← j p ( t ) ) z_{i} \mid x_{j} \sim \mathcal{N}\left(a_{i j} x_{j}+\widehat{p}_{i\leftarrow j}(t), \tau_{i\leftarrow j}^{p}(t)\right) zi​∣xj​∼N(aij​xj​+p ​i←j​(t),τi←jp​(t)) 其中， p ^ i → j ( t ) : = ∑ r ≠ j a i r x ^ i ← r ( t ) τ i → j p ( t ) : = ∑ r ≠ j ∣ a i r ∣ 2 τ r x ( t ) \begin{aligned} \widehat{p}_{i \rightarrow j}(t) &:=\sum_{r \neq j} a_{i r} \widehat{x}_{i \leftarrow r}(t) \\ \tau_{i \rightarrow j}^{p}(t) &:=\sum_{r \neq j}\left|a_{i r}\right|^{2} \tau_{r}^{x}(t) \end{aligned} p ​i→j​(t)τi→jp​(t)​:=r​=j∑​air​x i←r​(t):=r​=j∑​∣air​∣2τrx​(t)​ 那么， 我们进一步就有： Δ i → j ( t , x j ) ≈  const  + log ⁡ ∫ z i p Y ∣ Z ( y i ∣ z i ) N ( z i ; a i j x j + p ^ i ← j ( t ) , τ i ← j p ( t ) ) ⏟ ≜ H ( a i j x j + p ^ i ← j ( t ) , y i , τ i ← j p ( t ) ) (5) \Delta_{i \rightarrow j}\left(t, x_{j}\right) \approx \text { const }+\underbrace{\log \int_{z_{i}} p_{Y \mid Z}\left(y_{i} \mid z_{i}\right) \mathcal{N}\left(z_{i} ; a_{i j} x_{j}+\widehat{p}_{i\leftarrow j}(t), \tau_{i\leftarrow j}^{p}(t)\right)}_{\triangleq H\left(a_{i j} x_{j}+\widehat{p}_{i\leftarrow j}(t), y_{i}, \tau_{i\leftarrow j}^{p}(t)\right)}\tag{5} Δi→j​(t,xj​)≈ const +≜H(aij​xj​+p ​i←j​(t),yi​,τi←jp​(t)) log∫zi​​pY∣Z​(yi​∣zi​)N(zi​;aij​xj​+p ​i←j​(t),τi←jp​(t))​​(5) 这里，我们有： H ( p ^ , y , μ p ) ≜ log ⁡ ∫ p Y ∣ Z ( y ∣ z ) N ( z ; p ^ , μ p ) d z H\left(\widehat{p}, y, \mu^{p}\right) \triangleq \log \int p_{Y \mid Z}(y \mid z) \mathcal{N}\left(z ; \widehat{p}, \mu^{p}\right) d z H(p ​,y,μp)≜log∫pY∣Z​(y∣z)N(z;p ​,μp)dz ( μ p \mu^p μp 和 τ p \tau^p τp是一样的） 那么 类似 MAP版本的思路， 接下来我们要把所有 i ← j i\leftarrow j i←j 项替换掉， 定义变量： p ^ i ( t ) : = ∑ j a i j x ^ i ← j ( t ) = p ^ i → j + a i j x ^ i ← j ( t ) ， τ i p ( t ) = ∑ j ∣ a i j ∣ 2 τ j x ( t ) + a i j 2 τ j x ( t ) \widehat{p}_{i}(t):=\sum_{j} a_{i j} \widehat{x}_{i \leftarrow j}(t)=\widehat{p}_{i \rightarrow j}+a_{ij} \widehat{x}_{i \leftarrow j}(t)， \tau_{i}^{p}(t)=\sum_{j}\left|a_{i j}\right|^{2} \tau_{j}^{x}(t) + a_{ij}^2\tau_j^x(t) p ​i​(t):=j∑​aij​x i←j​(t)=p ​i→j​+aij​x i←j​(t)，τip​(t)=j∑​∣aij​∣2τjx​(t)+aij2​τjx​(t) 那么，(5)可以被写为 Δ i → j ( t , x j ) ≈ H ( p ^ i ( t ) + a i j ( x j − x ^ j ) , y i , τ i p ( t ) ) +  const.  \Delta_{i \rightarrow j}\left(t, x_{j}\right) \approx H\left(\widehat{p}_{i}(t)+a_{i j}\left(x_{j}-\widehat{x}_{j}\right), y_{i}, \tau_{i}^{p}(t)\right)+\text { const. } Δi→j​(t,xj​)≈H(p ​i​(t)+aij​(xj​−x j​),yi​,τip​(t))+ const.  可以通过泰勒展开得到： Δ i → j ( t , x j ) ≈  const  + s i ( t ) a i j ( x j − x ^ j ( t ) ) − τ i s ( t ) 2 a i j 2 ( x j − x ^ j ( t ) ) 2 = const ⁡ [ s i ( t ) a i j + a i j 2 τ i s ( t ) x ^ j ( t ) ] x j − τ i s ( t ) 2 a i j 2 x j 2 (6) \begin{aligned} \Delta_{i \rightarrow j}(&\left.t, x_{j}\right) \approx \text { const } \\ &+s_{i}(t) a_{i j}\left(x_{j}-\widehat{x}_{j}(t)\right)-\frac{\tau_{i}^{s}(t)}{2} a_{i j}^{2}\left(x_{j}-\widehat{x}_{j}(t)\right)^{2} \\ =& \operatorname{const}\left[s_{i}(t) a_{i j}+a_{i j}^{2} \tau_{i}^{s}(t) \widehat{x}_{j}(t)\right] x_{j} \\ &-\frac{\tau_{i}^{s}(t)}{2} a_{i j}^{2} x_{j}^{2} \end{aligned}\tag{6} Δi→j​(=​t,xj​)≈ const +si​(t)aij​(xj​−x j​(t))−2τis​(t)​aij2​(xj​−x j​(t))2const[si​(t)aij​+aij2​τis​(t)x j​(t)]xj​−2τis​(t)​aij2​xj2​​(6) 这两个表达式和之前一样，这里的近似是我们忽略了 O ( a i j 2 ) O\left(a_{i j}^{2}\right) O(aij2​)级的项。其中， s ^ i ( t ) = g out  ( t , p ^ i ( t ) , y i , τ i p ( t ) ) τ i s ( t ) = − ∂ ∂ p ^ g out  ( t , p ^ i ( t ) , y i , τ i p ( t ) ) \begin{aligned} \widehat{s}_{i}(t) &=g_{\text {out }}\left(t, \widehat{p}_{i}(t), y_{i}, \tau_{i}^{p}(t)\right) \\ \tau_{i}^{s}(t) &=-\frac{\partial}{\partial \widehat{p}} g_{\text {out }}\left(t, \widehat{p}_{i}(t), y_{i}, \tau_{i}^{p}(t)\right) \end{aligned} s i​(t)τis​(t)​=gout ​(t,p ​i​(t),yi​,τip​(t))=−∂p ​∂​gout ​(t,p ​i​(t),yi​,τip​(t))​ 这些都是和MAP版本的定义完全一致， 然而由于 H H H函数的不同， g o u t g_\mathrm{out} gout​的形式也截然不同。 因此，接下来我们要对 g o u t g_\mathrm{out} gout​进行推导， 也即推导 H H H的一阶导： H ′ ( p ^ , y , μ p ) = ∂ ∂ p ^ log ⁡ ∫ p Y ∣ Z ( y ∣ z ) 1 2 π μ p exp ⁡ ( − 1 2 μ p ( z − p ^ ) 2 ) d z = ∂ ∂ p ^ { log ⁡ 1 2 π μ p + log ⁡ ∫ z exp ⁡ ( log ⁡ p Y ∣ Z ( y ∣ z ) − 1 2 μ p ( z − p ^ ) 2 ) d z } = ∂ ∂ p ^ { − p ^ 2 2 μ p + log ⁡ ∫ exp ⁡ ( log ⁡ p Y ∣ Z ( y ∣ z ) − z 2 2 μ p + p ^ z μ p ) d z } = − p ^ μ p + ∂ ∂ p ^ log ⁡ [ μ p ∫ exp ⁡ ( ϕ ( u ) + p ^ u ) d u ]  via  u ≜ z μ p = − p ^ μ p + ∂ ∂ p ^ log ⁡ ∫ exp ⁡ ( ϕ ( u ) + p ^ u ) d u \begin{aligned} &H^{\prime}\left(\widehat{p}, y, \mu^{p}\right) \\ &\quad=\frac{\partial}{\partial \widehat{p}} \log \int p_{Y \mid Z}(y \mid z) \frac{1}{\sqrt{2 \pi \mu^{p}}} \exp \left(-\frac{1}{2 \mu^{p}}(z-\widehat{p})^{2}\right) d z \\ &=\frac{\partial}{\partial \widehat{p}}\left\{\log \frac{1}{\sqrt{2 \pi \mu^{p}}}+\log \int_{z} \exp \left(\log p_{Y \mid Z}(y \mid z)-\frac{1}{2 \mu^{p}}(z-\widehat{p})^{2}\right) d z\right\} \\ &=\frac{\partial}{\partial \widehat{p}}\left\{-\frac{\widehat{p}^{2}}{2 \mu^{p}}+\log \int \exp \left(\log p_{Y \mid Z}(y \mid z)-\frac{z^{2}}{2 \mu^{p}}+\frac{\widehat{p} z}{\mu^{p}}\right) d z\right\} \\ &=-\frac{\widehat{p}}{\mu^{p}}+\frac{\partial}{\partial \widehat{p}} \log \left[\mu^{p} \int \exp (\phi(u)+\widehat{p} u) d u\right] \text { via } u \triangleq \frac{z}{\mu^{p}} \\ &=-\frac{\widehat{p}}{\mu^{p}}+\frac{\partial}{\partial \widehat{p}} \log \int \exp (\phi(u)+\widehat{p} u) d u \end{aligned} ​H′(p ​,y,μp)=∂p ​∂​log∫pY∣Z​(y∣z)2πμp ​1​exp(−2μp1​(z−p ​)2)dz=∂p ​∂​{log2πμp ​1​+log∫z​exp(logpY∣Z​(y∣z)−2μp1​(z−p ​)2)dz}=∂p ​∂​{−2μpp ​2​+log∫exp(logpY∣Z​(y∣z)−2μpz2​+μpp ​z​)dz}=−μpp ​​+∂p ​∂​log[μp∫exp(ϕ(u)+p ​u)du] via u≜μpz​=−μpp ​​+∂p ​∂​log∫exp(ϕ(u)+p ​u)du​ 记 Z ( p ^ ) ≜ ∫ exp ⁡ ( ϕ ( u ) + p ^ u ) d u Z(\widehat{p}) \triangleq \int \exp (\phi(u)+\widehat{p} u) d u Z(p ​)≜∫exp(ϕ(u)+p ​u)du， 我们有如下数学公理： ∂ ∂ p ^ log ⁡ Z ( p ^ ) = E { u ∣ p ^ }  with  p U ∣ P ( u ∣ p ^ ) = exp ⁡ ( ϕ ( u ) + p ^ u ) Z ( p ^ ) ∂ 2 ∂ p ^ 2 log ⁡ Z ( p ^ ) = var ⁡ { u ∣ p ^ }  with  p U ∣ P ( u ∣ p ^ ) = exp ⁡ ( ϕ ( u ) + p ^ u Z ( p ^ ) . \begin{aligned} &\frac{\partial}{\partial \widehat{p}} \log Z(\widehat{p})=\mathrm{E}\{u \mid \widehat{p}\} \text { with } p_{U \mid P}(u \mid \widehat{p})=\frac{\exp (\phi(u)+\hat{p} u)}{Z(\hat{p})} \\ &\frac{\partial^{2}}{\partial \widehat{p}^{2}} \log Z(\widehat{p})=\operatorname{var}\{u \mid \widehat{p}\} \text { with } p_{U \mid P}(u \mid \widehat{p})=\frac{\exp (\phi(u)+\hat{p} u}{Z(\hat{p})} . \end{aligned} ​∂p ​∂​logZ(p ​)=E{u∣p ​} with pU∣P​(u∣p ​)=Z(p^​)exp(ϕ(u)+p^​u)​∂p ​2∂2​logZ(p ​)=var{u∣p ​} with pU∣P​(u∣p ​)=Z(p^​)exp(ϕ(u)+p^​u​.​ 可以通过简单的求导法则验证， 这是正确的。 因此， H ′ ( p ^ , y , μ p ) = − p ^ μ p + ∫ u exp ⁡ ( ϕ ( u ) + p ^ u ) Z ( p ^ ) d u = − p ^ μ p + ∫ z μ p exp ⁡ ( log ⁡ p Y ∣ Z ( y ∣ z ) − z 2 2 μ p + z p ^ μ p ) Z ( p ^ ) d z μ p  via  u ≜ z μ p = − p ^ μ p + 1 μ p ∫ z exp ⁡ ( log ⁡ p Y ∣ Z ( y ∣ z ) − 1 2 μ p ( z − p ^ ) 2 ) μ p Z ( p ^ ) exp ⁡ ( − p 2 2 μ p ) d z = − p ^ μ p + 1 μ p ∫ z p Y ∣ Z ( y ∣ z ) N ( z ; p ^ , μ p ) ∫ p Y ∣ Z ( y ∣ z ˉ ) N ( z ˉ ; p ^ , μ p ) d z ˉ d z = 1 μ p ( E { z ∣ y , p ^ ; μ p } − p ^ ) \begin{aligned} H^{\prime}\left(\widehat{p}, y, \mu^{p}\right) &=-\frac{\widehat{p}}{\mu^{p}}+\int u \frac{\exp (\phi(u)+\widehat{p} u)}{Z(\hat{p})} d u \\ &=-\frac{\widehat{p}}{\mu^{p}}+\int \frac{z}{\mu^{p}} \frac{\exp \left(\log p_{Y \mid Z}(y \mid z)-\frac{z^{2}}{2 \mu^{p}}+\frac{z \widehat{p}}{\mu^{p}}\right)}{Z(\widehat{p})} \frac{d z}{\mu^{p}} \text { via } u \triangleq \frac{z}{\mu^{p}} \\ &=-\frac{\widehat{p}}{\mu^{p}}+\frac{1}{\mu^{p}} \int z \frac{\exp \left(\log p_{Y \mid Z}(y \mid z)-\frac{1}{2 \mu^{p}}(z-\widehat{p})^{2}\right)}{\mu^{p} Z(\widehat{p}) \exp \left(-\frac{p^{2}}{2 \mu^{p}}\right)} d z \\ &=-\frac{\widehat{p}}{\mu^{p}}+\frac{1}{\mu^{p}} \int z \frac{p_{Y \mid Z}(y \mid z) \mathcal{N}\left(z ; \widehat{p}, \mu^{p}\right)}{\int p_{Y \mid Z}(y \mid \bar{z}) \mathcal{N}\left(\bar{z} ; \hat{p}, \mu^{p}\right) d \bar{z}} d z \\ &=\frac{1}{\mu^{p}}\left(\mathrm{E}\left\{z \mid y, \widehat{p} ; \mu^{p}\right\}-\widehat{p}\right) \end{aligned} H′(p ​,y,μp)​=−μpp ​​+∫uZ(p^​)exp(ϕ(u)+p ​u)​du=−μpp ​​+∫μpz​Z(p ​)exp(logpY∣Z​(y∣z)−2μpz2​+μpzp ​​)​μpdz​ via u≜μpz​=−μpp ​​+μp1​∫zμpZ(p ​)exp(−2μpp2​)exp(logpY∣Z​(y∣z)−2μp1​(z−p ​)2)​dz=−μpp ​​+μp1​∫z∫pY∣Z​(y∣zˉ)N(zˉ;p^​,μp)dzˉpY∣Z​(y∣z)N(z;p ​,μp)​dz=μp1​(E{z∣y,p ​;μp}−p ​)​ 因此，我们有： g out  ( p ^ , y , τ p ) : = 1 τ p ( z ^ 0 − p ^ ) , z ^ 0 : = E ( z ∣ p ^ , y , τ p ) g_{\text {out }}\left(\widehat{p}, y, \tau^{p}\right):=\frac{1}{\tau^{p}}\left(\widehat{z}^{0}-\widehat{p}\right), \quad \widehat{z}^{0}:=\mathbb{E}\left(z \mid \widehat{p}, y, \tau^{p}\right) gout ​(p ​,y,τp):=τp1​(z 0−p ​),z 0:=E(z∣p ​,y,τp) 这也是和MAP区别开的地方。 类似地可以得到： − ∂ ∂ p ^ g out  ( p ^ , y , τ p ) = 1 τ p ( 1 − var ⁡ ( z ∣ p ^ , y , τ p ) τ p ) -\frac{\partial}{\partial \widehat{p}} g_{\text {out }}\left(\widehat{p}, y, \tau^{p}\right)=\frac{1}{\tau^{p}}\left(1-\frac{\operatorname{var}\left(z \mid \widehat{p}, y, \tau^{p}\right)}{\tau^{p}}\right) −∂p ​∂​gout ​(p ​,y,τp)=τp1​(1−τpvar(z∣p ​,y,τp)​) 至此， Δ i → j \Delta_{i \rightarrow j} Δi→j​算是可以被近似得到了。 接下来估计 Δ i ← j \Delta_{i \leftarrow j} Δi←j​， 把(6)的结果代入可得到： Δ i ← j ( t + 1 , x j ) ≈ c o n s t + f i n ( x j , q j ) − 1 2 τ i ← j r ( t ) ( r ^ i ← j ( t ) − x j ) 2 \begin{aligned} &\Delta_{i \leftarrow j}\left(t+1, x_{j}\right) \approx \mathrm{const} \\ &\quad+\quad f_{\mathrm{in}}\left(x_{j}, q_{j}\right)-\frac{1}{2 \tau_{i \leftarrow j}^{r}(t)}\left(\widehat{r}_{i \leftarrow j}(t)-x_{j}\right)^{2} \end{aligned} ​Δi←j​(t+1,xj​)≈const+fin​(xj​,qj​)−2τi←jr​(t)1​(r i←j​(t)−xj​)2​ 其中， 1 τ i ← j r ( t ) = ∑ ℓ ≠ i a ℓ j 2 τ ℓ s ( t ) r ^ i ← j ( t ) = τ i ← j r ( t ) ∑ ℓ ≠ i [ s ℓ ( t ) a ℓ j + a ℓ j 2 τ ℓ s ( t ) x ^ j ( t ) ] = x ^ j ( t ) + τ i ← j r ( t ) ∑ ℓ ≠ i s ℓ ( t ) a ℓ j \begin{aligned} \frac{1}{\tau_{i \leftarrow j}^{r}(t)} &=\sum_{\ell \neq i} a_{\ell j}^{2} \tau_{\ell}^{s}(t) \\ \widehat{r}_{i \leftarrow j}(t) &=\tau_{i \leftarrow j}^{r}(t) \sum_{\ell \neq i}\left[s_{\ell}(t) a_{\ell j}+a_{\ell j}^{2} \tau_{\ell}^{s}(t) \widehat{x}_{j}(t)\right] \\ &=\widehat{x}_{j}(t)+\tau_{i \leftarrow j}^{r}(t) \sum_{\ell \neq i} s_{\ell}(t) a_{\ell j} \end{aligned} τi←jr​(t)1​r i←j​(t)​=ℓ​=i∑​aℓj2​τℓs​(t)=τi←jr​(t)ℓ​=i∑​[sℓ​(t)aℓj​+aℓj2​τℓs​(t)x j​(t)]=x j​(t)+τi←jr​(t)ℓ​=i∑​sℓ​(t)aℓj​​ 这步和MAP的步骤是完全一样。接下来，我们注意到： p Δ i ← j ( t , ⋅ ) ( x j ) ∝ exp ⁡ Δ i ← j ( t , x j ) ≈ 1 Z exp ⁡ F i n ( x j , r ^ i ← j ( t ) , q j , τ i ← j r ( t ) ) \begin{aligned} &p_{\Delta_{i \leftarrow j}(t, \cdot)}\left(x_{j}\right)\propto \exp \Delta_{i\leftarrow j}\left(t, x_{j}\right) \\ &\quad \approx \frac{1}{Z} \exp F_{\mathrm{in}}\left(x_{j}, \widehat{r}_{i \leftarrow j}(t), q_{j}, \tau_{i \leftarrow j}^{r}(t)\right) \end{aligned} ​pΔi←j​(t,⋅)​(xj​)∝expΔi←j​(t,xj​)≈Z1​expFin​(xj​,r i←j​(t),qj​,τi←jr​(t))​ 其中， F i n ( x , r ^ , q , τ r ) : = f i n ( x , q ) − 1 2 τ r ( r ^ − x ) 2 F_{\mathrm{in}}\left(x, \widehat{r}, q, \tau^{r}\right):=f_{\mathrm{in}}(x, q)-\frac{1}{2 \tau^{r}}(\widehat{r}-x)^{2} Fin​(x,r ,q,τr):=fin​(x,q)−2τr1​(r −x)2 注意到这又是熟悉的高斯分布变量的指数项， 因此，如果定义： g in  ( r ^ , q , τ r ) : = E [ X ∣ R ^ = r ^ , Q = q ] g_{\text {in }}\left(\widehat{r}, q, \tau^{r}\right):=\mathbb{E}[X \mid \widehat{R}=\widehat{r}, Q=q] gin ​(r ,q,τr):=E[X∣R =r ,Q=q] 这里的期望是对变量 R ^ = X + V , V ∼ N ( 0 , τ r ) \widehat{R}=X+V, \quad V \sim \mathcal{N}\left(0, \tau^{r}\right) R =X+V,V∼N(0,τr) 进行求取。 则我们有： x ^ i ← j ( t + 1 ) = E [ x j ∣ Δ i ← j ( t , ⋅ ) ] ≈ g in  ( r ^ i ← j ( t ) , q j , τ i ← j r ( t ) ) \widehat{x}_{i \leftarrow j}(t+1) =\mathbb{E}\left[x_{j} \mid \Delta_{i \leftarrow j}(t, \cdot)\right]\approx g_{\text {in }}\left(\widehat{r}_{i \leftarrow j}(t), q_{j}, \tau_{i \leftarrow j}^{r}(t)\right) x i←j​(t+1)=E[xj​∣Δi←j​(t,⋅)]≈gin ​(r i←j​(t),qj​,τi←jr​(t)) 至此，剩余步骤的推导都和MAP算法一致。 也就是说，两个版本的GAMP算法的不同仅仅体现在 g o u t g_\mathrm{out} gout​ 和 g i n g_\mathrm{in} gin​ 上。 作者做了一个表格用以对比： 由于MMSE 版本 和 MAP 版本高度相似， 这篇的推导写的较为简略。 没有关系， 我们重点聚焦在对GAMP算法的使用之上。 而从 GAMP 算法的 框图之中： 可以看到， 只有 s ^ i \hat{s}_i s^i​ 即 g o u t g_\mathrm{out} gout​, τ i s \tau_{i}^s τis​, x ^ j \hat{x}_j x^j​即 g i n g_\mathrm{in} gin​ 和 τ j x \tau_j^x τjx​ 是需要推导的， 其他的都是流水线无脑操作即可。 接下来我们以实例来展示。

AWGN 场景

我们考虑AWGN输出场景， 即 y = A x + n y = Ax + n y=Ax+n, 其中 n ∼ C N ( 0 , τ w ) n\sim\mathcal{CN}(0,\tau_w) n∼CN(0,τw​)为高斯白噪声。 此时我们有： f out  ( z , y ) : = log ⁡ p Y ∣ Z ( y ∣ z ) = c o n s t + 1 2 τ w ( z − y ) 2 f_{\text {out }}(z, y):=\log p_{Y \mid Z}(y \mid z)=const + \frac{1}{2\tau_w}(z-y)^2 fout ​(z,y):=logpY∣Z​(y∣z)=const+2τw​1​(z−y)2 那么推导MAP版本的 g o u t g_\mathrm{out} gout​ 为： g o u t = ( z ^ 0 − p ^ ) / τ p g_\mathrm{out}=\left(\widehat{z}^{0}-\widehat{p}\right) / \tau^{p} gout​=(z 0−p ​)/τp 其中 z ^ 0 : = arg ⁡ max ⁡ z F out  ( z , p ^ , y , τ p ) \widehat{z}^{0}:=\arg \max _{z} F_{\text {out }}\left(z, \widehat{p}, y, \tau^{p}\right) z 0:=argzmax​Fout ​(z,p ​,y,τp) 而 F out  ( z , p ^ , y , τ p ) : = f out  ( z , y ) − 1 2 τ p ( z − p ^ ) 2 = − 1 2 τ w ( z − y ) 2 − 1 2 τ p ( z − p ^ ) 2 F_{\text {out }}\left(z, \widehat{p}, y, \tau^{p}\right):=f_{\text {out }}(z, y)-\frac{1}{2 \tau^{p}}(z-\widehat{p})^{2}= -\frac{1}{2\tau_w}(z-y)^2 - \frac{1}{2 \tau^{p}}(z-\widehat{p})^{2} Fout ​(z,p ​,y,τp):=fout ​(z,y)−2τp1​(z−p ​)2=−2τw​1​(z−y)2−2τp1​(z−p ​)2 那么上式对 z z z求导， 得到： ( 1 τ p + 1 τ w ) z = 1 τ w y + 1 τ p p ^ , (\frac{1}{\tau^p} + \frac{1}{\tau_w})z= \frac{1}{\tau_w}y+ \frac{1}{\tau^p}\hat{p}, (τp1​+τw​1​)z=τw​1​y+τp1​p^​, 因此 z ^ 0 = τ p y + τ w p ^ τ p + τ w g o u t = y − p ^ τ p + τ w . \hat{z}^0 = \frac{\tau^py + \tau_w\hat{p}}{\tau^p + \tau_w}\\g_\mathrm{out}=\frac{y-\hat{p}}{\tau^p+\tau_w}. z^0=τp+τw​τpy+τw​p^​​gout​=τp+τw​y−p^​​. 那么很轻易又有： − ∂ ∂ p ^ g out  ( p ^ , y , τ p ) = 1 τ p + τ w -\frac{\partial}{\partial \widehat{p}} g_{\text {out }}\left(\widehat{p}, y, \tau^{p}\right)=\frac{1}{\tau^{p}+\tau^{w}} −∂p ​∂​gout ​(p ​,y,τp)=τp+τw1​

再考虑 MMSE 版本。 其 g o u t g_\mathrm{out} gout​ 函数定义为： g o u t ( p ^ , y , τ p ) : = 1 τ p ( z ^ 0 − p ^ ) , z ^ 0 : = E ( z ∣ p ^ , y , τ p ) g_{\mathrm{out}}\left(\widehat{p}, y, \tau^{p}\right):=\frac{1}{\tau^{p}}\left(\widehat{z}^{0}-\widehat{p}\right), \quad \widehat{z}^{0}:=\mathbb{E}\left(z \mid \widehat{p}, y, \tau^{p}\right) gout​(p ​,y,τp):=τp1​(z 0−p ​),z 0:=E(z∣p ​,y,τp) 其中， 变量 z z z服从分布： p ( z ∣ p ^ , y , τ p ) ∝ exp ⁡ F out  ( z , p ^ , y , τ p ) = − 1 2 τ w ( z − y ) 2 − 1 2 τ p ( z − p ^ ) 2 p\left(z \mid \widehat{p}, y, \tau^{p}\right) \propto \exp F_{\text {out }}\left(z, \widehat{p}, y, \tau^{p}\right)= -\frac{1}{2\tau_w}(z-y)^2 - \frac{1}{2 \tau^{p}}(z-\widehat{p})^{2} p(z∣p ​,y,τp)∝expFout ​(z,p ​,y,τp)=−2τw​1​(z−y)2−2τp1​(z−p ​)2 这里注意到， 事实上式可以看做是两个高斯变量PDF乘积的形式， 而这已有定理， 那就是 上式仍是一个高斯变量的PDF， 即： p ( z ∣ p ^ , y , τ p ) ∼ N ( z ^ 0 , τ z ) p\left(z \mid \widehat{p}, y, \tau^{p}\right) \sim \mathcal{N}\left(\widehat{z}^{0}, \tau^{z}\right) p(z∣p ​,y,τp)∼N(z 0,τz) 且均值和方差分别为： z ^ 0 : = p ^ + τ p τ w + τ p ( y − p ^ ) τ z : = τ w τ p τ w + τ p \begin{aligned} \widehat{z}^{0} &:=\widehat{p}+\frac{\tau^{p}}{\tau^{w}+\tau^{p}}(y-\widehat{p}) \\ \tau^{z} &:=\frac{\tau^{w} \tau^{p}}{\tau^{w}+\tau^{p}} \end{aligned} z 0τz​:=p ​+τw+τpτp​(y−p ​):=τw+τpτwτp​​ 具体的推导细节可以参考这篇博客 https://blog.csdn.net/chaosir1991/article/details/106910668/. 那么代入到 g o u t g_\mathrm{out} gout​ 中， 有： g o u t = y − p ^ τ p + τ w . g_\mathrm{out}=\frac{y-\hat{p}}{\tau^p+\tau_w}. gout​=τp+τw​y−p^​​. 因此， 在输出为AWGN的场景下， MAP 和 MMSE 完全等价。

接下来， 同样考虑 输入 也是 AWGN 场景， 即： p X ∣ Q ( x ∣ q ) = N ( q , τ x 0 ) p_{X \mid Q}(x \mid q)=\mathcal{N}\left(q, \tau^{x 0}\right) pX∣Q​(x∣q)=N(q,τx0) 那么根据定义， MAP 版本为： g in  ( r ^ , q , τ r ) : = arg ⁡ max ⁡ x F in  ( x , r ^ , q , τ r ) g_{\text {in }}\left(\widehat{r}, q, \tau^{r}\right):=\underset{x}{\arg \max } F_{\text {in }}\left(x, \widehat{r}, q, \tau^{r}\right) gin ​(r ,q,τr):=xargmax​Fin ​(x,r ,q,τr) 其中， F in  ( x , r ^ , q , τ r ) : = f in  ( x , q ) − 1 2 τ r ( r ^ − x ) 2 = − 1 2 τ x 0 ( q − x ) 2 − 1 2 τ r ( r ^ − x ) 2 F_{\text {in }}\left(x, \widehat{r}, q, \tau^{r}\right):=f_{\text {in }}(x, q)-\frac{1}{2 \tau^{r}}(\widehat{r}-x)^{2}=-\frac{1}{2 \tau^{x0}}(q-x)^{2}-\frac{1}{2 \tau^{r}}(\widehat{r}-x)^{2} Fin ​(x,r ,q,τr):=fin ​(x,q)−2τr1​(r −x)2=−2τx01​(q−x)2−2τr1​(r −x)2 同样， 对 x x x求导， 可得： g i n ( r ^ , q , τ r ) : = τ x 0 τ x 0 + τ r ( r ^ − q ) + q g_{\mathrm{in}}\left(\widehat{r}, q, \tau^{r}\right):=\frac{\tau^{x 0}}{\tau^{x 0}+\tau^{r}}(\widehat{r}-q)+q gin​(r ,q,τr):=τx0+τrτx0​(r −q)+q 也可轻易求得： τ r g i n ′ ( r ^ , q , τ r ) : = τ r 1 − τ r f i n ′ ′ ( x ^ , q ) = τ x 0 τ r τ x 0 + τ r . \tau^{r} g_{\mathrm{in}}^{\prime}\left(\widehat{r}, q, \tau^{r}\right) \quad:=\frac{\tau^{r}}{1-\tau^{r} f_{\mathrm{in}}^{\prime \prime}(\widehat{x}, q)}=\frac{\tau^{x 0} \tau^{r}}{\tau^{x 0}+\tau^{r}}. τrgin′​(r ,q,τr):=1−τrfin′′​(x ,q)τr​=τx0+τrτx0τr​. 而 MMSE 版本中则有： g in  ( r ^ , q , τ r ) : = E [ X ∣ R ^ = r ^ , Q = q ] g_{\text {in }}\left(\widehat{r}, q, \tau^{r}\right):=\mathbb{E}[X \mid \widehat{R}=\widehat{r}, Q=q] gin ​(r ,q,τr):=E[X∣R =r ,Q=q] 其中 R ^ = X + V , V ∼ N ( 0 , τ r ) \widehat{R}=X+V, \quad V \sim \mathcal{N}\left(0, \tau^{r}\right) R =X+V,V∼N(0,τr) 那么 根据贝叶斯公式我们有： p ( x ∣ r , q ) = p ( r ∣ x ) p ( x ∣ q ) p ( r ) ∝ p ( r ∣ x ) p ( x ∣ q ) p(x|r,q)=\frac{p(r|x)p(x|q)}{p(r)}\propto p(r|x)p(x|q) p(x∣r,q)=p(r)p(r∣x)p(x∣q)​∝p(r∣x)p(x∣q) 而右边这个，又正是刚刚说到的： 两个高斯分布PDF之积仍是高斯分布PDF， 其均值方差和上面用到的一样。 具体， 其均值为： g in  ( r ^ , q , τ r ) : = τ x 0 τ x 0 + τ r ( r ^ − q ) + q . g_{\text {in }}\left(\widehat{r}, q, \tau^{r}\right) \quad:=\frac{\tau^{x 0}}{\tau^{x 0}+\tau^{r}}(\widehat{r}-q)+q. gin ​(r ,q,τr):=τx0+τrτx0​(r −q)+q. 这再次和 MAP 版本一致。 那么也有： τ r g i n ′ ( r ^ , q , τ r ) : = τ x 0 τ r τ x 0 + τ r . \tau^{r} g_{\mathrm{in}}^{\prime}\left(\widehat{r}, q, \tau^{r}\right):=\frac{\tau^{x 0} \tau^{r}}{\tau^{x 0}+\tau^{r}}. τrgin′​(r ,q,τr):=τx0+τrτx0τr​.