OAMP的理解
Orthogonal Approximate Message Passing, OAMP
1 前言
本次博文主要介绍了OAMP论文，同时加了一些粗浅的理解。一方面，前两部分涉及了一些AMP相关的知识以及我自己给出的解释，另一方面，博文所推公式可能与论文稍有偏差，笔者才疏学浅，又是第一次写博客，这两个部分都不敢保证能理解到位和描述清晰，难免解释可能会有所偏差，甚至有歪曲原文的不当之处，还希望各位读者能够包涵。下一篇博文可能会是对VAMP的解读，希望这可以在年内完成，之后博客可能会推出一些与消息传递算法较相关的论文介绍，以及自认为有价值的积累。
2 绪论
简短回顾一下Approximate Message Passing (AMP) 的问题模型
y
=
A
x
+
n
(1a)
pmb y=pmb A pmb x + pmb n tag{1a}
yyy=AAAxxx+nnn(1a)
x
j
∼
P
X
(
x
)
,
∀
x
(1b)
mathit x_{j} sim P_{X}(x), forall x tag{1b}
xj∼PX(x),∀x(1b)
其中
A
∈
C
M
×
N
pmb A in mathbb C^{M times N}
AAA∈CM×N是感知矩阵，
n
∈
C
M
×
1
pmb n in mathbb C^{M times 1}
nnn∈CM×1是均值为
0
0
0，方差为
σ
2
sigma^{2}
σ2的高斯向量。AMP的一个重要性质就是其算法性能可以用state evolution来精确刻画，就是说给定真实的
x
pmb x
xxx，我们可以依据state evolution的结果
τ
t
tau^{t}
τt，预测AMP在第
t
t
t次迭代过程中估计得到的
x
^
t
{hat {pmb x}}^{t}
xxx^t与真实
x
pmb x
xxx的均方差（MSE）
E
[
∥
x
−
x
^
t
∥
2
2
]
mathbb E[{Vert pmb x-{hat {pmb x}}^{t} Vert}^2_{2} ]
E[∥xxx−xxx^t∥22]，表示为
τ
t
2
→
1
N
∥
x
−
x
^
t
∥
2
2
,
(
N
→
∞
)
(2)
{ tau_{t} }^{2} rightarrow frac{1}{N} {Vert pmb x-{hat {pmb x}}^{t} Vert}^2_{2}, (N rightarrow infty) tag{2}
τt2→N1∥xxx−xxx^t∥22,(N→∞)(2)
事实上，式(2)为趋近于而不是严格等于，因为
1
N
∥
x
−
x
^
t
∥
2
2
frac{1}{N} {Vert pmb x-{hat {pmb x}}^{t} Vert}^2_{2}
N1∥xxx−xxx^t∥22是关于
x
^
t
{hat {pmb x}}^{t}
xxx^t的二阶Lipschitz函数，满足如下收敛关系
1
N
∥
x
−
x
^
t
∥
2
2
→
E
[
∥
x
−
x
^
t
∥
2
2
]
,
(
N
→
∞
)
(3)
frac{1}{N} {Vert pmb x-{hat {pmb x}}^{t} Vert}^2_{2} rightarrow mathbb E[{Vert pmb x-{hat {pmb x}}^{t} Vert}^2_{2} ] , (N rightarrow infty) tag{3}
N1∥xxx−xxx^t∥22→E[∥xxx−xxx^t∥22],(N→∞)(3)
但遗憾的是，大部分情况下，只有当
A
pmb A
AAA为高斯矩阵或者次高斯矩阵时，state evolution才能与AMP估计的结果统一，如果感知矩阵的特征值分布与高斯矩阵的特征值分布相差较远，AMP的性能就不能保证，甚至可能会出现不收敛的情况。为了解决该问题，Junjie Ma和Li Ping提出了OAMP[1]。
AMP还有一个重要的点是其线性迭代过程中含有"onsager"这一项，它的作用是为了消除迭代过程中感知矩阵
A
pmb A
AAA与估计结果
x
^
t
{hat {pmb x}}^{t}
xxx^t之间的相关性。虽然OAMP把AMP中的“onsager”这一项给去掉了，但是为了补偿"onsager"原来的作用，OAMP在非线性估计中加入了divergence-free的约束（可能divergence-free这个概念有点抽象，简单理解为导数为0就行）。
备注：其实OAMP并没有严格意义的数学推导，作者先是给了两个独立性的假设（假设1和假设2），而OAMP-state evolution就是由该假设条件推出来的。让OAMP迭代开始(
t
=
0
t=0
t=0时刻)先保证独立性，满足两个假设条件之一，如果之后的每一次迭代假设1和假设2都能互相推出对方，那么我们可以认为每一次迭代的程始终可以保证独立条件成立，也就保证了state evolution。但略有遗憾，假设1和假设2只能“部分”地互相推出对方，可以保证不相关（由正交推出），但是不能保证独立。幸运的是，仿真结果表明，虽然迭代过程中的独立条件不能始终保持，但是state evolution还是能够和OAMP估计结果统一，论文作者将此描述为：假设1和假设2是只是state evolution的充分条件。
3 AMP
在开始OAMP之前，先回顾一下AMP。
3.1 AMP算法
假设矩阵
A
=
[
a
1
,
a
2
,
…
,
a
N
]
pmb A=[pmb a_{1}, pmb a_{2}, text{…},pmb a_{N}]
AAA=[aaa1,aaa2,…,aaaN]是列归一化的，即，
∀
i
∈
{
1
,
...
,
N
}
,
a
i
∈
C
M
×
1
{forall i} in {1,text{...} ,N}, pmb a_{i} in mathbb C^{M times 1}
∀i∈{1,...,N},aaai∈CM×1，满足
E
[
∥
a
i
∥
2
]
=
1
mathbb E[{Vert pmb a_{i} Vert}_2]=1
E[∥aaai∥2]=1。AMP的迭代过程如下
r
t
=
s
t
+
A
T
(
y
−
A
s
t
)
+
N
M
<
η
t
−
1
′
(
r
t
−
1
)
>
(
r
t
−
1
−
s
t
−
1
)
⏟
o
n
s
a
g
e
r
t
e
r
m
(4a)
pmb r^{t}=pmb s^{t} + pmb A^{T}(pmb y - pmb A pmb s^{t}) +\ underbrace {frac {N}{M} <{eta_{t-1}}^{'}( pmb r^{t-1})>( pmb r^{t-1}-pmb s^{t-1} )}_{onsager text{ } term} tag{4a}
rrrt=ssst+AAAT(yyy−AAAssst)+onsager term
MN<ηt−1′(rrrt−1)>(rrrt−1−ssst−1)(4a)
s
t
+
1
=
η
t
(
r
t
)
(4b)
pmb {s^{t+1}}={eta}_{t}(pmb r^t) tag{4b}
st+1st+1st+1=ηt(rrrt)(4b)
其中
η
t
{eta}_{t}
ηt是关于
r
t
pmb r^t
rrrt的一个Lipschitz连续函数(component-wise)，
s
t
pmb s^t
ssst是最后的估计，
r
t
−
s
t
pmb r^{t}-pmb s^{t}
rrrt−ssst其实就是一般AMP表达的“残差”(residual error)。式(4a)的最后一项就是前言所描述的"onsager"项，这一项是根据松弛信念传播算法（relaxed-Belief-Propagation, relaxed-BP）在极限条件下
M
,
N
→
∞
M,N rightarrow infty
M,N→∞时依据大数定律、中心极限定理和消息近似补偿（为了补偿经过近似的
O
(
1
N
)
O(frac {1}{sqrt N})
O(N
1)的消息），以及泰勒展开逐步推导出来的。
"onsager"项隐含的意义："onsager"项的存在确保了AMP-state evolution的正确性，但是在推导AMP之前，relaxed-BP里边的方差项其实跟state-evolution是对应的（只是对应，并不相等），区别在于当
M
,
N
<
∞
M,N lt infty
M,N<∞时，relaxed-BP的方差项（指消息传递过程中，消息分布的方差逐渐积累）一般记为
Σ
m
n
(
t
)
Sigma^n_{m}(t)
Σmn(t)，这里不去管具体的
n
,
m
n,m
n,m是什么含义了，
t
t
t是指迭代到第
t
t
t步，只是它隐式地强调了
A
pmb A
AAA和估计
s
t
pmb s^{t}
ssst具有相关性。而
M
,
N
→
∞
M,N rightarrow infty
M,N→∞时借助大数定理和中心极限定理进一步推导可以使得
lim
M
,
N
→
∞
Σ
m
n
(
t
)
→
Σ
(
t
)
lim_{M,Nto infty}Sigma^n_{m}(t) rightarrow Sigma(t)
limM,N→∞Σmn(t)→Σ(t)，即去掉了相关性。
3.2 AMP-state evolution与等效信号模型
为了方便之后OAMP的讨论，下面的叙述几乎都是基于OAMP原论文的，表达式跟一般的AMP表达式有些差异，是为了类比，在后面方便证明正交性，本质上并无差异。
定义两种误差，第一种误差是非线性估计结果
s
t
pmb s^t
ssst与真实值的差，第二种误差是线性估计结果
r
t
pmb r^t
rrrt与真实值的差，分别定义为
q
t
=
s
t
−
x
(5a)
pmb q^t = pmb s^t - pmb x tag{5a}
qqqt=ssst−xxx(5a)
h
t
=
r
t
−
x
(5b)
pmb h^t = pmb r^t - pmb x tag{5b}
hhht=rrrt−xxx(5b)
结合(5), (4)可以被展开为
h
t
=
(
I
−
A
T
A
)
q
t
+
A
T
n
+
N
M
<
η
t
−
1
′
(
x
+
h
t
−
1
)
>
(
h
t
−
1
−
q
t
−
1
)
(6a)
pmb h^t = (pmb I - pmb A^T pmb A)pmb q^t + pmb A^Tpmb n + \frac {N}{M}<eta^{'}_{t-1}(pmb x + pmb h^{t-1})>(pmb h^{t-1}-pmb q^{t-1}) tag{6a}
hhht=(III−AAATAAA)qqqt+AAATnnn+MN<ηt−1′(xxx+hhht−1)>(hhht−1−qqqt−1)(6a)
q
t
+
1
=
η
t
(
x
+
h
t
)
−
x
(6b)
pmb q^{t+1}=eta_{t}(pmb x + pmb h^t)-pmb x tag{6b}
qqqt+1=ηt(xxx+hhht)−xxx(6b)
式(6)并不是要给迭代算法，因为真实值
x
pmb x
xxx作为一个已知的量参与其中。
与AMP对应的evolution由下式直接给出，这里不展开描述，因为AMP-evolution的证明过程极其复杂，但是如果只是理解，有一个简单的思路，就是如上面所述，从relaxed-BP入手，着眼于其方差演变，可以在一定程度上帮助理解。
τ
t
2
=
N
M
v
t
2
+
σ
2
(7a)
tau^2_{t}=frac {N}{M} v^2_{t}+sigma^2 tag{7a}
τt2=MNvt2+σ2(7a)
v
t
+
1
2
=
E
{
[
η
t
(
X
+
τ
t
Z
)
−
X
]
2
}
(7b)
v^2_{t+1}=mathbb E{[eta_{t}(X+tau_{t}Z)-X]^2} tag{7b}
vt+12=E{[ηt(X+τtZ)−X]2}(7b)
注意式(7b)中的
X
,
Z
X,Z
X,Z表示随机变量，并且
Z
∼
N
(
0
,
1
)
Z sim mathcal N(0,1)
Z∼N(0,1)与
X
X
X独立。这里的理解对AMP以及AMP-state evolution是至关重要的，因为其独立性，我们可以将每次迭代估计得到的结果模型等效为
X
^
t
=
X
+
τ
t
2
Z
,
Z
∼
N
(
0
,
1
)
(8)
hat X^t=X+tau^2_{t}Z, Z sim mathcal N(0,1) tag{8}
X^t=X+τt2Z,Z∼N(0,1)(8)
这也意味着
X
^
t
∼
N
(
X
,
τ
t
2
)
(9)
hat X^t sim mathcal N(X,tau^2_{t}) tag{9}
X^t∼N(X,τt2)(9)
τ
t
2
tau^2_{t}
τt2指的是估计后的方差，结合式(7a)可以看出，方差项由高斯噪声的方差
σ
2
sigma^2
σ2和非线性估计结果与真实值之间的MSE构成。这反映了两点，一方面，高斯噪声与其他各个变量都保持独立，另一方面，AMP从始至终虽然都旨在接近真实值，因为
τ
t
tau^t
τt一直在变小，但是依然没有克服高斯噪声。
等效模型还有一个重要的作用就是，在部分
η
eta
η函数的推导过程中，使用等效模型的概念可以一定程度上简化推导过程，比如推导MMSE函数，或者退化为LMMSE等。以及，
τ
t
2
tau^2_{t}
τt2可能作为其中的一个参数，然而实际的仿真中我们并不知道真实
τ
t
2
tau^2_{t}
τt2，这个时候就需要将其近似为
τ
^
t
2
hat tau^2_{t}
τ^t2，具体表示为
τ
^
t
2
=
1
N
∥
r
t
−
s
t
∥
2
2
(10)
hat tau^2_{t}=frac {1}{N} {Vert pmb r^t - pmb s^t Vert}^2_{2} tag {10}
τ^t2=N1∥rrrt−ssst∥22(10)
4 OAMP
4.1 OAMP产生的动机
在阐述动机之前，先描述论文给出的一个定义
定义1： Divergence-free
对函数η
:
R
↦
R
eta: mathbb R mapsto mathbb R
η:R↦R，如果满足
E
[
η
′
(
R
)
]
=
0
mathbb E [eta^{'}(R)]=0
E[η′(R)]=0
那么就认为
η
eta
η是devergence-free
根据定义1，一个divergence-free的函数可以被构造成
η
(
r
)
=
C
⋅
(
η
^
(
r
)
−
E
[
η
^
(
R
)
]
⋅
r
)
(11)
eta(r)=C cdot (hat eta( r)-mathbb E[hat eta(R)] cdot r) tag{11}
η(r)=C⋅(η^(r)−E[η^(R)]⋅r)(11)
其中
η
^
hat eta
η^是任意一个一阶可导的函数，C是常数。
如果把AMP迭代公式(4)中的"onsager"项给去掉，
η
eta
η函数按照式(11)的方式给出，其中
η
^
hat eta
η^是软阈值函数（压缩感知常用来恢复稀疏信号的函数），在这样的设置下，作者发现即使感知矩阵
A
pmb A
AAA不是高斯矩阵，而是一个离散余弦变换矩阵，state evolution的结果却意外地与去"onsager"的AMP迭代结果一致。这个发现也就引出了OAMP。
4.2 去相关的线性估计
先给出两个定义
定义2： 酉不变矩阵(Unitarily-Invariant Matrix)
如果矩阵U
,
V
,
Σ
pmb U,pmb V,pmb Sigma
UUU,VVV,ΣΣΣ三者之间相互独立，并且
U
,
V
pmb U,pmb V
UUU,VVV满足Haar分布(随机各向同性)是正交阵，
Σ
pmb Sigma
ΣΣΣ是对角阵，那么认为
A
=
U
Σ
V
pmb A= pmb U pmb Sigma pmb V
AAA=UUUΣΣΣVVV是酉不变的。
定义3： 去相关矩阵
如果感知矩阵A
=
U
Σ
V
pmb A= pmb U pmb Sigma pmb V
AAA=UUUΣΣΣVVV是酉不变的，矩阵
W
pmb W
WWW如果满足
t
r
(
I
−
W
A
)
=
0
tr(pmb I - pmb W pmb A)=0
tr(III−WWWAAA)=0，就说
W
pmb W
WWW是关于
A
pmb A
AAA的一个去相关矩阵。指定准去相关矩阵
W
^
=
U
G
V
T
(12)
hat {pmb W}=pmb U pmb G pmb V^T tag{12}
WWW^=UUUGGGVVVT(12)
那么满足
t
r
(
I
−
W
A
)
tr(pmb I - pmb W pmb A)
tr(III−WWWAAA)的去相关矩阵
W
pmb W
WWW可以被构建为
W
=
N
t
r
(
W
^
A
)
W
^
(13)
pmb W= frac {N}{tr(hat {pmb W} pmb A)} hat {pmb W} tag{13}
WWW=tr(WWW^AAA)NWWW^(13)
其实定义3在论文中并没有直接给出，是我抽出来的，而且跟论文稍有出入，但是问题不大。我刚开始读到这里的去相关概念的时候非常不理解，如果矩阵A的映射是单射或者双射，可能稍微懂个大概，满射可能把
W
pmb W
WWW乘在右边，但是概念依然很模糊，后来请教了一下数院的学长，大概意思就说作者就是这么叫了而已。所以这里也就理解个大概吧，没有再细究。
下面给出三个常用的准去相关矩阵
（1）匹配滤波器(Matched Filter, MF)
W
^
M
F
=
A
T
(14a)
hat {pmb W}^{MF}=pmb A^T tag{14a}
WWW^MF=AAAT(14a)
（2）伪逆(Pseudo-inverse, PINV)
W
^
P
I
N
V
=
{
A
T
(
A
A
T
)
−
1
;
M
<
N
(
A
T
A
)
−
1
A
T
;
M
>
N
(14b)
hat {pmb W}^{PINV}= left{ begin{array}{lr} pmb A^T (pmb A pmb A^T)^{-1}; M<N \ (pmb A^T pmb A)^{-1} pmb A^T; M>N end{array} right tag{14b}.
WWW^PINV={AAAT(AAAAAAT)−1;M<N(AAATAAA)−1AAAT;M>N(14b)
（3）LMMSE
W
^
L
M
M
S
E
=
A
T
(
A
A
T
+
σ
2
v
2
I
)
−
1
(14c)
hat {pmb W}^{LMMSE}=pmb A^T(pmb A pmb A^T+frac {sigma^2}{v^2} pmb I)^{-1} tag{14c}
WWW^LMMSE=AAAT(AAAAAAT+v2σ2III)−1(14c)
其实式(14c)中
v
2
v^2
v2的含义有些微妙，将在后面合适的地方做更深的阐述。
LMMSE的简述
假设模型为y
=
A
x
+
n
pmb y = pmb A pmb x + pmb n
yyy=AAAxxx+nnn，假设随机向量
X
,
Y
X,Y
X,Y的均值都为0（实际当中不满足的话可以先减去均值），
n
∼
C
N
(
0
,
σ
2
I
)
pmb n sim mathcal C mathcal N(0,sigma^2 pmb I)
nnn∼CN(0,σ2III)，那么LMMSE的估计为
x
^
=
Σ
x
y
Σ
y
−
1
y
hat {pmb x} = Sigma_{xy} Sigma^{-1}_{y} pmb y
xxx^=ΣxyΣy−1yyy
其中
Σ
x
y
=
E
[
x
y
H
]
=
E
[
x
(
A
x
+
n
)
H
]
=
E
[
x
x
H
A
+
x
n
H
]
=
Σ
x
A
Sigma_{xy} = mathbb E[pmb x pmb y^{H}]=mathbb E[pmb x pmb (pmb A pmb x+pmb n)^{H}]=mathbb E[pmb x pmb x^Hpmb A+pmb xpmb n^H]=Sigma_{x} pmb A
Σxy=E[xxxyyyH]=E[xxx(((AAAxxx+nnn)H]=E[xxxxxxHAAA+xxxnnnH]=ΣxAAA
Σ
y
=
E
[
y
y
H
]
=
E
[
(
A
x
+
n
)
(
A
x
+
n
)
H
]
=
A
Σ
x
A
H
+
σ
2
I
begin{aligned} Sigma_{y} &= mathbb E[pmb y pmb y^{H}] \ &=mathbb E[(pmb A pmb x+pmb n) pmb (pmb A pmb x+pmb n)^{H}] \ &=pmb A Sigma_{x} pmb A^H+sigma^2 pmb I end{aligned}
Σy=E[yyyyyyH]=E[(AAAxxx+nnn)(((AAAxxx+nnn)H]=AAAΣxAAAH+σ2III
那么LMMSE的估计可以转化为
x
^
=
Σ
x
A
(
A
Σ
x
A
H
+
σ
2
I
)
−
1
y
hat {pmb x} = Sigma_{x} pmb A (pmb A Sigma_{x} pmb A^H+sigma^2 pmb I)^{-1}y
xxx^=ΣxAAA(AAAΣxAAAH+σ2III)−1y
该结果还可以继续延申，根据
(
E
+
B
C
D
)
−
1
=
E
−
1
−
E
−
1
B
(
C
−
1
+
D
E
−
1
B
)
−
1
D
E
−
1
(pmb E + pmb B pmb C pmb D)^{-1}=pmb E^{-1}- pmb E^{-1} pmb B (pmb C^{-1}+ D pmb E^{-1} pmb B)^{-1} pmb D pmb E^{-1}
(EEE+BBBCCCDDD)−1=EEE−1−EEE−1BBB(CCC−1+DEEE−1BBB)−1DDDEEE−1
将上述公式代入
Σ
y
Sigma_{y}
Σy项即可。
4.3 OAMP算法
OAMP的迭代公式
r
t
=
s
t
+
W
t
(
y
−
A
s
t
)
(15a)
pmb r^{t}=pmb s^{t} + pmb W_{t}(pmb y - pmb A pmb s^{t}) tag{15a}
rrrt=ssst+WWWt(yyy−AAAssst)(15a)
s
t
+
1
=
η
t
(
r
t
)
(15b)
pmb s^{t+1} = eta_{t}(pmb r^t) tag{15b}
ssst+1=ηt(rrrt)(15b)
其中
W
t
pmb W_{t}
WWWt是去相关矩阵，
η
t
eta_{t}
ηt是满足divergence-free的约束，即式(11)。将该式与AMP迭代式(4)做比较，可以发现线性估计中的矩阵
A
T
pmb A^T
AAAT变得更一般化，不再局限于匹配滤波，而且末尾缺少了"onsager"项，把“onsager”的作用加在了divergence-free约束里，这也跟OAMP动机部分所阐述的内容一致。
4.4 估计误差迭代与OAMP-state evolution
依然保持式(5a, 5b)所述的误差符号
h
t
,
q
t
pmb h^t, pmb q^t
hhht,qqqt，可以类比AMP中的式(6)写出OAMP的误差迭代，如下
h
t
=
B
t
q
t
+
W
t
n
(16a)
pmb h^t = pmb B_{t} pmb q^t + pmb W_{t} pmb n tag{16a}
hhht=BBBtqqqt+WWWtnnn(16a)
q
t
+
1
=
η
t
(
x
+
h
t
)
−
x
(16b)
pmb q^{t+1}=eta_{t}(pmb x+pmb h^t) - pmb x tag{16b}
qqqt+1=ηt(xxx+hhht)−xxx(16b)
其中
B
t
=
I
−
W
t
A
pmb B_{t} = pmb I-pmb W_{t} pmb A
BBBt=III−WWWtAAA，然后如AMP一样，指定
τ
t
2
=
1
N
E
[
∥
h
t
∥
2
2
]
(17a)
tau^2_{t}=frac {1}{N} mathbb E[{Vert pmb h^t Vert}^2_{2}] tag{17a}
τt2=N1E[∥hhht∥22](17a)
v
t
+
1
2
=
1
N
E
[
∥
q
t
+
1
∥
2
2
]
(17b)
v^2_{t+1}=frac {1}{N} mathbb E[{Vert pmb q^{t+1} Vert}^2_{2}] tag{17b}
vt+12=N1E[∥qqqt+1∥22](17b)
式(17)就是所谓的state evolution，可以对式(17a)做进一步推导
τ
t
2
=
1
N
E
[
∥
h
t
∥
2
2
]
=
1
N
E
[
∥
B
t
q
t
+
W
t
n
∥
2
2
]
=
1
N
{
E
[
∥
B
t
q
t
∥
2
2
]
+
E
[
∥
W
t
n
∥
2
2
]
}
=
1
N
{
E
[
t
r
(
(
q
t
)
T
B
t
T
B
t
q
t
)
)
]
+
E
[
t
r
(
n
T
W
t
T
W
t
n
)
]
}
=
1
N
{
E
[
t
r
(
B
t
T
B
t
)
]
⋅
E
[
∥
q
t
∥
2
2
]
+
E
[
t
r
(
W
t
T
W
t
)
]
⋅
E
[
∥
n
∥
2
2
]
}
=
1
N
{
E
[
t
r
(
B
t
T
B
t
)
]
⋅
N
v
t
2
+
E
[
t
r
(
W
t
T
W
t
)
]
⋅
M
σ
2
}
=
E
[
t
r
(
B
t
T
B
t
)
]
⋅
v
t
2
+
M
N
E
[
t
r
(
W
t
T
W
t
)
]
⋅
σ
2
begin{aligned} tau^2_{t}&= frac {1}{N} mathbb E[{Vert pmb h^t Vert}^2_{2}] \ &=frac {1}{N} mathbb E[{Vert pmb B_{t} pmb q^t + pmb W_{t} pmb n Vert}^2_{2}] \ &=frac {1}{N} { mathbb E [{Vert pmb B_t pmb q^t Vert}^2_{2}] + mathbb E [{Vert pmb W_t pmb n Vert}^2_{2}] } \ &=frac {1}{N} { mathbb E [tr ( (pmb q^t)^T pmb B^T_t pmb B_t pmb q^t) )] + mathbb E [tr(pmb n^T pmb W^T_t pmb W_t pmb n)] } \ &=frac {1}{N} { mathbb E [tr(pmb B^T_t pmb B_t )] cdot E[{Vert pmb q^{t} Vert}^2_{2}] + mathbb E [tr(pmb W^T_t pmb W_t )] cdot E[{Vert pmb n Vert}^2_{2}] } \ &=frac {1}{N} { mathbb E [tr(pmb B^T_t pmb B_t )] cdot N v^2_t + mathbb E [tr(pmb W^T_t pmb W_t )] cdot M sigma^2 } \ &=mathbb E [tr(pmb B^T_t pmb B_t )] cdot v^2_t + frac {M}{N} mathbb E [tr(pmb W^T_t pmb W_t )] cdot sigma^2 end{aligned}
τt2=N1E[∥hhht∥22]=N1E[∥BBBtqqqt+WWWtnnn∥22]=N1{E[∥BBBtqqqt∥22]+E[∥WWWtnnn∥22]}=N1{E[tr((qqqt)TBBBtTBBBtqqqt))]+E[tr(nnnTWWWtTWWWtnnn)]}=N1{E[tr(BBBtTBBBt)]⋅E[∥qqqt∥22]+E[tr(WWWtTWWWt)]⋅E[∥nnn∥22]}=N1{E[tr(BBBtTBBBt)]⋅Nvt2+E[tr(WWWtTWWWt)]⋅Mσ2}=E[tr(BBBtTBBBt)]⋅vt2+NME[tr(WWWtTWWWt)]⋅σ2
注意：上面推导出来的结果与OAMP论文式(23)给出的在系数上有差异，感觉应该上面是正确的，不管如何，思路应该没什么问题。那么，据此，就可以轻易地写出OAMP-state evolution
τ
t
2
=
E
[
t
r
(
B
t
T
B
t
)
]
⋅
v
t
2
+
M
N
E
[
t
r
(
W
t
T
W
t
)
]
⋅
σ
2
(18a)
tau^2_{t}=mathbb E [tr(pmb B^T_t pmb B_t )] cdot v^2_t + frac {M}{N} mathbb E [tr(pmb W^T_t pmb W_t )] cdot sigma^2 tag{18a}
τt2=E[tr(BBBtTBBBt)]⋅vt2+NME[tr(WWWtTWWWt)]⋅σ2(18a)
v
t
+
1
2
=
E
[
∣
η
t
(
X
+
τ
t
Z
)
−
X
∣
2
]
(18b)
v^2_{t+1}=mathbb E[{vert eta_t(X+tau_t Z) - X vert}^2] tag{18b}
vt+12=E[∣ηt(X+τtZ)−X∣2](18b)
其中
X
∼
P
X
(
x
)
X sim P_X(x)
X∼PX(x)，与
Z
∼
N
(
0
,
1
)
Z sim mathcal N(0,1)
Z∼N(0,1)独立。
4.5 关于OAMP的合理性以及两个重要假设
关于OAMP的合理性，前言部分已经做了简短的铺垫，这里再详细展开。论文作者提出了两个假设，虽然OAMP的证明并不严格，但是基于两个假设展开的讨论还是比较合理的。
假设1：式(16a)中的
h
t
∼
N
(
0
,
τ
t
2
)
pmb h^t sim mathcal N(0,tau^2_t)
hhht∼N(0,τt2)，并且独立于真实值
x
pmb x
xxx。
假设2：式(16b)中的q
t
+
1
pmb q^{t+1}
qqqt+1里的元素是独立同分布的(i.i.d)，并且独立于
A
,
n
pmb A, pmb n
AAA,nnn。
一般的条件会有
x
pmb x
xxx是i.i.d.的，独立于
A
,
n
pmb A, pmb n
AAA,nnn，在OAMP中，当迭代次数
t
=
−
1
t=-1
t=−1时，
q
0
=
−
x
pmb q^0=- pmb x
qqq0=−xxx，因此假设2在
t
=
−
1
t=-1
t=−1时是成立的。虽然OAMP的线性迭代式(15a)比AMP的线性迭代式(4a)少了"onsager"这一项，但是只有我们能证明假设2在迭代过程中一直成立，那么式(15a)便是合理的，因为"onsager"的作用就是去除迭代估计结果与
A
pmb A
AAA的相关性。接下来会提出两个推论来更直观地理解OAMP。
4.5.1 从假设2看假设1
推论1：如果假设2是成立的，矩阵
A
pmb A
AAA是酉不变的，
W
t
pmb W_t
WWWt是去相关矩阵，那么就有
h
t
pmb h^t
hhht的元素与
x
pmb x
xxx的元素不相关，以及，
h
t
pmb h^t
hhht的元素彼此之间不相关，而且它们拥有共同的方差，均值为0。
证明：从式(16b)
q
t
+
1
=
η
t
(
x
+
h
t
)
−
x
pmb q^{t+1}=eta_{t}(pmb x+pmb h^t) - pmb x
qqqt+1=ηt(xxx+hhht)−xxx可以看出，
q
t
pmb q^t
qqqt与
x
pmb x
xxx具有相关性，这可能会进一步导致
h
t
pmb h^t
hhht与
x
pmb x
xxx的相关，因为式(16a)中
h
t
pmb h^t
hhht由
q
t
pmb q^t
qqqt生成。但去相关矩阵
W
t
pmb W_t
WWWt的引入可以抑制此相关性，具体描述如下。
因为A
=
V
Σ
U
T
pmb A=pmb V Sigma pmb U^T
AAA=VVVΣUUUT，
W
t
=
U
G
t
V
T
pmb W_t=pmb U pmb G_t pmb V^T
WWWt=UUUGGGtVVVT，
B
=
I
−
W
t
A
=
U
(
I
−
G
t
Σ
)
U
T
pmb B = pmb I- pmb W_t pmb A = pmb U(pmb I-pmb G_t Sigma) pmb U^T
BBB=III−WWWtAAA=UUU(III−GGGtΣ)UUUT，那么
E
U
[
(
B
t
)
i
,
j
]
=
∑
m
=
1
N
E
[
U
i
,
m
U
j
,
m
]
⋅
(
1
−
g
m
λ
m
)
mathbb E_U[(pmb B_t)_{i,j}] = sum_{m=1}^N mathbb E[U_{i,m} U_{j,m}] cdot (1-g_m lambda_m)
EU[(BBBt)i,j]=m=1∑NE[Ui,mUj,m]⋅(1−gmλm)
其中
λ
m
,
g
m
lambda_m, g_m
λm,gm分别是矩阵
A
,
W
t
pmb A, pmb W_t
AAA,WWWt的奇异值，并且当
m
>
M
m>M
m>M时，
λ
m
=
g
m
=
0
lambda_m=g_m=0
λm=gm=0。
对于一个Haar分布的矩阵U
pmb U
UUU，有
E
[
U
i
,
m
U
j
,
m
]
=
{
0
;
i
≠
j
N
−
1
;
i
=
j
mathbb E[U_{i,m} U_{j,m}]= left{ begin{array}{lr} 0; i neq j \ N^{-1}; i=j end{array} right.
E[Ui,mUj,m]={0; i=jN−1; i=j
那么就有
E
U
[
(
B
t
)
i
,
j
]
=
{
0
;
i
≠
j
1
N
t
r
(
B
t
)
;
i
=
j
mathbb E_U[(pmb B_t)_{i,j}]= left{ begin{array}{lr} 0; i neq j \ frac {1}{N} tr ( pmb B_t ); i=j end{array} right.
EU[(BBBt)i,j]={0; i=jN1tr(BBBt); i=j
因为
W
t
pmb W_t
WWWt是去相关矩阵，根据定义3，有
t
r
(
B
t
)
=
t
r
(
I
−
W
t
A
)
=
0
tr(pmb B_t)=tr(pmb I - pmb W_t pmb A)=0
tr(BBBt)=tr(III−WWWtAAA)=0，所以
E
[
B
t
]
=
0
mathbb E[pmb B_t] = 0
E[BBBt]=0
假设2给出了条件，
q
t
pmb q^t
qqqt独立于
A
pmb A
AAA，显然也独立于
B
t
pmb B_t
BBBt，那么
E
[
h
t
]
=
E
[
B
t
q
t
]
+
E
[
W
t
n
]
=
E
[
B
t
]
E
[
q
t
]
+
E
[
W
t
]
E
[
n
]
=
0
N
(19)
begin{aligned} mathbb E[pmb h^t] &= mathbb E[pmb B_t pmb q^t] + mathbb E[pmb W_t pmb n] \ &=mathbb E[pmb B_t] mathbb E[pmb q^t] + mathbb E[pmb W_t] mathbb E[pmb n] \ &=pmb 0_N end{aligned} tag{19}
E[hhht]=E[BBBtqqqt]+E[WWWtnnn]=E[BBBt]E[qqqt]+E[WWWt]E[nnn]=000N(19)
又
h
t
=
B
t
q
t
+
W
t
n
pmb h^t = pmb B_{t} pmb q^t + pmb W_{t} pmb n
hhht=BBBtqqqt+WWWtnnn，要证
x
pmb x
xxx与
h
t
pmb h^t
hhht不相关，只需要证
x
pmb x
xxx与
B
t
q
t
pmb B_{t} pmb q^t
BBBtqqqt不相关（已经有
E
[
B
t
q
t
]
=
0
N
mathbb E[pmb B_t pmb q^t]=pmb 0_N
E[BBBtqqqt]=000N，所以满足正交性即可）
E
[
h
t
x
T
]
=
E
[
B
t
q
t
x
T
]
=
E
[
B
t
]
E
[
q
t
x
T
]
=
0
N
×
N
(20)
mathbb E[pmb h^t pmb x^T]= mathbb E[pmb B_t pmb q^t pmb x^T]=mathbb E[pmb B_t] mathbb E[ pmb q^t pmb x^T]=pmb 0_{N times N} tag{20}
E[hhhtxxxT]=E[BBBtqqqtxxxT]=E[BBBt]E[qqqtxxxT]=000N×N(20)
所以，根据
E
[
h
t
]
=
0
N
mathbb E[pmb h^t]=pmb 0_N
E[hhht]=000N和
E
[
h
t
x
T
]
=
0
N
×
N
mathbb E[pmb h^t pmb x^T]=pmb 0_{N times N}
E[hhhtxxxT]=000N×N（正交性），必然有
x
pmb x
xxx与
h
t
pmb h^t
hhht不相关。要证
h
t
pmb h^t
hhht的元素彼此之间不相关，而且它们拥有共同的方差，均值为0，只需证明
h
t
pmb h^t
hhht为对角阵的系数，这里省略。
事实上，因为
E
[
B
t
]
=
E
[
n
]
mathbb E[pmb B_t]=mathbb E[pmb n]
E[BBBt]=E[nnn]都为0，式(19)隐含了
h
t
,
q
t
pmb h^t ,pmb q^t
hhht,qqqt的正交性
E
[
h
t
(
q
t
)
T
]
=
0
N
×
N
(21)
mathbb E[pmb h^t (pmb q^t)^T]=pmb 0_{N times N} tag{21}
E[hhht(qqqt)T]=000N×N(21)
从推论1的证明过程中可以看出，去相关矩阵
W
t
pmb W_t
WWWt在里边起到了重要的去相关作用(间接地借助了
B
=
I
−
W
t
A
pmb B = pmb I- pmb W_t pmb A
BBB=III−WWWtAAA)。除此之外，还可以看出OAMP对矩阵
A
pmb A
AAA的特征值没有任何束缚，所以潜在的应用范围会更广一些相对于AMP。
4.5.2 从假设1看假设2
在这一部分，我们尝试基于假设1，来推出假设2，从式(16)可以看出，如果
q
t
+
1
pmb q^{t+1}
qqqt+1与
h
t
pmb h^t
hhht独立，那么
q
t
+
1
pmb q^{t+1}
qqqt+1与
A
,
n
pmb A, pmb n
AAA,nnn也独立，因为(16b)中
q
t
+
1
pmb q^{t+1}
qqqt+1只跟
h
t
pmb h^t
hhht和
x
pmb x
xxx有关系，那么就存在这样一条马尔可夫链(注意上标
t
t
t)
A
,
n
→
h
t
→
q
t
+
1
pmb A, pmb n rightarrow pmb h^t rightarrow pmb q^{t+1}
AAA,nnn→hhht→qqqt+1
也就是说，如果我们能够证明
q
t
+
1
pmb q^{t+1}
qqqt+1与
h
t
pmb h^t
hhht独立，那么假设2就自然而然成立。但遗憾的是我们并不能证明独立性，只能证明正交性，再推广到不相关（推论2将阐述）。
回到定义1里边所阐述的divergence-free函数，利用Lipchitz函数期望的收敛性，我们可以构建一个近似divergence-free函数，如下
η
t
(
r
t
)
=
C
⋅
{
η
^
t
(
r
t
)
−
(
1
N
∑
j
=
1
N
η
^
t
′
(
r
j
t
)
)
⋅
r
t
}
(22)
eta_t(pmb r^t)=C cdot { hat eta_t(pmb r^t)-mathbb (frac{1}{N} sum_{j=1}^N hat eta^{'}_t( r^t_j)) cdot pmb r^t } tag{22}
ηt(rrrt)=C⋅{η^t(rrrt)−(N1j=1∑Nη^t′(rjt))⋅rrrt}(22)
备注：式(22)的近似divergence-free函数与OAMP的迭代公式(15)结合，便是真正的实际所采用的OAMP算法。
推论2：如果
η
eta
η是divergence-free函数，那么
E
[
τ
t
Z
⋅
η
(
X
+
τ
t
Z
)
]
=
0
(23)
mathbb E[tau_t Z cdot eta(X+tau_t Z )]=0 tag{23}
E[τtZ⋅η(X+τtZ)]=0(23)
在证明上式之前，先引入一个重要引理
Stein引理：对∀
ψ
:
R
↦
R
forall psi: R mapsto R
∀ψ:R↦R，使得下式中的期望存在，有
E
[
Z
⋅
ψ
(
Z
)
]
=
E
[
ψ
′
(
Z
)
]
(24)
mathbb E[Z cdot psi(Z)]=mathbb E[psi^{'}(Z)] tag{24}
E[Z⋅ψ(Z)]=E[ψ′(Z)](24)
让
ϕ
(
Z
)
=
η
t
(
X
+
τ
t
Z
)
phi(Z)=eta_t(X+tau_t Z)
ϕ(Z)=ηt(X+τtZ)，根据Stein引理，
E
[
τ
t
Z
⋅
η
t
(
X
+
τ
t
Z
)
]
=
τ
t
⋅
E
X
{
E
Z
∣
X
[
Z
⋅
η
t
(
X
+
τ
t
Z
)
]
}
=
τ
t
2
⋅
E
X
{
E
Z
∣
X
[
η
t
′
(
X
+
τ
t
Z
)
]
}
=
τ
t
2
⋅
E
[
η
t
′
(
X
+
τ
t
Z
)
]
begin{aligned} mathbb E[tau_t Z cdot eta_t(X+tau_t Z )] &= tau_t cdot mathbb E_X { mathbb E_{Z|X} [Z cdot eta_t(X+tau_t Z )] } \ &=tau^2_t cdot mathbb E_X { mathbb E_{Z|X} [ eta^{'}_t(X+tau_t Z )] } \ &=tau^2_t cdot mathbb E [ eta^{'}_t(X+tau_t Z )] end{aligned}
E[τtZ⋅ηt(X+τtZ)]=τt⋅EX{EZ∣X[Z⋅ηt(X+τtZ)]}=τt2⋅EX{EZ∣X[ηt′(X+τtZ)]}=τt2⋅E[ηt′(X+τtZ)]
因为
η
t
eta_t
ηt是divergence-free函数，所以
E
[
η
t
′
(
X
+
τ
t
Z
)
]
=
0
mathbb E [ eta^{'}_t(X+tau_t Z )]=0
E[ηt′(X+τtZ)]=0，也就证明了(23)。而(23)等价于
E
[
(
R
t
−
X
)
⋅
(
η
t
(
R
t
)
−
X
)
]
=
0
(25)
mathbb E[(R^t-X) cdot (eta_t(R^t)-X)]=0 tag{25}
E[(Rt−X)⋅(ηt(Rt)−X)]=0(25)
其中
R
t
=
X
+
τ
t
Z
R^t= X+ tau_t Z
Rt=X+τtZ。
式(25)直接推广可以得到
E
[
h
t
(
q
t
+
1
)
T
]
=
0
N
×
N
(26)
mathbb E[pmb h^t (pmb q^{t+1})^T]=pmb 0_{N times N} tag{26}
E[hhht(qqqt+1)T]=000N×N(26)
强调为什么叫OAMP：式(21)说明了线性估计(16a)中，“input-error”
q
t
pmb q^t
qqqt和"output-error"
h
t
pmb h^t
hhht是正交的（由推论1得出）；式(26)说明了非线性估计(16b)中，"before-error"
h
t
pmb h^t
hhht和"after-error"
q
t
+
1
pmb q^{t+1}
qqqt+1是正交的（由推论2得出）。这也就是OAMP名字里正交的由来。
4.6 MSE估计和state evolution仿真
在开始这部分之前，先回顾一下式(14c)LMMSE中的参数
v
2
v^2
v2，OAMP迭代公式
r
t
=
s
t
+
W
t
(
y
−
A
s
t
)
=
s
t
+
W
t
(
A
(
x
−
s
t
)
+
n
)
pmb r^{t}=pmb s^{t} + pmb W_{t}(pmb y - pmb A pmb s^{t})=pmb s^{t} + pmb W_{t}(pmb A ( pmb x - pmb s^{t} )+ pmb n)
rrrt=ssst+WWWt(yyy−AAAssst)=ssst+WWWt(AAA(xxx−ssst)+nnn)，所以LMMSE中的参数
v
2
v^2
v2表示
E
[
(
x
−
s
t
)
(
x
−
s
t
)
T
]
=
v
2
I
mathbb E[( pmb x - pmb s^{t} ) ( pmb x - pmb s^{t} )^T]=v^2 pmb I
E[(xxx−ssst)(xxx−ssst)T]=v2III
两个MSE：
v
t
2
=
1
N
E
[
∥
q
t
∥
2
2
]
v^2_t=frac{1}{N}{mathbb E[{Vert pmb q^t Vert}^2_2]}
vt2=N1E[∥qqqt∥22]，
τ
t
2
=
1
N
E
[
∥
h
t
∥
2
2
]
tau^2_t=frac{1}{N}{mathbb E[{Vert pmb h^t Vert}^2_2]}
τt2=N1E[∥hhht∥22]，它们可以看作是去相关矩阵
W
t
pmb W_t
WWWt和divergence-free函数
η
t
eta_t
ηt的两个参数（这也就是为什么
W
t
pmb W_t
WWWt有下标
t
t
t的原因）
4.6.1 非线性均方误差的估计
非线性均方误差
v
t
2
v^2_t
vt2的估计表达式
v
^
t
2
=
1
N
∥
y
−
A
s
t
∥
2
2
−
M
⋅
σ
2
t
r
(
A
T
A
)
(27)
hat {v}^2_t=frac{1}{N} frac {{Vert {pmb y - pmb A pmb s^t} Vert}^2_2 - M cdot sigma^2} {tr(pmb A^T pmb A)} tag{27}
v^t2=N1tr(AAATAAA)∥yyy−AAAssst∥22−M⋅σ2(27)
(27)的理解：
E
[
∥
y
−
A
s
t
∥
2
2
]
=
E
[
∥
A
x
−
A
s
t
+
n
∥
2
2
]
=
E
[
∥
A
(
x
−
s
t
)
∥
2
2
]
+
E
[
∥
n
∥
2
2
]
=
E
[
∥
(
x
−
s
t
)
∥
2
2
]
⋅
E
[
t
r
(
A
T
A
)
]
+
M
⋅
σ
2
=
E
[
∥
q
t
∥
2
2
]
⋅
E
[
t
r
(
A
T
A
)
]
+
M
⋅
σ
2
begin{aligned} mathbb E[{Vert {pmb y - pmb A pmb s^t} Vert}^2_2] &= mathbb E[{Vert {pmb A pmb x - pmb A pmb s^t + pmb n} Vert}^2_2] \ &=mathbb E[{Vert {pmb A( pmb x - pmb s^t )} Vert}^2_2] + mathbb E[{Vert {pmb n} Vert}^2_2] \ &= mathbb E[{Vert {( pmb x - pmb s^t )} Vert}^2_2] cdot mathbb E[tr(pmb A^T pmb A)] + M cdot sigma^2 \ &=mathbb E[{Vert pmb q^t Vert}^2_2] cdot mathbb E[tr(pmb A^T pmb A)] + M cdot sigma^2 \ end{aligned}
E[∥yyy−AAAssst∥22]=E[∥AAAxxx−AAAssst+nnn∥22]=E[∥AAA(xxx−ssst)∥22]+E[∥nnn∥22]=E[∥(xxx−ssst)∥22]⋅E[tr(AAATAAA)]+M⋅σ2=E[∥qqqt∥22]⋅E[tr(AAATAAA)]+M⋅σ2
⇒
Rightarrow
⇒
v
^
t
2
=
1
N
E
[
∥
q
t
∥
2
2
]
=
1
N
E
[
∥
y
−
A
s
t
∥
2
2
]
−
M
⋅
σ
2
E
[
t
r
(
A
T
A
)
]
hat {v}^2_t=frac{1}{N} mathbb E[{Vert pmb q^t Vert}^2_2]=frac{1}{N} frac{mathbb E[{Vert {pmb y - pmb A pmb s^t} Vert}^2_2] - M cdot sigma^2}{mathbb E[tr(pmb A^T pmb A)]}
v^t2=N1E[∥qqqt∥22]=N1E[tr(AAATAAA)]E[∥yyy−AAAssst∥22]−M⋅σ2
式(27)跟论文差了一个系数
1
N
frac{1}{N}
N1，感觉(27)会合适一些。
4.6.2 线性均方误差的估计
非线性估计与式(18a)一致
τ
^
t
2
=
t
r
(
B
t
T
B
t
)
⋅
v
^
t
2
+
M
N
t
r
(
W
t
T
W
t
)
⋅
σ
2
(28)
hat {tau}^2_t=tr(pmb B^T_t pmb B_t ) cdot hat {v}^2_t + frac {M}{N} tr(pmb W^T_t pmb W_t ) cdot sigma^2 tag{28}
τ^t2=tr(BBBtTBBBt)⋅v^t2+NMtr(WWWtTWWWt)⋅σ2(28)
强调：如果仿真OAMP迭代过程中需要用到
v
^
t
2
hat {v}^2_t
v^t2和
τ
^
t
2
hat {tau}^2_t
τ^t2，那么就是由(27,28)确定的。
5 总结
如果线性估计中的
{
q
t
,
h
t
}
{ pmb q^t, pmb h^t }
{qqqt,hhht}相互独立，非线性非线性估计中的
{
q
t
+
1
,
h
t
}
{ pmb q^{t+1}, pmb h^t }
{qqqt+1,hhht}相互独立，那么假设1，2就自然而然成立。然而论文只能证明正交性，不能证明独立性，这也是OAMP里边正交的来源，虽然推论1，2弱于假设1，2，但是仿真结果表面OAMP-state evolution还是可靠的。即使对于一般的酉不变矩阵，对奇异值的分布没有严格的束缚，OAMP的性能依然可以被OAMP-state evolution表征，这也是AMP和AMP-state evolution所不能比拟的。因此，相对宽泛的感知矩阵使得OAMP的应用也更加广泛。