Orthogonal AMP
Junjie Ma, Li Ping
I Introduction
Consider the signal recovery problem for the following linear model:
Except when is Gaussian or for very small and , 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 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 while SE does not . When 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 is IID Gaussian, AMP is Bayes-optimal provided that the fixed-point of SE is unique.
The SE framework of AMP works with any . Such 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 represents a channel matrix). AMP has also been investigated for decoding sparse regression codes , which have theoretically capacity approaching performances.
The IID assumption for is crucial to the SE of AMP . When is not IID (especially when 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 has IID entries .
II-B State Evolution for AMP
Strictly speaking, (4) is not an algorithm since it involves 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 has IID Gaussian entries, SE can accurately characterize AMP, as shown in Theorem 1 below.
To see the implication of Theorem 1, let in (6). Then, Theorem 1 says that the empirical mean square error (MSE) of AMP defined by
converges to the predicted MSE (where is obtained using SE) defined by
II-C Limitation of AMP
The assumption that 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, in (9) is not a component-wise function as required in AMP. However, if Theorem 1 holds, will converge to a constant independent of each individual . In this case, is an approximate component-wise function and .
where is the thresholding function (which is commonly used in sparse signal recovery algorithms ) given in (47) with . A family of is obtained by changing . In particular, reduces to the soft-thresholding function when . We define a measure of the SE accuracy (after a sufficient number of iterations) as
By changing from 0 to 1, we obtain a family of . The solid line in Fig. 1 shows defined in (10) against for being IID Gaussian. We can see that SE is quite accurate in the whole range of shown (with ), 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 .
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 is still large for , which can be explained by the fact the Onsager term was ignored above. Interestingly, we can see that is very small at , where the Onsager term vanishes for the related 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
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 is the signal to be estimated and is independent of . For this additive Gaussian noise model, we define divergence-free estimator (or a divergence-free function of ) as follows.
A divergence-free function can be constructed as
where is an arbitrary function and an arbitrary constant.
III-C OAMP Algorithm
Starting with , OAMP proceeds as
where is de-correlated and is divergence-free. In the final stage, the output is
where is not necessarily divergence-free.
OAMP is different from the standard AMP in the following aspects:
In (19a), the function is restricted to be divergence-free. Consequently, the Onsager term vanishes.
We will show that, under certain assumptions, restricting to be de-correlated and 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 and . We can write an error recursion for OAMP (similar to that for AMP in (4)) as
where . Two error measures are introduced:
The SE for OAMP is defined by the following recursion
where is independent of . 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 .
in (21a) consists of IID zero-mean Gaussian entries independent of .
in (21b) consists of IID entries independent of and .
According to our earlier assumption below (1), is IID and independent of and . In OAMP, , so Assumption 2 holds for . 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 and a divergence-free , 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 from based on Assumption 2 (for ). 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 to be de-correlated, as shown below.
Suppose that Assumption 2 holds and is unitarily-invariant. If is de-correlated, then the entries of are uncorrelated with those of . Furthermore, the entries of 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 is unitarily invariant. Examples of such include the IID Gaussian matrix ensemble and the partial orthogonal ensemble . Note that there is no restriction on the eigenvalues of . Thus, OAMP is potentially applicable to a wider range of than AMP.
We can meet the de-correlated constraint using (14), in which 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 is independent of , then it is also independent of and , which can be seen from the Markov chain . Thus it is sufficient to ensure the independency between and . Similar to the discussion in Section III-F, we reduce our goal to ensuring orthogonality between and .
Suppose that Assumption 1 holds, we can construct an approximate divergence-free function 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 is a divergence-free function, then
where . Combining (28) with Definition 2, we arrive at (26). ∎
where . In (29), and 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 has zero mean.) Thus the divergence-free constrain on the NLE is to establish orthogonality between and .
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
Note that in (30) can be negative. We may use as a practical estimator for , where is a small positive constant. (Setting may cause a stability problem.)
Given , 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 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 and () denote the th diagonal entries of () and (), respectively. In (32), we define for ).
In (32), is for fixed and . Now, following , assume that the empirical cumulative distribution function (cdf) of , denoted by
converges to a limiting distribution when with a fixed ratio. Furthermore, assume that can be generated from as with a real-valued function. Then, (32) converges to
IV-B Optimal Structure of OAMP
The optimal and that minimize and 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 and in Lemma 1 are used, and are monotonically decreasing sequences. Furthermore, the stationary value of , denoted by , satisfies the following equation
Eqn. (41) is consistent with the fixed-point equation characterization of the MMSE performance for (1) (with 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 . Fig. 2 compares simulated MSE with SE prediction for OAMP and AMP. We first assume that the entries of are independently BPSK modulated, so 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 . 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 is s sparsity level and 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 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 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 since the former involves more iterations in this example.) OAMP-PINV has worse performance than AMP when but performs reasonably well for large . 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 consists of uniformly randomly selected rows of the identity matrix and 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 , 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 iterations. Even when the number of iterations of AMP is increased to , 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 is sufficiently large ( in this case). It should be noted that, when is larger, a smaller 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 generated by the NLE is generally correlated with , which may lead to the correlation between and . We will see below that a de-correlated LE can suppress this correlation.
where and denote the th diagonal entries of and , respectively. (We define for ). For a Haar distributed matrix , we have [43, Lemma 1.1 and Proposition 1.2]
From Assumption 1, is independent of (and so ). Then,
From (21a), to prove is uncorrelated with , we only need to prove is uncorrelated with since is independent of . This can be verified as
Following similar procedures, we can also verify that (i) the entries in are uncorrelated, and (ii) the entries of have identical variances. We omit the details here.
Appendix B Proof of Lemma 1
We can rewrite in (32) as
We now prove that in Lemma 1 is optimal for (54). To this end, define , . Applying the Cauchy-Schwarz inequality
where the right hand side of (56) is invariant to . The minimum in (56) is reached when
where 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 .
The following holds for any divergence-free function
Therefore, to prove Lemma 3, we only need to prove
Substituting into (65) yields
Since is a divergence-free function of , we have the following from (26)
Substituting (67) into (66), proving Lemma 3 becomes proving
We next prove the optimality of based on Lemma 3. Again, let be an arbitrary divergence-free function of . The estimation MSE of reads
Appendix C Proof of Lemma 2
We first verify the monotonicity of . From (39a) and after some manipulations, we obtain
To show the monotonicity of , we only need to show that
The derivative of can be computed based on the definition below (39). After some manipulations, the inequality in (76) becomes the inequality below
The monotonicity of can be proved in a similar way. Again, we only need to prove that
Note that . From [36, Proposition 9], we have
Appendix D Proof of Theorem 3
We first show that decrease monotonically. From (39b),
where (81d) is from the initialization of the SE. Since and is a monotonically increasing function, we have .
We now proceed by induction. Suppose that . Since both and are monotonically increasing, we have , which, together with the SE relationship , leads to . Hence, is a monotonically decreasing sequence.
The monotonicity of the sequence follows directly from the monotonicity of , the SE , and the fact that is a monotonically increasing function.
D-B Fixed Point Equation of SE
where denotes the -transform. For convenience, we further rewrite (83) as
where . Note the following relationship between the -transform and the -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