Orthogonal AMP

Junjie Ma, Li Ping

I Introduction

Consider the signal recovery problem for the following linear model:

Except when PX(x)P_{X}(x) is Gaussian or for very small MM and NN, finding the optimal solution to (1) (under, e.g., the minimum mean-squared error (MMSE) criterion ) can be computationally prohibitive. Approximate message passing (AMP) offers a computationally tractable option. AMP involves the iteration between two modules: one for linear estimation (LE) based on (1a) and the other for symbol-by-symbol non-linear estimation (NLE) based on (1b). An Onsager term is introduced to regulate the correlation problem during iterative processing.

When A\bm{A} contains zero-mean IID Gaussian (or sub-Gaussian) entries, the dynamical behavior of AMP can be characterized by a simple scalar recursion, referred to as state evolution (SE) . The latter bears similarity to density evolution (including EXIT analysis ) for message passing decoding algorithms. However, the underlying assumptions are different: density evolution requires sparsity in A\bm{A} while SE does not . When A\bm{A} is IID Gaussian, it is shown in that the fixed-point equation of the SE for AMP coincides with that of the MMSE performance for a large system. (The latter can be obtained using the replica method .) This implies that, when A\bm{A} is IID Gaussian, AMP is Bayes-optimal provided that the fixed-point of SE is unique.

The SE framework of AMP works with any PX(x)P_{X}(x). Such PX(x)P_{X}(x) can be the distribution of, e.g., amplitude or phase modulation that is widely used signal transmission. For this reason, AMP is also suitable for communication applications such as massive MIMO detection , and millimeter wave channel estimation (in which A\bm{A} represents a channel matrix). AMP has also been investigated for decoding sparse regression codes , which have theoretically capacity approaching performances.

The IID assumption for A\bm{A} is crucial to the SE of AMP . When A\bm{A} is not IID (especially when A\bm{A} is ill-conditioned), the accuracy of SE is not warranted and AMP may perform poorly . Various algorithms have been proposed to handle more general matrices , but most of the existing algorithms lack accurate SE characterization. An exception is the work in , which considers a closely related problem and uses a method different from this paper.

The work in this paper is motivated by our observation that, the SE for AMP is still relatively reliable for a wider family of matrices other than IID Gaussian ones when the Onsager term is small. Our contributions are summarized below.

We propose a modified AMP algorithm comprising of a de-correlated LE and a divergence-free NLEThe name is from , although the discussions therein are irrelevant to this paper.. The proposed algorithm allows LE structures beyond MF, such as pseudo-inverse (PINV) and linear MMSE (LMMSE). OAMP extends and provides new interpretations of our previous work in .

We derive an SE procedure for OAMP, which is accurate if the errors are independent during the iterative process. Independency, however, is a tricky condition. We will show that the use of a de-correlated LE and a divergence-free NLE makes the errors statistically orthogonal, hence the name orthogonal AMP (OAMP). Intuitively, such orthogonality partially satisfies the independency requirement. Our numerical results indicate that the SE predictions are reliable for various matrix ensembles (e.g., IID Gaussian, partial orthogonal and some ill-conditioned ones for which AMP does not work well) and also for various LE structures as mentioned above. Thus OAMP may have wider applications than AMP.

We derive optimal choices within the OAMP framework. We find that the fixed-point characterization of the SE is consistent with that of the optimal MMSE performance obtained by the replica method. This implies the potential optimality of OAMP. Compared with AMP, our result holds for the more general unitarily-invariant matrix ensemble.

We will provide numerical results to show that, compared with AMP, OAMP can achieve better MSE performance as well as faster convergence speed for ill-conditioned matrices. We will demonstrate the excellent performance of OAMP in communication systems with non-sparse binary phase shift keying (BPSK) signals as well as conventional sparse signals.

After we posted the preprint of this work , a proof was given for the state evolution of an OAMP related algorithm in systems involving unitarily-invariant matrices .

Part of the results in this paper have been published in . In this paper, we provide more detailed analysis and numerical results.

II AMP

The use of the Onsager term is the key to AMP. It regulates correlation during iterative processing and ensures the accuracy of SE when A\bm{A} has IID entries .

II-B State Evolution for AMP

Strictly speaking, (4) is not an algorithm since it involves x\bm{x} that is to be estimated. Nevertheless, (4) is convenient for the analysis of AMP discussed below.

The SE for AMP refers to the following recursion:

When A\bm{A} has IID Gaussian entries, SE can accurately characterize AMP, as shown in Theorem 1 below.

To see the implication of Theorem 1, let ψ(h,x)≡[ηt(x+h)−x]2\psi(h,x)\equiv\left[\eta_{t}(x+h)-x\right]^{2} in (6). Then, Theorem 1 says that the empirical mean square error (MSE) of AMP defined by

converges to the predicted MSE (where τt\tau_{t} is obtained using SE) defined by

II-C Limitation of AMP

The assumption that A\bm{A} contains IID entries is crucial to theorem 1. For other matrix ensembles, SE may become inaccurate. Here is an example. Consider the following function for the NLE in AMPStrictly speaking, ηt\eta_{t} in (9) is not a component-wise function as required in AMP. However, if Theorem 1 holds, ∑j=1Nη^t′(rjt)/N\sum_{j=1}^{N}\hat{\eta}_{t}^{\prime}(r_{j}^{t})/N will converge to a constant independent of each individual rjtr_{j}^{t}. In this case, ηt\eta_{t} is an approximate component-wise function and ∑j=1Nηt′(rjt)/N≈β⋅∑j=1Nη^t′(rjt)/N\sum_{j=1}^{N}{\eta}_{t}^{\prime}(r_{j}^{t})/N\approx\beta\cdot\sum_{j=1}^{N}\hat{\eta}_{t}^{\prime}(r_{j}^{t})/N.

where η^t\hat{\eta}_{t} is the thresholding function (which is commonly used in sparse signal recovery algorithms ) given in (47) with γt=1\gamma_{t}=1. A family of ηt\eta_{t} is obtained by changing β\beta. In particular, ηt\eta_{t} reduces to the soft-thresholding function η^t\hat{\eta}_{t} when β=1\beta=1. We define a measure of the SE accuracy (after a sufficient number of iterations) as

By changing β\beta from 0 to 1, we obtain a family of ηt\eta_{t}. The solid line in Fig. 1 shows EE defined in (10) against β\beta for A\bm{A} being IID Gaussian. We can see that SE is quite accurate in the whole range of β\beta shown (with E<10−2E<10^{-2}), which is consistent with the result in Theorem 1.

Clearly, this is inconsistent with the SE in (5a). The problem is caused by the discrepancy in eigenvalue distributions: (11) above is derived from the eigenvalue distribution of a partial DCT matrix while (5a) from that of an IID Gaussian A\bm{A}.

How about replacing (5a) by (11) for the partial DCT matrix? This is shown by the solid line with triangle markers in Fig. 1. We can see that EE is still large for β>0\beta>0, which can be explained by the fact the Onsager term was ignored above. Interestingly, we can see that EE is very small at β=0\beta=0, where the Onsager term vanishes for the related ηt\eta_{t} in (9). This observation motivates the work presented below.

III Orthogonal AMP

In this section, we first introduce the concepts for de-correlated and divergence-free structures for the LE and NLE. We then discuss the OAMP algorithm and its properties.

The following are some common examples of such W^\hat{\bm{W}}

We will discuss the properties of de-correlated LE in Section III-F later.

III-B Divergence-free Estimator

Consider signal estimation from an observation corrupted by additive Gaussian noise

where X∼PX(x)X\sim P_{X}(x) is the signal to be estimated and is independent of Z∼N(0,1)Z\sim\mathcal{N}(0,1). For this additive Gaussian noise model, we define divergence-free estimator (or a divergence-free function of RR) as follows.

A divergence-free function η\eta can be constructed as

where η^\hat{\eta} is an arbitrary function and CC an arbitrary constant.

III-C OAMP Algorithm

Starting with s0=0\bm{s}^{0}=\mathbf{0}, OAMP proceeds as

where Wt\bm{W}_{t} is de-correlated and ηt\eta_{t} is divergence-free. In the final stage, the output is

where ηtout\eta_{t}^{{\rm{out}}} is not necessarily divergence-free.

OAMP is different from the standard AMP in the following aspects:

In (19a), the function ηt\eta_{t} is restricted to be divergence-free. Consequently, the Onsager term vanishes.

We will show that, under certain assumptions, restricting Wt\bm{W}_{t} to be de-correlated and ηt\eta_{t} to be divergence-fee ensure the orthogonality between the input and output “error” terms for both LE and NLE. The name “orthogonal AMP” comes from this fact.

III-D OAMP Error Recursion and SE

Similar to (3), define the error terms as ht≡rt−x\bm{h}^{t}\equiv\bm{r}^{t}-\bm{x} and qt≡st−x\bm{q}^{t}\equiv\bm{s}^{t}-\bm{x}. We can write an error recursion for OAMP (similar to that for AMP in (4)) as

where Bt≡I−WtA\bm{B}_{t}\equiv\bm{I}-\bm{W}_{t}\bm{A}. Two error measures are introduced:

The SE for OAMP is defined by the following recursion

where X∼PX(x)X\sim P_{X}(x) is independent of Z∼N(0,1)Z\sim\mathcal{N}(0,1). Also, at the final stage, the MSE is predicted as

III-E Rationales for OAMP

It is straightforward to verify that the SE in (23) is consistent with the error recursion in (21), provided that the following two assumptions hold for every tt.

ht\bm{h}^{t} in (21a) consists of IID zero-mean Gaussian entries independent of x\bm{x}.

qt+1\bm{q}^{t+1} in (21b) consists of IID entries independent of A\bm{A} and n\bm{n}.

According to our earlier assumption below (1), x\bm{x} is IID and independent of A\bm{A} and n\bm{n}. In OAMP, q0=−x\bm{q}^{0}=-\bm{x}, so Assumption 2 holds for t=−1t=-1. Thus the two Assumptions will hold if we can prove that they imply each other in the iterative process. Unfortunately, so far, we cannot.

Assumptions 1 and 2 are only sufficient conditions for the SE. Even if they do not hold exactly, the SE may still be valid. In Section V, we will show that the SE for OAMP is accurate for a wide range of sensing matrices using simulation results. In the following two subsections, we will see that, with a de-correlated Wt\bm{W}_{t} and a divergence-free ηt\eta_{t}, Assumptions 1 and 2 can partially imply each other. We emphasize that the discussions below are to provide intuitions for OAMP, which are by no means rigorous.

III-F Intuitions for the LE Structure

Eqn. (19a) performs linear estimation of x\bm{x} from y\bm{y} based on Assumption 2 (for qt\bm{q}^{t}). We first consider ensuring Assumption 1 based on Assumption 2. The independence requirements in Assumption 1 are difficult to handle. We reduce our goal to remove the correlation among the variables involved. This is achieved by restricting Wt\bm{W}_{t} to be de-correlated, as shown below.

Suppose that Assumption 2 holds and A\bm{A} is unitarily-invariant. If Wt\bm{W}_{t} is de-correlated, then the entries of ht\bm{h}^{t} are uncorrelated with those of x\bm{x}. Furthermore, the entries of ht\bm{h}^{t} in (21a) are mutually uncorrelated with zero-mean and identical variances.

The name “de-correlated” LE comes from Proposition 1.

A key condition to Proposition 1 is that the sensing matrix A\bm{A} is unitarily invariant. Examples of such A\bm{A} include the IID Gaussian matrix ensemble and the partial orthogonal ensemble . Note that there is no restriction on the eigenvalues of A\bm{A}. Thus, OAMP is potentially applicable to a wider range of A\bm{A} than AMP.

We can meet the de-correlated constraint using (14), in which W^t\hat{\bm{W}}_{t} can be chosen from those in (15). Thus OAMP has more choices for the LE than AMP, which makes the former potentially more efficient.

III-G Intuitions for the NLE Structure

We next consider ensuring Assumption 2 based on Assumption 1. From (21), if qt+1\bm{q}^{t+1} is independent of ht\bm{h}^{t}, then it is also independent of A\bm{A} and n\bm{n}, which can be seen from the Markov chain A,n→ht→qt+1\bm{A},\bm{n}\to\bm{h}^{t}\to\bm{q}^{t+1}. Thus it is sufficient to ensure the independency between qt+1\bm{q}^{t+1} and ht\bm{h}^{t}. Similar to the discussion in Section III-F, we reduce our goal to ensuring orthogonality between qt+1\bm{q}^{t+1} and ht\bm{h}^{t}.

Suppose that Assumption 1 holds, we can construct an approximate divergence-free function ηt\eta_{t} according to (18):

All the numerical results about OAMP shown in Section V are based on (19) and (25).

There is an inherent orthogonality property associated with divergence-free functions.

If η\eta is a divergence-free function, then

where ηt′(X+τtZ)≡ηt′(R)∣R=X+τtZ\eta_{t}^{\prime}(X+\tau_{t}Z)\equiv\eta_{t}^{\prime}(R)|_{R=X+\tau_{t}Z}. Combining (28) with Definition 2, we arrive at (26). ∎

where Rt≡X+τtZR^{t}\equiv X+\tau_{t}Z. In (29), Rt−XR^{t}-X and ηt(Rt)−X\eta_{t}(R^{t})-X represent, respectively, the error terms before and after the estimation. Eqn. (29) indicates that these two error terms are orthogonal. (They are also uncorrelated as Rt−XR^{t}-X has zero mean.) Thus the divergence-free constrain on the NLE is to establish orthogonality between qt+1\bm{q}^{t+1} and ht\bm{h}^{t}.

III-H Brief Summary

If the input and output errors of the LE and NLE are independent of each other, Assumptions 1 and 2 naturally hold. However, independency is generally a tricky issue. We thus turn to orthogonality instead. The name “orthogonal AMP” came from this fact. Propositions 1 and 2 are weaker than Assumptions 1 and 2. Nevertheless, our extensive numerical study (see Section V) indicates that the SE in (23) is indeed reliable for OAMP.

Also note that each of Propositions 1 and 2 depends on one assumption, so they do not ensure orthogonality in the overall process. Nevertheless, we observed from numerical results that the orthogonality property is accurate for with unitarily-invariant matrices.

III-I MSE Estimation

We can adopt the following estimator [34, Eqn. (71)] for vt2v_{t}^{2}

Note that v^t2\hat{v}_{t}^{2} in (30) can be negative. We may use max⁡(v^t2,ϵ)\max(\hat{v}_{t}^{2},\epsilon) as a practical estimator for vt2v_{t}^{2}, where ϵ\epsilon is a small positive constant. (Setting ϵ=0\epsilon=0 may cause a stability problem.)

Given v^t2\hat{v}_{t}^{2}, τt2\tau_{t}^{2} can be estimated using (23a):

In certain cases, Eqn. (31) can be simplified to more concise formulas. For example, (31) simplifies to \hat{\tau}_{t}^{2}=\left({N-M}\right)/M\cdot\hat{v}_{t}^{2}+N/M^{2}\cdot{\rm{tr}}\big{\{}{\left({{\bm{AA}}^{\rm{T}}}\right)^{-1}}\big{\}}\cdot\sigma^{2} when Wt\bm{W}_{t} is given by the PINV estimator in (15b) together with (14). Also, simple closed-form asymptotic expression exists for (31) for certain matrix ensembles. For example, (23a) converges to (42a), (42b) and (42c) for IID Gaussian matrices with MF, PINV and LMMSE linear estimators, respectively.

The numerical results presented in Section V are obtained based on approximations in (30) and (31).

IV Optimization Structures for OAMP

In this section, we derive the optimal LE and NLE structures for OAMP based on SE. We show that OAMP can potentially achieve optimal performance, provided that its SE is reliable.

where λi\lambda_{i} and g^i\hat{g}_{i} (i=1,…,Mi=1,\ldots,M) denote the iith diagonal entries of Σ\bm{\Sigma} (M×NM\times N) and G^t\hat{\bm{G}}_{t} (N×MN\times M), respectively. In (32), we define λi=g^i=0\lambda_{i}=\hat{g}_{i}=0 for i=M+1,…,Ni=M+1,\ldots,N).

In (32), Φt(vt2)\Phi_{t}(v_{t}^{2}) is for fixed {λi}\{\lambda_{i}\} and {g^i}\{\hat{g}_{i}\}. Now, following , assume that the empirical cumulative distribution function (cdf) of {λ12,…,λN2}\{\lambda_{1}^{2},\ldots,\lambda_{N}^{2}\}, denoted by

converges to a limiting distribution when M,N→∞M,N\to\infty with a fixed ratio. Furthermore, assume that g^i\hat{g}_{i} can be generated from λi\lambda_{i} as g^i=g^t(vt2,λi)\hat{g}_{i}=\hat{g}_{t}(v_{t}^{2},\lambda_{i}) with g^t\hat{g}_{t} a real-valued function. Then, (32) converges to

IV-B Optimal Structure of OAMP

The optimal Wt\bm{W}_{t} and ηt\eta_{t} that minimize Φt\Phi_{t} and Ψt\Psi_{t} in (32) and (35) are given by

According to the state evolution process, the final MSE can be expressed as

IV-C Potential Optimality of OAMP

When the optimal {Wt⋆}\{\bm{W}_{t}^{\star}\} and {ηt⋆}\{\eta_{t}^{\star}\} in Lemma 1 are used, {vt2}\{v_{t}^{2}\} and {τt2}\{\tau_{t}^{2}\} are monotonically decreasing sequences. Furthermore, the stationary value of τt2\tau_{t}^{2}, denoted by τ∞2\tau_{\infty}^{2}, satisfies the following equation

Eqn. (41) is consistent with the fixed-point equation characterization of the MMSE performance for (1) (with A\bm{A} being unitarily-invariant) via the replica method [10, Eqn. (17)][21, Eqn. (30)]. This implies that OAMP can potentially achieve the optimal MSE performance. We can see that the de-correlated and divergence-free constraints on LE and NLE, though restrictive, do not affect the potential optimality of OAMP.

V Numerical Study

We start from an IID Gaussian matrix where Ai,j∼N(0,1/M)A_{i,j}\sim\mathcal{N}(0,1/M). Fig. 2 compares simulated MSE with SE prediction for OAMP and AMP. We first assume that the entries of x\bm{x} are independently BPSK modulated, so x\bm{x} is not sparse. This is a typical detection problem in massive MIMO applications. Fig. 2 compares simulated MSEs with SE prediction for OAMP and AMP. In Fig. 2, OAMP-MF, OAMP-PINV and OAMP-LMMSE refer to, respectively, OAMP algorithms with the MF, PINV and LMMSE estimators given in (15) and the normalization in (14). The asymptotic SE formula in (34) becomes, respectively,

where c≡(N−M)/Mc\equiv(N-M)/M. Comparing (42a) and (42b), we see that OAMP-PINV has better interference cancellation property than OAMP-MF (but less robust to noise). This is consistent with the observation in Fig. 2 (which represents a high SNR scenario) that OAMP-PINV can outperform OAMP-MF.

From Fig. 2, we observe good agreement between the simulated and predicted MSE for all curves. Furthermore, we see that AMP has the same convergent value as OAMP-LMMSE for IID Gaussian matrices, while the latter converges faster. Following the approach in , we can prove this observation but the details are omitted due to space limitation.

V-B General Unitarily-invariant Matrix

where ρ∈(0,1]\rho\in(0,1] is s sparsity level and δ(⋅)\delta(\cdot) is the Dirac delta function.

Fig. 3 shows the simulated and predicted MSEs for OAMP for the above ill-conditioned sensing matrix. The SE of OAMP is based on the empirical form in (32) as {λi}\{\lambda_{i}\} are fixed in this example. We can make the following observations.

The performances of AMP and OAMP-MF deteriorate in this case. The SE prediction for AMP is not shown in Fig. 3 since it is noticeably different from the simulation result. (See Fig. 1 for a similar issue.)

The performance of OAMP is strongly affected by the LE structure. OAMP-PINV and OAMP-LMMSE significantly outperform OAMP-MF.

The most interesting point is that the SE in (36) can accurately predict the OAMP simulation results for all the LE structures in Fig. 3. We observed in simulations that such good agreement also holds for LEs beyond the three options shown in Fig. 3.

Fig. 4 compares the MSE performances of AMP, OAMP and genie-aided MMSE (where the positions of the non-zero entries are known) as the condition number of A\bm{A} varies. AMP with adaptive damping (AMP-damping) (based on the Matlab code released by its authorsAvailable at http://sourceforge.net/projects/gampmatlab/ and the parameters used in [17, Fig. 1]) and GAMP-ADMM are also shown. From Fig. 4, we can see that the performance of OAMP-LMMSE is significantly better than those of AMP, AMP-damping and ADMM-GAMP for highly ill-conditioned scenarios. (ADMM-GAMP slightly outperforms OAMP-LMMSE for κ≤100\kappa\leq 100 since the former involves more iterations in this example.) OAMP-PINV has worse performance than AMP when κ≥10\kappa\geq 10 but performs reasonably well for large κ\kappa. OAMP-MF does not work well and thus not included.

For the schemes shown in Fig. 4, AMP have the lowest complexity. OAMP-PINV requires one additional matrix inversion, but it can be pre-computed as it remains unchanged during the iterations. Both OAMP-LMMSE and ADMM-GAMP require matrix inversions in each iteration. As pointed out in , it may be possible to replace the matrix inversion in ADMM-GAMP using an iterative method such as conjugate gradient . Similar approximation should be possible for OAMP as well.

V-C Partial Orthogonal Matrix

Therefore, the complexity of OAMP-LMMSE is the same as AMP.

Unitarily invariant matrices with the partial orthogonality constraint becomes partial Haar-distributed matrices (i.e., uniformly distributed among all partial orthogonal matrices). We next consider the following partial orthogonal matrix

where S\bm{S} consists of MM uniformly randomly selected rows of the identity matrix and U\bm{U} is an Haar-distributed orthogonal matrix. We will also consider deterministic orthogonal matrices, which are important in compressed sensing and found applications in, e.g., MRI . For a partial orthogonal A\bm{A}, the three approaches in Fig. 2, i.e., OAMP-MF, OAMP-PINV and OAMP-LMMSE, become identical. The related complexity is the same as AMP. In this case, the SE equation in (32) becomes

Fig. 5 compares OAMP with AMP in recovering Bernoulli-Gaussian signals with a partial DCT matrix. Following , we will use the empirical phase transition curve (PTC) to characterize the sparsity-undersampling tradeoff. A recovery algorithm “succeeds” with high probability below the PTC and “fails” above it. The empirical PTCs are generated according to [34, Section IV-A]. We see that OAMP considerably outperforms AMP when both algorithms are fixed to 5050 iterations. Even when the number of iterations of AMP is increased to 500500, OAMP still slightly outperforms AMP at relatively high sparsity levels.

Fig. 6 shows the accuracy of SE for OAMP with partial orthogonal matrices. Three matrices are considered: a partial Haar matrix, a partial DCT matrix and a partial Hadamard matrix. From Fig. 6, we see that the simulated MSE performances agree well with state evolution predictions for all the three types of partial orthogonal matrices when NN is sufficiently large (N=8192N=8192 in this case). It should be noted that, when M/NM/N is larger, a smaller NN will suffice to guarantee good agreement between simulation and SE prediction.

The NLEs used in Figs. 2-6 are based on the optimized structure given in Lemma 1. Fig. 7 shows the OAMP SE accuracy with the following soft-thresholding function :

VI Conclusions

AMP performs excellently for IID Gaussian transform matrices. The performance of AMP can be characterized by SE in this case. However, for other matrix ensembles, the SE for AMP is not directly applicable and its performance is not warranted.

In this paper, we proposed an OAMP algorithm based on a de-correlated LE and a divergence-free NLE. Our numerical results indicate that OAMP could be characterized by SE for general unitarily-invariant matrices with much relaxed requirements on the eigenvalue distribution and LE structure. This makes OAMP suitable for a wider range of applications than AMP, especially for applications with ill-conditioned transform matrices and partial orthogonal matrices. We also derived the optimal structures for OAMP and showed that the corresponding SE fixed point potentially coincides with that of the Bayes-optimal performance obtained by the replica method.

VII Acknowledgement

The authors would like to thank Dr. Ulugbek Kamilov and Prof. Phil Schniter for generously sharing their Matlab code for ADMM-GAMP.

Appendix A Proof of Proposition 1

It is seen from (21b) that qt\bm{q}^{t} generated by the NLE is generally correlated with x\bm{x}, which may lead to the correlation between x\bm{x} and ht\bm{h}^{t}. We will see below that a de-correlated LE can suppress this correlation.

where gmg_{m} and λm\lambda_{m} denote the (m,m)(m,m)th diagonal entries of Gt\bm{G}_{t} and Σ\bm{\Sigma}, respectively. (We define gm=λm=0g_{m}=\lambda_{m}=0 for m=M+1,…,Nm=M+1,\ldots,N). For a Haar distributed matrix U\bm{U}, we have [43, Lemma 1.1 and Proposition 1.2]

From Assumption 1, qt\bm{q}^{t} is independent of A\bm{A} (and so Bt\bm{B}_{t}). Then,

From (21a), to prove x\bm{x} is uncorrelated with ht\bm{h}^{t}, we only need to prove x\bm{x} is uncorrelated with Btqt\bm{B}_{t}\bm{q}^{t} since Wtn\bm{W}_{t}\bm{n} is independent of x\bm{x}. This can be verified as

Following similar procedures, we can also verify that (i) the entries in ht\bm{h}^{t} are uncorrelated, and (ii) the entries of ht\bm{h}^{t} have identical variances. We omit the details here.

Appendix B Proof of Lemma 1

We can rewrite Φt(vt2)\Phi_{t}(v_{t}^{2}) in (32) as

We now prove that Wt⋆\bm{W}_{t}^{\star} in Lemma 1 is optimal for (54). To this end, define ai≡g^ivt2λi2+σ2a_{i}\equiv\hat{g}_{i}\sqrt{{v_{t}^{2}\lambda_{i}^{2}+\sigma^{2}}}, bi≡λi/vt2λi2+σ2b_{i}\equiv\lambda_{i}/\sqrt{{v_{t}^{2}\lambda_{i}^{2}+\sigma^{2}}}. Applying the Cauchy-Schwarz inequality

where the right hand side of (56) is invariant to {g^i}\{\hat{g}_{i}\}. The minimum in (56) is reached when

where CC is an arbitrary constant. From (57),

The SE equation in (35) are obtained based on the following signal model

The following identity is from [44, Eqn. (123)]

Lemma 3 below is the key to prove the optimality of ηt⋆\eta_{t}^{\star}.

The following holds for any divergence-free function ηt\eta_{t}

Therefore, to prove Lemma 3, we only need to prove

Substituting Rt=X+τtZR^{t}=X+\tau_{t}Z into (65) yields

Since ηt\eta_{t} is a divergence-free function of RtR^{t}, we have the following from (26)

Substituting (67) into (66), proving Lemma 3 becomes proving

We next prove the optimality of ηt⋆\eta_{t}^{\star} based on Lemma 3. Again, let ηt\eta_{t} be an arbitrary divergence-free function of RtR^{t}. The estimation MSE of ηt\eta_{t} reads

Appendix C Proof of Lemma 2

We first verify the monotonicity of Φ⋆\Phi^{\star}. From (39a) and after some manipulations, we obtain

To show the monotonicity of Φ⋆\Phi^{\star}, we only need to show that

The derivative of mmseA(vt2)mmse_{A}\left({v_{t}^{2}}\right) can be computed based on the definition below (39). After some manipulations, the inequality in (76) becomes the inequality below

The monotonicity of Ψ⋆\Psi^{\star} can be proved in a similar way. Again, we only need to prove that

Note that mmseB(τt2)=E{[X−E{X∣Rt=X+τtZ}]2}mmse_{B}\left({\tau_{t}^{2}}\right)={\rm{E}}\left\{{\left[{X-{\rm{E}}\left\{{X|R^{t}=X{\rm{+}}\tau_{t}Z}\right\}}\right]^{2}}\right\}. From [36, Proposition 9], we have

Appendix D Proof of Theorem 3

We first show that {vt2}\{v_{t}^{2}\} decrease monotonically. From (39b),

where (81d) is from the initialization of the SE. Since Φ⋆(v02)<∞\Phi^{\star}(v_{0}^{2})<\infty and Ψ⋆\Psi^{\star} is a monotonically increasing function, we have v12=Ψ⋆(Φ⋆(v02))<v02v_{1}^{2}=\Psi^{\star}\left(\Phi^{\star}(v_{0}^{2})\right)<v_{0}^{2}.

We now proceed by induction. Suppose that vt2<vt−12v_{t}^{2}<v_{t-1}^{2}. Since both Φ⋆\Phi^{\star} and Ψ⋆\Psi^{\star} are monotonically increasing, we have Ψ⋆(Φ⋆(vt2))<Ψ⋆(Φ⋆(vt−12))\Psi^{\star}\left(\Phi^{\star}(v_{t}^{2})\right)<\Psi^{\star}\left(\Phi^{\star}(v_{t-1}^{2})\right), which, together with the SE relationship vt+12=Ψ⋆(Φ⋆(vt2))v_{t+1}^{2}=\Psi^{\star}\left(\Phi^{\star}(v_{t}^{2})\right), leads to vt+12<vt2v_{t+1}^{2}<v_{t}^{2}. Hence, {vt2}\{v_{t}^{2}\} is a monotonically decreasing sequence.

The monotonicity of the sequence {τt2}\{\tau_{t}^{2}\} follows directly from the monotonicity of {vt2}\{v_{t}^{2}\}, the SE τt2=Φ⋆(vt2)\tau_{t}^{2}=\Phi^{\star}(v_{t}^{2}), and the fact that Φ⋆\Phi^{\star} is a monotonically increasing function.

D-B Fixed Point Equation of SE

where ηATA\eta_{{\bm{A}}^{\rm{T}}{\bm{A}}} denotes the η\eta-transform. For convenience, we further rewrite (83) as

where γ≡vt2/σ2\gamma\equiv v_{t}^{2}/\sigma^{2}. Note the following relationship between the η\eta-transform and the RR-transform [32, Eqn. (2.74)]

where the second equality in (86) is from (36a) and (39a). We can rewrite the SE equations in (39a) and (39b) as follows

Substituting (88) into (86), we get the desired fixed point equation

References