前言

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

模型背景

我们旨在解决上图这样的问题，已知输入q\mathbf{q}q （先验信息），已知输出 y\mathbf{y}y（后验信息)，已知变换矩阵A\mathbf{A}A，反推出变量x\mathbf{x}x。以AWGN信道举例：
y=z+w=Ax+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∈RnF(x,z,q,y),z^=Ax^\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}} xmap :=x∈RnargmaxF(x,z,q,y),z=Ax
即最大化后验概率。我举一个例子，比如经过计算得到， x=0.5x=0.5x=0.5 的概率是0.20.20.2，而x=0.1x=0.1x=0.1的概率是0.190.190.19， x=0.11x=0.11x=0.11的概率是0.180.180.18， x=0.09x=0.09x=0.09的概率是0.170.170.17。此时，如果使用最大后验 MAP 估计，得到的xxx估计结果就是x=0.5x=0.5x=0.5。但在这种情况下这就有失偏颇，因为他完全没有考虑其他情况的概率。这就有点像通信中的硬判决，一刀切，因此，综合了所有可能情况的软判决，可能更显客观，这也是我们要介绍的MMSE估计，即：
x^mmse:=E[x∣y,q]\widehat{\mathrm{x}}^{\mathrm{mmse}}:=\mathbb{E}[\mathrm{x} \mid \mathrm{y}, \mathbf{q}] xmmse:=E[x∣y,q]
以期望作为估计值。下面，我们就推导以MMSE为估计准则的GAMP版本。

回顾下MAP中的消息传递：

在每次迭代中，实际要计算两个消息算子：
Δi→j(t,xj)=const+max⁡x\xjfout (zi,yi)+∑r≠jΔi←r(t,xr)(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\xjmaxfout (zi,yi)+r=j∑Δi←r(t,xr)(1)
和
Δi←j(t+1,xj)=const+fin(xj,qj)+∑ℓ≠iΔℓ→j(t,xj)(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⁡\maxmax后的结果。然而在MMSE准则下，我们要传递的消息就变成了：
Δi→j(t,xj)=log⁡E(pY∣Z(yi,zi)∣xj)+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(xj)=const +log⁡pX∣Q(xj∣qj)+∑l≠iΔl→j(xj)(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)
注意到，在这个消息传递中，传递的都是似然信息了，也就是说有：
pi←r(xr)∝exp⁡Δi←r(t,xr)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)中的期望是对随机变量 zi=aiTxz_{i}=\mathbf{a}_{i}^{T} \mathbf{x}zi=aiTx 进行求取，此时 xjx_jxj 是固定的，而x\mathbf{x}x 中的其余项，则服从独立分布 pi←r(xr)∝exp⁡Δi←r(t,xr)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估计一样，我们最后要得到的估计值也是由下式得出：
pj(xj)∝exp⁡Δj(t,xj)p_{j}\left(x_{j}\right) \propto \exp \Delta_{j}\left(t, x_{j}\right) pj(xj)∝expΔj(t,xj)
其中
Δj(t+1,xj)=fin(xj,qj)+∑iΔi→j(t,xj)\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[xj∣Δj(t,⋅)]x^i←j(t):=E[xj∣Δi←j(t,⋅)]τjx(t):=var⁡[xj∣Δj(t,⋅)]τi←jx(t):=var⁡[xj∣Δ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} xj(t)xi←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}_jx^j 就是我们最后要得到的估计量。

我们先对 (3) 进行近似，这可以写为：
Δi→j(t,xj)=const +log⁡∫{xr}r≠jpY∣Z(yi∣aijxj+∑r≠jairxr⏟≜zi)∏r≠jeΔi←r(t,xr).\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=jpY∣Z(yi∣≜ziaijxj+r=j∑airxr)r=j∏eΔi←r(t,xr).
而当矩阵维度较大时，那么，根据中心极限定理，相互独立的随机变量，其和服从高斯分布，且均值就是各自均值之和，而方差则是各自方差之和，因此，有：
zi∣xj∼N(aijxj+p^i←j(t),τi←jp(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(aijxj+pi←j(t),τi←jp(t))
其中，
p^i→j(t):=∑r≠jairx^i←r(t)τi→jp(t):=∑r≠j∣air∣2τrx(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} pi→j(t)τi→jp(t):=r=j∑airxi←r(t):=r=j∑∣air∣2τrx(t)
那么，我们进一步就有：
Δi→j(t,xj)≈const +log⁡∫zipY∣Z(yi∣zi)N(zi;aijxj+p^i←j(t),τi←jp(t))⏟≜H(aijxj+p^i←j(t),yi,τi←jp(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(aijxj+pi←j(t),yi,τi←jp(t))log∫zipY∣Z(yi∣zi)N(zi;aijxj+pi←j(t),τi←jp(t))(5)
这里，我们有：
H(p^,y,μp)≜log⁡∫pY∣Z(y∣z)N(z;p^,μp)dzH\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←ji\leftarrow ji←j 项替换掉，定义变量：
p^i(t):=∑jaijx^i←j(t)=p^i→j+aijx^i←j(t)，τip(t)=∑j∣aij∣2τjx(t)+aij2τjx(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) pi(t):=j∑aijxi←j(t)=pi→j+aijxi←j(t)，τip(t)=j∑∣aij∣2τjx(t)+aij2τjx(t)
那么，(5)可以被写为
Δi→j(t,xj)≈H(p^i(t)+aij(xj−x^j),yi,τip(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(pi(t)+aij(xj−xj),yi,τip(t))+ const.
可以通过泰勒展开得到：
Δi→j(t,xj)≈const +si(t)aij(xj−x^j(t))−τis(t)2aij2(xj−x^j(t))2=const⁡[si(t)aij+aij2τis(t)x^j(t)]xj−τis(t)2aij2xj2(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−xj(t))−2τis(t)aij2(xj−xj(t))2const[si(t)aij+aij2τis(t)xj(t)]xj−2τis(t)aij2xj2(6)
这两个表达式和之前一样，这里的近似是我们忽略了 O(aij2)O\left(a_{i j}^{2}\right)O(aij2)级的项。其中，
s^i(t)=gout (t,p^i(t),yi,τip(t))τis(t)=−∂∂p^gout (t,p^i(t),yi,τip(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} si(t)τis(t)=gout (t,pi(t),yi,τip(t))=−∂p∂gout (t,pi(t),yi,τip(t))
这些都是和MAP版本的定义完全一致，然而由于HHH函数的不同， goutg_\mathrm{out}gout的形式也截然不同。因此，接下来我们要对 goutg_\mathrm{out}gout进行推导，也即推导 HHH的一阶导：
H′(p^,y,μp)=∂∂p^log⁡∫pY∣Z(y∣z)12πμpexp⁡(−12μp(z−p^)2)dz=∂∂p^{log⁡12πμp+log⁡∫zexp⁡(log⁡pY∣Z(y∣z)−12μp(z−p^)2)dz}=∂∂p^{−p^22μp+log⁡∫exp⁡(log⁡pY∣Z(y∣z)−z22μp+p^zμp)dz}=−p^μp+∂∂p^log⁡[μp∫exp⁡(ϕ(u)+p^u)du]via u≜zμp=−p^μp+∂∂p^log⁡∫exp⁡(ϕ(u)+p^u)du\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πμp1exp(−2μp1(z−p)2)dz=∂p∂{log2πμp1+log∫zexp(logpY∣Z(y∣z)−2μp1(z−p)2)dz}=∂p∂{−2μpp2+log∫exp(logpY∣Z(y∣z)−2μpz2+μppz)dz}=−μpp+∂p∂log[μp∫exp(ϕ(u)+pu)du] via u≜μpz=−μpp+∂p∂log∫exp(ϕ(u)+pu)du
记 Z(p^)≜∫exp⁡(ϕ(u)+p^u)duZ(\widehat{p}) \triangleq \int \exp (\phi(u)+\widehat{p} u) d uZ(p)≜∫exp(ϕ(u)+pu)du，我们有如下数学公理：
∂∂p^log⁡Z(p^)=E{u∣p^}with pU∣P(u∣p^)=exp⁡(ϕ(u)+p^u)Z(p^)∂2∂p^2log⁡Z(p^)=var⁡{u∣p^}with pU∣P(u∣p^)=exp⁡(ϕ(u)+p^uZ(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)∂p2∂2logZ(p)=var{u∣p} with pU∣P(u∣p)=Z(p^)exp(ϕ(u)+p^u.
可以通过简单的求导法则验证，这是正确的。因此，
H′(p^,y,μp)=−p^μp+∫uexp⁡(ϕ(u)+p^u)Z(p^)du=−p^μp+∫zμpexp⁡(log⁡pY∣Z(y∣z)−z22μp+zp^μp)Z(p^)dzμpvia u≜zμp=−p^μp+1μp∫zexp⁡(log⁡pY∣Z(y∣z)−12μp(z−p^)2)μpZ(p^)exp⁡(−p22μp)dz=−p^μp+1μp∫zpY∣Z(y∣z)N(z;p^,μp)∫pY∣Z(y∣zˉ)N(zˉ;p^,μp)dzˉdz=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)+pu)du=−μpp+∫μpzZ(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)
因此，我们有：
gout (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(z0−p),z0:=E(z∣p,y,τp)
这也是和MAP区别开的地方。类似地可以得到：
−∂∂p^gout (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,xj)≈const+fin(xj,qj)−12τi←jr(t)(r^i←j(t)−xj)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(ri←j(t)−xj)2
其中，
1τi←jr(t)=∑ℓ≠iaℓj2τℓs(t)r^i←j(t)=τi←jr(t)∑ℓ≠i[sℓ(t)aℓj+aℓj2τℓs(t)x^j(t)]=x^j(t)+τi←jr(t)∑ℓ≠isℓ(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)1ri←j(t)=ℓ=i∑aℓj2τℓs(t)=τi←jr(t)ℓ=i∑[sℓ(t)aℓj+aℓj2τℓs(t)xj(t)]=xj(t)+τi←jr(t)ℓ=i∑sℓ(t)aℓj
这步和MAP的步骤是完全一样。接下来，我们注意到：
pΔi←j(t,⋅)(xj)∝exp⁡Δi←j(t,xj)≈1Zexp⁡Fin(xj,r^i←j(t),qj,τi←jr(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)≈Z1expFin(xj,ri←j(t),qj,τi←jr(t))
其中，
Fin(x,r^,q,τr):=fin(x,q)−12τr(r^−x)2F_{\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
注意到这又是熟悉的高斯分布变量的指数项，因此，如果定义：
gin (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[xj∣Δi←j(t,⋅)]≈gin (r^i←j(t),qj,τi←jr(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) xi←j(t+1)=E[xj∣Δi←j(t,⋅)]≈gin (ri←j(t),qj,τi←jr(t))
至此，剩余步骤的推导都和MAP算法一致。也就是说，两个版本的GAMP算法的不同仅仅体现在goutg_\mathrm{out}gout 和 ging_\mathrm{in}gin 上。作者做了一个表格用以对比：

由于MMSE 版本和 MAP 版本高度相似，这篇的推导写的较为简略。没有关系，我们重点聚焦在对GAMP算法的使用之上。而从 GAMP 算法的框图之中：

可以看到，只有 s^i\hat{s}_is^i 即 goutg_\mathrm{out}gout, τis\tau_{i}^sτis, x^j\hat{x}_jx^j即 ging_\mathrm{in}gin 和 τjx\tau_j^xτjx 是需要推导的，其他的都是流水线无脑操作即可。接下来我们以实例来展示。

AWGN 场景

我们考虑AWGN输出场景，即 y=Ax+ny = Ax + ny=Ax+n, 其中 n∼CN(0,τw)n\sim\mathcal{CN}(0,\tau_w)n∼CN(0,τw)为高斯白噪声。此时我们有：
fout (z,y):=log⁡pY∣Z(y∣z)=const+12τw(z−y)2f_{\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τw1(z−y)2
那么推导MAP版本的 goutg_\mathrm{out}gout 为：
gout=(z^0−p^)/τpg_\mathrm{out}=\left(\widehat{z}^{0}-\widehat{p}\right) / \tau^{p} gout=(z0−p)/τp
其中
z^0:=arg⁡max⁡zFout (z,p^,y,τp)\widehat{z}^{0}:=\arg \max _{z} F_{\text {out }}\left(z, \widehat{p}, y, \tau^{p}\right) z0:=argzmaxFout (z,p,y,τp)
而
Fout (z,p^,y,τp):=fout (z,y)−12τp(z−p^)2=−12τw(z−y)2−12τp(z−p^)2F_{\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τw1(z−y)2−2τp1(z−p)2
那么上式对zzz求导，得到：
(1τp+1τw)z=1τwy+1τpp^,(\frac{1}{\tau^p} + \frac{1}{\tau_w})z= \frac{1}{\tau_w}y+ \frac{1}{\tau^p}\hat{p},(τp1+τw1)z=τw1y+τp1p^,
因此
z^0=τpy+τwp^τp+τwgout=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+τwp^gout=τp+τwy−p^.
那么很轻易又有：
−∂∂p^gout (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 版本。其goutg_\mathrm{out}gout 函数定义为：
gout(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(z0−p),z0:=E(z∣p,y,τp)
其中，变量zzz服从分布：
p(z∣p^,y,τp)∝exp⁡Fout (z,p^,y,τp)=−12τw(z−y)2−12τp(z−p^)2p\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τw1(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(z0,τ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} z0τz:=p+τw+τpτp(y−p):=τw+τpτwτp
具体的推导细节可以参考这篇博客 https://blog.csdn.net/chaosir1991/article/details/106910668/.
那么代入到 goutg_\mathrm{out}gout 中，有：
gout=y−p^τp+τw.g_\mathrm{out}=\frac{y-\hat{p}}{\tau^p+\tau_w}.gout=τp+τwy−p^.
因此，在输出为AWGN的场景下， MAP 和 MMSE 完全等价。

接下来，同样考虑输入也是 AWGN 场景，即：
pX∣Q(x∣q)=N(q,τx0)p_{X \mid Q}(x \mid q)=\mathcal{N}\left(q, \tau^{x 0}\right) pX∣Q(x∣q)=N(q,τx0)
那么根据定义， MAP 版本为：
gin (r^,q,τr):=arg⁡max⁡xFin (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):=xargmaxFin (x,r,q,τr)
其中，
Fin (x,r^,q,τr):=fin (x,q)−12τr(r^−x)2=−12τx0(q−x)2−12τr(r^−x)2F_{\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
同样，对xxx求导，可得：
gin(r^,q,τr):=τx0τx0+τr(r^−q)+qg_{\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
也可轻易求得：
τrgin′(r^,q,τr):=τr1−τrfin′′(x^,q)=τx0τrτx0+τ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 版本中则有：
gin (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，其均值方差和上面用到的一样。具体，其均值为：
gin (r^,q,τr):=τx0τx0+τ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 版本一致。那么也有：
τrgin′(r^,q,τr):=τx0τrτx0+τ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.

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

深入浅出K-Means算法
深入浅出K-Means算法摘要:在数据挖掘中,K-Means算法是一种 cluster analysis 的算法,其主要是来计算数据聚集的算法,主要通过不断地取离种子点最近均值的算法. 在数据挖掘中 ...
【测试算法】深入浅出Pairwise 算法
深入浅出Pairwise 算法作者:王勇软件测试是软件开发中很重要的一环,在软件成本中也占着很大的比重.本文在介绍pairwise算法的基础上,提出了针对某一类问题的扩展算法并加以实现. 本文的组 ...
深入浅出Pairwise 算法
深入浅出Pairwise 算法作者:王勇软件测试是软件开发中很重要的一环,在软件成本中也占着很大的比重.本文在介绍pairwise算法的基础上,提出了针对某一类问题的扩展算法并加以实现. 本文的组 ...
## ***电池SOC仿真系列-基于扩展卡尔曼(EKF)算法的SOC估计（内含代码等资料）***
## ***电池SOC仿真系列-基于扩展卡尔曼(EKF)算法的SOC估计(内含代码等资料)*** ## 1 研究背景电池的荷电状态(SOC)代表的是电池当前的剩余容量,数值定义是电池剩余电量与电池额 ...
1325: 深入浅出学算法020-阶乘和(sum)
1325: 深入浅出学算法020-阶乘和(sum) 欢迎使用Markdown编辑器 #include<bits/stdc++.h> using namespace std; int b[3 ...
深入浅出 pairwise 算法
深入浅出Pairwise 算法作者:王勇软件测试是软件开发中很重要的一环,在软件成本中也占着很大的比重.本文在介绍pairwise算法的基础上,提出了针对某一类问题的扩展算法并加以实现. 本文的组 ...
有限算法下的技术实现路线
有限算法下的技术实现路线从信息技术本身来说,它最大的优势是可以快速处理复杂的运算.我们企业做信息化要善于利用这样的优势,而不是因为它有这样的优势就考虑怎么用它,那样非常容易作茧自缚的.在不考虑业务逻 ...
计算机操作系统实验银行家算法,实验六银行家算法(下)
实验六银行家算法(下) 一.实验说明实验说明:本次实验主要是对银行家算法进行进一步的实践学习,掌握银行家算法的整体流程,理解程序测试时每一步的当前状态,能对当前的资源分配进行预判断. 二.实验要求 ...
Problem E: 深入浅出学算法019-求n的阶乘
Problem E: 深入浅出学算法019-求n的阶乘 Time Limit: 1 Sec Memory Limit: 64 MB Submit: 5077 Solved: 3148 Descri ...
深入浅出学算法007-统计求和
4006: 深入浅出学算法007-统计求和 Time Limit: 1 Sec Memory Limit: 64 MB Submit: 4335 Solved: 2014 Description 求含 ...

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

前言

模型背景

MMSE

AWGN 场景

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

最新文章

热门文章