Iterative Reconstruction of Rank-One Matrices in Noise

Introduction

where $\mathbf{W}$ represents some unknown noise and $\sqrt{m}$ is a normalization factor. The problem can be considered as a rank-one special case of finding a low-rank matrix in the presence of noise. Such low-rank estimation problems arise in a range of applications including blind channel estimation , antenna array processing , subspace system identification , and principal component or factor analysis .

When the noise term $\mathbf{W}$ is zero, the vector pair $(\mathbf{u}_{0},\mathbf{v}_{0})$ can be recovered exactly, up to a scaling, from the maximal left and right singular vectors of $\mathbf{A}$ . However, in the presence of noise, the rank-one matrix can in general only be estimated approximately. In this case, a priori information or constraints on $(\mathbf{u}_{0},\mathbf{v}_{0})$ may improve the estimation. Such constraints arise, for example, in factor analysis in statistics, where one of the factors is often constrained to be either positive or sparse . Similar sparsity constraints occur in the problem of dictionary learning . Also, in digital communications, one of the factors could come from a discrete QAM constellation.

In this paper, we impose the constraints on $\mathbf{u}_{0}$ and $\mathbf{v}_{0}$ in a Bayesian setting where $\mathbf{u}_{0}$ and $\mathbf{v}_{0}$ are assumed to be independent from one another with i.i.d. densities of the form

for some scalar random variables $U_{0}$ and $V_{0}$. The noise $\mathbf{W}$ in (1) is assumed to have i.i.d. Gaussian components where $W_{ij}\sim{\mathcal{N}}(0,\tau_{w})$ for some variance $\tau_{w}>0$. Under this model, the posterior density of $\mathbf{u}_{0}$ and $\mathbf{v}_{0}$ is given by

where $\|\cdot\|_{F}$ is the Frobenius norm. Given this posterior density, we consider two estimation problems:

MAP estimation: In this problem, we wish to find the maximum a posteriori estimates

MMSE estimation: In this problem, we wish to find the posterior mean or, equivalently, the minimum mean squared error estimate,

We may also be interested in the posterior variances and posterior marginals.

Exact computation of either of these estimates is generally intractable. Even though the components $\mathbf{u}_{0}$ and $\mathbf{v}_{0}$ are assumed to have independent components, the posterior density (3) is not, in general, separable due to the term $\mathbf{u}_{0}\mathbf{v}_{0}^{T}$. Thus, the MAP estimate requires a search over $m$ and $n$-dimensional vectors, $(\mathbf{u}_{0},\mathbf{v}_{0})$ and the MMSE estimate requires integration over this $m+n$-dimensional space. However, due to the separability assumption on the priors (2), it is computationally simpler to alternately estimate the factors $\mathbf{u}_{0}$ and $\mathbf{v}_{0}$. This form of alternating estimation is used, for example, in the classic alternating power method for finding maximal singular values and some alternating methods in sparse or non-negative dictionary learning .

The approach in this work also uses a form of alternating estimation, but based on a recently-developed powerful class of algorithms known as Approximate Message Passing or AMP. AMP methods are derived from Gaussian and quadratic approximations of loopy belief propagation in graphical models and were originally used for problems in compressed sensing . The AMP methodology has been successfully applied in a range of applications . In recent years, there has been growing interest of AMP for matrix factorization problems as well . The problem considered in this paper can be considered as a simple rank one case of these problems.

Our main contribution in this work is to show that, for the rank one MMSE and MAP estimation problems, the proposed IterFac algorithm admits an asymptotically-exact characterization in the limit when $m,n\rightarrow\infty$ with $m/n$ constant and the components of the true vectors $\mathbf{u}_{0}$ and $\mathbf{v}_{0}$ have limiting empirical distributions. In this scenario, we show that the empirical joint distribution of the components of $\mathbf{u}_{0}$ and $\mathbf{v}_{0}$ and the corresponding estimates from the IterFac algorithm are described by a simple scalar equivalent model where the IterFac component estimates are identically distributed to scalar estimates of the variables corrupted by Gaussian noise. Moreover, the effective Gaussian noise level in this model is described by a simple set of scalar state evolution (SE) equations. This scalar equivalent model is common in analyses of AMP methods as well as replica analyses of related estimation problems . From the scalar equivalent model, one can compute the asymptotic value of almost any component-separable metric including mean-squared error or correlation. Thus, in addition to being computationally simple and general, the IterFac algorithm admits a precise analysis in the case of Gaussian noise. Moreover, since fixed points of the IterFac algorithm correspond, under suitable circumstances, to local maxima of objectives such as (4), the analysis can be used to characterize the behavior of such minima—even if alternate algorithms to IterFac are used.

The main analytical tool is a recently-developed technique by Bayati and Montanari used in the analysis of AMP algorithms. This work proved that, in the limit for large Gaussian mixing matrices, the behavior of AMP estimates can be described by a scalar equivalent model with effective noise levels defined by certain scalar state evolution (SE) equations. Similar SE analyses have appeared in related contexts . To prove the SE equations for the IterFac algorithm, we apply a key theorem from with a simple change of variables and a slight modification to account for parameter adaptation.

A conference version of this paper appeared in . This paper provides all the proofs along with more detailed discussions and simulations. Since the original publication of the conference version in , several other research groups have extended and built on the work. Importantly, has shown that in the case of certain discrete priors, the IterFac algorithm is provably Bayes optimal in a large system limit. It was shown in , that an algorithm closely related to IterFac could provide the best known scaling laws for the hidden clique problem. More recently, and have used related AMP-type methods for matrix factorization problems with rank greater than one.

Iterative Rank-One Factorization

for some scalar functions $G_{u}(\cdot)$ and $G_{v}(\cdot)$. The choice of the factor selection function will depend on whether IterFac is used for MAP or MMSE estimation.

To describe the factor selection functions for the MAP estimation problem (4), Let $\lambda_{u}(t)$ and $\lambda_{v}(t)$ be sequences of pairs of parameters

Then, for each $t$, consider random vectors $\mathbf{p}(t)$ and $\mathbf{q}(t)$ given by,

The random variables $\mathbf{p}(t)$ and $\mathbf{q}(t)$ in (8) are simply scaled versions of the true vectors $\mathbf{u}_{0}$ and $\mathbf{v}_{0}$ with additive white Gaussian noise (AWGN). Then, for the MAP problem we take the factor selection functions to be the MAP estimates of $\mathbf{u}_{0}$ and $\mathbf{v}_{0}$ given the observations $\mathbf{p}(t)$ and $\mathbf{q}(t)$:

Importantly, due to the separability assumption on the priors (2), these MAP estimates can be computed componentwise:

Hence, the IterFac algorithm replaces the vector-valued MAP estimation of $\mathbf{u}_{0},\mathbf{v}_{0}$ from the joint density (3), with a sequence of scalar MAP estimation problems along with multiplications by $\mathbf{A}$ and $\mathbf{A}^{T}$.

MMSE Estimation:

For the MMSE estimation problem, we simply take the factor selection functions the MMSE estimates with respect to the random variables (8),

Again, due to the separability assumption (2), these MMSE estimation problems can be performed componentwise. Hence MMSE IterFac reduces the vector MMSE problem to a sequence of scalar problems in Gaussian noise.

Intuitive Algorithm Derivation and Analysis

Before we formally analyze the algorithm, it is instructive to understand the intuition behind the method. For both the MAP and MMSE versions of IterFac, first observe that

where (a) follows from line 6 in Algorithm 1; (b) follows from the assumptions on the measurements (1) and in (c), we have used the definitions

For the first iteration, $\mathbf{W}$ is a large zero mean matrix independent of the initial condition $\mathbf{v}(0)$. Also, $\mu_{u}(0)=0$ and hence the components of $\mathbf{z}_{u}(0)$ will be asymptotically Gaussian zero mean variables due to the Central Limit Theorem. Hence, the vector $\mathbf{p}(0)$ will be distributed as the true vector $\mathbf{u}_{0}$ with Gaussian noise as in the model (8). Therefore, we can take an initial MAP or MMSE estimate of $\mathbf{u}_{0}$ from the scalar estimation functions in (9) or (11).

Unfortunately, in subsequent iterations with $t>0$, $\mathbf{v}(t)$ is no longer independent of $\mathbf{W}$ and hence $\mathbf{W}\mathbf{v}(t)$ will not, in general, be a Gaussian zero mean random vector. However, remarkably, we will show that with the specific choice of $\mu_{u}(t)$ in Algorithm 1, the addition of the term $\mu_{u}(t)\mathbf{u}(t)$ in (13) “debiases” the noise term so that $\mathbf{z}_{u}(t)$ is asymptotically zero mean Gaussian. Thus, $\mathbf{p}(t)$ continues to be distributed as in (8) and we can construct MAP or MMSE estimates of $\mathbf{u}_{0}$ from the vector $\mathbf{p}(t)$. For example, using the MMSE estimation function (11) with appropriate selection of the parameters $\lambda_{u}(t)$, we obtain that the estimate for $\mathbf{u}(t\!+\!1)$ is given by the MMSE estimate

For the MAP estimation problem, the MAP selection function (9) computes the MAP estimate for $\mathbf{u}_{0}$.

Similarly, for $\mathbf{q}(t)$, we see that

for $\gamma_{v}(t)=\mathbf{u}_{0}^{T}\mathbf{u}(t\!+\!1)/m$ and $\mathbf{z}_{v}(t)=(1/\sqrt{m})\mathbf{W}^{T}\mathbf{u}(t\!+\!1)+\mu_{v}(t)\mathbf{v}(t)$. Then, $\mathbf{v}(t\!+\!1)$ is taken as the MMSE or MAP estimate of $\mathbf{v}_{0}$ from the vector $\mathbf{q}(t)$.

Thus, we see that the algorithm attempts to produce a sequence of unbiased estimates $\mathbf{p}(t)$ and $\mathbf{q}(t)$ for scaled versions of $\mathbf{u}_{0}$ and $\mathbf{v}_{0}$. Then, it constructs the MAP or MMSE estimates $\mathbf{u}(t)$ and $\mathbf{v}(t\!+\!1)$ from these unbiased estimates.

2 State Evolution Analysis

The above “derivation” of the algorithm suggests a possible method to analyze the IterFac algorithm. Consider a sequence of problems indexed by $n$ and suppose that the dimension $m=m(n)$ grows linearly with $n$ in that

for some $\beta>0$. For each $n$ and iteration number $t$, define the sets

which are the set of the components of the true vectors $u_{0i}$ and $v_{0j}$ and their corresponding estimates $u_{i}(t)$ and $v_{j}(t)$. To characterize the quality of the estimates, we would like to describe the distributions of the pairs in $\theta_{u}(t)$ and $\theta_{v}(t)$.

Our formal analysis below will provide an exact characterization of these distributions. Specifically, we will show that the sets have empirical limits of the form

for some random variable pairs $(U_{0},U(t)$ and $(V_{0},V(t))$. The precise definition of empirical limits is given in Appendix A.1, but loosely, it means that the empirical distribution of the components in $\theta_{u}(t)$ and $\theta_{v}(t)$ converge to certain random variable pairs.

In addition, we can inductively derive what the distributions are for the limits in (17). To make matters simple, suppose that the parameters $\lambda_{u}(t)$ and $\lambda_{v}(t)$ are not selected adaptively but instead given by a fixed sequences $\overline{\lambda}_{u}(t)$ and $\overline{\lambda}_{v}(t)$. Also, given random variables in the limits of (17), define the deterministic constants

Now suppose that the second limit in (17) holds for some $t$. Then, $\gamma_{u}(t)$ in (13) would have the following limit:

Also, if we ignore the debiasing term with $\mu_{u}(t)$ and assume that $\mathbf{W}$ is independent of $\mathbf{v}(t)$, the variance of the components of $\mathbf{z}_{u}(t)$ in (13) would be

Thus, we would expect a typical component of $p_{i}(t)$ in (12) to be distributed as

Due to the separability assumption (6), each component $u_{i}(t\!+\!1)=G_{u}(p_{i}(t),\overline{\lambda}_{u}(t))$. So, we would expect the components to follow the distribution

This provides an exact description of the expected joint density of $(U_{0},U(t\!+\!1))$. From this density we can then compute $\overline{\alpha}_{u0}(t\!+\!1)$ and $\overline{\alpha}_{u1}(t\!+\!1)$ in (18).

A similar calculation, again assuming the “debiasing” and Gaussianity assumptions are valid, shows that the limiting empirical distribution of $\theta_{v}(t\!+\!1)$ in (17) should follow

This provides the joint density $(V_{0},V(t\!+\!1))$ from which we can compute $\overline{\alpha}_{v0}(t\!+\!1)$ and $\overline{\alpha}_{v1}(t\!+\!1)$ in (18). Thus, we have provided a simple method to recursively compute the joint densities of the limits in (17) and their second-order statistics (18).

Asymptotic Analysis under Gaussian Noise

In the above intuitive analysis, we did not formally describe the sense of convergence nor offer any formal justification for the Gaussianity assumptions. In addition, we assumed that the parameters $\lambda_{u}(t)$ and $\lambda_{v}(t)$ were fixed sequences. In reality, we will need them to be data adaptive. We will make the above arguments rigorous under the following assumptions. Note that although we will be interested in MAP or MMSE estimation functions, our analysis will apply to arbitrary factor selection functions $G_{u}(\cdot)$ and $G_{v}(\cdot)$.

Consider a sequence of random realizations of the estimation problem in Section 1 indexed by the dimension $n$. The matrix $\mathbf{A}$ and the parameters in Algorithm 1 satisfy the following:

For each $n$, the output dimension $m=m(n)$ is deterministic and scales linearly with $n$ as (15).

The factor selection functions $G_{u}(\mathbf{p},\lambda_{u})$ and $G_{v}(\mathbf{q},\lambda_{v})$ in lines 7 and 12 are componentwise separable in that for all component indices $i$ and $j$,

for some scalar functions $G_{u}(p,\lambda_{u})$ and $G_{v}(q,\lambda_{v})$. The scalar functions must be differentiable in $p$ and $q$. Moreover, for every $t$, the functions $G_{u}(p,\lambda_{u})$ and $\partial G_{u}(p,\lambda_{u})/\partial p$ must be Lipschitz continuous in $p$ with a Lipschitz constant that is continuous in $\lambda_{u}$, and continuous in $\lambda_{u}$ uniformly over $p$. Similarly, for every $t$, the functions $G_{v}(q,\lambda_{v})$ and $\partial G_{v}(q,\lambda_{v})/\partial q$ must be Lipschitz continuous in $q$ with a Lipschitz constant that is continuous in $\lambda_{v}$, and continuous in $\lambda_{v}$ uniformly over $q$.

The parameters $\lambda_{u}(t)$ and $\lambda_{v}(t)$ are computed via

for (possibly vector-valued) functions $\phi_{\lambda u}(\cdot)$ and $\phi_{\lambda v}(\cdot)$ that are pseudo-Lipschitz continuous of order $p=2$.

For $t=0$, the sets (17) empirically converge with bounded moments of order $2$ to the limits

for some random variable pairs $(U_{0},U(0)$ and $(V_{0},V(0))$. See Appendix A.1 for a precise definition of the empirical convergence used here.

The assumptions need some explanations. Assumptions 4.1(a) and (b) simply state that we are considering an asymptotic analysis for certain large matrices $\mathbf{A}$ consisting of a random rank one matrix plus Gaussian noise. The analysis of Algorithm 1 for higher ranks is still not known, but we provide some possible ideas later. Assumption 4.1(c) is a mild condition on the factor selection functions. In particular, the separability assumption holds for the MAP or MMSE functions (9) and (11) under separable priors.

Assumption (d) allows for the parameters $\lambda_{u}(t)$ and $\lambda_{v}(t)$ in the factor selection functions to be data dependent, provided that they can each be determined via empirical averages of some function of the most recent data. Assumption (e) is the initial induction hypothesis.

Under these assumptions, we can recursively define the sequence of random variables $(U_{0},U(t))$ and $(V_{0},V(t))$ as described above. For the parameters $\overline{\lambda}_{u}(t)$ and $\overline{\lambda}_{v}(t)$, define them from the expectations

These are the limiting values we would expect given the adaptation rules (22).

Under Assumption 4.1, the sets $\theta_{u}(t)$ and $\theta_{v}(t)$ in (16) converge empirically with bounded moments of order $p=2$ with the limits in (17).

2 Scalar Equivalent Model

The main contribution of Theorem 4.1 is that it provides a simple scalar equivalent model for the asymptotic behavior of the algorithm. The sets $\theta_{u}(t)=\{(u_{0i},u_{i}(t))\}$ and $\theta_{v}(t)=\{(v_{0j},v_{j}(t))\}$ in (16) are the components of true vectors $\mathbf{u}_{0}$ and $\mathbf{v}_{0}$ and their estimates $\mathbf{u}(t)$ and $\mathbf{v}(t)$. The theorem shows that empirical distribution of these components are asymptotically equivalent to simple random variable pairs $(U_{0},U(t))$ and $(V_{0},V(t))$ given by (19) and (20). In this scalar system, the variable $U(t\!+\!1)$ is the output of the factor selection function $G_{u}(\cdot)$ applied to a scaled and Gaussian noise-corrupted version of the true variable $U_{0}$. Similarly, $V(t\!+\!1)$ is the output of the factor selection function $G_{v}(\cdot)$ applied to a scaled and Gaussian noise-corrupted version of the true variable $V_{0}$. Following , we can thus call the result a single-letter characterization of the algorithm.

From this single-letter characterization, one can exactly compute a large class of performance metrics of the algorithm. Specifically, the empirical convergence of $\theta_{u}(t)$ shows that for any pseudo-Lipschitz function $\phi(u_{0},u)$ of order $p$, the following limit exists almost surely:

where the expectation on the right-hand side is over the variables $(U_{0},U(t))$ with $U_{0}$ identical to the variable in the limit in (23) and $U(t)$ given by (19). This expectation can thus be explicitly evaluated by a simple two-dimensional integral and consequently any component-separable performance metric based on a suitably continuous loss function $\phi(u_{0},u)$ can be exactly computed.

For example, if we take $\phi(u_{0},u)=(u-u_{0})^{2}$ we can compute the asymptotic mean squared error of the estimate,

Also, for each $t$, define the empirical second-order statistics

From the second order statistics, we can compute the asymptotic correlation coefficient between $\mathbf{u}_{0}$ and its estimate $\mathbf{u}$ given by

Similarly, the asymptotic correlation coefficient between $\mathbf{v}_{0}$ and $\mathbf{v}$ has a simple expression

The correlation coefficient is useful, since we know that, without additional constraints, the terms $\mathbf{u}_{0}$ and $\mathbf{v}_{0}$ can only be estimated up to a scalar. The correlation coefficient is scale invariant.

Examples

The SE analysis can be used to exactly predict the asymptotic behavior of the IterFac algorithm for any smooth scalar estimation functions, including the MAP or MMSE functions (9) and (11). There are, however, two cases, where the SE equations have particularly simple and interesting solutions: linear estimation functions and MMSE functions.

Suppose we use linear selection functions of the form

where the parameters $\lambda_{u}$ and $\lambda_{v}$ allow for normalization or other scalings of the outputs. Linear selection functions of the form (31) arise when one selects $G_{u}(\cdot)$ and $G_{v}(\cdot)$ from the MAP or MMSE functions (9) or (11) with Gaussian priors.

With Gaussian priors, the correct solution to the MAP estimate (4) is for $(\widehat{\mathbf{u}},\widehat{\mathbf{v}})$ to be the (appropriately scaled) left and right maximal singular vectors of $\mathbf{A}$. We will thus call the estimates $(\mathbf{u}(t),\mathbf{v}(t))$ of Algorithm 1 and linear selection functions (31) the estimated maximal singular vectors.

Consider the state evolution equation (18), (19), and (20) with the linear selection functions (31). Then:

The asymptotic correlation coefficients (29) and (30) satisfy the following recursive rules:

For any positive initial condition, $\rho_{v}(0)>0$, the asymptotic correlation coefficients converge to the limits

The theorem provides a set of recursive equations for the asymptotic correlation coefficients $\rho_{u}(t)$ and $\rho_{v}(t)$ along with simple expressions for the limiting values as $t\rightarrow\infty$. We thus obtain exactly how correlated the estimated maximal singular vectors of a matrix $\mathbf{A}$ of the form (1) are to the rank one factors $(\mathbf{u}_{0},\mathbf{v}_{0})$. The proof of the theorem also provides expressions for the second-order statistics in (18) to be used in the scalar equivalent model.

The fixed point expressions (33) agree with the more general results in that derive the correlations for ranks greater than one and low-rank recovery with missing entries. Similar results can also be found in . An interesting consequence of the expressions in (33) is that unless

the asymptotic correlation coefficients are exactly zero. The ratio $\tau_{u}\tau_{v}/\tau_{w}$ can be interpreted as a scaled SNR.

2 Minimum Mean-Squared Error Estimation

Next suppose we use the MMSE selection functions (11). Using the scalar equivalent models (19) and (20), we take the scaling factor and noise parameters as

where $D\sim\mathcal{N}(0,1)$ is independent of $U_{0}$ and $V_{0}$. That is, ${\mathcal{E}}_{u}(\eta_{u})$ and ${\mathcal{E}}_{v}(\eta_{v})$ are the mean-squared errors of estimating $U_{0}$ and $V_{0}$ from observations $Y$ with SNRs of $\eta_{u}$ and $\eta_{v}$. The functions ${\mathcal{E}}_{u}(\cdot)$ and ${\mathcal{E}}_{v}(\cdot)$ arise in a range of estimation problems and the analytic and functional properties of these functions can be found in .

For all $t$, the asymptotic correlation coefficients (29) and (30) satisfy the recursive relationships

If, in addition, ${\mathcal{E}}_{u}(\eta_{u})$ and ${\mathcal{E}}_{v}(\eta_{v})$ are continuous, then for any positive initial condition, $\rho_{v}(0)>0$, as $t\rightarrow\infty$, the asymptotic correlation coefficients $(\rho_{u}(t),\rho_{v}(t))$ increase monotonically to fixed points $(\rho_{u}^{*},\rho_{u}^{*})$ of (37) with $\rho_{v}^{*}>0$.

Again, we see that we can obtain simple, explicit recursions for the asymptotic correlations. Moreover, the asymptotic correlations provably converge to fixed points of the SE equations. The proof of the theorem also provides expressions for the second-order statistics in (18) to be used in the scalar equivalent model.

3 Zero Initial Conditions

Of course, with zero mean random variables, a more sensible initial condition is to take $\mathbf{v}(0)$ to be some non-zero random vector, as is commonly done in power algorithm recursions for computing maximal singular vectors. To understand the behavior of the algorithm under this random initial condition, let

where we have explicitly denoted the dependence on the problem dimension $n$. From (30), we have that $\lim_{n\rightarrow\infty}\rho_{v}(t,n)=\rho_{v}(t)$ for all $t$. Also, with a random initial condition $\mathbf{v}(0)$ independent of $\mathbf{v}_{0}$, it can be checked that $\rho_{v}(0,n)=O(1/n)$ so that

Hence, from the SE equations (37), $\rho_{v}(t)=\rho_{u}(t)=0$ for all $t$. That is,

This limit suggests that, even with random initial condition, the IterFac algorithm will not produce a useful estimate.

However, it is still possible that the limit in the opposite order of (39) may be non-zero:

That is, for each $n$, it may be possible to obtain a non-zero correlation, but the number of iterations for convergence increases with $n$ since the algorithm starts from a decreasingly small initial correlation. Unfortunately, our SE analysis cannot make predictions on limits in the order of (40).

We can however analyze the following limit:

If $\sqrt{\beta}\tau_{u}\tau_{v}>\tau_{w}$,

Conversely, if $\sqrt{\beta}\tau_{u}\tau_{v}<\tau_{w}$,

The result of the lemma is somewhat disappointing. The lemma shows that $\sqrt{\beta}\tau_{u}\tau_{v}>\tau_{w}$ is essentially necessary and sufficient for the IterFac algorithm with MMSE estimates to be able to overcome arbitrarily small initial conditions and obtain an estimate with a non-zero correlation to the true vector. Unfortunately, this is the identical to the condition (34) for the linear estimator to obtain a non-zero correlation. Thus, the IterFac algorithm with MMSE estimates performs no better than simple linear estimation in the initial iterations when the priors have zero means. Since linear estimation is equivalent to finding maximal singular vectors without any particular constraints, we could interpret Lemma 5.5 as saying that the IterFac algorithm under MMSE estimation cannot exploit structure in the components in the initial iterations. As a result, in low SNRs it may be necessary to use other algorithms as an initial condition for IterFac – such procedures, however, require further study.

Numerical Simulation

which provides a simple model for a sparse, positive vector. The parameter $\lambda$ is the fraction of nonzero components and is set in this simulation to $\lambda=0.1$. Note that these components have a non-zero mean so the difficulties of Section 5.3 are avoided. The dimensions are $(m,n)=(1000,500)$, and the noise level $\tau_{w}$ is set according to the scaled SNR defined as

The results of the simulation are shown in Fig. 1, which shows the simulated and SE-predicted performance of the IterFac algorithm with both the linear and MMSE selection functions for the priors on $\mathbf{u}$ and $\mathbf{v}$. The algorithm is run for $t=10$ iterations and the plot shows the median of the final correlation coefficient $\rho_{v}(t)$ over 50 Monte Carlo trials at each value of SNR. It can be seen that the performance of the IterFac algorithm for both the linear and MMSE estimates are in excellent agreement with the SE predictions. The correlation coefficient of the linear estimator also matches the correlation of the estimates produced from the maximal singular vectors of $\mathbf{A}$. This is not surprising since, with linear selection functions, the IterFac algorithm is essentially an iterative method to find the maximal singular vectors. The figure also shows the benefit of exploiting the prior on $\mathbf{v}_{0}$, which is evident from the superior performance of the MMSE estimate over the linear reconstruction.

Fig. 2 shows the correlation coefficient as a function of the iteration number for the MMSE and linear methods for two values of the SNR. Again, we see that the SE predictions are closely matched to the median simulated values. In addition, we see that we get good convergence within 4 to 8 iterations. Based on the SE predictions, this number will not scale with dimension and hence the simulation suggests that only a small number of iterations will be necessary for even very large problems. All code for the simulation can be found in the public GAMP sourceforge project .

Limitations and Extensions

There are several potential lines for future work – some of which have already been explored in other works made since the original publication of this paper in .

The algorithm presented in this paper considers only rank one matrices. The works have proposed AMP-type algorithms for more general classes of matrix factorization problems and provided analyses of these methods based on replica methods and other ideas from statistical physics.

Unknown priors

The MMSE estimator in Section 5.2 requires exact knowledge of the priors on $U_{0}$ and $V_{0}$ as well as the Gaussian noise level $\tau_{w}$. In many problems in statistics, these are not known. There are two possible solutions that may be investigated in the future. One method is to parameterize the distributions of $U_{0}$ and $V_{0}$ and estimate these parameters in the MMSE selection functions (11) – potentially through an EM type procedure as in . This EM type approach with hidden hyperparameters has been recently successfully used in a related approximate message passing method in . The analysis of the such parameter learning could possibly be accounted for through the adaptation parameters $\lambda_{u}(t)$ and $\lambda_{v}(t)$. A second approach is to assume that the distributions of $U_{0}$ and $V_{0}$ belong to a known family of distributions and then find a min-max solution. Such a min-max technique was proposed for AMP recovery of sparse vectors in . See also .

Optimality

While the current paper characterizes the performance of the IterFac algorithm, it remains open how far that performance is to optimal estimation such as the joint MMSE estimates of $\mathbf{u}_{0}$ and $\mathbf{v}_{0}$. AMP methods have been shown to be provably optimal in a wide range of scenarios . More recently, has proven that for sparse binary priors, the IterFac algorithm was provably Bayes optimal for certain sparsity levels. A potential line of future work is to see if these results can be extended to more general settings.

Conclusions

We have presented a computationally-efficient method for estimating rank-one matrices in noise. The estimation problem is reduced to a sequence of scalar AWGN estimation problems which can be performed easily for a large class of priors or regularization functions on the coefficients. In the case of Gaussian noise, the asymptotic performance of the algorithm is exactly characterized by a set of scalar state evolution equations which appear to match the performance at moderate dimensions well. Thus, the methodology is computationally simple, general and admits a precise analysis in certain asymptotic, random settings. Future work include extensions to higher rank matrices and handling of the cases where the priors are not known.

Acknowledgments

The authors thank Vivek Goyal, Ulugbek Kamilov, Andrea Montanari and Phil Schniter for detailed comments on an earlier draft.

Appendix A Appendices: Proofs

Now suppose that for each $n=1,2,\ldots$, we have a set of vectors

When the nature of convergence is clear, we may write (with some abuse of notation)

A.2 Bayati–Montanari Recursions with Adaptation

and $\xi_{v}(t)$ and $\xi_{u}(t)$ are scalar step sizes given by

The recursions (46) to (48) are identical to the recursions analyzed in , except for the introduction of the parameters $\nu_{u}(t)$ and $\nu_{v}(t)$. We will call these parameters adaptation parameters since they enable the functions $H_{u}(\cdot)$ and $H_{v}(\cdot)$ to have some data dependence, not explicitly considered in . Similar to the selection of the parameters $\lambda_{u}(t)$ and $\lambda_{v}(t)$ in (22), we assume that, in each iteration $t$, the adaptation parameters are selected by functions of the form,

where $\phi_{u}(\cdot)$ and $\phi_{v}(\cdot)$ are (possibly vector-valued) pseudo-Lipschitz continuous of order $p$ for some $p>1$. Thus, the values of $\nu_{u}(t)$ and $\nu_{v}(t)$ depend on the outputs $\mathbf{v}(t)$ and $\mathbf{u}(t\!+\!1)$. Note that in equations (47) to (49), $d_{i}$, $u_{0i}$, $b_{j}$ and $v_{0j}$ are the components of the vectors $\mathbf{d}$, $\mathbf{u}_{0}$, $\mathbf{b}$ and $\mathbf{v}_{0}$, respectively. The algorithm is initialized with $t=0$, $\xi_{u}(0)=0$ and some vector $\mathbf{v}(0)$.

Now, similar to Section 4, consider a sequence of random realizations of the parameters indexed by the input dimension $n$. For each $n$, we assume that the output dimension $m=m(n)$ is deterministic and scales linearly as in (15) for some $\beta\geq 0$. Assume that the transform matrix $\mathbf{S}$ has i.i.d. Gaussian components $s_{ij}\sim{\cal N}(0,1/m)$. Also assume that the empirical limits in (23) hold with bounded moments of order $2p-2$ for some limiting random variables $(U_{0},U(0))$ and $(V_{0},V(0))$. We will also assume the following continuity assumptions on $H_{u}(\cdot)$ and $H_{v}(\cdot)$:

The function $H_{u}(t,b,u_{0},\nu_{u})$ satisfies the following continuity conditions:

For every $\nu_{u}$ and $t$, $H_{u}(t,b,u_{0},\nu_{u})$ and its derivative $\partial H_{u}(t,b,u_{0},\nu_{u})/\partial b$ are Lipschitz continuous in $b$ and $u_{0}$ for some Lipschitz constant that is continuous in $\nu_{u}$; and

For every $\nu_{u}$ and $t$, $H_{u}(t,b,u_{0},\nu_{u})$ and $\partial H_{u}(t,b,u_{0},\nu_{u})/\partial b$ are is continuous at $\nu_{u}$ uniformly over $(b,u_{0})$.

The function $H_{v}(t,d,v_{0},\nu_{v})$ satisfies the analogous continuity assumptions in $d$, $v_{0}$ and $\nu_{v}$.

Under these assumption, we will show, as in Section 4, that for any fixed iteration $t$, the sets $\theta_{u}(t)$ and $\theta_{v}(t)$ in (16) converge empirically to the limits (17) for some random variable pairs $(U_{0},U(t))$ and $(V_{0},V(t))$. The random variable $U_{0}$ is identical to the variable in the limit (23) and, for $t\geq 0$, $U(t)$ is given by

for some deterministic constants $\overline{\nu}_{v}(t)$ and $\tau^{b}(t)$ that will be defined below. Similarly, the random variable $V_{0}$ is identical to the variable in the limit (23) and, for $t\geq 0$, $V(t)$ is given by

for some constants $\overline{\nu}_{v}(t)$ and $\tau^{d}(t)$, also defined below.

The constants $\tau^{b}(t)$, $\tau^{d}(t)$, $\overline{\nu}_{u}(t)$ and $\overline{\nu}_{v}(t)$ can be computed recursively through the following state evolution equations

where the expectations are over the random variables $U(t)$ and $V(t)$ above and initialized with

With these definitions, we can now state the adaptive version of the result from Bayati and Montanari . Although the full proof requires that $p=2$, much of the proof is valid for $p\geq 2$. Hence, where possible, we provide the steps for the general $p$ case.

Consider the recursion in (46) to (49) satisfying the above assumptions for the case when $p=2$. Then, for any fixed iteration number $t$, the sets $\theta_{u}(t)$ and $\theta_{v}(t)$ in (16) converge empirically to the limits (17) with bounded moments of order $p=2$ to the random variable pairs $(U_{0},U(t))$ and $(V_{0},V(t))$ described above.

Now, Bayati and Montanari’s result in shows that this non-adaptive algorithm satisfies the required properties. That is, the following limits hold with bounded moments of order $p$,

So, the limits (17) will be shown if we can prove the following limits hold almost surely for all $t$:

where $\|\cdot\|_{p}$ is the $p$-norm. In the course of proving (57), we will also show the following limits hold almost surely,

The proof of the limits (57) and (58) can be demonstrated via induction on $t$ with the following straightforward (but somewhat tedious) continuity argument:

To begin the induction argument, first note that the non-adaptive algorithm has the same initial condition as the adaptive algorithm. That is, $\mathbf{v}^{*}(0)=\mathbf{v}(0)$ and $\xi_{u}(0)=\xi_{u}^{*}(0)=0$. Also, since $\xi_{u}(0)=\xi_{u}^{*}(0)=0$, from (46a), the initial value of $\mathbf{u}(t)$ does not matter. So, without loss of generality, we can assume that the initial condition satisfies $\mathbf{u}(0)=\mathbf{u}^{*}(0)$. Thus, the limits (57) and (58c) hold for $t=0$.

We now proceed by induction. Suppose that the limits (57) and (58c) hold almost surely for some $t\geq 0$. Since $\mathbf{S}$ has i.i.d. components with zero mean and variance $1/m$, by the Marceko-Pastur Theorem , the maximum singular value of $\mathbf{S}$ is bounded. For $p=2$, the maximum singular value is the $p$-norm operator norm, and therefore, there exists a constant $C_{S}>0$ such that

Substituting the bound (59) into (46a), we obtain

Now, since $p\geq 1$, we have that for any positive numbers $a$ and $b$,

Applying (61) into (60), and the fact that $\lim_{n}m/n=\beta$, we obtain that

for some other constant $C^{\prime}_{S}>0$. Now, since $\mathbf{u}^{*}(t)$ is the output of the non-adaptive algorithm it satisfies the limit

Substituting the bound (63) along with induction hypotheses, (57) and (58c) into (62) shows (58a).

Next, to prove the limit (58e), first observe that since $\phi_{u}(\cdot)$ is pseudo-Lipschitz continuous of order $p$, we have that $\overline{\nu}_{u}(t)$ in (52b) can be replaced by the limit of the empirical means

where the limit holds almost surely. Combining (64) with the (49),

Applying the fact that $\phi_{u}(\cdot)$ is pseudo-Lipschitz continuous of order $p$ to (65), we obtain that there exists a constant $L_{v}>0$ such that

where the last step is due to Hölder’s inequality with the exponents

Now, similar to the proof of (63) one can show that the non-adaptive output satisfies the limit

Also, from the induction hypothesis (57), it follows that the non-adaptive output must satisfy the same limit

Applying the bounds (67) and (68) and the limit (66) shows that (58e) holds almost surely.

Now the limit in (58e) and the second limit in (57) together with the continuity conditions on $H_{u}(\cdot)$ in Assumption A.1 show that the first limit in (57) holds almost surely for $t+1$ and (58d) holds almost surely for $t$. Using (46b), the proof of the limit (58b) is similar to the proof of (58a). These limits in turn show that the second limit in (57) and the limits in (58c) hold almost surely for $t+1$. We have thus shown that if (57) and (58c) hold almost surely for some $t$, they hold for $t+1$. Thus, by induction they hold for all $t$. Finally, applying the limits (56), (57) and a continuity argument shows that the desired limits (17) hold almost surely.

A.3 Proof of Theorem 4.1

The theorem directly follows from the adaptive Bayati–Montanari recursion theorem, Theorem A.2 above, with some change of variables. Specifically, let

where $\mathbf{W}$ is the Gaussian noise in the rank one model in Assumption 4.1(b). Since $\mathbf{W}$ has i.i.d. components with Gaussian distributions ${\mathcal{N}}(0,\tau_{w})$, the components of $\mathbf{S}$ will be i.i.d. with distributions ${\mathcal{N}}(0,1/m)$.

Now, using the rank one model for $\mathbf{A}$ in (1)

where the last step is from the definition of $\alpha_{v1}(t)$ in (26). Substituting (70) into the the update rule for $\mathbf{p}(t)$ in line 6 of Algorithm 1, we obtain

Note that we have used the fact that $\beta=n/m$. Hence, if we define

Similarly, one can show that if we define

which are the adaptation functions in (22) with additional components for the second-order statistics $uu_{0}$ and $vv_{0}$. Since $\phi_{\lambda u}(t,\cdot)$ and $\phi_{\lambda v}(t,\cdot)$ are pseudo-Lipschitz of order $p$, so are $\phi_{u}(t,\cdot)$ and $\phi_{v}(t,\cdot)$. Taking the empirical means over each of the two components of $\phi_{u}(\cdot)$ and $\phi_{v}(\cdot)$, and applying (22) and (26), we see that if $\nu_{u}(t)$ and $\nu_{v}(t)$ are defined as in (49),

Therefore, $\nu_{u}(t)$ and $\nu_{v}(t)$ are vectors containing the parameters $\lambda_{u}(t)$ and $\lambda_{v}(t)$ for the factor selection functions in lines 7 and 12 of Algorithm 1 as well as the second-order statistics $\alpha_{u1}(t)$ and $\alpha_{v1}(t)$. Now, for $\nu_{u}=(\alpha_{v1},\lambda_{u})$ and $\nu_{v}=(\alpha_{u1},\lambda_{v})$ define the scalar functions

Since $G_{u}(p,\lambda_{u})$ and $G_{v}(q,\lambda_{v})$ satisfy the continuity conditions in Assumption 4.1(c), $H_{u}(t,b,u_{0},\lambda_{u})$ and $H_{v}(t,d,v_{0},\lambda_{v})$ satisfy Assumption A.1. In addition, the componentwise separability assumption in (21) implies that the updates in lines 7 and 12 of Algorithm 1 can be rewritten as

Thus, combining (78) and (79) with the definitions of $\mathbf{b}(t)$ and $\mathbf{d}(t)$ in (72) and (74), we obtain

where (a) follows from the definition of $\xi_{u}(t)$ in (73); (b) is the setting for $\mu_{u}(t\!+\!1)$ in line 8; and (c) follows from the definition of $H_{v}(t,d)$ in (78b). Similarly, one can show that

Equations (73), (75), (77), (78), (81) and (82) exactly match the recursions in equations (46) to (49). Therefore, Theorem A.2 shows that the limits (17) hold in a sense that the sets $\theta_{u}(t)$ and $\theta_{v}(t)$ converge empirically with bounded moments of order $p$.

We next show that the limits $U(t)$ and $V(t)$ on the right-hand side of (17) match the descriptions in (19) and (20). First, define $\overline{\alpha}_{u1}(t)$ and $\overline{\alpha}_{v1}(t)$ in (18) and $\overline{\lambda}_{u}(t)$ and $\overline{\lambda}_{v}(t)$ as in (24). Then, from (76), the expectations $\overline{\nu}_{u}(t)$ and $\overline{\nu}_{v}(t)$ in (52b) are given by

where $B(t)\sim{\mathcal{N}}(0,\tau^{b}(t))$. Therefore, if we let $Z_{u}(t)=\sqrt{\tau_{w}}B(t)$, then $Z_{u}(t)$ is zero mean Gaussian with variance

where follows from (52a) and (b) follows from the definition of $\overline{\alpha}_{v0}(t)$ in (18). Substituting $Z_{u}(t)=\sqrt{\tau_{w}}B(t)$ into (84) we obtain the model for $U(t\!+\!1)$ in (19). Similarly, using (51) and (78b), we can obtain the model $V(t\!+\!1)$ in (20). Thus, we have proven that the random variables $(U_{0},U(t))$ and $(V_{0},V(t))$ are described by (19) and (20), and this completes the proof.

A.4 Proof of Theorem 5.1

The theorem is proven by simply evaluating the second order statistics. We begin with $\overline{\alpha}_{u1}(t\!+\!1)$:

where (a) is the definition in (18); (b) follows from (19) and (31); (c) follows from (19); and (d) follows from the independence of $Z_{u}(t)$ and $U_{0}$ and the definition of $\tau_{u}$ in (28). Similarly, one can show that

Substituting (85) and (86) into (29), we obtain the asymptotic correlation coefficient

where the last step follows from (30). This proves the first equation in (32).

Applying these equations into (29) and (30), we obtain the recursion (32). Hence, we have proven part (a) of the theorem.

For part (b), we need the following simple lemma.

Suppose that $H:\rightarrow$ is continuous and monotonically increasing, and $x(t)$ is a sequence satisfying the recursive relation

for some initial condition $x(0)\in$. Then, either $x(t)$ monotonically increases or decreases to some $x^{*}=H(x^{*})$.

This can be proven similar to [42, lemma 7].

To apply Lemma A.4, observe that the recursions (32) show that

then it follows from (32) that $\rho_{u}(t\!+\!1)=H(\rho_{u}(t))$ for all $t$. By taking the derivative, it can be checked that $H(\rho_{u})$ is monotonically increasing. It follows from Lemma A.4 that $\rho_{u}(t)\rightarrow\rho_{u}^{*}$ for some fixed point $\rho_{u}^{*}=H(\rho_{u}^{*})$ with $\rho_{u}^{*}\in$.

Now, there are only two fixed point solutions to $\rho_{u}^{*}=H(\rho_{u}^{*})$: $\rho_{u}^{*}=0$ and

then $\rho_{u}^{*}$ in (90) is not positive, so the only fixed solution in $$is$\rho_{u}^{*}=0$. Therefore, when (91) is not satisfied,$\rho_{u}(t)$must converge to the zero fixed point:$\rho_{u}(t)\rightarrow 0$.

Now, suppose that (91) is satisfied. In this case, we claim that $\rho_{u}(t)\rightarrow\rho_{u}^{*}$ where $\rho_{u}^{*}$ is in (90). We prove this claim by contradiction and suppose, instead, that $\rho_{u}(t)$ converges to the other fixed point: $\rho_{u}(t)\rightarrow 0$. Since Lemma A.4 shows that $\rho_{u}(t)$ must be either monotonically increasing or decreasing, the only way $\rho_{u}(t)\rightarrow 0$ is that $\rho_{u}(t)$ monotonically decreases to zero. But, when (91) is satisfied, it can be checked that for $\rho_{u}(t)$ sufficiently small and positive, $\rho_{u}(t\!+\!1)>H(\rho_{u}(t))$. This contradicts the fact that $\rho_{u}(t)$ is monotonically decreasing, and therefore, $\rho_{u}(t)$ must converge to the other fixed point $\rho_{u}^{*}$ in (90).

Hence, we have shown that when (91) is not satisfied, $\rho_{u}(t)\rightarrow 0$, and when (91) is satisfied $\rho_{u}(t)\rightarrow\rho_{u}^{*}$ in (90). This is equivalent to the first limit in (33). The second limit (33) is proved similarly.

A.5 Proof of Theorem 5.3

and therefore $\overline{\alpha}_{u0}(t)=\overline{\alpha}_{u1}(t)$. Now, consider the measurement $P(t)$ in (19). The SNR in this channel is

Since $U(t\!+\!1)$ is the conditional expectation of $U_{0}$ given $P(t)$, the mean-squared error is given by ${\mathcal{E}}_{u}(\eta_{u}(t))$ defined in (36). Therefore,

where (a) follows from expanding the square and substituting in the definitions in (18) and (28); and (b) follows from the fact that $\overline{\alpha}_{u0}(t\!+\!1)=\overline{\alpha}_{u1}(t\!+\!1)$ proven above. We have thus proven that

Therefore, the asymptotic correlation coefficient is given by

where (a) follows from (29) and (b) follows from (95). Substituting in (93) into (96) proves the first equation in (37). The second recursion in (37) can be proven similarly.

Hence, from (30), the initial correlation coefficient is

which agrees with the statement in the theorem. This proves part (a).

To prove part (b), we again use Lemma A.4. Define the functions

From (37), it follows that $\rho_{v}(t\!+\!1)=H(\rho_{v}(t))$. Now, ${\mathcal{E}}_{u}(\eta_{u})$ and ${\mathcal{E}}_{v}(\eta_{v})$ defined in (36) are the mean-squared errors of $U_{0}$ and $V_{0}$ under AWGN estimation measurements with SNRs $\eta_{u}$ and $\eta_{v}$. Therefore, ${\mathcal{E}}_{u}(\eta_{u})$ and ${\mathcal{E}}_{v}(\eta_{v})$ must be monotonically decreasing in $\eta_{u}$ and $\eta_{v}$. Therefore, $H_{u}(\rho_{v})$ and $H_{v}(\rho_{u})$ in (97) are monotonically increasing functions and thus so is the concatenated function $H(\rho_{v})$ in (98). Also, since the assumption of part (b) is that ${\mathcal{E}}_{u}(\eta_{u})$ and ${\mathcal{E}}_{v}(\eta_{v})$ are continuous, $H(\rho_{v})$ is also continuous. It follows from Lemma A.4 that $\rho_{v}(t)\rightarrow\rho_{v}^{*}$ where $\rho_{v}^{*}$ is a fixed point of (37).

It remains to show $\rho_{v}^{*}>0$. Observe that

Therefore, the limit point $\rho_{v}^{*}$ of $\rho_{v}(t)$ must satisfy $\rho_{v}^{*}\geq\rho_{v}(0)>0$.

A.6 Proof of Lemma 5.5

Towards this end, we first use a standard result (see, for example, ) that, for any prior on $U_{0}$ with bounded second moments, the mean-squared error in (36) satisfies

The term $\tau_{u}/(1+\eta_{u}\tau_{u})$ is the minimum mean-squared error with linear estimation of $U_{0}$ from an AWGN noise-corrupted measurement $Y=\sqrt{\eta_{u}}U_{0}+D$, $D\sim{\mathcal{N}}(0,1)$. Equation (99) arises from the fact that linear estimation is optimal in low SNRs – see for details. Using (99), we can compute the derivative of $H_{u}$,

We now apply a standard linearization analysis of the nonlinear system $\rho_{v}(t\!+\!1)=H(\rho_{v}(t))$ around the fixed point $\rho_{v}=0$. See, for example, . If $\sqrt{\beta}\tau_{u}\tau_{v}<\tau_{w}$ then $H^{\prime}(0)<1$ and the fixed point is stable. Thus, for any $\rho_{v}(0)$ sufficiently small $\rho_{v}(t)\rightarrow 0$. This proves part (b) of the lemma.

On the other hand, if $\sqrt{\beta}\tau_{u}\tau_{v}>\tau_{w}$ then $H^{\prime}(0)>1$ and the fixed point is unstable. This will imply that for any $\rho_{v}(0)>0$, $\rho_{v}(t)$ will diverge from zero. But, we know from Theorem 5.3 that $\rho_{v}(t)$ must converge to some fixed point of $\rho_{v}=H(\rho_{v})$. So, the limit point must be positive. This proves part (a) of the lemma.

References