AMP State Evolution (SE)的计算

$t=1$时，${\mathcal{E}}^{(t)}=\mathbb{E}[X^2]$，SE的迭代式为
$$\begin{array}{rl}{\tau }_r^{(t)}& ={\sigma }^2+\frac{1}{\delta }{\mathcal{E}}^{(t)}\\ {\mathcal{E}}^{(t+1)}& =\mathbb{E}{\left|{\eta }^{(t)}\left(X+Z\right)-X\right|}^2,Z\sim \mathcal{N}(0,{\tau }_r^{(t)})\end{array}$$

撰写的时候存在一定的符号乱用，之后的$\tau$即指${\tau }_r^{(t)}$。

注意到，${\mathcal{E}}^{(t+1)}=\mathbb{E}{\left|{\eta }^{(t)}\left(X+Z\right)-X\right|}^2$是关于随机变量$X,Z$求期望，这里我们认为$X,Z$之间相互独立，因此
$${\mathcal{E}}^{(t+1)}=\int \int p_{X,Z}(X,Z){\left|{\eta }^{(t)}\left(X+Z\right)-X\right|}^{2}dXdZ$$

令$R=X+Z$，不难得到$p_{X,R}(X,R)$等价于$p_{X,Z}(X,Z)$（因为$p_{X,Z}(X=x_0,Z=z_0)=p_{X,R}(X=x_0,R=x_0+z_0)$）,因此我们可以把${\mathcal{E}}^{(t+1)}$写为（这里我们考虑${\eta }^{(t)}$为MMSE函数）
$$\begin{array}{rl}{\mathcal{E}}^{(t+1)}& =\int \int p_{X,R}(X,R){\left|{\eta }^{(t)}\left(R\right)-X\right|}^{2}dXdR\\ & =\int \int p_R(R)p_{X|R}(X|R){\left|\mathbb{E}\left[X|R\right]-X\right|}^{2}dXdR\\ & =\int p_R(R)\mathrm{var}[X|R]dR\end{array}$$

因为$p_{X,Z}(X,Z)=p_X(X)\cdot p_X(Z)=p_X(X)\mathcal{N}(z;0,{\tau }_r^{(t)})$，因此可以得到
$$\begin{array}{rl}{p}_{X,R}(X,R)& =p_X(X)\mathcal{N}(R-X;0,{\tau }_r^{(t)})\\ & =p_X(X)\mathcal{N}(R;X,{\tau }_r^{(t)})\end{array}$$

X的先验为为伯努利高斯分布时

$$p_X(X)\equiv (1-\rho )\delta (X)+\rho \mathcal{N}(\mu ,\nu )$$

那么，$X,R$的联合分布可写为
$$\begin{array}{rl}{p}_{X,R}(X,R)& =\left((1-\rho )\delta (X)+\rho \mathcal{N}(X;\mu ,\nu )\right)\cdot \mathcal{N}(X;R,{\tau }_r^{(t)})\\ & =(1-\rho )\delta (X)\mathcal{N}(X;R,{\tau }_r^{(t)})+\rho \mathcal{N}(X;\mu ,\nu )\mathcal{N}(X;R,{\tau }_r^{(t)})\end{array}$$

进一步，我们可以得到关于$R$的边缘分布
$$\begin{array}{rl}{p}_R(R)& =\int p_{X,R}(X,R)dX\\ & =(1-\rho )\mathcal{N}(0;R,{\tau }_r^{(t)})+\rho \int \mathcal{N}(X;\mu ,\nu )\mathcal{N}(X;R,{\tau }_r^{(t)})dX\\ & =(1-\rho )\mathcal{N}(0;R,{\tau }_r^{(t)})+\rho \frac{1}{\sqrt{2\pi (\nu +{\tau }_r^{(t)})}}\exp \left\{-\frac{(R-\mu {)}^2}{2(\nu +{\tau }_r^{(t)})}\right\}\\ & =(1-\rho )\mathcal{N}(R;0,{\tau }_r^{(t)})+\rho \mathcal{N}(R;\mu ,\nu +{\tau }_r^{(t)})\end{array}$$

我们计算后验均值$\mathbb{E}[X|R]$
$$\begin{array}{rl}\mathbb{E}[X|R]& =\int X{p}_{X|R}(X|R)dX\\ & =\frac{1}{p_R(R)}\int X{p}_{X,R}(X,R)dX\\ & =\frac{1}{p_R(R)}\cdot \rho \cdot \int X\mathcal{N}(X;\mu ,\nu )\mathcal{N}(X;R,{\tau }_r^{(t)})dX\\ & =\frac{1}{p_R(R)}\cdot \rho \cdot \mathcal{N}(R;\mu ,\nu +{\tau }_r^{(t)})\cdot \frac{\mu {\tau }_r^{(t)}+R\nu }{\nu +{\tau }_r^{(t)}}\\ & =\frac{\rho \cdot \mathcal{N}(R;\mu ,\nu +{\tau }_r^{(t)})}{(1-\rho )\mathcal{N}(R;0,{\tau }_r^{(t)})+\rho \mathcal{N}(R;\mu ,\nu +{\tau }_r^{(t)})}\cdot \frac{\mu {\tau }_r^{(t)}+R\nu }{\nu +{\tau }_r^{(t)}}\\ & =\frac{1}{1+\frac{1-\rho }{\rho }\cdot \frac{\mathcal{N}(R;0,{\tau }_r^{(t)})}{\mathcal{N}(R;\mu ,\nu +{\tau }_r^{(t)})}}\cdot \frac{\mu {\tau }_r^{(t)}+R\nu }{\nu +{\tau }_r^{(t)}}\end{array}$$

对于大多数问题，我们有$\mu =0$，令$\varphi (R)=\frac{1-\rho }{\rho }\cdot \frac{\mathcal{N}(R;0,{\tau }_r^{(t)})}{\mathcal{N}(R;\mu ,\nu +{\tau }_r^{(t)})}$, 当$\mu =0$时，我们有
$$\begin{array}{r}\mathbb{E}[X|R]=\frac{1}{1+\varphi (R)}\cdot \frac{R\nu }{\nu +{\tau }_r^{(t)}}\end{array}$$

进一步，可以得到
$$\frac{\partial }{\partial R}\mathbb{E}[X|R]=\frac{\nu }{\nu +{\tau }_r^{(t)}}\cdot \frac{1+\varphi (R)+R^2(\frac{1}{\tau }-\frac{1}{\tau +\nu })\varphi (R)}{(1+\varphi (R){)}^2}$$

根据AWGN后验均值与方差的关系，
$$\mathrm{var}[X|R]={\tau }_r^{(t)}\cdot \frac{\partial }{\partial R}\mathbb{E}[X|R]=\frac{\nu {\tau }_r^{(t)}}{\nu +{\tau }_r^{(t)}}\cdot \frac{1+\varphi (R)+R^2(\frac{1}{\tau }-\frac{1}{\tau +\nu })\varphi (R)}{(1+\varphi (R){)}^2}$$

不难看出，$p_R(R)$可以表示为$p_R(R)=(1+\varphi (R))\cdot \rho \mathcal{N}(R;0,\nu +{\tau }_r^{(t)})$，因此
$$\begin{array}{rl}{\mathcal{E}}^{(t+1)}& =\int p_R(R)\mathrm{var}[X|R]dR\\ & =\int (1+\varphi (R))\cdot \rho \mathcal{N}(R;0,\nu +{\tau }_r^{(t)})\cdot \frac{\nu {\tau }_r^{(t)}}{\nu +{\tau }_r^{(t)}}\cdot \frac{1+\varphi (R)+R^2(\frac{1}{\tau }-\frac{1}{\tau +\nu })\varphi (R)}{(1+\varphi (R){)}^2}\\ & =\rho \frac{\nu {\tau }_r^{(t)}}{\nu +{\tau }_r^{(t)}}\cdot \int \mathcal{N}(R;0,\nu +{\tau }_r^{(t)})\left\{1+R^2\frac{\nu }{\tau (\tau +\nu )}\left(1-\frac{1}{1+\varphi (R)}\right)\right\}dR\\ & =\rho \frac{\nu {\tau }_r^{(t)}}{\nu +{\tau }_r^{(t)}}\cdot \left\{1+\frac{\nu }{\tau }-\frac{\nu }{\tau (\tau +\nu )}\int \frac{R^2}{1+\varphi (R)}\mathcal{N}(R;0,\nu +{\tau }_r^{(t)})dR\right\}\\ & =\rho \cdot \nu -\rho \cdot \frac{{\nu }^2}{(\nu +\tau {)}^2}\int \frac{R^2}{1+\varphi (R)}\mathcal{N}(R;0,\nu +{\tau }_r^{(t)})dR\end{array}$$

其中
$$\frac{R^2}{1+\varphi (R)}=\frac{R^2}{1+\frac{1-\rho }{\rho }\cdot \sqrt{\frac{\nu +\tau }{\tau }}\exp \left(-\frac{\nu }{2\tau (\tau +\nu )}{R}^2\right)}$$

总结

$t=1$时，${\mathcal{E}}^{(t)}=\mathbb{E}[X^2]$，SE的迭代式为
$$\begin{array}{rl}{\tau }_r^{(t)}& ={\sigma }^2+\frac{1}{\delta }{\mathcal{E}}^{(t)}\\ {\mathcal{E}}^{(t+1)}& =\mathbb{E}{\left|{\eta }^{(t)}\left(X+Z\right)-X\right|}^2,Z\sim \mathcal{N}(0,{\tau }_r^{(t)})\end{array}$$

当$X$的先验是伯努利高斯时：$X\sim (1-\rho )\delta (X)+\rho \mathcal{N}(0,\nu )$，可以把SE表征为
$t=1$时，${\mathcal{E}}^{(t)}=\mathbb{E}[X^2]=\rho \nu$，SE的迭代式为
$$\begin{array}{rl}{\tau }_r^{(t)}& ={\sigma }^2+\frac{1}{\delta }{\mathcal{E}}^{(t)}\\ {\mathcal{E}}^{(t+1)}& =\rho \cdot \nu -\rho \cdot \frac{{\nu }^2}{(\nu +{\tau }_r^{(t)}{)}^2}\int \frac{R^2}{1+\varphi (R)}\mathcal{N}(R;0,\nu +{\tau }_r^{(t)})dR\end{array}$$

其中
$$\frac{R^2}{1+\varphi (R)}=\frac{R^2}{1+\frac{1-\rho }{\rho }\cdot \sqrt{\frac{\nu +{\tau }_r^{(t)}}{{\tau }_r^{(t)}}}\exp \left(-\frac{\nu }{2{\tau }_r^{(t)}({\tau }_r^{(t)}+\nu )}{R}^2\right)}$$

SE部分的MATLAB代码

sigma2 = 0.0165;
rho = 0.1;
nu = 1 / rho;
lim = inf;
Iteration = 50;
delta = 0.6;   % underdetermined ratio
SE_MSE = zeros(Iteration, 1);
SE_tau2 = zeros(Iteration, 1);

SE_MSE(1) = rho * nu;
SE_tau2(1) = sigma2 + 1/delta * SE_MSE(1);
for it = 2: Iteration
    tau = SE_tau2( it - 1 );
    f = @(r) r.*r ./ ...
        (    1   +   (1-rho)./rho .* sqrt( 1 + nu./tau )   .*   exp( max(-60, -0.5 * r .* r .* nu./ tau./(nu + tau) ) )     ) ...
        .* normpdf(r, 0, sqrt(nu + tau))  ;  % -60 is a bound for numerical robustness

    SE_MSE(it) = rho * nu - rho * nu^2 / (nu + tau)^2 * integral(f,-lim,lim);
    SE_tau2(it) = sigma2 + 1/delta * SE_MSE(it);
end

当N=20000时，AMP的MSE性能与SE一致

附录：两个高斯随机变量相乘

$$\begin{array}{rl}& \mathcal{N}(x;{\mu }_1,{\nu }_1)\mathcal{N}(x;{\mu }_2,{\nu }_2)\\ & =\frac{1}{\sqrt{2\pi {\nu }_1}}\exp \left(-\frac{(x-{\mu }_1{)}^2}{2{\nu }_1}\right)\cdot \frac{1}{\sqrt{2\pi {\nu }_2}}\exp \left(-\frac{(x-{\mu }_2{)}^2}{2{\nu }_2}\right)\\ & =\frac{1}{2\pi \sqrt{{\nu }_1{\nu }_2}}\exp \left\{-\frac{{\nu }_2(x^2-2{\mu }_{1}x+{\mu }_1^2)}{2{\nu }_1{\nu }_2}-\frac{{\nu }_1(x^2-2{\mu }_{2}x+{\mu }_2^2)}{2{\nu }_1{\nu }_2}\right\}\\ & =\frac{1}{2\pi \sqrt{{\nu }_1{\nu }_2}}\exp \left\{-\frac{({\nu }_1+{\nu }_2)x^2-2({\mu }_1{\nu }_2+{\mu }_2{\nu }_1)x+{\mu }_1^2{\nu }_2+{\mu }_2^2{\nu }_1}{2{\nu }_1{\nu }_2}\right\}\\ & =\frac{1}{2\pi \sqrt{{\nu }_1{\nu }_2}}\exp \left\{-\frac{x^2-2\frac{({\mu }_1{\nu }_2+{\mu }_2{\nu }_1)}{{\nu }_1+{\nu }_2}x+{\left(\frac{({\mu }_1{\nu }_2+{\mu }_2{\nu }_1)}{{\nu }_1+{\nu }_2}\right)}^2+\frac{{\mu }_1^2{\nu }_2+{\mu }_2^2{\nu }_1}{{\nu }_1+{\nu }_2}-{\left(\frac{({\mu }_1{\nu }_2+{\mu }_2{\nu }_1)}{{\nu }_1+{\nu }_2}\right)}^2}{2\frac{{\nu }_1{\nu }_2}{{\nu }_1+{\nu }_2}}\right\}\\ & =\frac{1}{\sqrt{2\pi ({\nu }_1+{\nu }_2)}}\exp \left\{-\frac{{\mu }_1^2{\nu }_2+{\mu }_2^2{\nu }_1-\frac{({\mu }_1{\nu }_2+{\mu }_2{\nu }_1{)}^2}{{\nu }_1+{\nu }_2}}{2{\nu }_1{\nu }_2}\right\}\cdot \frac{1}{\sqrt{2\pi \frac{{\nu }_1{\nu }_2}{{\nu }_1+{\nu }_2}}}\exp \left\{-\frac{{\left(x-\frac{{\mu }_1{\nu }_2+{\mu }_2{\nu }_1}{{\nu }_1+{\nu }_2}\right)}^2}{2\frac{{\nu }_1{\nu }_2}{{\nu }_1+{\nu }_2}}\right\}\\ & =\frac{1}{\sqrt{2\pi ({\nu }_1+{\nu }_2)}}\exp \left\{-\frac{({\mu }_1-{\mu }_2{)}^2}{2({\nu }_1+{\nu }_2)}\right\}\cdot \frac{1}{\sqrt{2\pi \frac{{\nu }_1{\nu }_2}{{\nu }_1+{\nu }_2}}}\exp \left\{-\frac{{\left(x-\frac{{\mu }_1{\nu }_2+{\mu }_2{\nu }_1}{{\nu }_1+{\nu }_2}\right)}^2}{2\frac{{\nu }_1{\nu }_2}{{\nu }_1+{\nu }_2}}\right\}\end{array}$$

求导时，我们借助了在线求导网站: https://mathdf.com/der/cn/