Bayes-Optimal Convolutional AMP
Keigo Takeuchi
cs.IT math.IT




















































IEEE transactions on information theory, Vol. , No. ,
Takeuchi: Bayes-Optimal Convolutional AMP












0000–0000/00$00.00  2020 IEEE












Bayes-Optimal Convolutional AMP
    Keigo Takeuchi, Member, IEEE


The author was in part supported by the Grant-in-Aid 
for Scientific Research (B) (JSPS KAKENHI Grant Numbers 18H01441 and 
21H01326), Japan. 
The material in this paper was presented in part at 2019 IEEE 
International Symposium on Information Theory and submitted in part 
to 2021 IEEE International Symposium on Information Theory. 

K. Takeuchi is with the Department of Electrical and Electronic Information Engineering, Toyohashi University of Technology, Toyohashi 441-8580, Japan (e-mail: takeuchi@ee.tut.ac.jp).

    December 30, 2023
============================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================






This paper proposes Bayes-optimal convolutional approximate 
message-passing (CAMP) for signal recovery in compressed sensing. 
CAMP uses the same low-complexity matched filter (MF) for interference 
suppression as approximate message-passing (AMP). To improve the 
convergence property of AMP for ill-conditioned sensing matrices, 
the so-called Onsager correction term in AMP is replaced by a convolution 
of all preceding messages. The tap coefficients in the convolution 
are determined so as to realize asymptotic Gaussianity of estimation 
errors via state evolution (SE) under the assumption of orthogonally 
invariant sensing matrices. An SE equation is derived to optimize the 
sequence of denoisers in CAMP. The optimized CAMP is proved 
to be Bayes-optimal for all orthogonally invariant sensing matrices 
if the SE equation converges to a fixed-point and 
if the fixed-point is unique. 
For sensing matrices with low-to-moderate condition 
numbers, CAMP can achieve the same performance as high-complexity 
orthogonal/vector AMP that requires the linear minimum mean-square error 
(LMMSE) filter instead of the MF.




Compressed sensing, approximate message-passing (AMP), 
orthogonal/vector AMP, convolutional AMP, large system limit, 
state evolution.    
















§ INTRODUCTION



 §.§ Compressed Sensing

Compressed sensing (CS) <cit.> is 
a powerful technique for recovering sparse signals from compressed 
measurements. Under the assumption of linear measurements, CS is formulated 
as estimation of a sparse signal vector x∈ℝ^N 
from a compressed measurement vector y∈ℝ^M 
(M≤ N) and a sensing matrix A∈ℝ^M× N, 
given by 

    y = Ax + w,
   
where w∈ℝ^M is an unknown additive noise vector. 

For simplicity in information-theoretic discussion <cit.>, suppose that 
the signal vector x has independent and identically 
distributed (i.i.d.) elements. Sparsity of signals is measured with the 
Rényi information dimension <cit.> of each signal element. 
When each signal takes non-zero real values with probability ρ∈[0, 1], 
the information dimension is equal to ρ. In the noiseless case 
w=0, Wu and Verdú <cit.> proved that, 
if and only if the compression rate δ=M/N is equal to or larger than 
the information dimension, there are some sensing matrix A and 
method for signal recovery such that the signal vector x can be  
recovered with negligibly small error probability in the large 
system limit, where M and N tend to infinity with the compression rate 
δ kept constant. Thus, an important issue in CS is a construction of 
practical sensing matrices and a low-complexity algorithm for signal recovery 
achieving the information-theoretic compression limit. 

Important examples of sensing matrices are zero-mean i.i.d. 
sensing matrices <cit.> and random sensing matrices with 
orthogonal rows <cit.>. The information-theoretic compression 
limit of zero-mean i.i.d. sensing matrices was analyzed with the non-rigorous 
replica method <cit.>—a tool developed in statistical 
mechanics <cit.>. The compression limit is characterized 
via a potential function called free energy. The results themselves were 
rigorously justified in <cit.> while 
the justification of the replica method is still open. It is a simple 
exercise to prove that the compression limit for zero-mean i.i.d. sensing 
matrices is equal to the Rényi information dimension in the noiseless case, 
by using a relationship between the information dimension and 
mutual information <cit.>.  



Random sensing matrices with orthogonal rows can be constructed 
efficiently in terms of both time and space complexity while zero-mean 
i.i.d. sensing matrices require O(MN) time and memory 
for matrix-vector multiplication.  
When the fast Fourier transform or fast Walsh-Hadamard transform is used, 
the matrix-vector multiplication needs O(Nlog N) time and 
O(N) memory. Thus, random sensing matrices with orthogonal rows 
are preferable from a practical point of view.  
 
The class of orthogonally invariant matrices includes zero-mean i.i.d. 
Gaussian matrices and Haar orthogonal matrices <cit.>, of 
which the latter is regarded as an idealized model of random matrices with 
orthogonal rows. The class allows us to analyze the information-theoretic 
compression limit in signal recovery. 
The replica method <cit.> was used 
to analyze the compression limit for orthogonally invariant sensing matrices. 
The replica results themselves were justified in <cit.>. 
In particular, Haar orthogonal matrices achieve the Welch lower 
bound <cit.> and were proved to be optimal for 
Gaussian <cit.> and general <cit.> signals. 
In the noiseless case, of course, Haar orthogonal sensing matrices achieve 
the compression rate that is equal to the Rényi information dimension. 

In practical systems, the measurement vector is subject not only to 
additive noise but also to multiplicative noise. A typical example is 
fading in wireless communication systems <cit.>. 
The effective sensing matrix containing fading influence may be 
ill-conditioned even if a Haar orthogonal sensing matrix is used. 
Such effective sensing matrices can be modeled as orthogonally invariant 
matrices. Thus, an ultimate algorithm for signal recovery is required to be low 
complexity and Bayes-optimal for all orthogonally invariant sensing matrices. 




 §.§ Message-Passing

A promising solution to signal recovery is message-passing (MP). 
Approximate message-passing (AMP) <cit.> is a low-complexity and 
powerful algorithm for signal recovery from zero-mean i.i.d. sub-Gaussian 
measurements. Bayes-optimal AMP is regarded as an exact large-system 
approximation of loopy belief propagation (BP) <cit.>. The main 
feature of AMP is the so-called Onsager correction to realize asymptotic 
Gaussianity of the estimation errors before denoising. The Onsager 
correction originates from that in the Thouless-Anderson-Palmer (TAP) 
equation <cit.> for a solvable  spin glass model with i.i.d. 
interaction between all spins <cit.>. The Onsager correction 
cancels intractable dependencies of the current estimation error on past 
estimation errors due to i.i.d. dense sensing matrices. 

The convergence property of AMP was analyzed rigorously via state 
evolution (SE) <cit.>, inspired by Bolthausen's conditioning 
technique <cit.>. SE is a dense counterpart of density 
evolution <cit.> in sparse systems. SE tracks a few 
state variables to describe rigorous dynamics of MP in the large system 
limit. SE analysis in <cit.> implies that AMP is 
Bayes-optimal for zero-mean i.i.d. sub-Gaussian sensing matrices when the 
compression rate δ is larger than a certain value called BP 
threshold <cit.>. Spatial 
coupling <cit.> is needed to 
realize the optimality of AMP for any compression rate. However, this paper 
does not consider spatial coupling since spatial coupling is a universal 
technique <cit.> to improve the performance of MP. 

A disadvantage of AMP is that AMP fails to converge when the sensing matrix 
is non-zero mean <cit.> or ill-conditioned <cit.>.  
To solve this issue, orthogonal AMP (OAMP) <cit.> and vector 
AMP <cit.> were proposed. The two MP algorithms are 
equivalent to each other. Bayes-optimal OAMP/VAMP can be regarded as an exact 
large-system approximation of expectation propagation 
(EP) <cit.>. 
Rigorous SE analysis <cit.>   
proved that OAMP/VAMP is Bayes-optimal for orthogonally invariant sensing 
matrices when the compression rate is larger than BP threshold. While 
non-zero mean matrices are outside the class of orthogonally invariant 
matrices, numerical simulations in <cit.> indicated that OAMP/VAMP 
can treat the non-zero mean case. 

A prototype of OAMP/VAMP was originally proposed by Opper and 
Winther <cit.>. Historically, they <cit.> 
generalized the Onsager correction in the TAP equation <cit.> from   
zero-mean i.i.d. spin interaction to orthogonally invariant interaction.  
Their method was formulated as the expectation-consistency (EC) 
approximation <cit.>. The EC approximation itself does not produce MP 
algorithms but a potential function of which a local minimum should be solved 
with some MP algorithm. OAMP/VAMP can be derived from an EP-type 
iteration–called a single loop algorithm <cit.>—to solve a local 
minimum of the EC potential. See <cit.> 
for the derivation of OAMP/VAMP via the EC approximation. 
 
The main weakness of OAMP/VAMP is a per-iteration requirement of the linear 
minimum mean-square error (LMMSE) filter, of which the time complexity is 
O(M^3+M^2N) per iteration. The singular-value decomposition (SVD) 
of the sensing matrix allows us to circumvent the use of the LMMSE 
filter <cit.>. However, the complexity of the SVD itself is high 
in general. The performance of OAMP/VAMP degrades 
significantly when the LMMSE filter is replaced by the low-complexity 
matched filter (MF) <cit.> used in AMP. Thus, OAMP/VAMP can be applied 
only to limited problems in which the SVD of the sensing matrix is 
computed efficiently.  

In summary, it is still open to construct a low-complexity and Bayes-optimal 
MP algorithm for all orthogonally invariant sensing matrices. The purpose of 
this paper is to tackle the design issue of such ultimate MP algorithms. 



 §.§ Methodology

The main idea of this paper is to extend the class of MP algorithms.  
Conventional MP algorithms use update rules that depend only on messages 
in the latest iteration. Long-memory MP algorithms considered in this paper 
are allowed to depend on messages in all preceding iterations. 

This class of long-memory MP algorithms was motivated by SE  
analysis of AMP for orthogonally invariant sensing matrices <cit.>. 
When the asymptotic singular-value distribution of the sensing matrix is equal 
to that of zero-mean i.i.d. Gaussian matrices, the error model of AMP was 
proved to be an instance of a general error model <cit.>, 
in which each error depends on errors in all preceding iterations.   
This result implies that the Onsager correction in AMP uses messages in all 
preceding iterations to realize the asymptotic Gaussianity of the current 
estimation error while the representation itself of the correction term looks 
as if only messages in the latest iteration are utilized. Inspired by this 
observation, we consider long-memory MP algorithms as a starting point.  

The idea of long-memory MP was originally proposed in Opper, Çakmak, 
and Winther's paper <cit.> to solve the TAP 
equations for spin glass models with orthogonally invariant interaction.    
Their methodology was based on non-rigorous dynamical functional theory. 
After the initial submission of this paper, their results were rigorously 
justified via SE in <cit.>. 

The proposed design of long-memory MP consists of three steps: A first step 
is an establishment of rigorous SE for analyzing the dynamics of 
long-memory MP algorithms for orthogonally invariant sensing 
matrices. This step has been already established in <cit.> 
by generalizing conventional SE analysis <cit.> to 
the long-memory case. The SE analysis provides a sufficient condition for 
a long-memory MP algorithm to have Gaussian-distributed estimation errors 
in the large system limit.  
The main advantage in the SE analysis is to provide a systematic design 
of long-memory MP that satisfies the asymptotic Gaussianity in estimation 
errors while the class of long-memory MP is slightly smaller than in 
<cit.>.    

A second step is to modify the Onsager correction in AMP so as to satisfy the 
sufficient condition for the asymptotic Gaussianity. A solvable class of 
long-memory MP was proposed in <cit.>, where the Onsager 
correction was defined as a convolution of messages in all preceding 
iterations. The tap coefficients in the convolution were determined so as to 
satisfy the sufficient condition. Thus, long-memory MP proposed 
in <cit.> was called convolutional AMP (CAMP) and 
is the main object of this paper. 

This paper generalizes CAMP in <cit.>, motivated by an 
implementation of OAMP/VAMP based on conjugate gradient 
(CG) <cit.>. OAMP/VAMP applies the LMMSE filter to a message 
z∈ℝ^M 
after interference subtraction. The LMMSE filter is decomposed into a 
noise-whitening filter and MF. In principle, CG approximates the output of 
the noise-whitening filter with a vector in the Krylov subspace spanned by 
{z, AA^Tz, 
(AA^T)^2z, 
…} i.e. a finite weighted sum of 
{(AA^T)^jz}. 
On the other hand, messages in the original CAMP <cit.> are in 
the 0th Krylov subspace {αz: α∈ℝ} 
since only the MF is used. To fill this gap, we generalize a convolution 
of all preceding messages in the original CAMP <cit.> to that 
of affine transforms of the preceding messages. 

The last step is to optimize the sequence of denoisers in CAMP. 
The last step is a new contribution of this paper, submitted 
to <cit.>. The optimization requires information on 
the distribution of the estimation errors before denoising in each 
iteration. Since the estimation errors  are asymptotically 
Gaussian-distributed, we need to track the dynamics of the variance of 
the estimation errors. To analyze this dynamics, we utilize the SE analysis 
established in the first step. 




 §.§ Contributions

The contributions of this paper are sixfold: A first contribution 
(Theorem <ref> in Section <ref>) is to propose a general error 
model for long-memory MP and prove the asymptotic Gaussianity of estimation 
errors in the general error model via rigorous SE under the assumption of 
orthogonally invariant sensing matrices. The general error model contains 
both error models of AMP and OAMP/VAMP.  

A second contribution (Section <ref>) is the addition of a convolution 
proportional to AA^T to the Onsager 
correction in <cit.>, according to the above-mentioned 
argument on the Krylov subspace. This addition improves the convergence 
property of CAMP. 

A third contribution (Theorem <ref> in Section <ref>) is 
to design tap coefficients in the convolution so as to guarantee the 
asymptotic Gaussianity of estimation errors for all orthogonally invariant 
sensing matrices. Part of the tap coefficients are used to realize the 
asymptotic Gaussianity. The remaining coefficients can be utilized to improve 
the convergence property of CAMP.   

A fourth contribution (Theorem <ref> in Section <ref>) 
is to present the designed tap coefficients in closed-form. 
This closed-form representation circumvents numerical instability in solving 
the tap coefficients numerically. The third and fourth contributions are 
based on the same proof strategy as in <cit.>. 
 
A fifth contribution (Theorems <ref> and 
<ref> in Section <ref>) 
is to optimize the sequence of denoisers in CAMP. An SE equation is derived to 
describe the dynamics of the variance of the estimation errors before 
denoising in CAMP. The SE equation is a two-dimensional 
nonlinear difference equation. By analyzing the fixed-point of the SE equation, 
we prove that optimized CAMP is Bayes-optimal for all orthogonally invariant 
sensing matrices if the SE equation converges to a fixed-point and 
if the fixed-point is unique. 

The last contribution (Section <ref>) is numerical evaluation of CAMP. 
The remaining parameters in the Bayes-optimal CAMP are optimized numerically 
to improve the convergence property. Numerical simulations show that the CAMP 
can converge for sensing matrices with larger condition numbers than the 
original CAMP <cit.> when the design parameters are optimized.  
The CAMP can achieve the same performance as OAMP/VAMP for sensing matrices 
with low-to-moderate condition numbers while it is inferior to OAMP/VAMP for 
high condition numbers. 



 §.§ Organization

The remainder of this paper is organized as follows: After summarizing the 
notation used in this paper, we present a unified SE framework for analyzing 
long-memory MP under the assumption of orthogonally invariant sensing 
matrices in Section <ref>. This section corresponds to the first step 
for proposing Bayes-optimal CAMP. 

In Section <ref>, we propose CAMP with design parameters. 
This section corresponds to the remaining two steps for establishing 
Bayes-optimal CAMP. The proposed 
CAMP is more general than in <cit.>. We utilize the SE framework 
established in Section <ref> to determine the tap coefficients in 
CAMP that guarantee the asymptotic Gaussianity of estimation errors. To 
design the remaining design parameters, we derive an SE equation to optimize 
the performance of signal recovery. 

Section <ref> presents numerical results. The remaining design parameters 
in CAMP are optimized via numerical simulations. 
The optimized CAMP is compared to conventional AMP and OAMP/VAMP via the 
SE equation and numerical simulations. Section <ref> concludes this 
paper. The details for the proofs of the main theorems are presented in  
appendices. 



 §.§ Notation

For a matrix M, the transpose of M is 
denoted by M^T. The notation 
Tr(A) represents the trace of a square matrix 
A. For a symmetric matrix A, the minimum 
eigenvalue of A is written as 
λ_min(A). 
The notation 𝒪_M× N denotes 
the space of all possible M× N matrices with orthonormal columns for 
M≥ N and orthonormal rows for M<N. In particular, 
𝒪_N× N reduces to the space 𝒪_N of 
all possible N× N orthogonal matrices. 

For a vector v, 
the notation diag(v) denotes the diagonal matrix 
of which the nth diagonal element is equal to v_n=[v]_n. 
The norm v=√(v^Tv) 
represents the Euclidean norm. 
For a matrix M_i with an index i, the tth column 
of M_i is denoted by m_i,t. Furthermore, 
we write the nth element of m_i,t as m_i,t,n. 

The Kronecker delta is denoted by δ_τ,t while the Dirac delta 
function is represented as δ(·). We write the Gaussian 
distribution with mean μ and covariance Σ 
as 𝒩(μ,Σ). The notations 
 and  denote almost sure convergence and equivalence, 
respectively.   

We use the notational convention ∑_t=t_1^t_2⋯=0 and 
∏_t=t_1^t_2⋯=1 for t_1>t_2. 
For any multivariate function ϕ:ℝ^t→ℝ, the 
notation ∂_t'ϕ for t'=0,…,t-1 denotes the partial 
derivative of ϕ with respect to the t'th variable x_t', 

    ∂_t'ϕ = ∂ϕ/∂ x_t'
    (x_0,…,x_t-1).

For any vector v∈ℝ^N, 
the notation ⟨v⟩=N^-1∑_n=1^Nv_n represents 
the arithmetic mean of the elements. 
For any scalar function f:∈ℝ→ℝ, the notation 
f(v) means the element-wise application of f to a vector 
v, i.e. [f(v)]_n=f(v_n). 

For a sequence {p_t}_t=0^∞, we define the Z-transform of 
{p_t} as 

    P(z) = ∑_t=0^∞p_tz^-t.

For two sequences {p_t, q_t}_t=0^∞, 
we define the convolution operator * as 

    p_t+i*q_t+j = ∑_τ=0^tp_τ+iq_t-τ+j,

with p_t=0 and q_t=0 for t<0. For finite-length sequences 
{p_t}_t=0^T of length T+1, we transform them into infinite-length 
sequences by adding p_t=0 and q_t=0 for all t>T. 

For two arrays {a_t',t, b_t',t: t', t=0,…, ∞}, 
we write the two-dimensional convolution as 

    a_t'+i,t+j*b_t'+k,t+l 
    = ∑_τ'=0^t'∑_τ=0^ta_τ'+i,τ+jb_t'-τ'+k,t-τ+l,

where a_t',t=0 and b_t',t=0 are defined for t'<0 or t<0. 

Whether a convolution is one-dimensional can be distinguished as follows: 
A convolution is one-dimensional, such as a_t+i*b_t+j, when both operands 
contain only one identical subscript. On the other hand, a convolution is 
two-dimensional, such as (a_t'a_t+i)*b_t'+j,t, when both operands 
include two identical subscripts. 
 




§ UNIFIED FRAMEWORK
 


 §.§ Definitions and Assumptions

We define the statistical properties of the random variables in the 
measurement model (<ref>). 
The performance of MP is commonly measured in terms of the mean-square error 
(MSE). Nonetheless, we follow <cit.> to consider 
a general performance measure in terms of separable and pseudo-Lipschitz 
functions while we assume the separability and Lipschitz-continuity 
for denoisers. 


A vector-valued function f=(f_1,…,f_N)^T
:ℝ^N× t→ℝ^N is said to be 
separable if [f(x_1,…,x_t)]_n
=f_n(x_1,n,…,x_t,n) holds for all 
x_i∈ℝ^N. 
  


A function f:ℝ^t→ℝ is said to be pseudo-Lipschitz 
of order k <cit.> if there are some Lipschitz constant L>0 and 
some order k∈ℕ such that 
for all x∈ℝ^t and y∈ℝ^t 

    |f(x) - f(y)| 
    ≤ L(1 + x^k-1 + y^k-1)
    x - y.



By definition, any pseudo-Lipschitz function of order k=1 is 
Lipschitz-continuous. A vector-valued function 
f=(f_1,…,f_N)^T is pseudo-Lipschitz if all 
element functions {f_n} are pseudo-Lipschitz. 


A separable pseudo-Lipschitz function 
f:ℝ^N× t→ℝ^N is said to be proper 
if the Lipschitz constant L_n>0 of the nth function f_n satisfies 

    lim sup_N→∞1/N∑_n=1^NL_n^j<∞

for any j∈ℕ. 


A proper pseudo-Lipschitz function allows us apply a proof strategy for 
pseudo-Lipschitz functions with n-independent Lipschitz constant L_n=L 
to the n-dependent case straightforwardly.    
The space of all possible separable and proper pseudo-Lipschitz functions of 
order k is denoted by 𝒫ℒ(k). We have the inclusion relation 
𝒫ℒ(k)⊂𝒫ℒ(k') for all k<k' since 
x^k≤x^k' holds for 
x≫1. 












We assume statistical properties of the signal vector associated 
with separable and proper pseudo-Lipschitz functions of order k≥ 2. 
Note that the integer k in the following assumptions is an identical 
parameter that is equal to the order of separable and proper 
pseudo-Lipschitz functions used in SE to measure the performance of MP. 
If the MSE is considered, the integer k is set to 2.  

 
The signal vector x satisfies the following strong law of 
large numbers:   

    ⟨f(x)⟩
    - 𝔼[
     ⟨f(x)⟩]0

as N→∞ for any separable and proper pseudo-Lipschitz function 
f:ℝ^N→ℝ^N of order k≥2. Furthermore, 
x has zero-mean and bounded (2k-2+ϵ)th moments 
for some ϵ>0.


Assumption <ref> follows from the classical strong law of large 
numbers when x has i.i.d. elements. 


An orthogonal matrix V∈𝒪_N is said to be 
Haar-distributed <cit.> if V is orthogonally invariant, 
i.e. V∼ΦVΨ 
for all orthogonal matrices Φ, 
Ψ∈𝒪_N independent of V. 

  
The sensing matrix A is right-orthogonally invariant, i.e. 
A∼AΨ for any orthogonal 
matrix Ψ∈𝒪_N independent of A. 
More precisely, the orthogonal matrix V∈𝒪_N 
in the SVD A=UΣV^T is Haar-distributed and independent of 
UΣ. Furthermore, the empirical eigenvalue 
distribution of A^TA converges almost 
surely to a compactly supported deterministic distribution with unit first 
moment in the large system limit. 


The assumption of unit first moment implies the almost sure convergence 
N^-1Tr(A^TA)1 
in the large system limit. Assumption <ref> holds when 
A has zero-mean i.i.d. Gaussian elements with variance 
M^-1. As shown in SE, the asymptotic Gaussianity of 
estimation errors in MP depends heavily on the Haar assumption of 
V. Intuitively, the orthogonal transform Va of a vector a∈ℝ^N is distributed as 
N^-1/2az in which 
z∼𝒩(0,I_N) is a 
standard Gaussian vector and independent of a. When  
the amplitude N^-1/2a tends to a constant as N→∞, 
the vector Va looks like a Gaussian vector. 
This is a rough intuition on the asymptotic Gaussianity of estimation 
errors. 

 
The noise vector w is orthogonally 
invariant, i.e. w∼Φw for 
any orthogonal matrix Φ∈𝒪_M independent of 
w. Furthermore, w has zero-mean,  
lim_M→∞M^-1w^2σ^2>0, and 
bounded (2k-2+ϵ)th moments for some ϵ>0.


Assumption <ref> holds when w∼𝒩(0,σ^2I_M) is an additive white 
Gaussian noise (AWGN) vector. It holds for 
U^Tw when the sensing matrix 
A is 
left-orthogonally invariant, i.e. A∼ΦA for any orthogonal matrix Φ∈𝒪_M 
independent of A.  




 §.§ General Error Model
 
We propose a unified framework of SE for analyzing MP algorithms 
that have asymptotically Gaussian-distributed estimation errors for 
orthogonally invariant sensing matrices. Instead of starting with concrete 
MP algorithms, we consider a general class of error models. The proposed class 
does not necessarily contain the error models of all possible long-memory 
MP algorithms. However, it is a natural class of error models that allows 
us to prove the asymptotic Gaussianity of estimation errors for orthogonally 
invariant sensing matrices via a generalization of conventional 
SE <cit.>.  

Let h_t∈ℝ^N and 
q_t+1∈ℝ^N denote error vectors in iteration t 
before and after denoising, respectively. We assume that the error vectors 
are recursively given by 

    b_t = 
    V^Tq̃_t,   q̃_t 
    = q_t - ∑_t'=0^t-1⟨∂_t'ψ_t-1⟩h_t',


    m_t 
    = ϕ_t(B_t+1, w̃; 
    λ),
 

    h_t 
    = Vm̃_t, 
      m̃_t 
    = m_t - ∑_t'=0^t⟨∂_t'ϕ_t⟩b_t',


    q_t+1 
    = ψ_t
    (H_t+1, x),

with q_0=-x. In (<ref>), 
the orthogonal matrix V∈𝒪_N consists of the 
right-singular vectors in the SVD 
A=UΣV^T, with U∈𝒪_M. 
In (<ref>) and (<ref>), we have defined B_t+1
=(b_0,…,b_t) and 
H_t+1=(h_0,…,h_t). 
Furthermore, λ∈ℝ^N is the vector of 
eigenvalues of A^TA. 
The vector w̃∈ℝ^N is given by  

    w̃
    =[ U^Tw;    0 ],
 
where w is the additive noise vector in (<ref>). 

The vector-valued functions 
ϕ_t:ℝ^N×(t+3)→ℝ^N and 
ψ_t:ℝ^N×(t+2)→ℝ^N are assumed 
to be separable, nonlinear, and proper Lipschitz-continuous. 
 
The functions ϕ_t and ψ_t are 
separable. The nonlinearities ϕ_t≠∑_t'=0^tD_t'b_t' and 
ψ_t≠∑_t'=0^tD̃_̃t̃'̃h_t' 
hold for all diagonal matrices {D_t', 
D̃_t'}. The function ϕ_t is 
Lipschitz-continuous with respect to the first t+2 variables and proper 
while ψ_t is proper Lipschitz-continuous with respect to 
all variables.  


It might be possible to relax Assumption <ref> to 
the non-separable case <cit.>. For simplicity, 
however, this paper postulates separable denoisers. 
The nonlinearity is a technical condition for circumventing the zero norm  
N^-1q̃_t^2=0 or 
N^-1m̃_t^2=0, which implies error-free 
estimation N^-1b_t^2=0 or 
N^-1h_t^2=0.  

By definition, the nth function ϕ_t,n has a λ_n-dependent 
Lipschitz constant L_n=L_n(λ_n). Thus, the proper assumption for 
ϕ_t may be regarded as a 
condition on the asymptotic eigenvalue distribution of 
A^TA, as well as a condition on the denoiser 
ϕ_t. For example, ϕ_t is proper 
when the asymptotic eigenvalue distribution has a compact support and when 
the Lipschitz constant L_n(λ_n) itself is a pseudo-Lipschitz 
function of λ_n.  

The main feature of the general error model is in the definitions of 
q̃_t and m̃_t. The 
second terms on the right-hand sides (RHSs) of (<ref>) and (<ref>) 
are correction terms to realize the asymptotic Gaussianity of 
{b_t} and {h_t}. The correction terms 
are a modification of conventional correction that allows us to prove  
the asymptotic Gaussianity via a natural generalization <cit.> 
of Stein's lemma used in conventional SE <cit.>. See 
Lemma <ref> in Appendix <ref> for the details.  

The following examples imply that the general error model (<ref>)–(<ref>) 
contains those of OAMP/VAMP and AMP. 


Consider OAMP/VAMP <cit.> with 
a sequence of scalar denoisers f_t:ℝ→ℝ: 

    x_A→B,t = x_B→A,t 
    + γ_tA^TW_t^-1
    (y - Ax_B→A,t),


    v_A→B,t = γ_t - v_B→A,t,


    W_t 
    = σ^2I_M 
    + v_B→A,tAA^T,


    γ_t^-1 
    = 1/NTr(
     W_t^-1AA^T),


    x_B→A,t+1 
    = v_B→A,t+1(
     f_t(x_A→B,t)/ξ_t v_A→B,t 
     - x_A→B,t/v_A→B,t),


    1/v_B→A,t+1 
    = 1/ξ_tv_A→B,t 
    - 1/v_A→B,t,

with ξ_t =⟨ f_t'(x_A→B^t)⟩.  

It is an exercise to prove that the error model of the OAMP/VAMP is an 
instance of the general error model with

    [ϕ_t(b_t, w̃; 
    λ)]_n 
    = b_t,n 
    - γ_tλ_nb_t,n
    -γ_t√(λ_n)w̃_n/σ^2+v_B→A,tλ_n,


    ψ_t(h_t, x) 
    = f_t(x+h_t) - x/1 - ξ_t,

by using the fact that ξ_t converges almost surely to a constant in the 
large system limit <cit.>.  
The two separable functions ψ_t and ϕ_t 
for the OAMP/VAMP depend only on the vectors b_t and 
h_t in the latest iteration. 



Consider AMP <cit.> with 
a sequence of scalar denoisers f_t:ℝ→ℝ: 

    x_t+1 
    = f_t(x_t + A^Tz_t),


    z_t 
    = y - Ax_t 
    + ξ_t-1/δz_t-1.


Suppose that the empirical eigenvalue distribution of A^TA is equal to that for zero-mean i.i.d. Gaussian matrix 
A in the large system limit. Then, 
the error model of the AMP was proved in <cit.> to be an 
instance of the general error model with 
rl
ϕ_t
   = (I_N-Λ)b_t 
- ξ_t-1/δb_t-1 
+ diag({√(λ_n)})w̃ 


+    ξ_t-1{
 (
  1+1/δ
 )I_N - Λ
}ϕ_t-1
- ξ_t-1ξ_t-2/δϕ_t-2, 
 

    ψ_t(h_t, x) 
    = f_t(x+h_t) - x,
 
with Λ=diag(λ) and 
ξ_t =⟨ f_t'(x+h_t)⟩.  
Note that ϕ_t is a function of B_t+1 while 
ψ_t is a function of h_t. 





 §.§ State Evolution

A rigorous SE result for the general error model (<ref>)–(<ref>) is 
presented in the large system limit. 

 
Suppose that Assumptions <ref>–<ref> 
hold. Then, the following properties hold for all t=0,… and 
t'=0,…, t in the large system limit: 


  *  
The inner products 
N^-1m̃_t^Tm̃_t' and 
N^-1q̃_t^Tq̃_t' converge 
almost surely to some constants π_t,t'∈ℝ and 
κ_t,t'∈ℝ, respectively. 


  * Suppose that 
ψ̃_t(H_t+1,x):
ℝ^N×(t+2)→ℝ^N is a separable and proper  
pseudo-Lipschitz function of order k, 
that ϕ̃_t(B_t+1, 
w̃; λ):
ℝ^N×(t+3)→ℝ^N is separable, pseudo-Lipschitz 
of order k with respect to the first t+2 variables, and proper, and that 
Z_t+1=(z_0,…,z_t)
∈ℝ^N×(t+1) denotes a zero-mean Gaussian random matrix with 
covariance 𝔼[z_τz_τ'^T]
=π_τ,τ'I_N for all τ, τ'=0,…,t, while 
a zero-mean Gaussian random matrix 
Z̃_t+1=(z̃_0,…,
z̃_t)∈ℝ^N×(t+1) has covariance 
𝔼[z̃_τz̃_τ'^T]
=κ_τ,τ'I_N. Then, 

    ⟨ψ̃_t(H_t+1, x)⟩
    - 𝔼[
     ⟨ψ̃_t(Z_t+1, 
     x)⟩] 0,


    ⟨ϕ̃_t(B_t+1, 
    w̃;λ)⟩
    - 𝔼[
     ⟨ϕ̃_t(Z̃_t+1, 
     w̃;λ)⟩] 0.

In evaluating the expectation in (<ref>), 
U^Tw in (<ref>) 
follows the zero-mean Gaussian distribution with covariance 
σ^2I_M.  
In particular, for k=1   

    ⟨∂_t'ψ̃_t(H_t+1, 
    x)⟩
    - 𝔼[
     ⟨∂_t'ψ̃_t(Z_t+1, 
     x)⟩] 0,


    ⟨∂_t'ϕ̃_t(B_t+1, 
    w̃;λ)⟩
    - 𝔼[
     ⟨∂_t'ϕ̃_t(Z̃_t+1, 
     w̃;λ)⟩] 0.



  * Suppose that 
ψ̃_t(H_t+1,x):
ℝ^N×(t+2)→ℝ^N is separable and 
proper Lipschitz-continuous, and that  
ϕ̃_t(B_t+1, 
w̃;λ):
ℝ^N×(t+3)→ℝ^N is separable, 
Lipschitz-continuous with respect to the first t+2 variables, and proper. 
Then,   

    1/Nh_t'^T(
     ψ̃_t
     - ∑_τ=0^t⟨∂_τψ̃_t⟩h_τ)
    0,


    1/Nb_t'^T(
     ϕ̃_t
     - ∑_τ=0^t⟨∂_τϕ̃_t⟩b_τ)
    0.




See Appendix <ref>. 


Properties (<ref>) and (<ref>) are used to evaluate 
the performance of MP by specifying the functions 
ψ̃_t and ϕ̃_t according 
to a performance measure. An important observation is the asymptotic 
Gaussianity of 
H_t+1 and B_t+1. In evaluating the performance 
of MP, we can replace them with tractable Gaussian random matrices 
Z_t+1 and Z̃_t+1.  

The asymptotic Gaussianity originates from the definitions of 
q̃_t and m̃_t in 
(<ref>) and (<ref>). Properties (<ref>) and 
(<ref>) imply the asymptotic orthogonality 
N^-1h_t'^Tq̃_t+10 
and N^-1b_t'^Tm̃_t0.  
This orthogonality is used to prove that the distributions of 
H_t+1 and B_t+1 are asymptotically Gaussian.  

Properties (<ref>) and (<ref>) can be 
regarded as computation formulas to evaluate 
N^-1h_t'^Tψ̃_t and 
N^-1b_t'^Tϕ̃_t. They 
can be computed via linear combinations of 
{N^-1h_t'^Th_τ}_τ=0^t  
and {N^-1b_t'^Tb_τ}_τ=0^t.   
In particular, (<ref>), (<ref>), and Property <ref>) in 
Theorem <ref> imply 
N^-1h_t'^Th_τπ_t',τ and 
N^-1b_t'^Tb_τκ_t',τ. 
Furthermore, the coefficients in the linear combinations can be computed 
with (<ref>) and (<ref>). From these 
observations, the SE equations of the general error model are given as 
dynamical systems with respect to {π_t,t', κ_t,t'} in general.  

We do not derive SE equations with respect to {π_t,t', 
κ_t,t'} in a general form. Instead, we derive SE equations 
after specifying MP. The usefulness of Theorem <ref> is clarified 
in deriving SE equations. 



§ SIGNAL RECOVERY
 


 §.§ Convolutional Approximate Message-Passing
 
Let x_t∈ℝ^N denote an estimator of 
the signal vector x in iteration t. CAMP computes the 
estimator x_t recursively as 

    x_t+1 
    = f_t(x_t + A^Tz_t),


    z_t 
    = y - Ax_t 
    + ∑_τ=0^t-1ξ_τ^(t-1)(
    θ_t-τAA^T
    - g_t-τI_M)z_τ,

with the initial condition x_0=0, 
where ξ_τ^(t-1)=∏_t'=τ^(t-1)ξ_t' is the product of 
{ξ_t'} given by 

    ξ_t = 
    ⟨
     f_t'(x_t + A^Tz_t) 
    ⟩.


In (<ref>) and (<ref>), A and y are the 
sensing matrix and the measurement vector in (<ref>), respectively. 
The functions {f_t: ℝ→ℝ} are a sequence of 
Lipschitz-continuous denoisers. 
The tap coefficients {g_τ∈ℝ} and 
{θ_τ∈ℝ} in the convolution are design parameters. 
The parameters {θ_τ} are optimized to improve the performance of 
the CAMP while {g_τ} are determined so as to realize the asymptotic 
Gaussianity of the estimation errors via Theorem <ref>. 

To impose the initial condition x_0=0, 
it is convenient to introduce the notational convention 
f_-1(·)=0, which is used throughout this paper. 


The CAMP is a generalization of AMP <cit.> and reduces to AMP when 
g_1=-δ^-1, g_τ=0 for τ>1, and θ_τ=0 hold. 
Also, as a generalization of CAMP in <cit.>, the affine transform 
(θ_t-τAA^T
-g_t-τI_M)
z_τ has been applied before the convolution. 
Nonetheless, the proposed MP is called CAMP simply. 
In particular, the MP algorithm reduces to the original 
CAMP <cit.> when θ_τ=0 is assumed. 


 
The design parameters {θ_τ} are not required and can be set to 
zero for sensing matrices with identical non-zero singular values since 
AA^T reduces 
to the identity matrix with the exception of a constant factor. 
Thus, non-zero parameters {θ_τ} should be introduced only for the 
case of non-identical singular values. 





 §.§ Error Model
 
To design the parameters g_τ and θ_τ via 
Theorem <ref>, we derive an error model of the CAMP. 
Let h_t=x_t 
+ A^Tz_t-x and 
q_t+1=x_t+1-x denote the error 
vectors before and after denoising f_t, respectively. Then, we have 

    q_t+1 = f_t(x + h_t) 
    - x≡ψ_t(h_t, x),


    q̃_t+1=q_t+1-ξ_th_t.

Using the notational convention f_-1(·)=0, we obtain the initial 
condition q_0=-x imposed in the general 
error model.  

We define m_t=V^Th_t  
and b_t=V^Tq̃_t 
to formulate the error model of the CAMP in a form corresponding to 
the general error model (<ref>)–(<ref>). Substituting the definition  
h_t=x_t 
+ A^Tz_t-x into 
m_t=V^Th_t yields 

    m_t 
    = V^Tq_t 
    + Σ^TU^Tz_t,

where we have used the definition 
q_t=x_t-x and 
the SVD A=UΣV^T. We utilize the definitions (<ref>), 
b_t=V^Tq̃_t, and 
m_t=V^Th_t to obtain  

    V^Tq_t  
    = b_t + ξ_t-1m_t-1.

Combining these two equations yields 

    Σ^TU^Tz_t
    = m_t 
    - b_t - ξ_t-1m_t-1.


To obtain a closed-form equation with respect to m_t, 
we left-multiply (<ref>) by Σ^TU^T and use (<ref>) to have   
rl
Σ^TU^Tz_t 
   = - ΛV^Tq_t 
+ Σ^TU^Tw 


+ ∑_τ=0^t-1   ξ_τ^(t-1)(θ_t-τΛ
- g_t-τI_M)
Σ^TU^Tz_τ, 

with Λ=
Σ^TΣ. Substituting 
(<ref>) and 
(<ref>) into this expression, we arrive at 
rl
m_t 
= (I_N - Λ)   
(b_t + ξ_t-1m_t-1) 
+ Σ^TU^Tw 


+ ∑_τ=0^t-1   ξ_τ^(t-1)(θ_t-τΛ 
- g_t-τI_M) 

   ·(m_τ 
- b_τ - ξ_τ-1m_τ-1), 
  

where any vector with a negative index is set to zero. 
This expression implies that ϕ_t for the CAMP depends on 
all messages B_t+1. 

We note that Assumption <ref> holds under 
Assumption <ref> since the denoiser f_t has been 
assumed to be Lipschitz-continuous. 



 §.§ Asymptotic Gaussianity
 
We compare the obtained error model with the general error 
model (<ref>)–(<ref>). The only difference is in (<ref>): 
The correction m̃_t of m_t is 
used to define h_t in the general error model while  
no correction is performed in the error model of the CAMP. 
Thus, the general error model contains the error model of the CAMP when 
⟨∂_t'm_t⟩=0 holds for all 
t'=0,…,t. In the CAMP, the parameters {g_τ} are determined 
so as to guarantee ⟨∂_t'm_t⟩=0 
in the large system limit.  

Let μ_j denote the jth moment of the asymptotic eigenvalue 
distribution of A^TA, given by 

    μ_j = lim_M=δ N→∞1/NTr(Λ^j).

Assumption <ref> implies μ_1=1. 
We define a coupled dynamical system {g_τ^(j)} determined via 
the tap coefficients {g_τ} and {θ_τ} as  

    g_0^(j) = μ_j+1 - μ_j,

rl
g_1^(j)
=    g_0^(j) - g_0^(j+1)
- g_1(g_0^(j) + μ_j)  


   + θ_1(g_0^(j+1) + μ_j+1),   


rl
g_τ^(j)
=    g_τ-1^(j) - g_τ-1^(j+1) 
- g_τμ_j + θ_τμ_j+1


   + ∑_τ'=0^τ-1(θ_τ-τ'g_τ'^(j+1) 
- g_τ-τ'g_τ'^(j))


   - ∑_τ'=1^τ-1(θ_τ-τ'g_τ'-1^(j+1) 
- g_τ-τ'g_τ'-1^(j))  

for τ>1.  

  
Suppose that Assumptions <ref>–<ref> 
hold, that the denoiser f_t is Lipschitz-continuous, and that 
the tap coefficients {g_τ} and {θ_τ} in the CAMP satisfy  

    g_1 = θ_1(g_0^(1) + 1) - g_0^(1) ,


    g_τ
    = θ_τ - g_τ-1^(1) 
    + ∑_τ'=0^τ-1θ_τ-τ'g_τ'^(1)
    - ∑_τ'=1^τ-1θ_τ-τ'g_τ'-1^(1),

where {g_τ^(1)} is governed by the dynamical 
system (<ref>)–(<ref>). 
Then, ⟨∂_t'm_t⟩→0 holds in the large 
system limit, i.e. the error model of the CAMP is included into the general 
error model. 


Let 

    g_t',t^(j) 
    = -lim_M=δ N→∞⟨Λ^j∂_t'm_t⟩.

It is sufficient to prove g_t',t^(j)ξ_t'^(t-1)g_t-t'^(j)+o(1) 
and g_τ^(0)=0 under the notational convention ξ_t'^(t)=1 
for t'>t. The latter property g_τ^(0)=0 follows from 
(<ref>) for τ=0, (<ref>) and 
(<ref>) for τ=1, and from (<ref>) and (<ref>). 
See Appendix <ref> for the proof of the former property. 


Throughout this paper, we assume that the tap coefficients {g_τ} and 
{θ_τ} satisfy (<ref>) and (<ref>). Thus, 
Theorem <ref> implies that  
the asymptotic Gaussianity is guaranteed for the CAMP. In principle, it is 
possible to compute the tap coefficients by solving the coupled dynamical 
system (<ref>)–(<ref>) numerically for a given moment sequence 
{μ_j}. However, numerical evaluation indicated that the 
dynamical system is unstable against numerical errors when the moment 
sequence {μ_j} is a diverging sequence. Thus, we need 
a closed-form solution to the tap coefficients.  

To present the closed-form solution, we define the η-transform of the 
asymptotic eigenvalue distribution of 
A^TA <cit.> as 

    η(x) = lim_M=δ N→∞1/NTr{(
      I_N + xA^TA)^-1}.

By definition, we have the power-series expansion 

    η(x) = lim_M=δ N→∞1/N∑_n=1^N1/1+xλ_n
    = ∑_j=0^∞μ_j(-x)^j

for |x|<1/max{λ_n}. 
Let G(z) denote the generating function of the tap coefficients 
{g_τ} given by  

    G(z) = ∑_τ=0^∞g_τz^-τ, 
       g_0=1.

Similarly, we write the generating function of 
{θ_τ} with θ_0=1 as Θ(z). 

 
Suppose that the tap coefficients {g_τ} and {θ_τ} 
satisfy (<ref>) and (<ref>). Then, 
the generating functions G(z) and Θ(z) of   
{g_τ} and {θ_τ} satisfy 

    η(
     1 - (1-z^-1)Θ(z)/(1-z^-1)G(z)) 
    = (1-z^-1)Θ(z),

where η denotes the η-transform of the asymptotic eigenvalue 
distribution of A^TA. 


See Appendix <ref>. 


Suppose that the η-transform is given. Since the η-transform has 
the inverse function, from Theorem <ref> we have 
(1-z^-1)G(z)=[1-(1-z^-1)Θ(z)]/η^-1((1-z^-1)Θ(z)) for a 
fixed generating function Θ(z). Each tap coefficient g_τ can be 
computed by evaluating the coefficient of the τth-order term in 
G(z). 

 
Suppose that the sensing matrix A has independent Gaussian 
elements with mean √(γ/M) and variance (1-γ)/M for any 
γ∈[0, 1). Then, the tap coefficient g_t is given by 

    g_t = (
     1 - 1/δ)θ_t 
    + 1/δ∑_τ=0^t(θ_τ-θ_τ-1)θ_t-τ

for fixed tap coefficients {θ_t}. 


We shall evaluate the generating function G(z). 
The R-transform R(x) <cit.> of the asymptotic 
eigenvalue distribution of A^TA is given by 

    R(x) = δ/δ - x.

Using Theorem <ref> and the relationship between the 
R-transform and the η-transform <cit.> 

    η(x) = 1/1 + xR(-xη(x)),

we obtain 

    G(z)  
    = [
     1 - 1/δ + (1-z^-1)/δΘ(z)
    ]Θ(z),

which implies the time-domain expression (<ref>). 


In particular, we consider the original CAMP θ_τ=0 for τ>0. 
In this case, we have g_1=-δ^-1 and g_τ=0. As remarked in 
<cit.>, 
the original CAMP reduces to the AMP for the i.i.d. Gaussian sensing matrix.  

 
Suppose that the sensing matrix A has M identical non-zero 
singular values for M≤ N, i.e. 
AA^T=δ^-1I_M. 
Then, the tap coefficient g_t in the original CAMP θ_t=0 for t>0 
is given by g_τ=1-δ^-1 for all τ≥1. 


We evaluate the generating function G(z). 
By definition, the η-transform is given by 

    η(x) = 1/N(
     M/1+xδ^-1 + N-M
    ) = 1 - δ + δ^2/δ+x.

Using Theorem <ref> and Θ(z)=1 yields 

    G(z)
    = 1-δ^-1z^-1/1-z^-1
    = 1 + ∑_j=1^∞(
     1 - 1/δ)z^-j,

which implies g_τ=1-δ^-1 for all τ≥1. 
 

 
Suppose that the sensing matrix A has non-zero singular 
values σ_0≥⋯≥σ_M-1>0 satisfying condition number 
κ=σ_0/σ_M-1>1, σ_m/σ_m-1=κ^-1/(M-1), 
and σ_0^2=N(1-κ^-2/(M-1))/(1-κ^-2M/(M-1)). 
Assume θ_t=0 for all t>t_1 for some t_1∈ℕ. 
Let α_0^(j)=1 and 

    α_t^(j) 
    = {[    C^t/j/(t/j)!θ̅_j^t/j if t is divisible by j,;                       0               otherwise ].

for t∈ℕ and j∈{1,…,t_1}, 
with θ̅_t=θ_t-1-θ_t and  
C=2δ^-1lnκ. 
Define p_0=q̅_0=1 and 

    p_t
    = - β_t^(t_1)/κ^2-1,


    q̅_t 
    = 1/θ̅_1(
     β_t+1^(t_1)/C 
     - ∑_τ=1^t_1θ̅_τ+1q̅_t-τ)

for t>0, 
with β_t^(t_1)=α_t^(1)*α_t^(2)*⋯*α_t^(t_1). 
Then, the tap coefficient g_t is recursively given by

    g_t = p_t - ∑_τ=1^tq_τg_t-τ,

with 

    q_t = q̅_t - q̅_t-1.



We first evaluate the inverse of the η-transform. By definition, 
σ_m^2=κ^-2m/(M-1)σ_0^2 holds. Thus, we have 
rl
μ_j 
=    1/N∑_m=0^M-1σ_m^2j
= σ_0^2j1 - κ^-2jM/(M-1)/N(1 - κ^-2j/(M-1))


→    (
 C/1-κ^-2
)^j
1-κ^-2j/Cj

in the large system limit, where we have used the convergence 
N(1-κ^-a/(M-1))→δ^-1alnκ for any a∈ℝ. 
We note the series-expansion ln(1+x)=∑_j=1^∞(-1)^j-1j^-1x^j 
for |x|<1 to obtain 

    η(x) = 1 + ∑_j=1^∞(-x)^jμ_j 
    = 1 - 1/Cln(
     κ^2-1 + κ^2Cx/κ^2-1 + Cx),

which implies the inverse function 

    η^-1(x)
    = (κ^2-1){e^C(1 - x) - 1}/C{κ^2 - e^C(1 - x)}.


We next evaluate the generating function G(z). 
Using Theorem <ref> yields G(z)=P(z)/Q(z), with 

    P(z) 
    = κ^2 - e^CΘ̅(z)/κ^2-1,


    Q(z)
    = (1-z^-1)Q̅(z),   Q̅(z) 
    = e^CΘ̅(z) - 1/CΘ̅(z),


    Θ̅(z) = ∑_t=1^∞θ̅_tz^-t.


Finally, we derive a time-domain expression of G(z). 
It is an exercise to 
confirm that the series-expansions of P(z) 
and Q̅(z) have the coefficients p_t and q̅_t 
for the tth-order terms, respectively. Then, the Z-transform of 
(<ref>) is equal to P(z)/Q(z). 


The sequences {p_τ} and {q_τ} in Corollary <ref> 
define the generating functions P(z) and Q(z) with p_0=q_0=1, 
respectively, which satisfy 
G(z)=P(z)/Q(z). Thus, we derive an SE equation in time domain in terms of 
{p_τ, q_τ}, rather than {g_τ}.  




 §.§ SE Equation
 
We design the tap coefficients {θ_τ} so as to 
minimize the MSE N^-1x_t-x^2 for 
the CAMP estimator x_t in the large system limit. 
For that purpose, we derive an SE equation that describes  
the dynamics of the MSE. For simplicity, we assume i.i.d. signals. 

The CAMP has no closed-form SE equation with respect to the MSEs  
N^-1x_t-x^2 in general. Instead, 
it has a closed-form SE equation with respect to the correlations 

    d_t'+1,t+1= 𝔼[
     {f_t'(x_1+z_t') - x_1}{f_t(x_1+z_t) - x_1}],

where {z_t} denote zero-mean Gaussian random variables with 
covariance a_t',t=𝔼[z_t'z_t]. In particular, 
d_t+1,t+1 corresponds to the MSE of the CAMP estimator in iteration t.  

As an asymptotic alternative to ξ_t, we use the following 
quantity:  

    ξ̅_t = 𝔼[
     f_t'(x_1 + z_t)
    ],

which is a function of a_t,t. 
The notation ξ̅_t'^(t) is defined in the same manner as in 
ξ_t'^(t). 

  
Assume that Assumptions <ref>–<ref> hold, 
that the denoiser f_t is Lipschitz-continuous, and 
that the signal vector x has i.i.d. elements. 
Suppose that the generating functions G and Θ 
for the tap coefficients {g_τ} 
and {θ_τ}—given in (<ref>)—satisfy the 
condition (<ref>) in Theorem <ref>. 


  * Define generating functions A(y,z), D(y,z), and Σ(y,z) as 

    A(y,z) 
    = ∑_t',t=0^∞a_t',t/ξ̅_0^(t'-1)ξ̅_0^(t-1)
    y^-t'z^-t,


    D(y,z) 
    = ∑_t',t=0^∞d_t',t/ξ̅_0^(t'-1)ξ̅_0^(t-1)
    y^-t'z^-t,


    Σ(y,z) 
    = ∑_t',t=0^∞σ^2/ξ̅_0^(t'-1)ξ̅_0^(t-1)
    y^-t'z^-t.

Then, the correlation N^-1(x_t'-x)^T
(x_t-x) converges almost 
surely to d_t',t in the large system limit, which satisfies the following 
SE equation in terms of the generating functions:
rl 
F_G,Θ(y,z)A(y,z) 
=    {
 G(z)Δ_Θ - Θ(z)Δ_G 
}D(y,z) 


   + (
 Δ_Θ_1  - Δ_Θ   
)Σ(y,z), 


with 
rl
F_G,Θ(y,z) 
=    (y^-1+z^-1-1)[G(z)Δ_Θ - Θ(z)Δ_G]


   + Δ_G_1 - Δ_G,  

where the notations S_1(z)=z^-1S(z) and Δ_S
=[S(y) - S(z)]/(y^-1-z^-1) have been used for 
any generating function S(z). 


  * Suppose that G(z) is represented as G(z)=P(z)/Q(z) for the 
generating functions P(z) and Q(z) of some sequences {p_τ} and 
{q_τ} with p_0=1 and q_0=1. Let r_t=q_t*θ_t.  
Then, the SE equation (<ref>) reduces to 
l
∑_τ'=0^t'∑_τ=0^tξ̅_t'-τ'^(t'-1)ξ̅_t-τ^(t-1)
{ 𝔇_τ',τa_t'-τ',t-τ  


- (p_τ*r_τ'+τ+1 - r_τ*p_τ'+τ+1)d_t'-τ',t-τ


- σ^2[
 (q_τ'q_τ)*(θ_τ'+τ - θ_τ'+τ+1)
]} 
=0,   

where all variables with negative indices are set to zero, with 
rl
𝔇_τ',τ
   = (p_τ'+τ - p_τ'+τ+1)*q_τ + (p_τ - p_τ-1)*q_τ'+τ+1 


+    (p_τ-1-p_τ)*r_τ'+τ+1 + (r_τ - r_τ-1)*p_τ'+τ+1
 

+    p_τ*(r_τ'+τ - δ_τ',0r_τ) 
- r_τ*(p_τ'+τ - δ_τ',0p_τ). 

In solving the SE equation (<ref>), we impose 
the initial condition d_0,0=1 and boundary conditions 
d_0,τ+1=d_τ+1,0=-𝔼[x_1{f_τ(x_1+z_τ)-x_1}] 
for any τ.  



See Appendix <ref>. 


The SE equation (<ref>) in time domain is useful for numerical 
evaluation of {a_t',t} while the generating-function 
representation (<ref>) is used in fixed-point 
analysis. To apply Corollary <ref>, we have represented the 
generating function G(z) as G(z)=P(z)/Q(z). If G(z) is given directly, 
the functions P(z)=G(z) and Q(z)=1 can be used. In this case, we have 
p_τ=g_τ, q_τ=δ_τ,0, and r_τ=θ_τ.   

Note that d_t'+1,t+1 given in (<ref>) is a function of 
{a_t',t, a_t',t', a_t,t}, 
so that the SE equation (<ref>) in time domain is a nonlinear 
difference equation with respect to {a_t',t} for given tap 
coefficients {g_τ} and {θ_τ}.  
Theorem <ref> allows us to compute the MSEs 
a_t,t and d_t+1,t+1 before and after denoising.  

The SE equation (<ref>) in time domain can be solved recursively 
by extracting the term 𝔇_0,0a_t',t for τ'=τ=0 
in the sum and moving the other terms to the RHS. More precisely, 
we can solve the SE equation (<ref>) as follows:   



  * Let t=0 and solve a_0,0 with the SE equation (<ref>) 
and the initial condition d_0,0=1. 


  * 
Suppose that {a_τ',τ, d_τ',τ} have been 
obtained for all τ',τ=0,…,t. Use the boundary condition 
d_0,t+1 in Theorem <ref> and compute d_τ,t+1 with
the definition (<ref>) for all τ=1,…,t+1 
while the symmetry d_t+1,τ=d_τ,t+1 is used in the lower triangular 
elements. 


  * Compute a_τ,t+1 with the SE equation (<ref>) in the   
order τ=0,…,t+1 while the symmetry a_t+1,τ=a_τ,t+1 is used 
in the upper triangular elements. 


  * If some termination conditions are satisfied, output {a_τ',τ, 
d_τ',τ}. Otherwise, update 
t:=t+1 and go back to Step <ref>). 


We can define the Bayes-optimal denoiser f_t via the 
MSE d_t+1,t+1 in the large system limit. A denoiser f_t is 
said to be Bayes-optimal if f_t=𝔼[x_1|x_1+z_t] is    
the posterior mean of x_1 given an AWGN observation x_1 + z_t with 
z_t∼𝒩(0,a_t,t). We write the Bayes-optimal denoiser 
as f_t(·)=f_opt(·; a_t,t). 

The boundary condition d_0,τ+1 in Theorem <ref> 
has a simple representation for 
the Bayes-optimal denoiser f_opt. Since the posterior mean estimator 
f_opt(x_1+z_τ; a_τ,τ) is uncorrelated with the 
estimation error f_opt(x_1+z_τ; a_τ,τ)-x_1, we 
obtain 
rl
d_0,τ+1 
=    𝔼[{f_opt(x_1+z_τ; a_τ,τ) - x_1 
- f_opt(x_1+z_τ; a_τ,τ)} 

   ·{f_opt(x_1+z_τ; a_τ,τ)-x_1}] 


=    𝔼[{f_opt(x_1+z_τ; a_τ,τ)-x_1}^2]
= d_τ+1,τ+1. 



  
Consider the Bayes-optimal denoiser under the same assumptions as in 
Theorem <ref>. Suppose 
that the SE equation (<ref>) in time domain converges to a 
fixed-point {a_s, d_s}, 
i.e. lim_t',t→∞a_t',t=a_s and 
lim_t',t→∞d_t',t=d_s. If Θ(ξ_s^-1)=1 
and 1+(ξ_s-1)dΘ(ξ_s^-1)/(dz^-1)≠0 
hold for ξ_s=d_s/a_s, then 
the fixed-point {a_s, d_s} of the SE 
equation (<ref>) satisfies 

    a_s
    = σ^2/R(-d_s/σ^2), 
       
    d_s=𝔼[
     {f_opt(x_1+z_s; a_s) - x_1}^2],

with z_s∼𝒩(0,a_s), 
where R(x) denotes the R-transform of the asymptotic eigenvalue distribution 
of A^TA. 


See Appendix <ref>. 


The fixed-point equations given in (<ref>) coincide with those for 
describing the asymptotic performance of the posterior mean estimator 
of the signal vector x <cit.>.  
This coincidence implies that the CAMP with Bayes-optimal denoisers 
is Bayes-optimal if the SE equation (<ref>) converges 
toward a fixed-point and if the fixed-point is unique. 
Thus, we refer to CAMP with Bayes-optimal denoisers as Bayes-optimal CAMP. 



 §.§ Implementation
 
We summarize the implementation of the Bayes-optimal CAMP. We need to 
specify the sequence of denoisers {f_t} and 
the tap coefficients {g_τ, θ_τ} in (<ref>) and 
(<ref>). For simplicity, assume θ_τ=0 for all τ>2. 
To impose the condition Θ(a_s/d_s)=1 in 
Theorem <ref>, we use θ_0=1, 
θ_1=-θ d_s/a_s, and 
θ_2=θ∈ℝ, in which 
a_s and d_s are a solution to the fixed-point 
equations (<ref>). In particular, the CAMP reduces to the 
original one in <cit.> for θ=0.  

For a given parameter θ, the tap coefficients {g_τ} are 
determined via Theorem <ref>. More precisely, we use 
the coefficients {p_τ, q_τ} in the rational generating function 
G(z)=P(z)/Q(z). See Corollaries <ref>–<ref> for 
examples of the coefficients. 

For given parameters {θ, p_τ, q_τ}, we can solve the 
SE equation (<ref>) numerically. The obtained parameter 
a_t,t is used to determine the Bayes-optimal denoiser 
f_t(·)=f_opt(·; a_t,t). 

Damping <cit.> is a well-known technique to improve the convergence 
property in finite-sized systems. In damped CAMP, the update 
rule (<ref>) is replaced by 

    x_t+1 
    = ζ f_t(x_t+A^Tz_t)
    + (1-ζ)x_t,

with damping factor ζ∈[0, 1]. In solving the SE 
equation (<ref>), the associated parameters d_t'+1,t+1 and 
ξ̅_t in (<ref>) and (<ref>) are damped as 
follows: 
rl
d_t'+1,t+1 
=    ζ𝔼[
 { f_t'(x_1+z_t') - x_1}{f_t(x_1+z_t) - x_1}
] 

   + (1-ζ)d_t',t,  
 

    ξ̅_t 
    = ζ𝔼[
     f_t'(x_1 + z_t)
    ] + (1-ζ)ξ̅_t-1.

In particular, no damping is applied for ζ=1. 



Table <ref> lists time and space complexity of the CAMP, AMP, and 
OAMP/VAMP. Let t denote the number of iterations. 
We assume that the scalar parameters in the CAMP can be computed in 
O(t^4) time. In particular, computation of {a_t,t} via the 
SE equation (<ref>) is dominant.  

To compute the update rule (<ref>) in the CAMP efficiently,  
the vectors z_t∈ℝ^M and 
AA^Tz_t∈ℝ^M 
are computed and stored in iteration t. We need O(MN) space 
complexity to store the sensing matrix A, which is dominant 
for the case t≪ N. Furthermore, the time complexity is dominated by 
matrix-vector multiplications. 

In the OAMP/VAMP, the SVD of A requires dominant complexity 
unless the sensing matrix has a special structure that enables efficient SVD 
computation. As a result, the OAMP/VAMP has higher complexity than the AMP 
and CAMP while the CAMP has comparable complexity to the AMP for t≪ N. 



§ NUMERICAL RESULTS
 


 §.§ Simulation Conditions
 
The Bayes-optimal CAMP—called CAMP simply—is compared to the AMP and 
OAMP/VAMP. In all numerical results, 10^5 independent trials were 
simulated. We assumed the AWGN noise w∼𝒩(
0,σ^2I_M) and 
 i.i.d. Bernoulli-Gaussian signals with signal 
density ρ∈[0, 1] in the measurement model (<ref>). 
The probability density function (pdf) of x_n is given by

    p(x_n) = (1-ρ)δ(x_n) + ρ/√(2π/ρ)
    e^-x_n^2/2/ρ.
 
Since x_n has zero mean and unit variance, the signal-to-noise ratio (SNR) 
is equal to 1/σ^2. See Appendix <ref> for evaluation 
of the correlation d_t'+1,t+1 given in (<ref>).  

Corollary <ref> was used to simulate ill-conditioned sensing 
matrices A. The non-zero singular values 
{σ_m} of A are uniquely determined via the 
condition number κ. To reduce the complexity of the OAMP/VAMP, 
we assumed the SVD structure 
A=diag{σ_0,…,σ_M-1,
0}V^T. Note that the CAMP does not 
require this SVD structure. The CAMP only needs the right-orthogonal 
invariance of A. For a further reduction in the complexity, 
we used the Hadamard matrix V^T∈𝒪_N with 
the rows permuted uniformly and randomly. This matrix A 
is a practical alternative of right-orthogonally invariant matrices.  



We simulated damped AMP <cit.> with the same Bayes-optimal 
denoiser f_t(·)=f_opt(·; v_t) as in the 
CAMP. The variance parameter v_t was computed via the SE equation 

    v_t = σ^2 + 1/δMMSE(v_t-1), 
      MMSE(v_-1)=1,
 
with 

    MMSE(v) 
    = 𝔼[
     {f_opt(x_1+√(v)z; v) - x_1}^2],

where z∼𝒩(0,1) denotes the standard Gaussian random variable 
independent of x_1. The SE equation (<ref>) was derived in 
<cit.> 
under the assumption of zero-mean i.i.d. Gaussian sensing matrix with 
compression rate δ=M/N. Furthermore, ξ_t in (<ref>) 
was replaced by the asymptotic value 
ξ̅_t=MMSE(v_t)/v_t <cit.>. 
To improve the convergence 
property of the AMP, we replaced the update rule (<ref>) with 
the damped rule  

    x_t+1 
    = ζ f_t(x_t + A^Tz_t)
    + (1-ζ)x_t.

Note that SE cannot describe the exact 
dynamics of AMP when damping is employed. 


For the OAMP/VAMP <cit.>, we used the Bayes-optimal 
denoiser f_t(·)=f_opt(·;
v̅_A→B,t) computed via the SE 
equations <cit.>

    v̅_A→B,t 
    = γ̅_t - v̅_B→A,t,   v̅_B→A,0=1,
 

    1/v̅_B→A,t+1 
    = 1/MMSE(v̅_A→B,t) 
    - 1/v̅_A→B,t,

with 

    γ̅_t^-1 
    = lim_M=δ N→∞1/N∑_m=0^M-1σ_m^2/σ^2 + v̅_B→A,tσ_m^2.

To improve the convergence property, we applied the damping technique: 
The messages x_B→A,t+1 and 
v_B→A,t+1 in (<ref>) 
were replaced by the damped messages  
ζx_B→A,t+1
+(1-ζ)x_B→A,t and 
ζv̅_B→A,t+1
+(1-ζ)v̅_B→A,t, respectively.  
Note that damped EP cannot be described via SE. 




 §.§ Ill-Conditioned Sensing Matrices

We first consider the parameter θ in the CAMP defined in 
Section <ref>. From Theorem <ref>, we know that 
the CAMP is Bayes-optimal for any θ if it converges. Thus, the 
parameter θ only affects the convergence property of the CAMP.  

Figure <ref> shows the MSEs of the CAMP for a sensing matrix 
with condition number κ=5 defined in Corollary <ref>. 
As a baseline, we plotted the asymptotic MSE of the Bayes-optimal signal 
recovery <cit.>. 
The CAMP with θ=2 and ζ=0.9 converges to the Bayes-optimal 
performance more slowly than that with θ=0 and ζ=0.85. This 
observation does not necessarily imply that θ=0 is the best option. 
When the damping factor ζ=0.9 is used, the CAMP converges for θ=2 
in the finite-sized system while it diverges for θ=0. Thus, we 
conclude that using non-zero θ≠0 improves the stability of the 
CAMP in finite-sized systems. 



The CAMP is compared to the AMP and OAMP/VAMP for sensing matrices with 
unit condition number, i.e. orthogonal rows. As noted in Remark <ref>, 
without loss of generality, we can use θ=0 for this case. 
In this case, the 
OAMP/VAMP has comparable complexity to the AMP since the SVD of the sensing 
matrix is not required. Figure <ref> shows that the OAMP/VAMP is the best 
in terms of the convergence speed among the three MP algorithms. 




We next consider a sensing matrix with condition number κ=10. 
As shown in Fig. <ref>, the AMP cannot approach the Bayes-optimal 
performance. The CAMP converges to the Bayes-optimal performance more slowly 
than the OAMP/VAMP while the CAMP does not require high-complexity SVD of 
the sensing matrix. Especially in large systems, thus, the CAMP should need 
lower complexity to achieve the Bayes-optimal performance than the OAMP/VAMP. 



We investigate the influence of the condition number κ 
shown in Fig. <ref>. In evaluating the SE of the CAMP as a baseline, 
the parameter θ was optimized for each condition number while no 
damping was employed. 
In particular, the parameter θ was set to -0.7 for 
κ≥17. Otherwise, θ=0 was used. 
See Table <ref> for the parameters used in the three algorithms, 
which were numerically optimized for each condition number. More precisely, 
the parameters were selected so as to achieve the fastest convergence among 
all possible parameters that approach the best MSE in the last iteration. 



The AMP has poor performance with the exception of small condition numbers. 
The CAMP achieves the Bayes-optimal performance for low-to-moderate condition 
numbers. However, it is inferior to the high-complexity OAMP/VAMP for large 
condition numbers. These observations are consistent with the SE results of 
the CAMP. The SE prediction of the MSE changes rapidly from the Bayes-optimal 
performance to a large value around a condition number κ≈18 while 
the OAMP/VAMP still achieves the Bayes-optimal performance for κ>18.   
This is because the CAMP fails to converge for κ>18. As a result, 
we cannot use Theorem <ref> to 
claim the Bayes-optimality of the CAMP. Thus, we conclude that the CAMP is 
Bayes-optimal in a strictly smaller class of sensing matrices than the 
OAMP/VAMP.  



Finally, we investigate the convergence properties of the CAMP for high 
condition numbers. Figure <ref> shows the correlation d_t',t in the 
CAMP for t'=0,…,t. For the condition number κ=16, 
the correlation d_t',t converges toward the Bayes-optimal MSE for all t' 
as t increases. This provides numerical evidence for the assumption in 
Theorem <ref>: the convergence of the CAMP toward 
a fixed-point. 

The results for κ=17 imply that the CAMP fails to converge. 
A soliton-like quasi-steady wave propagates as t 
grows, while the CAMP does not diverge. As implied from Fig. <ref>, 
using non-zero θ≠0 allows us to avoid the occurrence of such a 
wave for κ=17. However, such waves occur for any θ when the 
condition number is larger than κ≈18.  

Intuitively, the occurrence of soliton-like waves can be understood as 
follows: The SE equation (<ref>) in time domain becomes unstable 
for high condition numbers, so that a_t',t increases as t grows.  
However, larger a_t',t results in a geometrically smaller forgetting factor 
ξ̅_t-τ^(t-1) in (<ref>), which suppresses the 
divergence of a_t',t. As a result, a soliton-like quasi-steady wave 
occurs for high condition numbers. 



§ CONCLUSIONS
  
The Bayes-optimal CAMP solves the disadvantages of AMP and OAMP/VAMP,  
and realizes their advantages for orthogonally invariant sensing matrices 
with low-to-moderate condition numbers: The Bayes-optimal CAMP is an efficient 
MP algorithm that has comparable complexity to AMP. Furthermore, the CAMP has 
been proved to be Bayes-optimal for all orthogonally invariant sensing 
matrices if it converges. High-complexity OAMP/VAMP is Bayes-optimal 
for this class of sensing matrices while AMP is not. The CAMP converges 
for sensing matrices with low-to-moderate condition numbers while it fails 
to converge for high condition numbers. 


A disadvantage of CAMP is to need all moments of the asymptotic 
singular-value distribution of the sensing matrix. In general, computation 
of the moments requires high complexity unless their closed-form is available. 
To circumvent this issue, deep unfolding <cit.> might  
be utilized to learn the tap coefficients in the Onsager correction without 
using the asymptotic singular-value distribution.  

The CAMP has a room for improvement especially in finite-sized and 
ill-conditioned sensing matrices. One option is a replacement of scalar  
parameters determined via the SE equation with empirical estimators that 
depend on the measurements, as considered in AMP and OAMP/VAMP. 

Another option is a damping technique that keeps the asymptotic Gaussianity 
of estimation errors. This paper used a heuristic damping technique to improve 
the convergence property of the CAMP. However, the heuristic damping breaks 
the asymptotic Gaussianity. Damped CAMP should be designed via 
Theorem <ref> to guarantee the asymptotic Gaussianity. 
A recent paper <cit.> proposed long-memory damping in the MF-based 
interference cancellation to improve the convergence property of long-memory 
MP. A possible future work is to design CAMP with long-memory damping. 






§ PROOF OF THEOREM <REF>

. 


 §.§ Formulation
 
We use Bolthausen's conditioning technique <cit.> to prove 
Theorem <ref>. In the technique, the random variables are 
classified into three groups: V, 
𝔉={λ,w̃,x}, 
and 𝔈_t,t'={B_t', M̃_t', 
H_t, Q̃_t+1} with 
Q̃_t+1=(q̃_0,…,
q̃_t) and M̃_t
=(m̃_0,…,m̃_t-1).  
The random variables in 𝔉 are fixed throughout the proof of 
Theorem <ref> while V is averaged out.  

The set 𝔈_t,t contains all messages just before updating 
b_t=V^Tq̃_t  
while 𝔈_t,t+1 includes all messages just before updating  
h_t=Vm̃_t. The main 
part in the conditioning technique is evaluation of  
the conditional distribution of b_t given 
𝔈_t,t and 𝔉 via that of V. 

Theorem <ref> is proved by induction. More precisely, we prove 
a theorem obtained by adding several technical results to 
Theorem <ref>. Before presenting the theorem, 
we first define several notations. 

The notation o(1) denotes a finite-dimensional vector 
with vanishing norm. 
For a tall matrix M∈ℝ^N× t with rank r≤ t, 
the SVD of M is denoted by 
M=Φ_MΣ_MΨ_M^T, 
with Φ_M=
(Φ_M^∥, 
Φ_M^⊥). The matrix  
Φ_M^∥∈𝒪_N× r consists 
of all left-singular vectors corresponding to r non-zero singular values 
while Φ_M^⊥∈𝒪_N×(N-r) 
is composed of left-singular vectors corresponding to N-r zero singular 
values. The matrix P_M^∥
=M(M^TM)^-1M^T is the projection to the space spanned by the 
columns of M while P_M^⊥
=I-P_M^∥ is the projection 
to the orthogonal complement. Note that P_M^∥
=Φ_M^∥
(Φ_M^∥)^T and 
P_M^⊥
=Φ_M^⊥
(Φ_M^⊥)^T hold. 

In the following theorem, we call the system with respect to 
{B_t, M̃_t} module A while 
we refer to that for {H_t, Q̃_t+1} 
as module B. 
  
Suppose that Assumptions <ref>–<ref> 
hold. Then, the following properties in module A hold  
for all τ=0,1,… in the large system limit.  


  * 
Let β_τ=(Q̃_τ^TQ̃_τ)^-1Q̃_τ^Tq̃_τ, q̃_τ^⊥
=P_Q̃_τ^⊥q̃_τ, and 

    ω̃_τ 
    = Ṽ^T
    (Φ_(Q̃_τ, H_τ)^⊥)^Tq̃_τ,

where Ṽ∈𝒪_N-2τ is a Haar orthogonal 
matrix and independent of 𝔉 and 𝔈_τ,τ. Then, 
for τ>0 

    b_τ∼B_τβ_τ 
    + M̃_τo(1) 
    + B_τo(1) 
    + Φ_(B_τ, M̃_τ)^⊥ω̃_τ

conditioned on 𝔉 and 𝔈_τ,τ  
in the large system limit, with 

    lim_M=δ N→∞1/N{ω̃_τ^2 
     - q̃_τ^⊥^2} 0.



  * 
Suppose that ϕ̃_τ(B_τ+1, 
w̃, λ):
ℝ^N×(τ+3)→ℝ^N 
is separable, pseudo-Lipschitz of order k with respect to 
the first τ+2 variables, and proper. 
If N^-1q̃_t^Tq̃_t' 
converges almost surely to some constant κ_t,t'∈ℝ 
in the large system limit for all t, t'=0,…,τ, then   

    ⟨ϕ̃_τ(B_τ+1, 
    w̃;λ)⟩
    - 𝔼[
     ⟨ϕ̃_τ(Z̃_τ+1, 
     w̃,λ)⟩] 0.

In (<ref>), 
Z̃_τ+1=(z̃_0,…,
z̃_τ)∈ℝ^N×(τ+1) denotes a zero-mean 
Gaussian random matrix with covariance 𝔼[z̃_tz̃_t'^T]
=κ_t,t'I_N for all t, t'=0,…,τ. In evaluating 
the expectation in (<ref>), 
U^Tw in (<ref>) follows 
the zero-mean Gaussian distribution with covariance 
σ^2I_M.  
In particular, for k=1 we have 

    ⟨∂_τ'ϕ̃_τ(B_τ+1, 
    w̃;λ)⟩
    - 𝔼[
     ⟨∂_τ'ϕ̃_τ
     (Z̃_τ+1, 
     w̃;λ)⟩] 0

for all τ'=0,…,τ. 


  *  
Suppose that ϕ̃_τ(B_τ+1,
w̃;λ):
ℝ^N×(τ+3)→ℝ^N is separable, 
Lipschitz-continuous with respect to the first τ+2 variables, 
and proper. Then, 

    1/Nb_τ'^T(
     ϕ̃_τ
     - ∑_t'=0^τ⟨∂_t'ϕ̃_τ⟩b_t')
    0

for all τ'=0,…,τ. 


  * 
The inner product 
N^-1m̃_τ'^Tm̃_τ 
converges almost surely to some constant π_τ',τ∈ℝ for 
all τ'=0,…,τ. 


  * 
For some ϵ>0 and C>0, 

    lim_M=δ N→∞𝔼[
     |m̃_τ,n|^2k-2+ϵ]<∞,
 

    lim inf_M=δ N→∞λ_min(
     1/NM̃_τ+1^TM̃_τ+1)
     C.



The following properties in module B hold for all τ=0,1,… 
in the large system limit. 


  *  
Let α_τ=(M̃_τ^TM̃_τ)^-1M̃_τ^Tm̃_τ, m̃_0^⊥
=m̃_0, m̃_τ^⊥
=P_M̃_τ^⊥m̃_τ, and  

    ω_τ 
    = {[             Ṽ(Φ_b_0^⊥)^Tm̃_0                     for τ=0,; Ṽ
    (Φ_(M̃_τ, 
    B_τ+1)^⊥)^Tm̃_τ                     for τ>0,;                              ].

where Ṽ∈𝒪_N-(2τ+1) is a Haar orthogonal 
matrix and independent of 𝔉 and 𝔈_τ,τ+1. Then, 
we have 

    h_0∼ o(1)q̃_0 
    + Φ_q̃_0^⊥ω_τ,

conditioned on 𝔉 and 
𝔈_0,1={b_0, m̃_0, 
q̃_0} in the large system limit. For τ>0 

    h_τ∼H_τα_τ 
    + Q̃_τ+1o(1) 
    + H_τo(1)  
    + Φ_(H_τ, Q̃_τ+1)^⊥ω_τ,

conditioned on 𝔉 and 𝔈_τ,τ+1 
in the large system limit, with 

    lim_M=δ N→∞1/N{ω_τ^2  
     - m̃_τ^⊥^2} 0.



  * 
Suppose that ψ̃_τ(H_τ+1, 
x):ℝ^N×(τ+2)→ℝ^N 
is a separable and proper pseudo-Lipschitz function of order k.  
If N^-1m̃_t^Tm̃_t' 
converges almost surely to some constant π_t,t'∈ℝ 
in the large system limit for all t, t'=0,…,τ, then   

    ⟨ψ̃_τ(H_τ+1, 
    x)⟩
    - 𝔼[
     ⟨ψ̃_τ(Z_τ+1, 
     x)⟩] 0,

where Z_τ+1=(z_0,…,z_τ)
∈ℝ^N×(τ+1) denotes a zero-mean Gaussian random matrix with 
covariance 𝔼[z_tz_t'^T]
=π_t,t'I_N for all t, t'=0,…,τ. 
In particular, for k=1 we have 

    ⟨∂_τ'ψ̃_τ(H_τ+1, 
    x)⟩
    - 𝔼[
     ⟨∂_τ'ψ̃_τ(Z_τ+1, 
     x)⟩] 0

for all τ'=0,…,τ. 


  * 
Suppose that ψ̃_τ(H_τ+1,
x):ℝ^N×(τ+2)→ℝ^N is a 
separable and proper Lipschitz-continuous function. Then, 

    1/Nh_τ'^T(
     ψ̃_τ
     - ∑_t'=0^τ⟨∂_t'ψ̃_τ⟩h_t')
    0

for all τ'=0,…,τ. 


  *  
The inner product 
N^-1q̃_τ'^Tq̃_τ+1 
converges almost surely to some constant π_τ',τ+1∈ℝ for 
all τ'=0,…,τ+1. 


  * 
For some ϵ>0 and C>0, 

    lim_M=δ N→∞𝔼[
     |q̃_τ+1,n|^2+ϵ]<∞,
 

    lim inf_M=δ N→∞λ_min(
     1/NQ̃_τ+2^TQ̃_τ+2)
     C.




We summarize useful lemmas used in the proof of Theorem <ref> 
by induction. 
 
 
Suppose that X∈ℝ^N× t has full rank for 
0<t<N, and consider the noiseless and compressed observation 
Y∈ℝ^N× t of V given by 

    Y = VX.

Then, the conditional distribution of the Haar orthogonal matrix 
V given X and Y satisfies 

    V|_X, Y∼Y(Y^TY)^-1X^T + Φ_Y^⊥Ṽ(Φ_X^⊥)^T,

where Ṽ∈𝒪_N-t is a Haar orthogonal matrix 
independent of X and Y. 


The following lemma is a generalization of Stein's lemma.  
The lemma is proved under a different assumption from in <cit.>. 

 
Let z=(z_1,…,z_t)^T∼𝒩(0,Σ) for any positive 
definite covariance matrix Σ. 
If f:ℝ^t→ℝ is Lipschitz-continuous, 
then we have 

    𝔼[z_1f(z)] 
    = ∑_t'=1^t𝔼[z_1z_t']𝔼[∂_t'f(z)].



We first confirm that both sides of (<ref>) are bounded. 
For the left-hand side (LHS), we find 
f(z)= O(z) as 
z→∞ since f is Lipschitz-continuous. 
Thus, 𝔼[z_1f(z)] is 
bounded for z∼𝒩(0,Σ). 

For the RHS, we use the Lipschitz-continuity of f to find that 
there is some Lipschitz-constant L>0 such that 

    |
     f(z+Δe_t') - f(z)/Δ|
    ≤ L

holds for any Δ≠0, 
where e_t'∈ℝ^t 
is the t'th column of I_t. This implies that each 
partial derivative ∂_t'f is bounded almost everywhere since the 
partial derivatives of any Lipschitz-continuous function exist almost 
everywhere. Thus, 𝔼[∂_t'f(z)] is 
bounded. These observations indicate the boundedness of both sides 
in (<ref>). 


For the eigen-decomposition 
Σ=ΦΛΦ^T, we use the change of variables 
z̃=Φ^Tz to obtain 

    𝔼[z_1f(z)] 
    = ∑_τ=1^t[Φ]_1,τ𝔼[z̃_τf(Φz̃)] 
    = ∑_τ=1^t[Φ]_1,τ𝔼[z̃_τg(z̃_τ)],

with g(z̃_τ)
=𝔼[f(Φz̃)|z̃_τ]. 

We prove that g is Lipschitz-continuous. 
Let z̃_x denote the vector obtained by replacing 
z̃_τ in z̃ with x. 
Since z̃∼𝒩(0,
Λ) has independent elements, we have 
rl
|g(x) - g(y)| 
≤    𝔼[
 |
  f(Φz̃_x) 
  - f(Φz̃_y) 
 |
] 

≤    L𝔼[
 Φ(z̃_x - z̃_y)
] 

=    L𝔼[
 z̃_x - z̃_y
]
= L|x-y|, 

where the second inequality follows from the Lipschitz-continuity of f with 
a Lipschitz-constant L>0. 
Thus, g(z̃_τ) is Lipschitz-continuous, so that 
g(z̃_τ) is differentiable almost everywhere. 

Since z̃∼𝒩(0,
Λ) holds, Stein's lemma <cit.> yields 
rl
𝔼[z_1f(z)] 
   = ∑_τ=1^t[Φ]_1,τ
𝔼[z̃_τ^2]𝔼[
 g'(z̃_τ)
] 

=   ∑_τ=1^t[Φ]_1,τ
[Λ]_τ,τ𝔼[
 ∑_t'=1^t[Φ]_t',τ∂_t'
 f(z)
]. 

Using the identity 

    ∑_τ=1^t[Φ]_1,τ
    [Λ]_τ,τ[Φ]_t',τ
    = [ΦΛΦ^T
    ]_1,t'
    = 𝔼[z_1z_t'],

we arrive at Lemma <ref>. 


 
For t∈ℕ, 
suppose that f:ℝ^N×(t+1)→ℝ^N is 
separable and pseudo-Lipschitz of order k. Let L_n>0 denote 
a Lipschitz constant of the nth element [f]_n. The 
sequence of Lipschitz constants is assumed to satisfy 

    lim sup_N→∞1/N∑_n=1^NL_n^2<∞.

Let ϵ=(ϵ_1,…,ϵ_N)^T∈ℝ^N denote a vector that satisfies  

    lim_N→∞1/N∑_n=1^NL_nϵ_n^20,


    lim sup_N→∞1/N∑_n=1^NL_nϵ_n^2k-2∞.
 
Suppose that A_t+1=(a_0,…,a_t)
∈ℝ^N×(t+1)  satisfies 

    lim sup_N→∞1/N∑_n=1^NL_n^ia_t',n^2k-2∞  for i=1, 2.
 
For t'>0, let E=(e_1^T,…, 
e_N^T)^T∈ℝ^N× t' denote 
a matrix that satisfies   

    lim sup_N→∞1/N∑_n=1^NL_ne_n^max{2, 2k-2}∞,


    lim inf_N→∞λ_min(
     1/NE^HE)  C

for some constant C>0. 
Suppose that ω∈ℝ^N-t' is an orthogonally 
invariant random vector conditioned on 
ϵ, A_t+1, and E. 
For some v>0, postulate the following: 

    lim_N→∞1/Nω^2 v>0.

Let z∼𝒩(0,vI_N) 
denote a standard Gaussian random vector independent of the other 
random variables. Then,  

    lim_N→∞⟨f(A_t,
     a_t+ϵ
     +Φ_E^⊥ω)
     - 𝔼_z[f(
     A_t, 
     a_t+z)]
    ⟩ 0.
  





 §.§ Module A for τ=0

[Proof of Property <ref> for τ=0]
The latter property (<ref>) follows from the 
former property (<ref>) and a technical result proved 
in <cit.>. Thus, we only prove the former  
property for τ=0. 

Property (<ref>) follows from Lemma <ref> 
for f(w̃,b̃_0)
=ϕ̃_0(b̃_0,
w̃;λ) with 
a_0=w̃, 
a_1+ϵ=0, 
Φ_E^⊥=I_N, and  
ω=b̃_0. 
We confirm all conditions in Lemma <ref>. 
Applying Hölder's inequality for any ϵ>0, we have    

    1/N∑_n=1^NL_n^iw̃_n^2k-2≤(
     1/N∑_n=1^NL_n^ip)^1/p(
     1/N∑_n=1^Nw̃_n^2k-2+ϵ)^1/q

for i=1, 2, with q=1+ϵ/(2k-2) and p^-1=1-q^-1, which is bounded 
because of Assumption <ref>. Furthermore, the definition 
b_0=-V^Tx implies the 
orthogonal invariance and N^-1b_0^2 1. 
Thus, all conditions in Lemma <ref> hold. 
Using Lemma <ref>, we obtain 

    ⟨ϕ̃_0(b̃_0,
    w̃; λ)⟩ 
    - 𝔼_z̃_0[
     ⟨ϕ̃_0(z̃_0,
     w̃;λ)⟩] 0,

with z̃_0∼𝒩(0,
I_N). 

We repeat the use of Lemma <ref> for 
f(z̃_0,w̃)
=ϕ̃_0(z̃_0,
w̃; λ) with 
a_0=z̃_0 and 
ω=w̃. Using Lemma <ref> 
from Assumption <ref> and applying Assumption <ref>, 
we obtain 

    ⟨ϕ̃_0(z̃_0,
    w̃;λ)⟩ 
    - 𝔼[
     ⟨ϕ̃_0(z̃_0,
     w̃;λ)⟩] 0.

In evaluating the expectation over w̃, 
the first M elements U^Tw 
in (<ref>) follow 
𝒩(0,σ^2I_M). 
Combining these results, we arrive at (<ref>) for τ=0. 


[Proof of <ref> for τ=0] 
The LHS of (<ref>) is a separable and 
proper pseudo-Lipschitz function of order 2. We can use 
(<ref>) for τ=0 to find that the LHS of 
(<ref>) converges almost surely to its expectation in which 
b_0 and ⟨∂_0ϕ̃_0⟩ are replaced by 
z̃_0∼𝒩(0,
I_N) and the expected one, 
respectively. Thus, it is sufficient to evaluate the expectation.  

The function f(z̃_0; 
w̃, λ)
=ϕ̃_0(z̃_0,
w̃;λ) 
- 𝔼[ ⟨∂_0ϕ̃_0⟩]
z̃_0 is a separable Lipschitz-continuous function of 
z̃_0. 
Thus, we can use Lemma <ref> to obtain 
rl
   1/N𝔼[
 z̃_0^T(
  ϕ̃_0
  - 𝔼[
   ⟨∂_0ϕ̃_0
   ⟩]z̃_0
 )
] 

=    1/N∑_n=1^N𝔼[
 z̃_0,n^2
]𝔼[
 ∂_0ϕ̃_0,n
]
- 𝔼[
 ⟨∂_0ϕ̃_0
 ⟩]=0. 

Thus, (<ref>) holds for τ=0. 


[Proof of <ref> for τ=0] 
From the definition (<ref>) of m̃_0 and 
(<ref>), we find the orthogonality 
N^-1b_0^Tm̃_00. 
Using this orthogonality and (<ref>) for τ=0 yields 
rl
   1/Nm̃_0^2
1/Nm_0^Tm̃_0 + o(1) 


=    1/Nm_0^Tm_0
- 𝔼[
 ⟨∂_0ϕ_0⟩]m_0^Tb_0/N
+ o(1). 


The first and second terms are separable and proper pseudo-Lipschitz functions 
of order 2. From (<ref>) for τ=0, they converge  
almost surely to their expected terms. Thus, 
N^-1m̃_0^2 converges almost surely to a 
constant.  


[Proof of Property <ref> for τ=0]
The latter property (<ref>) for τ=0 follows from the nonlinearity 
of ϕ_0 in Assumption <ref>. 
Thus, we only prove the former property (<ref>) for τ=0. 

The proper Lipschitz-continuity in Assumption <ref> 
implies the upper bound |m̃_0,n|≤ C_n(1+|b_0,n|
+|w̃_0,n|) for some λ_n-dependent constant C_n. 
From Assumptions <ref> and <ref>, we find that 
b_0 and w̃ have 
bounded (2k-2+ϵ)th moments for some ϵ>0. Thus, 
we obtain the former property (<ref>) for τ=0.  




 §.§ Module B for τ=0
 
[Proof of Property <ref> for τ=0]
Lemma <ref> for the constraint 
Vb_0=q̃_0 implies  

    V∼q̃_0b_0^T/q̃_0^2 
    + Φ_q̃_0^⊥Ṽ
    (Φ_b_0^⊥)^T
 
conditioned on 𝔉 and 𝔈_0,0, where 
Ṽ∈𝒪_N-1 is Haar orthogonal and independent 
of b_0 and q̃_0. 
Using the definition (<ref>) of h_0 and 
the orthogonality N^-1b_0^Tm̃_00 obtained from Property <ref> 
for τ=0, we obtain (<ref>).  

To complete the proof of Property <ref> for τ=0, we prove 
(<ref>) for τ=0. By definition, 

    1/Nω̃_0^2 
    = 1/Nm̃_0^TP_b_0^⊥m̃_01/Nm̃_0^2,

where the last equality follows from the orthogonality 
N^-1b_0^Tm̃_00. 
Thus, (<ref>) holds for τ=0, because of 
the notational convention m̃_0^⊥
=m̃_0. 


[Proof of Property <ref> for τ=0] 
Since the latter property (<ref>) follows from the 
former property (<ref>), we only prove the former  
property for τ=0. Using Property <ref> for τ=0 and 
Lemma <ref> for f(x,h_0)
=ψ̃_0(h_0,x) with 
a_0=x, a_1=0, 
ϵ=o(1)q̃_0, 
E=q̃_0, and 
ω=ω_0, we obtain 

    ⟨ψ̃_0(h_0, 
    x)⟩
    - 𝔼_z_0[
     ⟨ψ̃_0(z_0, 
     x)⟩] 0,
 
with z_0∼𝒩(0,
π_0,0I_N). Applying Assumption <ref> to 
the second term, we arrive at (<ref>) for τ=0. 


[Proof of Properties <ref> and 
<ref> for τ=0]
Repeat the proofs of Properties <ref> and 
<ref> for τ=0. 


[Proof of Property <ref> for τ=0]
The former property (<ref>) for τ=0 is obtained by repeating 
the proof of (<ref>) for τ=0. See <cit.> 
for the proof of the latter property (<ref>) for τ=0. 




 §.§ Proof by Induction
 
Suppose that Theorem <ref> is correct for all 
τ<t. In a proof by induction we need to prove all properties in 
modules A and B for τ=t. Since the properties for module B can be proved 
by repeating the proofs for module A, we only prove the properties for 
module A. 

[Proof of Property <ref> for τ=t]
The matrix (B_t, M̃_t) has full rank 
from the induction hypotheses (<ref>) and (<ref>) for 
τ=t-1, as well as the orthogonality N^-1b_τ^Tm̃_τ'0 for all τ, τ'<t.  
Using Lemma <ref> for the constraint 
(Q̃_t, H_t) 
=V(B_t, M̃_t), we obtain  
rl
V 
=    (Q̃_t,H_t)
[ Q̃_t^TQ̃_t  Q̃_t^TH_t;  H_t^TQ̃_t   H_t^TH_t ]^-1
[  B_t^T; M̃_t^T ] 

   + Φ_(Q̃_t, H_t)^⊥
Ṽ(Φ_(B_t, 
M̃_t)^⊥)^T 

conditioned on 𝔉 and 𝔈_t,t. 
Applying the orthogonality N^-1b_τ^Tm̃_τ'0 and N^-1h_τ^Tq̃_τ'0 obtained from the induction 
hypotheses <ref> and <ref> for τ<t, 
as well as the definition (<ref>) of b_t, we have 
rl
b_t
∼    B_t(Q̃_t^T
Q̃_t)^-1Q̃_t^T
q̃_t  
+ B_to(1) 
+ M̃_to(1) 


   + Φ_(B_t, 
M̃_t)^⊥Ṽ^T
(Φ_(Q̃_t, H_t)^⊥)^T
q̃_t 

conditioned on 𝔉 and 𝔈_t,t, which is 
equivalent to (<ref>) for τ=t. 

To complete the proof of Property <ref> for τ=t, 
we shall prove (<ref>). By definition, 

    ω̃_t^2/N 
    = q̃_t^TP_(Q̃_t,H_t)^⊥q̃_t/Nq̃_t^TP_Q̃_t^⊥q̃_t/N+ o(1),

where the last equality follows from the orthogonality 
N^-1h_τ^Tq̃_τ'0. 
Thus, (<ref>) holds for τ=t.  


[Proof of Property <ref> for τ=t]
Since the latter property (<ref>) follows from the 
former property (<ref>), we only prove the former  
property for τ=t. 

We use Property <ref> for τ=t and 
Lemma <ref> for the function f(w̃, 
B_t, b_t) 
=ϕ̃_t(B_t+1,w̃;
λ) with A_t+1=(w̃, 
B_t), a_t+1
=B_tβ_t, 
ϵ=M̃_to(1) 
+ B_to(1), 
E=(B_t,M̃_t), and 
ω=ω̃. Then, 

    ⟨ϕ̃_t(B_t+1, 
    w̃;λ)⟩
    - 𝔼_z̃_t[
     ⟨ϕ̃_t(B_t, 
     B_tβ_t+z̃_t, 
     w̃,λ)⟩] 0,

where z̃_t has independent zero-mean Gaussian 
elements with variance μ_t N^-1q̃_t^⊥^2. Repeating this argument yields 

    ⟨ϕ̃_t(B_t+1, 
    w̃;λ)⟩
    - 𝔼[
     ⟨ϕ̃_t(Z̃_t+1, 
     w̃,λ)⟩] 0,

where Z̃_t+1 is a zero-mean Gaussian random matrix 
having independent elements. In evaluating the expectation over 
w̃,  
U^Tw in (<ref>) follows 
the zero-mean Gaussian distribution with covariance 
σ^2I_M. 

To complete the proof of (<ref>) for τ=t, we evaluate 
the covariance of Z_t+1. By construction, we have 
N^-1𝔼[z_τ^Tz_τ']
= N^-1b_τ^Tb_τ'κ_τ,τ'+o(1). Thus, the former 
property (<ref>) is correct for τ=t. 


[Proof of Property <ref> for τ=t]
The LHS of (<ref>) is a separable and proper 
pseudo-Lipschitz function of order 2. We can use 
(<ref>) for τ=t to find that the LHS of 
(<ref>) converges almost surely to its expectation in which 
B_t+1 and ⟨∂_t'ϕ̃_t⟩ are replaced by Z̃_t+1 and the expected one, 
respectively. Thus, it is sufficient to evaluate the expectation.  

Since the function f(Z̃_t+1; 
w̃, λ)
=ϕ̃_t(Z̃_t+1,
w̃;λ) - ∑_t'=0^t𝔼[ ⟨∂_t'ϕ̃_t⟩]
z̃_t' is separable and Lipschitz-continuous with 
respect to Z̃_t+1, 
we can use Lemma <ref> to obtain 
rl
   1/N𝔼[
 z̃_τ'^T(
  ϕ̃_t
  - ∑_t'=0^t𝔼[
   ⟨∂_t'ϕ̃_t
   ⟩]z̃_t'
 )
] 

=    1/N∑_n=1^N∑_t'=0^t
𝔼[z̃_τ',nz̃_t,n]
𝔼[
 ∂_t'ϕ̃_t,n
] 

   - ∑_t'=0^t𝔼[
 ⟨∂_t'ϕ̃_t
 ⟩]𝔼[z̃_τ'^T
z̃_t']/N=0. 

Thus, (<ref>) holds for τ=t. 


[Proof of Properties <ref> and 
<ref> for τ=t]
Repeat the proofs of Properties <ref> and 
<ref> for τ=0. In particular, see <cit.> 
for the proof of (<ref>) for τ=t. 




§ PROOF OF THEOREM <REF>

 
In evaluating the derivative in g_t',t^(j), the parameter ξ_t 
requires a careful treatment since it depends on B_t+1 via 
h_t. If the general error model contained the error model of 
the CAMP, we could use (<ref>) in Theorem <ref> to 
prove that ξ_t converges almost surely to a 
B_t+1-independent constant ξ̅_t in the large system 
limit. To use Theorem <ref>, however, we have to prove the 
inclusion of the CAMP error model into the general error model. 
To circumvent this dilemma, we prove g_t-τ,t^(j)ξ_t-τ^(t-1)g_τ^(j)+o(1) for all t and τ=0,…,t 
by induction. 

We consider the case τ=0, in which the expression (<ref>) 
requires no special treatments in computing the derivative. 
Differentiating (<ref>) with respect to the tth variable yields 

    g_t,t^(j) = μ_j+1 - μ_j,

where μ_j denotes the jth moment (<ref>) of the asymptotic 
eigenvalue distribution of A^TA. 
Comparing (<ref>) and (<ref>), we have 
g_t,t^(j)=g_0^(j) for all t. 

Suppose that there is some t>0 such that 
g_t'-τ,t'^(j)ξ_t'-τ^(t'-1)g_τ^(j)+o(1) is correct 
for all t'<t and τ=0,…,t'. Then, (<ref>) in 
Theorem <ref> implies that ξ_t' converges almost 
surely to a constant ξ̅_t' for any t'<t. 
We need to prove g_t-τ,t^(j)ξ_t-τ^(t-1)g_τ^(j)+o(1) 
for all τ=0,…,t. 

We first consider the case τ=1 since we have already proved the case 
τ=0. Differentiating (<ref>) with respect to the (t-1)th 
variable yields  
rl
g_t-1,t^(j)
=    ξ̅_t-1(g_t-1,t-1^(j) - g_t-1,t-1^(j+1))
- ξ̅_t-1g_1(g_t-1,t-1^(j) + μ_j)


   + ξ̅_t-1θ_1(g_t-1,t-1^(j+1) + μ_j+1).  

Using g_t,t^(j)=g_0^(j) and (<ref>), we arrive at 
g_t-1,t^(j)ξ_t-1g_1^(j)+o(1).  

We next consider the case τ>1. 
Differentiating (<ref>) with respect to the (t-τ)th variable, 
we have 
l
g_t-τ,t^(j) 
= ξ̅_t-1(g_t-τ,t-1^(j) - g_t-τ,t-1^(j+1))


+ ∑_τ'=t-τ^t-1ξ̅_τ'^(t-1)(θ_t-τ'g_t-τ,τ'^(j+1)
- g_t-τ'g_t-τ,τ'^(j))


- ∑_τ'=t-τ+1^t-1ξ̅_τ'-1^(t-1)
(θ_t-τ'g_t-τ,τ'-1^(j+1) - g_t-τ'g_t-τ,τ'-1^(j)) 

 
+ ξ̅_t-τ^(t-1)(θ_τμ_j+1 - g_τμ_j). 

Using (<ref>) and the induction hypothesis 
g_t'-τ,t'^(j)ξ_t'-τ^(t'-1)
g_τ^(j)+o(1) for all t'<t and τ=0,…,t', 
we find g_t-τ,t^(j)ξ_t-τ^(t-1)g_τ^(j)+o(1).  



§ PROOF OF THEOREM <REF>


Let G(x,z) denote the generating function of {g_τ^(j)} given by  

    G(x,z) = ∑_j=0^∞G_j(z)x^j,

with 

    G_j(z) = ∑_τ=0^∞g_τ^(j)z^-τ.

It is possible to prove that G(x,z) is given by 

    G(x,z) = {Θ(z) - xG(z)}η(-x) - Θ(z)/xG̃(z) + 1 - Θ̃(z),

with G̃(z)=(1-z^-1)G(z) and Θ̃(z)=(1-z^-1)Θ(z). 
Let -x^* denote a pole of the generating function, i.e. 
x^* = [1 - Θ̃(z)]/G̃(z). Since the generating function 
is analytical, the numerator of (<ref>) at x=-x^* must be zero. 

    {Θ(z) + x^*G(z)}η(x^*) - Θ(z) = 0,

which is equivalent to (<ref>). 
 
To complete the proof of Theorem <ref>, 
we prove (<ref>). The proof is a simple exercise of the Z-transform. 
We first compute G_j(z) given by 

    G_j(z) 
    = g_0^(j) + g_1^(j)z^-1
    + ∑_τ=2^∞g_τ^(j)z^-τ.

To evaluate the last term with (<ref>), we note 

    ∑_τ=2^∞g_τ-1^(j)z^-τ
    = z^-1∑_τ=1^∞g_τ^(j)z^-τ
    = z^-1{
     G_j(z) - g_0^(j)},

rl
   ∑_τ=2^∞∑_τ'=0^τ-1g_τ-τ'g_τ'^(j)z^-τ 


=    g_0^(j)∑_τ=2^∞g_τz^-τ 
+ ∑_τ'=1^∞∑_τ=τ'+1^∞
g_τ-τ'g_τ'^(j)z^-τ 

=    [
 G(z) - 1
]G_j(z) - g_1g_0^(j)z^-1, 

rl
   ∑_τ=2^∞∑_τ'=1^τ-1g_τ-τ'g_τ'-1^(j)z^-τ 


=    ∑_τ'=1^∞∑_τ=τ'+1^∞
g_τ-τ'g_τ'-1^(j)z^-τ 

=    [
 G(z) - 1
]z^-1G_j(z). 

Combining (<ref>), (<ref>), (<ref>), and these results, 
we arrive at  
rl
G_j(z) 
=    [1-G̃(z)]G_j(z) - [1-Θ̃(z)]G_j+1(z) 

   - μ_jG(z) + μ_j+1Θ(z). 


We next evaluate G(x,z). Substituting (<ref>) into the definition 
of G(x,z) yields 
rl
G(x,z)   
=    [1-G̃(z)]G(x,z) - [1 - Θ̃(z)]G(x,z)/x 

   - η(-x)G(z) + η(-x)-1/xΘ(z),   

where we have used the definition (<ref>) and 
the identity G_0(z)=0 obtained from 
Theorem <ref>. Solving this equation with respect to 
G(x,z), we obtain (<ref>). 




§ PROOF OF THEOREM <REF>

 


 §.§ SE Equations

The proof of Theorem <ref> consists of four steps: 
A first step is a derivation of the SE equations, which is a dynamical 
system that describes the dynamics of five variables with three indices.  
A second step is evaluation of the generating functions for the five 
variables. The step is a simple exercise of the Z-transform. 
In a third step, we evaluate the obtained generating functions at poles 
to prove the SE equation (<ref>) in terms of 
the generating functions. The last step is a derivation of the SE 
equation (<ref>) in time domain via the inverse Z-transform. 

Let a_t',t^(j)=N^-1m_t'^TΛ^jm_t, 
b_t',t^(j)=N^-1b_t'^TΛ^jm_t, 
c_t',t=N^-1q̃_t'^Tq̃_t, d_t',t=N^-1q_t'^Tq_t, and e_t^(j)=N^-1w^TUΣΛ^jm_t. 
Theorem <ref> implies the asymptotic orthogonality between 
b_t' and m_t.  
We use the definition (<ref>) to obtain  
l
a_t',t^(j) 
b_t,t'^(j) - b_t,t'^(j+1) 
+ ξ̅_t-1(a_t',t-1^(j) - a_t',t-1^(j+1))
+ e_t'^(j) 


+ ∑_τ=0^t-1ξ̅_τ^(t-1)θ_t-τ
(a_t',τ^(j+1) - b_τ,t'^(j+1) - ξ̅_τ-1a_t',τ-1^(j+1))

 
- ∑_τ=0^t-1ξ̅_τ^(t-1)g_t-τ
(a_t',τ^(j) - b_τ,t'^(j) - ξ̅_τ-1a_t',τ-1^(j))
+ o(1), 


where we have replaced ξ_t with the asymptotic value ξ̅_t. 
Applying (<ref>) in Theorem <ref> and 
(<ref>) yields 
rl
b_t',t^(j) 
    (μ_j - μ_j+1)c_t',t  
+ ξ̅_t-1(b_t',t-1^(j) - b_t',t-1^(j+1)) + o(1)


+    ∑_τ=0^t-1ξ̅_τ^(t-1)θ_t-τ
(b_t',τ^(j+1) - μ_j+1c_t',τ - ξ̅_τ-1b_t',τ-1^(j+1))


-    ∑_τ=0^t-1ξ̅_τ^(t-1)g_t-τ
(b_t',τ^(j) - μ_jc_t',τ - ξ̅_τ-1b_t',τ-1^(j)).  

Using (<ref>) in Theorem <ref>, 
(<ref>), and (<ref>), we have  

    c_t'+1,t+1q_t'+1^Tq̃_t+1/N
    + o(1) 
     d_t'+1,t+1 - ξ̅_tξ̅_t'a_t',t^(0)
    + o(1).

Applying (<ref>) in Theorem <ref> yields 

    d_t'+1,t+1𝔼[
     {f_t'(x_1+z_t') - x_1}{f_t(x_1+z_t) - x_1}],

where {z_t} are zero-mean Gaussian random variables with 
covariance 𝔼[z_t'z_t]=a_t',t^(0). 
Finally, we use (<ref>) in Theorem <ref> to 
obtain  
rl
e_t^(j) 
    ξ̅_t-1(e_t-1^(j) - e_t-1^(j+1))
+ σ^2μ_j+1 + o(1)


   + ∑_τ=0^t-1ξ̅_τ^(t-1)θ_t-τ
(e_τ^(j+1) - ξ̅_τ-1e_τ-1^(j+1))


   - ∑_τ=0^t-1ξ̅_τ^(t-1)g_t-τ
(e_τ^(j) - ξ̅_τ-1e_τ-1^(j)).  


To transform the summations in these equations to convolution, 
we use the change of variables 
a_t',t^(j)=ξ̅_0^(t'-1)ξ̅_0^(t-1)ã_t',t^(j). 
Similarly, we define b̃_t',t^(j), c̃_t',t, and 
d̃_t',t while we use 
e_t'^(j)=ξ̅_0^(t'-1)ξ̅_0^(t-1)ẽ_t',t^(j). 
Then, the SE equations (<ref>)–(<ref>) reduce to 
rl
ã_t',t^(j) 
    b̃_t,t'^(j) - b̃_t,t'^(j+1) 
+ ã_t',t-1^(j) - ã_t',t-1^(j+1)
+ ẽ_t',t^(j) 


   + ∑_τ=0^t-1θ_t-τ
(ã_t',τ^(j+1) - b̃_τ,t'^(j+1) 
- ã_t',τ-1^(j+1)) 
 
   - ∑_τ=0^t-1g_t-τ
(ã_t',τ^(j) - b̃_τ,t'^(j) 
- ã_t',τ-1^(j)) + o(1),  

rl
b̃_t',t^(j) 
    (μ_j - μ_j+1)c̃_t',t  
+ b̃_t',t-1^(j) - b̃_t',t-1^(j+1) + o(1) 


   + ∑_τ=0^t-1θ_t-τ
(b̃_t',τ^(j+1) - μ_j+1c̃_t',τ 
- b̃_t',τ-1^(j+1))


   - ∑_τ=0^t-1g_t-τ
(b̃_t',τ^(j) - μ_jc̃_t',τ 
- b̃_t',τ-1^(j)), 


    c̃_t'+1,t+1d̃_t'+1,t+1 - ã_t',t^(0) + o(1),

rl
ẽ_t',t^(j) 
    ẽ_t'-1,t^(j) - ẽ_t'-1,t^(j+1)
+ μ_j+1σ_t',t^2 + o(1)


   + ∑_τ=0^t'-1θ_t'-τ
(ẽ_τ,t^(j+1) - ẽ_τ-1,t^(j+1))


   - ∑_τ=0^t'-1g_t'-τ
(ẽ_τ,t^(j) - ẽ_τ-1,t^(j)), 


with 

    σ_t',t^2 = σ^2/ξ̅_0^(t'-1)ξ̅_0^(t-1).


In principle, it is possible to solve the coupled dynamical 
system (<ref>), (<ref>)–(<ref>) numerically. 
However, numerical evaluation is a challenging task due to instability 
against numerical errors.  



 §.§ Generating Functions

We solve the coupled dynamical system via the Z-transform. 
Define the generating function of ã_t',t^(j) as  

    A(x,y,z) = ∑_j=0^∞x^jA_j(y,z),

with 

    A_j(y,z) = ∑_t',t=0^∞ã_t',t^(j)y^-t'z^-t.

Similarly, we write the generating functions of {b̃_t',t^(j)}, 
{c̃_t',t}, {d̃_t',t}, {ẽ_t',t^(j)}, and 
{σ_t',t^2} as B(x,y,z), C(y,z), D(y,z), E(x,y,z), and 
Σ(y,z), respectively. 

To evaluate the generating function A_j(y,z), we utilize  
rl
   ∑_t'=0^∞y^-t'∑_t=1^∞z^-t∑_τ=0^t-1g_t-τ
ã_t',τ-k^(j) 

=    ∑_t'=0^∞y^-t'∑_τ=0^∞∑_t=τ+1^∞
z^-tg_t-τã_t',τ-k^(j) 
 
=   z^-k[
 G(z) - 1
]A_j(y,z)   

for any integer k, where we have used the definition (<ref>) 
of G(z). From (<ref>), we have 
l
A_j(y,z) 
B_j(z,y) - B_j+1(z,y) + A_j(y,z)/z - A_j+1(y,z)/z  


+ [
 Θ(z) - 1
]{
 A_j+1(y,z) - B_j+1(z,y) - A_j+1(y,z)/z 
} 

- [
 G(z) - 1
]{
 A_j(y,z) - B_j(z,y) - A_j(y,z)/z
}

 
+ E_j(y,z).


Similarly, we can derive 
l
B_j(y,z) 
(μ_j - μ_j+1)C(y,z)  
+ B_j(y,z)/z - B_j+1(y,z)/z 


+ [
 Θ(z) - 1
]{
 B_j+1(y,z) - μ_j+1C(y,z) - B_j+1(y,z)/z 
}


- [
 G(z) - 1
]{
 B_j(y,z) - μ_jC(y,z) - B_j(y,z)/z 
} + o(1), 
 



    C(y,z)  D(y,z) - (yz)^-1A_0(y,z) + o(1),

rl
E_j(y,z)
    E_j(y,z)/y - E_j+1(y,z)/y 
+ μ_j+1Σ(y,z) + o(1)


   + (1-y^-1)[Θ(y) - 1]E_j+1(y,z) 
 
   - (1-y^-1)[G(y) - 1]E_j(y,z). 


We next substitute (<ref>) into (<ref>) to obtain 
l
{
 xG̃(z) + 1 - Θ̃(z) 
}A(x,y,z) 
[1 - Θ̃(z)]A_0(y,z) 

+ {
 xG̃(z) - Θ̃(z) 
}B(x,z,y)/1-z^-1   
+ xE(x,y,z) + o(1), 

with G̃(z)=(1-z^-1)G(z) and Θ̃(z)=(1-z^-1)Θ(z), 
where we have used the identity B_0(y,z) o(1) obtained from the 
asymptotic orthogonality between b_t' and 
m_t. Similarly, we use (<ref>) and 
(<ref>) to obtain 

    B(x,y,z) 
    [xG̃(z) - Θ̃(z)]η(-x) + Θ̃(z)/xG̃(z) + 1 - Θ̃(z)C(y,z)/1-z^-1
    + o(1).

Furthermore, we have 
rl
E(x,y,z)
    1-Θ̃(y)/xG̃(y) + 1 - Θ̃(y)E_0(y,z) 

 
   + η(-x)-1/xG̃(y) + 1 - Θ̃(y)Σ(y,z) 
+ o(1). 




 §.§ Evaluation at Poles

The equations (<ref>), (<ref>), (<ref>), and (<ref>) 
provide all information about the generating functions. However, 
we are interested only in those at x=0. To extract this information,   
we focus on the poles of A(x,y,z) and E(x,y,z). Let -x^* denote 
the pole of A(x,y,z) given by 

    x^* = 1-Θ̃(z)/G̃(z).

Since A(x,y,z) is analytical, the RHS of (<ref>) has to be 
zero at x=-x^*. 

    B(-x^*,z,y)/1-z^-1 [1-Θ̃(z)]A_0(y,z) - x^*E(-x^*,y,z)
    + o(1).


Similarly, we use (<ref>) and Theorem <ref> 
to obtain  

    E_0(y,z) Σ(y,z) + o(1).

Thus, (<ref>) reduces to  
rl
   E(-x^*,y,z) 

    [Θ̃(z)  - Θ̃(y)]G̃(z)Σ(y,z)/G̃(y)Θ̃(z) - Θ̃(y)G̃(z)
+ G̃(z) - G̃(y)  + o(1).  

 

Evaluating B(x,z,y) given via (<ref>) at x=-x^* yields 
rl
B(-x^*,z,y)/1-z^-1 
    Θ(y)G(z) - G(y)Θ(z)/G̃(y)Θ̃(z) - Θ̃(y)G̃(z) 
+ G̃(z) - G̃(y) 


   ·[1-Θ̃(z)]C(y,z) + o(1),  

where we have used Θ̃(z)=(1-z^-1)Θ(z), 
G̃(z)=(1-z^-1)G(z), and the symmetry C(z,y)=C(y,z).  
Substituting (<ref>), (<ref>), and (<ref>) into (<ref>),  
we obtain 
l 
F_G,Θ(y,z)A_0(y,z) 
Θ(y)G(z) - G(y)Θ(z)/y^-1 - z^-1D(y,z)  

+ (1-z^-1)Θ(z) - (1-y^-1)Θ(y)/y^-1 - z^-1Σ(y,z) 
+ o(1),  

with 
rl
F_G,Θ(y,z) 
=    (y^-1+z^-1-1)[Θ(y)G(z) - G(y)Θ(z)]/y^-1 - z^-1 


   + (1-z^-1)G(z) - (1-y^-1)G(y)/y^-1 - z^-1. 


We transform the SE equation (<ref>) into another 
generating-function representation that is suited for deriving time-domain 
representation. 
Let S denote the generating function of some sequence {s_t}. 
We use the notations S_1(z)=z^-1S(z), Δ_S, and Δ_S_1, 
given by 

    Δ_S 
    = S(y) - S(z)/y^-1-z^-1,

which is a function of y and z. The inverse Z-transform of these 
generating functions can be evaluated straightforwardly, as shown shortly. 
We use these notations to re-write the SE 
equation (<ref>) as  
rl
F_G,Θ(y,z)A_0(y,z) 
    {
 G(z)Δ_Θ - Θ(z)Δ_G 
}D(y,z) 


+    (
 Δ_Θ_1  - Δ_Θ   
)Σ(y,z) + o(1), 


with  
rl
F_G,Θ(y,z)
=    (y^-1+z^-1-1)[G(z)Δ_Θ - Θ(z)Δ_G]

 
   + Δ_G_1 - Δ_G,  

where G_1(z)=z^-1G(z) and Θ_1(z)=z^-1Θ(z) 
are defined in the same manner as in S_1(z). The SE 
equation (<ref>) is equivalent to the former 
statement in Theorem <ref>. 



 §.§ Time-Domain Representation
 
We transform the SE equation (<ref>) into 
a time-domain representation that is suitable for numerical evaluation. 
Suppose that G(z) is represented as G(z)=P(z)/Q(z). 
Let R(z) denote the generating function of {r_t}, i.e. 
R(z)=Q(z)Θ(z). 
We multiply both sides of the SE equation (<ref>) 
by Q(y)Q(z) to obtain  
rl
   F_P,Q,Θ(y,z)A_0(y,z) 
 
    {
 P(z)Δ_R - R(z)Δ_P 
}D(y,z) 


+    Q(y)Q(z)(
 Δ_Θ_1  - Δ_Θ   
)Σ(y,z) + o(1), 


with  
rl
F_P,Q,Θ(y,z)
=    [Δ_P_1 - Δ_P]Q(z) + (1 - z^-1)P(z)Δ_Q 

 
   + (z^-1-1)[P(z)Δ_R - R(z)Δ_P] 

   + y^-1[P(z)Δ_R - R(z)Δ_P].  


It is possible to evaluate the inverse Z-transform of S_1(z), 
Δ_S, Δ_S_1, and z^-1Δ_S 
for any generating function S(z). 
By definition, we have 

    S_1(z) = ∑_t=0^∞s_tz^-(t+1) 
    = ∑_t=0^∞s_t-1z^-t,

where the convention s_-1=0 has been used. Thus, S_1(z) is the 
generating function of the sequence {s_t-1}. 

For Δ_S, we obtain  
rl
   Δ_S 
= ∑_τ=1^∞s_τy^-τ - z^-τ/y^-1 - z^-1
= ∑_τ=1^∞∑_τ'=0^τ-1s_τy^-τ'z^-(τ-τ'-1) 


=    ∑_τ'=0^∞∑_τ=τ'+1^∞s_τy^-τ'z^-(τ-τ'-1)
= ∑_τ'=0^∞∑_τ=0^∞s_τ'+τ+1y^-τ'z^-τ, 


 
which implies that Δ_S is the generating function of 
the two-dimensional array s_t',t=s_t'+t+1. 

We combine these results to evaluate the inverse Z-transform of the remaining  
generating functions. For S_1(z), Δ_S_1 is 
the generating function of {s_t'+t}. 
Since y^-1 is the generating function of 
δ_t',1δ_t,0 and since Δ_S is the generating function 
of s_t',t=s_t'+t+1, y^-1Δ_S is the generating function of 
the two-dimensional convolution:

    (δ_t',1δ_t,0)*s_t',t 
    = s_t'-1,t
    = s_t'+t - δ_t',0s_t,
  
where the last expression is due to the convention s_-1,t=0. 
See Table <ref> for a summary of these results. 



We evaluate the inverse Z-transform of (<ref>). 
It is a simple exercise to confirm that 
(<ref>) is equal to the Z-transform of 
the following difference equation:  
rl
𝔇_t',t*ã_t',t^(0) 
   (p_t*r_t'+t+1 - r_t*p_t'+t+1)*d̃_t',t
 

+    (q_t'q_t)*(θ_t'+t - θ_t'+t+1)*σ_t',t^2
+ o(1),  
   

with 
rl
𝔇_t',t
=    (p_t'+t - p_t'+t+1)*q_t + (p_t - p_t-1)*q_t'+t+1 


+    (p_t-1-p_t)*r_t'+t+1 + (r_t - r_t-1)*p_t'+t+1
 

+    p_t*(r_t'+t - δ_t',0r_t) 
- r_t*(p_t'+t - δ_t',0p_t), 

where all variables with negative indices are set to zero. 
Multiplying (<ref>) by 
ξ̅_0^(t'-1)ξ̅_0^(t-1) and 
using the definitions ã_τ',τ^(0)=a_τ',τ^(0)/
(ξ̅_0^(τ'-1)ξ̅_0^(τ-1)), 
d̃_τ',τ^(0)=d_τ',τ/
(ξ̅_0^(τ'-1)ξ̅_0^(τ-1)), and 
σ_τ',τ^2=σ^2/(ξ̅_0^(τ'-1)ξ̅_0^(τ-1)), 
we arrive at the SE equation (<ref>) in time domain, with 
the superscript in a_τ',τ^(0) omitted. 

Finally, we use the notational convention f_-1(·)=0 to obtain initial 
and boundary conditions. From the definition (<ref>) of 
d_t'+1,t+1, we have the initial condition 
d_0,0=𝔼[x_1^2]=1. Similarly, we use 
(<ref>) to obtain the boundary condition  
d_0,τ+1= -𝔼[x_1{f_τ(x_1+z_τ) - x_1}]. 
The boundary condition d_τ+1,0=d_0,τ+1 follows from the symmetry. 




§ PROOF OF THEOREM <REF>
 
 
Without the loss of generality, we assume p_t=g_t and q_t=δ_t,0. 
Then, the SE equation (<ref>) in time domain reduces to 
rl
∑_τ'=0^t'∑_τ=0^tξ̅_t'-τ'^(t'-1)ξ̅_t-τ^(t-1)
{ 𝔇_τ',τa_t'-τ',t-τ     


- (g_τ*θ_τ'+τ+1 - θ_τ*g_τ'+τ+1)d_t'-τ',t-τ   


- σ^2(
 θ_τ'+τ - θ_τ'+τ+1
)}    
=0,   

with
rl
𝔇_τ',τ
=    g_τ'+τ - g_τ'+τ+1 +(g_τ-1-g_τ)*θ_τ'+τ+1  


   + (θ_τ - θ_τ-1)*g_τ'+τ+1
+ g_τ*(θ_τ'+τ - δ_τ',0θ_τ) 
 

   - θ_τ*(g_τ'+τ - δ_τ',0g_τ). 

 
We evaluate a fixed-point of the reduced SE 
equation (<ref>) for the Bayes-optimal denoiser 
f_opt. 
Suppose that lim_t',t→∞a_t',t=a_s, 
lim_t',t→∞d_t',t=d_s, and 
lim_t→∞ξ̅_t=ξ_s hold. 
The main feature of the Bayes-optimal denoiser is the identity 
ξ_s=d_s/a_s <cit.>. 
We use this identity and the assumptions in Theorem <ref> 
to prove the fixed-point (<ref>). 

Taking the limits t',t→∞ in (<ref>) yields 
rl
   a_s∑_τ',τ=0^∞
𝔇_τ',τ(ξ_s^-1)^-τ'-τ  

=    d_s∑_τ',τ=0^∞
(g_τ*θ_τ'+τ+1 - θ_τ*g_τ'+τ+1)
(ξ_s^-1)^-τ'-τ 

   + σ^2∑_τ',τ=0^∞
(θ_τ'+τ - θ_τ'+τ+1)
(ξ_s^-1)^-τ'-τ. 

We use the properties of the Z-transform in Table <ref> 
and the identity ξ_s=d_s/a_s 
to find 

    F_G,Θ(y,z)d_s/ξ_s
    ={G(z)Δ_Θ - Θ(z)Δ_G}d_s  
    + (Δ_Θ_1 - Δ_Θ)σ^2

in the limit y, z→ξ_s^-1, where F_G,Θ is given by 
(<ref>). 

Series-expanding Δ_S with respect to z^-1 at z=y up to the 
first order yields  

    lim_y,z→ξ_s^-1Δ_S
    = dS/dz^-1(ξ_s^-1).

Similarly, we have 

    lim_y,z→ξ_s^-1Δ_S_1
    = S(ξ_s^-1) + ξ_sdS/dz^-1(ξ_s^-1),

Applying these results to (<ref>) with (<ref>) yields 

    {
     1 + (ξ_s-1)
     dΘ/dz^-1(ξ_s^-1) 
    }{G(ξ_s^-1)d_s/ξ_s
     - σ^2} 
    = 0,

where we have used the assumption Θ(ξ_s^-1)=1. 
Since 1+(ξ_s-1)dΘ(ξ_s^-1)/(dz^-1)≠0 
has been assumed, we arrive at 

    G(ξ_s^-1)/ξ_s 
    = σ^2/d_s.


To prove the fixed-point (<ref>), we use the 
relationship (<ref>) between the η-transform and 
the R-transform. Evaluating (<ref>) at 
x=x^* given in (<ref>) and using Theorem <ref>, 
we obtain 

    G(z) 
    = Θ(z) R(
     -1 - (1-z^-1)Θ(z)/G(z)Θ(z)
    ).

Letting z=ξ_s^-1 and applying the assumption 
Θ(ξ_s^-1)=1 yields 

    G(ξ_s^-1) 
    = R(
     - ξ_s/G(ξ_s^-1)).

Substituting (<ref>) into this identity and using 
ξ_s=d_s/a_s, we arrive at 

    a_s
    = σ^2/R(- d_s/σ^2).




§ EVALUATION OF (<REF>) FOR BERNOULLI-GAUSSIAN SIGNALS
 



 §.§ Summary

We evaluate the correlation (<ref>) for the Bernoulli-Gaussian 
signals. This appendix is organized as an independent section of the other 
parts. Thus, we use different notations from the other parts. 

Let A∈{0, 1} denote a Bernoulli random variable taking 1 with 
probability ρ∈[0, 1]. Suppose that Z∼𝒩(0,ρ^-1) is 
independent of A and a zero-mean Gaussian random variable with 
variance ρ^-1. We consider estimation of a Bernoulli-Gaussian signal 
X=AZ on the basis of two dependent noisy observations,

    Y_t'=X+W_t',    Y_t=X+W_t,
 
with 

    [ W_t';  W_t ]∼𝒩(0,Σ), 
      Σ
    =[ a_t',t'  a_t',t;  a_t',t   a_t,t ],

where Σ is positive definite. 
The goal of this appendix is to evaluate the correlation d_t'+1,t+1 of 
the estimation errors for the Bayes-optimal denoiser 
f_opt(Y_t;a_t,t)=𝔼[X|Y_t],   

    d_t'+1,t+1=𝔼[{f_opt(Y_t';a_t',t') - X}{f_opt(Y_t;a_t,t) - X}].


Before presenting the derived expression of the 
correlation (<ref>), we first introduce 
several definitions. We write the pdf of a zero-mean Gaussian random variable 
Y with variance σ^2 as p_G(y;σ^2), with 

    p_G(y; σ^2) 
    = 1/√(2πσ^2)exp(
     -y^2/2σ^2).

The pdf of a Gaussian mixture is defined as 

    p_GM(y;a_t,t) 
    = ρ p_G(y;ρ^-1+a_t,t)
    + (1-ρ)p_G(y;a_t,t),

which is used to represent the marginal pdf of Y_t. As proved in 
Appendix <ref>, the probability of 
A=1 given Y_t is given by 
Pr(A=1|Y_t=y)=π(y;a_t,t), with 

    π(y,a_t,t) = ρ p_G(y;ρ^-1+a_t,t)/p_GM(y;a_t,t).

The Bayes-optimal denoiser f_opt(Y_t;a_t,t) 
is derived in the same appendix: 

    f_opt(y;a_t,t) 
    = y/1+ρ a_t,tπ(y,a_t,t),

where the conditional probability π(y,a_t,t) is 
given by (<ref>). 

We write the MSE function MSE(a_t,t) as 
rl
   MSE(a_t,t) 

=    a_t,t/1+ρa_t,t 
+ 𝔼[{1-π(Y_t,a_t,t)}{f_opt(Y_t;a_t,t)}^2]  


+    𝔼[
 π(Y_t,a_t,t)
 {
  Y_t/1+ρa_t,t 
  - f_opt(Y_t;a_t,t)
 }^2
],  

where the Bayes-optimal denoiser f_opt is given in 
(<ref>). 
In (<ref>), the expectation  is over 
Y_t∼ p_GM(y;a_t,t) given in (<ref>).  


The joint pdf of {Y_t', Y_t} is represented as  

    p(Y_t', Y_t)=ρ p(Y_t', Y_t|A=1) + (1-ρ)p(Y_t', Y_t|A=0).

As proved in Appendix <ref>, 
the conditional pdf p(Y_t', Y_t|A) is given by
rl
   p(Y_t', Y_t|A=a) 
= p_G(Y_t';ρ^-1a+a_t',t') 

   ·p_G(
 Y_t - a+ρa_t',t/a+ρa_t',t'Y_t'; 
 a+ρa_t,t/ρ - (a+ρa_t',t)^2/ρ(a+ρa_t',t')
) 



for a=0, 1. 

 


  * Let MSE(a_t,t) denote the MSE function (<ref>). Then,  

    d_t+1,t+1 = MSE(a_t,t).



  * For t'≠ t, let 

    v_t',t = a_t',t'a_t,t-a_t',t^2/a_t',t' + a_t,t - 2a_t',t.

Then, the correlation d_t'+1,t+1 for t'≠ t is given by 
l
d_t'+1,t+1
= 𝔼[f_opt(Y_t';a_t',t')f_opt(Y_t;a_t,t)] 


+ 𝔼[
 π(Y_t',t;v_t',t){
  (
   Y_t',t/1 + ρv_t',t
  )^2
  + ρ^-1v_t',t/ρ^-1+v_t',t 
 .
. 

.
 .
  - Y_t',t[f_opt(Y_t';a_t',t')+f_opt(Y_t;a_t,t)]/1 + ρv_t',t 
 }
],  


with 

    Y_t',t = (a_t,t - a_t',t)Y_t' + (a_t',t' - a_t',t)Y_t/a_t',t' + a_t,t - 2a_t',t,

where the expectation in (<ref>) over {Y_t',Y_t} is evaluated 
via the joint pdf (<ref>).  



See from Appendix <ref> to Appendix <ref>. 


Proposition <ref> implies that d_t'+1,t+1 for t'≠ t 
requires numerical computation of the double integrals.  



 §.§ Bayes-Optimal Denoiser
  
We compute the Bayes-optimal denoiser 
f_opt(Y_t;a_t,t)=𝔼[X|Y_t], given by 
rl
f_opt(Y_t;a_t,t) 
=    𝔼[
 .
  𝔼[AZ|Y_t, A]
 | Y_t 
] 

=    𝔼[Z|Y_t, A=1]Pr(A=1 | Y_t). 

Note that f_opt is different from the true posterior mean estimator 
(PME) 𝔼[X|Y_t', Y_t]. 

We first evaluate the former factor 𝔼[Z|Y_t, A=1]. 
Since Y_t=Z+W_t given A=1 is the AWGN observation of 
Z∼𝒩(0,ρ^-1), we obtain the well-known LMMSE estimator 

    𝔼[Z|Y_t, A=1] 
    = ρ^-1Y_t/ρ^-1+a_t,t,

which implies the Bayes-optimal denoiser (<ref>). 

We next prove that the latter factor Pr(A=1 | Y_t) is equal to 
π(Y_t;a_t,t) given in (<ref>). By definition,  

    Pr(A=1 | Y_t)
    = ρ p(Y_t|A=1)/p(Y_t).

For the numerator, we have  
rl
   p(Y_t|A=1) 
= 𝔼_Z[p(Y_t|A=1, Z)]

=    𝔼_Z[p_G(Y_t-Z;a_t,t)]  
= p_G(Y_t;ρ^-1+a_t,t), 

where the last equality follows from the fact that 
Z+W_t is a zero-mean Gaussian random variable with 
variance ρ^-1+a_t,t. 

The denominator p(Y_t) is computed in the same manner,   
rl
   p(Y_t) = ρp(Y_t | A=1) + (1-ρ)p(Y_t | A=0) 

=    ρp_G(Y_t;ρ^-1+a_t,t)
+ (1-ρ)p_G(Y_t;a_t,t), 

which is equal to p_GM(Y_t;a_t,t) given in (<ref>). 
Combining these results, we arrive at Pr(A=1|Y_t)
=π(Y_t;a_t,t) given in (<ref>). 



 §.§ MSE
 
To evaluate the MSE 
d_t+1,t+1=𝔼[{X-f_opt(Y_t;a_t,t)}^2], we focus on 
the posterior variance  
𝔼[{X-f_opt(Y_t;a_t,t)}^2| Y_t]. 
By definition, 
rl
   𝔼[{X-f_opt(Y_t;a_t,t)}^2| Y_t] 

=    Pr(A=1|Y_t)
𝔼[{Z-f_opt(Y_t;a_t,t)}^2| Y_t, A=1] 

   + {1 - Pr(A=1|Y_t)}{f_opt(Y_t;a_t,t)}^2, 

with Pr(A=1|Y_t)=π(Y_t,a_t,t) given in (<ref>). 

Let 𝔼[Z|Y_t, A=1] denote the PME of Z conditioned on 
Y_t and A=1, given in (<ref>). 
The conditional expectation in the first term can be evaluated as follows: 
rl
   𝔼[{Z-f_opt(Y_t;a_t,t)}^2| Y_t, A=1] 
= 𝔼[{Z-𝔼[Z| Y_t, A=1] 

   .
  .
   .
   + 𝔼[Z|Y_t, A=1 ] -f_opt(Y_t;a_t,t)
   }^2
 | Y_t, A=1 
] 

   = ρ^-1a_t,t/ρ^-1+a_t,t 
+ {
 Y_t/1+ρa_t,t - f_opt(Y_t;a_t,t)
}^2. 

Combining these results and taking the expectation over 
Y_t∼ p(Y_t)=p_GM(Y_t;a_t,t) given in (<ref>), 
we arrive at the MSE (<ref>).  



 §.§ Sufficient Statistic

As a preliminary step for computing the correlation (<ref>) 
for t'≠ t, 
we derive a sufficient statistic of X based on 
the two correlated observations {Y_t', Y_t}. 

Let Σ^-1/2 denote a square root of 
Σ^-1, 
i.e. (Σ^-1/2)^2=Σ^-1. Applying 
the noise whitening filter Σ^-1/2 to 
the observation vector (Y_t', Y_t)^T yields 

    Σ^-1/2[ Y_t';  Y_t ] 
    = Σ^-1/21_2X
    + Σ^-1/2[ W_t';  W_t ],

with 1_2=(1, 1)^T. 
Note that the effective noise vector—the second term on the RHS—follows 
the standard Gaussian distribution. It is well-known that the MF output is 
a sufficient statistic of X when the effective noise vector has zero-mean 
i.i.d. Gaussian elements. Applying the MF  
(Σ^-1/21_2)^T
/1_2^TΣ^-11_2 
to (<ref>), 
we arrive at a sufficient statistic Y_t',t, given by 

    Y_t',t = 1_2^TΣ^-1/1_2^TΣ^-11_2[ Y_t';  Y_t ] 
    = X + W_t',t,

with 

    W_t',t = 1_2^TΣ^-1/1_2^TΣ^-11_2[ W_t';  W_t ].

It is straightforward to confirm that 
the sufficient statistic (<ref>) reduces to 
(<ref>). Furthermore, we find W_t',t∼𝒩(0,v_t',t), 
with v_t',t=(1_2^TΣ^-11_2)^-1, which reduces to 
(<ref>). 



 §.§ Correlation
  
To evaluate the correlation (<ref>) for t'≠ t, 
we first derive a few quantities associated with the sufficient 
statistic (<ref>). 

The probability of A=1 given Y_t' and Y_t is equal to that of A=1 
given the sufficient statistic (<ref>). Thus, 
repeating the derivation of Pr(A=1|Y_t)=π(Y_t;a_t,t) given 
in (<ref>), we have 

    Pr(A=1|Y_t',Y_t) = π(Y_t',t; v_t',t),
 
where Y_t',t and v_t',t are given by (<ref>) and (<ref>). 
Similarly, repeating the derivation of (<ref>) implies that 
the PME 𝔼[Z|Y_t', Y_t, A=1] reduces to 

    𝔼[Z|Y_t', Y_t, A=1] 
    = Y_t',t/1 + ρ v_t',t.

Furthermore, the true PME 𝔼[X|Y_t',Y_y] is given by 

    𝔼[X|Y_t',Y_t]
    = f_opt(Y_t',t; v_t',t).


We next evaluate the posterior covariance 
rl
    𝔼[{f_opt(Y_t';a_t',t') - X}
{f_opt(Y_t;a_t,t) - X} | Y_t', Y_t] 

=    Pr(A=0|Y_t', Y_t)f_opt(Y_t';a_t',t')
f_opt(Y_t;a_t,t) 

   + Pr(A=1|Y_t', Y_t)𝔼[{f_opt(Y_t';a_t',t') - Z}


   ·{f_opt(Y_t;a_t,t) - Z} | Y_t', Y_t, A=1 ].  
 

Substituting (<ref>) into (<ref>) and 
using f_opt(Y_τ;a_τ,τ) - Z
= {f_opt(Y_τ;a_τ,τ) - 𝔼[Z|Y_t',Y_t,A=1]} 
+ {𝔼[Z|Y_t',Y_t,A=1] - Z} with (<ref>) for 
τ=t', t, we have 
rl
   𝔼[{f_opt(Y_t';a_t',t') - X}
{f_opt(Y_t;a_t,t) - X} | Y_t', Y_t] 

=    {1 - π(Y_t',t;v_t',t)}f_opt(Y_t';a_t',t')
f_opt(Y_t;a_t,t) 

   + π(Y_t',t;v_t',t){
 [
  f_opt(Y_t';a_t',t') - Y_t',t/1 + ρv_t',t
 ]
. 

   .
 ·[
  f_opt(Y_t;a_t,t) - Y_t',t/1 + ρv_t',t
 ]
 + ρ^-1v_t',t/ρ^-1+v_t',t
},   

where Y_t',t is computed with {Y_t', Y_t}, as given in (<ref>). 

Finally, we derive the correlation (<ref>). 
Taking the expectation of the posterior 
covariance (<ref>) over Y_t' and Y_t,  
we arrive at (<ref>).  



 §.§ Joint pdf
 
To compute the expectation in (<ref>), 
we need the conditional pdf p(Y_t', Y_t|A) 
in the joint pdf (<ref>) of {Y_t', Y_t}.  

We first evaluate the conditional distribution of W_t given W_t'. 
Let 

    W_t = α W_t' + √(β)W̃,

with some constants α∈ℝ and β>0,  
where W̃ is a standard Gaussian random variable independent of 
W_t'. Computing the correlation 𝔼[W_t'W_t] and 
variance 𝔼[W_t^2], we obtain 

    𝔼[W_t'W_t]
    = α𝔼[W_t'^2],


    𝔼[W_t^2] 
    = α^2𝔼[W_t'^2] + β.

We use the definitions 𝔼[W_τ^2]=a_τ,τ for τ=t', t 
and 𝔼[W_t'W_t]=a_t',t to have  
α=a_t',t/a_t',t' and β= a_t,t - a_t',t^2/a_t',t'. 
Thus, (<ref>) implies  

    W_t conditioned on W_t'∼𝒩(
     a_t',tW_t'/a_t',t', 
     a_t,t - a_t',t^2/a_t',t').


We next evaluate the conditional pdf p(Y_t', Y_t|A) for A=0.  
Since Y_τ=W_τ holds for A=0, we have 
rl
   p(Y_t', Y_t|A=0)
= p(W_t'=Y_t', W_t=Y_t)  

=    p_G(
 Y_t - a_t',t/a_t',t'Y_t';  
 a_t,t - a_t',t^2/a_t',t'
)p_G(Y_t'; a_t',t'). 

For A=1, we use Y_τ=Z+W_τ to find that 
{Y_t', Y_t} given A=1 are zero-mean Gaussian random variables 
with covariance, 

    𝔼[Y_τ^2 | A=1] = ρ^-1 + a_τ,τ  for τ=t', t,
  

    𝔼[Y_t'Y_t | A=1] = ρ^-1 + a_t',t.

Repeating the derivation of (<ref>), we obtain 
rl
   p(Y_t', Y_t|A=1) 
= p_G(Y_t';ρ^-1+a_t',t') 

   ·p_G(
 Y_t - ρ^-1+a_t',t/ρ^-1+a_t',t'Y_t'; 
 ρ^-1+a_t,t - (ρ^-1+a_t',t)^2/ρ^-1+a_t',t'
). 


Combining these results, we arrive at the conditional 
pdf (<ref>). 




§ ACKNOWLEDGMENT

The author thanks the anonymous reviewers for their suggestions that have 
improved the quality of the manuscript greatly. 























IEEEtran
























































