Near-optimal matrix recovery from random linear measurements

TL;DR

This paper introduces a fast, near-optimal matrix recovery algorithm using approximate message passing and singular value shrinkage, outperforming previous methods and approaching theoretical limits.

Contribution

The authors develop a novel AMP-based matrix recovery algorithm with exponential convergence and near-optimal measurement efficiency, surpassing existing methods like NNM.

Findings

Algorithm converges exponentially fast.

Phase transition curve is better than NNM.

Approaches the information-theoretic lower bound.

Abstract

In matrix recovery from random linear measurements, one is interested in recovering an unknown $M$-by-$N$ matrix $X_0$ from $n<MN$ measurements $y_i=Tr(A_i^T X_0)$ where each $A_i$ is an $M$-by-$N$ measurement matrix with i.i.d random entries, $i=1,\ldots,n$. We present a novel matrix recovery algorithm, based on approximate message passing, which iteratively applies an optimal singular value shrinker -- a nonconvex nonlinearity tailored specifically for matrix estimation. Our algorithm typically converges exponentially fast, offering a significant speedup over previously suggested matrix recovery algorithms, such as iterative solvers for Nuclear Norm Minimization (NNM). It is well known that there is a recovery tradeoff between the information content of the object $X_0$ to be recovered (specifically, its matrix rank $r$) and the number of linear measurements $n$ from which recovery is to be attempted. The precise tradeoff between $r$ and $n$, beyond which recovery by a given algorithm becomes possible, traces the so-called phase transition curve of that algorithm in the $(r,n)$ plane. The phase transition curve of our algorithm is noticeably better than that of NNM. Interestingly, it is close to the information-theoretic lower bound for the minimal number of measurements needed for matrix recovery, making it not only state-of-the-art in terms of convergence rate, but also near-optimal in terms of the matrices it successfully recovers.

Equations

$$X_{t+1}$$

$$z_t$$

$$X_{t+1}$$

$$z_t$$

$$b_t=\frac{1}{n}\nabla \cdot {\eta }_{t-1}\left(X_{t-1}+{\mathcal{A}}^{\ast }{z}_{t-1}\right).$$

$${\eta }_t(W)={\hat{\sigma }}_t\eta (W/{\hat{\sigma }}_t)$$

$${\eta }_{\lambda }^{soft}(y)=(y-\lambda {)}_{+}$$

$${\text{min}}_W{\left|\left|W\right|\right|}_{\ast }$$

$$\text{subject to}\mathcal{A}(W)=y,$$

$$\frac{M}{N}\to \beta \in (0,1],\frac{r}{M}\to \rho \in (0,1],\frac{n}{NM}\to \delta \in (0,1],$$

$${\delta }^{\ast }(\rho ,\beta )=\underset{N\to \infty }{\lim }{\delta }^{\ast }(N,M,r).$$

$${\delta }^{nnm}(\rho ,\beta )=\underset{N\to \infty }{\lim }\underset{\lambda }{\min }\underset{r(X)\le \rho M}{\max }\frac{1}{MN}\mathbb{E}{\left|\left|{\eta }_{\lambda }^{soft}(Y)-X\right|\right|}_F^2.$$

$${\delta }^{\ast }(\rho )=\mathcal{M}(\rho |\eta ).$$

$${\eta }^{spiked}(y)=\{\begin{array}{ll}\frac{1}{y}\sqrt{(y^2-\beta -1{)}^2-4\beta }& y\ge 1+\sqrt{\beta }\\ 0& y\le 1+\sqrt{\beta }\end{array}.$$

$$\underset{rank(X)\le r}{\max }\frac{1}{NM}\mathbb{E}{\left|\left|X-\hat{X}(X+Z)\right|\right|}_F^2\ge \frac{r}{M}+\frac{r}{N}-\frac{r^2+r}{MN}.$$

$$\underset{N\to \infty }{\lim }\frac{r}{M}+\frac{r}{N}-\frac{r^2+r}{MN}=\rho (1+\beta -\beta \rho ),$$

$${\delta }^{it}(\rho ,\beta )=\underset{N\to \infty }{\lim }\frac{r(N+M-r)}{NM}=\rho (1+\beta -\beta \rho ).$$

$$\widehat{{\sigma }_t}=\frac{1}{{\Phi }^{-1}(0.75)}median(\left|z_t\right|),$$

$$\begin{array}{rl}\nabla \cdot \eta (X)& =\sum _{i=1}^M\frac{d\eta }{d\sigma }({\sigma }_i)+2\sum _{i<j}\frac{\eta ({\sigma }_i){\sigma }_i-\eta ({\sigma }_j){\sigma }_j}{{\sigma }_i^2-{\sigma }_j^2}\\ & +(N-M)\sum _{i=1}^M\frac{\eta ({\sigma }_i)}{{\sigma }_i}\end{array}$$

$${\eta }^{opt}(\sigma )=\alpha \sqrt{N}\cdot {\eta }^{spiked}\left(\frac{\sigma }{\alpha \sqrt{N}}\right),$$

$$\displaystyle\frac{d\eta^{spiked}}{d\sigma}(\sigma)=\Bigg{(}-\frac{1}{\sigma^{2}}\sqrt{(\sigma^{2}-\beta-1)^{2}-4\beta}+\frac{2\sigma(\sigma^{2}-\beta-1)}{\sqrt{(\sigma^{2}-\beta-1)^{2}-4\beta}}\Bigg{)}$$

Full text

Modern datasets often take the form of large matrices. When the full dataset is not directly observable, the scientist can only obtain measurements, from which she hopes to recover the dataset. Let $X$ be our $M$-by-$N$ data matrix, and consider $n$ linear measurements of $X$, $y_{i}=Tr(A_{i}^{\top}X)$ ($i=1,\ldots,n$), where $A_{1},\ldots,A_{n}$ are measurement matrices. Basic linear algebra tells us that, in order to reconstruct $X$ from these measurements using a linear algorithm, successful recovery is only possible if $n\geq MN$. In recent years, intense research in applied mathematics, optimization, and information theory has shown that, when the rank $r=rank(X)$ is low, nonlinear algorithms based on convex optimization allow exact recovery of $X$ from only $O(rM+rN)$ measurements, up to logarithmic factors, whereby solving a severely underdetermined system of linear equations. In recent years, research in matrix recovery has flourished spanning both theory [5, 9, 19, 23, 39, 7] and efficient algorithms [4, 27, 46, 21, 45, 22]. Applications have been developed in fields ranging widely, from video and image processing [6, 49] and system identification [26, 39] to quantum state tomography [20] and collaborative filtering [23, 5]; see also the recent survey [10].

It is convenient to denote the linear measurement scheme by a linear operator $\mathcal{A}:\mathbb{R}^{M\times N}\to\mathbb{R}^{n}$, so that $y=\mathcal{A}(X)$ are the $n$ linear measurements available to the scientist. As a simple model for random linear measurements, one takes $A_{i}$ to be white Gaussian vectors, or equivalently, assumes that the entries of $\mathcal{A}$ are independent and identically distributed $\mathcal{N}(0,1/n)$ random variables. These are known as Gaussian measurements; more generally we will assume that the entries of $\mathcal{A}$ are i.i.d from some zero-mean distribution with variance $1/n$, and refer to this model simply as random linear measurements.

As matrix recovery problems contain $MN$ degrees of freedom, their size grows quickly with $M$ and $N$, making direct solvers infeasible; nearly all existing algorithms approach the problem iteratively. In fact, most matrix recovery algorithms proceed by some variation of iterative singular value shrinkage

[TABLE]

where $X_{t}$ is the current estimate for $X$, $\mu_{t}$ is a step size and $z_{t}\in\mathbb{R}^{n}$ is the current residual vector. Here, $\eta_{t}:\mathbb{R}^{M\times N}\to\mathbb{R}^{M\times N}$ is a singular value shrinker, namely, a matrix function that applies the same univariate nonlinearity to each of the singular values of its matrix argument. Abusing notation, here and below we write $\eta_{t}$ both for the univariate nonlinearity and for the corresponding matrix function.

In this paper we propose two new matrix recovery algorithms, based on the Approximate Message Passing (AMP) framework [15]. These algorithms offer significant improvement over variations of (1). Starting at $X_{0}=0$, we propose the Matrix AMP iteration

[TABLE]

where $\eta_{t}$ are a sequence of singular value denoisers and

[TABLE]

Here, $\nabla\cdot\eta_{t}$ is the divergence of the matrix function $\eta_{t}$. These formulas are a very natural extension of the AMP framework, originally developed for recovery of sparse vectors, to the matrix recovery problem. Any single singular value shrinker $\eta$ can yield a full-blown AMP matrix recovery algorithm by setting

[TABLE]

in (2), where $\hat{\sigma}_{t}$ is the current noise level estimate rigorously defined below. Comparison of Eqs. (1) and (2) reveals a subtle yet crucial difference between the two algorithms, known as the Onsager correction term, which is extensively discussed in the AMP literature [15, 16, 14, 32, 30].

The first algorithm we present, AMP-SVST, is a Matrix AMP iteration based on the famous soft thresholding nonlinearity

[TABLE]

with a tuning parameter that depends on $M$,$N$ and $r$ in a manner described below. The second algorithm we present, AMP-OPT, is a Matrix AMP iteration based on a variation of the asymptotically optimal singular value shrinker for low-rank matrix denoising [18]. The explicit form of shrinker, which we denote by $\eta^{opt}$, is based on (10) and (15) below. Figure 1 compares the two shrinkers with the hard thresholding nonlinearity.

The main discovery reported here is that by plugging the soft thresholding singular value shrinker $\eta^{soft}_{\lambda}$ and the optimal singular value shrinker $\eta^{opt}$ into the Matrix AMP iteration (2), *one obtains matrix recovery algorithms which, to the best of our knowledge, meet or exceed the state-of-the-art * in the case of random linear measurements, in the following two aspects:

(i) Convergence rate. Whenever recovery is possible, our algorithms typically converge exponentially fast in the number of iterations, a phenomenon that has been documented in AMP iterations for vector recovery [15, 14]. We bring substantial evidence that the convergence rate of our algorithms compares very favorably with state-of-the-art matrix recovery methods, including first order algorithms based on variations of (1).

(ii) Number of measurements required. Figure 2, which summarizes results of massive computer simulations described below, shows that AMP-OPT  requires a near-optimal number of measurements for successful recovery, as we now elaborate. To give a concrete example, when $r=rank(X)$ is low, AMP-OPT  requires $\approx r(N+M)$ measurements, which compares favorably with the $\approx 3r(N+M)$ measurements required by the popular Nuclear Norm Minimization algorithm, and agrees with the information-theoretic lower bound.

Phase transitions in matrix recovery

How many measurements $n$ are required for a given matrix recovery algorithm to correctly recover a matrix $X$ with $rank(X)=r$? Certainly as $r$ grows (so that the information content of the object to be recovered grows), $n$ must grow as well. Many variations of (1) solve the Nuclear Norm Minimization (NNM) convex program

[TABLE]

where $\left|\left|\cdot\right|\right|_{*}$ denotes the nuclear norm of a matrix, namely, the sum of its singular values. In the Gaussian measurements model, several authors have noted a phase transition phenomenon for NNM, whereby for a given value of $r=rank(X)$ there is a fairly precise number of samples $n(r,M,N)$ such that NNM from $n$ measurements typically succeeds if $n>n(r,M,N)$ and typically fails otherwise [39, 13, 33, 2]. Following these authors, we consider a sequence of recovery problems obeying a proportional growth model as $N\to\infty$,

[TABLE]

where $\beta$, $\rho$ and $\delta$ are known as the aspect ratio, rank fraction and undersampling ratio, respectively. Letting $\delta^{*}(r,M,N)=n(r,M,N)/(MN)$ be the location of the phase transition we consider the formal limit

[TABLE]

While most results in the literature focus on NNM, it is tempting to study any matrix recovery algorithm using the same lens, and compare competing algorithms using their phase transition curves, an approach that has proven very successful in compressed sensing [17, 12, 16, 28].

Phase transitions and minimax MSE

It has been empirically [13] and theoretically [2, 33] demonstrated that NNM (and any iterative singular value shrinkage algorithm converging to the NNM solution) exhibits a phase transition, which we will denote by $\delta^{{nnm}}(\rho,\beta)$. These authors show that $\delta^{{nnm}}(\rho,\beta)$ is given exactly by the minimax Mean Square Error (MSE) of matrix denoising by singular value soft thresholding over matrices of asymptotic rank fraction $\rho$ and aspect ratio $\beta$. Formally,

[TABLE]

Here, $Y=X+Z$ with $Z$ a matrix with i.i.d $\mathcal{N}(0,1)$ entries, and $\left|\left|\cdot\right|\right|_{F}^{2}$ is the squared Frobenius norm, namely the sum of squares of the matrix entries. The quantity on the right hand side has been calculated explicitly as a function of $\beta$ and $\rho$, see [11, 13] and the SI Appendix, section 5.

While (8) connects the phase transition of NNM, a convex optimization algorithm, to the minimax MSE of a matrix denoiser, a similar connection has been observed and verified for many different AMP iterations. Specifically, under mild conditions on a shrinker $\eta$, and under the assumption known as AMP State Evolution, [14, 3] found that the phase transition $\delta^{*}(\rho)$ of the AMP iteration, corresponding to $\eta$, matches $\mathcal{M}(\rho|\eta)$, the asymptotic minimax MSE of denoising a signal with rank-fraction $\rho$ using $\eta$. Formally,

[TABLE]

This fundamental connection has been observed for multiple vector recovery problems [14], with the rank fraction $\rho$ replaced by some other measure of sparsity or information content.

In this paper we observe that the fundamental connection (9) holds for the Matrix AMP iteration (2) with two very different choices of shrinker $\eta$. When $\eta=\eta^{soft}_{\lambda}$ (with $\lambda$ optimally tuned), one naturally conjectures that the phase transition of the corresponding Matrix AMP iteration will exactly match the asymptotic minimax MSE of matrix denoising by $\eta^{soft}_{\lambda}$, namely, will have the same phase transition curve as NNM, given by (8). While the literature points to a deep connection between decision theory (specifically, minimax MSE) and compressed sensing phase transitions in matrix recovery [13, 2, 33], this observation is, to the best of our knowledge, the first instance where such connection is made direct and explicit, in the same sense that [15] has identified a similar connection in sparse vector recovery problems.

Matrix denoising and optimal singular value shrinkage

The identity (9) as it holds for AMP-based matrix recovery algorithms suggests that by designing a singular value shrinker with lower asymptotic minimax MSE, one might obtain a matrix recovery algorithm with an improved phase transition. This can hypothetically be achieved by plugging in the improved shrinker into (4) and using the AMP iteration (2).

Interestingly, in an asymptotic model of low-rank matrices observed in white noise, simple, closed-form formulae are available for the optimal singular value shrinker [18, 44]. The optimal singular value shrinker can be shown to be

[TABLE]

Adapting the shrinker $\eta^{spiked}$ to the proportional growth model (6) requires technical detail and is deferred to the SI Appendix, see section 4.2. We denote the adapted version of (10), formally defined in (15) below, by $\eta^{opt}$. The shrinkers $\eta^{opt}$ and $\eta^{soft}_{\lambda}$ are shown in Figure 1 above.

Now, if indeed the adapted version of (10) has an appealing worst-case (hence minimax) MSE, one could hope, according to (9), that the corresponding AMP iteration would have an appealing phase transition.

Unification of lower bounds

There is reason to believe that the minimax MSE of the $\eta^{opt}$ is not only appealing, but in fact near-optimal. Indeed, [11] show that for any measurable matrix denoiser $\hat{X}:\mathbb{R}^{M\times N}\to\mathbb{R}^{M\times N}$ (not necessarily based on singular value shrinkage) the minimax MSE of $\eta$ is lower-bounded by

[TABLE]

Moreover, [11] show that this lower bound, if achieved, is achieved by a singular value shrinker. Empirical evidence, included in the SI Appendix (see section 4.3), suggest that the minimax MSE of $\eta^{opt}$ is close to the lower bound (11).

The connection (9) suggests that no Matrix AMP iteration can achieve a phase transition lower than the limit of (11), namely

[TABLE]

and that the Matrix AMP iteration based on $\eta^{opt}$ would have a phase transition close to the best phase transition achievable by any matrix AMP iteration.

Surprisingly, simple dimension considerations imply that no matrix recovery (AMP or other) can have better uniform guarantees: there is a simple information-theoretic lower bound on the phase transition one could possibly hope for in a matrix recovery algorithm, which we denote by $\delta^{{it}}(\rho,\beta)$. Indeed, as the set of rank-$r$ $M$-by-$N$ matrices is a smooth manifold of dimension $r(M+N-r)$ embedded in $\mathbb{R}^{M\times N}$ (see, e.g, [47]), faithful recovery of any matrix in the manifold requires at least $n\geq r(N+M-r)$ linear measurements. It follows that

[TABLE]

This lower bound agrees, to order $1/N$, with (11), which stems from an altogether different consideration. In other words, if a Matrix AMP iteration based on $\eta^{opt}$ is shown to have a phase transition close to the curve $\delta(\rho)=\rho(1+\beta-\beta\rho)$, as predicted by (9), then it is in fact near-optimal among any matrix recovery algorithm whatsoever, in terms of the number of measurements required for successful recovery.

Detailed description of the proposed algorithms

Complete specification of the proposed algorithms requires full specification of $\nabla\cdot\eta_{t}$ from (2), as well as $\hat{\sigma}_{t}$ and $\eta$ from (4). To estimate the noise level $\hat{\sigma}_{t}$ we can use [15, 14]

[TABLE]

where $z_{t}$ is defined in (2). The divergence term in (2) for a general singular value shrinker is given by [8, 11]

[TABLE]

where $\sigma_{1}>\ldots>\sigma_{M}>0$ are the singular values of the matrix argument $X$, whose spectrum is assumed to be non-degenerate. (Indeed, degenerate spectrum occurs with probability zero in any matrix model of interest). The derivative $d\eta/d\sigma$ is shown next for each choice of shrinker.

Algorithm 1: AMP-SVST.

We use $\eta^{soft}_{\lambda}$ (5) with a tuning parameter $\lambda=\sqrt{N}\lambda^{*}(\rho,\beta)$ achieving the asymptotic minimax MSE in the right hand side of (8). Details of calculating $\lambda^{*}(\rho,\beta)$ are deferred to the SI Appendix (see section 5). The (weak) derivative of $\eta^{soft}_{\lambda}$ is simply $d\eta^{soft}_{\lambda}/d\sigma=\mathbf{1}_{\{\sigma>\lambda\}}$ .

Algorithm 2: AMP-OPT.

We use a properly calibrated version of $\eta^{spiked}$ from (10), defined by

[TABLE]

where $\alpha=(\sqrt{1-\beta\rho}+\sqrt{\beta-\beta\rho})/(1+\sqrt{\beta})$. The full derivation of (15) and $\alpha$ is deferred to the SI Appendix, section 4.2. Now, the derivative of $\eta^{spiked}$ from (10) is given by

[TABLE]

when $\sigma>1+\sqrt{\beta}$, so that the derivative of $\eta^{opt}$ is given by $(d\eta^{spiked}/d\sigma)\left(\sigma/(\alpha\sqrt{N})\right)$.

Main hypotheses

This brief announcement tests two hypotheses regarding the merits of our proposed algorithms, by conducting substantial computer experiments generating large numbers of random problem instances.

Phase transitions.

Based on (9) we hypothesize that the phase transition of the Matrix AMP algorithm AMP-SVST  matches (8) exactly, and that the phase transition of AMP-OPT  is close to the lower bound (12).

Convergence rates.

One of the appealing properties of AMP algorithms for sparse vector recovery is their exponential rate of convergence [15, 14] – see SI Appendix section 6.3 for more details. We hypothesize that when $NM$ is large both AMP-SVST  and AMP-OPT  similarly demonstrate an exponential rate of convergence whenever recovery is possible.

Methods

We use statistical methods to check for agreement between the hypothesis and the predicted phase transition. We further compare both phase transitions and convergence rates of our proposed algorithms with two different state-of-the-art algorithms for matrix recovery.

Comparison with state-of-the-art.

The Matrix AMP framework we propose, and in particular the proposed algorithms AMP-SVST  and AMP-OPT, are essentially iterative singular value thresholding methods. We therefore compare them with well-known algorithms based on variations of iterative singular value shrinkage, (1). (SI Appendix section 7.5 also compares AMP-OPT  with an involved, state-of-the-art AMP algorithm, developed from first principles, which has a completely different form.) In order to compare convergence rates of the algorithms close to their respective phase transitions, we match AMP-SVST  (resp. AMP-OPT) with an algorithm whose phase transition is expected to be similar.

Accelerated proximal gradient singular value thresholding (APG) [46]: An accelerated proximal gradient descent algorithm for NNM, based on soft thresholding of singular values. As an NNM solver, APG  is expected to have the same phase transition as AMP-SVST. Being a simple, general-purpose solver, it serves as a reasonable baseline for our evaluation111Some related gradient descent algorithms have been shown to converge exponentially fast, e.g. [1]. However they sometimes fail to converge to the NNM solution [34] and are therefore not suitable candidates for comparison to AMP-SVST. . 2. 2.

Normalized iterated hard thresholding (NIHT) [45]: An alternating projection algorithm, with an adaptive stepsize, based on hard thresholding of singular values. In [45] the authors bring some evidence that this algorithm exhibits a phase transition which is not far from the information-theoretical lower bound. To the best of our knowledge, its convergence rate is also state-of-the-art.

Both algorithms are described in detail in the SI Appendix, section 3. We studied both the computational complexity and the phase transition of all four algorithms, as follows.

Evaluating convergence rates.

All four algorithms under study are iterative, with roughly the same per-iteration complexity. To compare their overall complexity, we measured the rate of convergence of each algorithm. We recorded the relative error $\Delta_{t}:=\left|\left|X_{t}-X\right|\right|_{F}/\left|\left|X\right|\right|_{F}$ of each algorithm under study over the number of iterations $t$, with $X$ being the true (unknown) matrix to be recovered. Due to space considerations, a detailed account of our methodology is deferred to the SI Appendix, see section 2.

Evaluating phase transitions.

Our analysis follows the methods of [12, 31, 13]. For each of the proposed algorithms we assigned a suspected phase transition $\delta^{*}$. The suspected phase transition for AMP-SVST  was $\delta^{*}=\delta^{{nnm}}$ from (8), and for AMP-OPT  and NIHT  we used $\delta^{*}=\delta^{{it}}$ from (12). For several problem dimensions $N$ and parameter combinations $M=\beta N,r=\rho M$ we performed many random recovery experiments, using undersampling ratios $\delta$ above and below $\delta^{*}$. We used matrices with fully degenerate spectrum, a configuration known to be least favorable for AMP-SVST  but possibly not for AMP-OPT  (see important discussion in the SI Appendix, section 6.2). At each point in parameter space, we chronicled the empirical probability of success $\pi(r|n,N,M)$. We then fitted a logit curve $\log(\hat{\pi}/(1-\hat{\pi}))=a+b(\delta-\delta^{*})$. Setting $\pi=1/2$, we obtain an estimate $\hat{\delta}=\delta^{*}-a/b$ of the real phase transition at $\rho,\beta$. Due to space consideration, we defer further details of our methodology to the SI Appendix (see section 2).

Simulation software platform.

Empirical evaluation of matrix recovery phase transitions requires computer simulation on a massive scale, which poses a software development challenge. To efficiently conduct the simulation reported in this paper, we developed a dedicated software platform inspired by [13, 31]. Our framework is written in Python and is made available on the Data and Code Supplement [40]. We used the Spark parallelization framework [48] to orchestrate the parallel computations required, and executed the code on large Amazon Web Services machines using roughly $15,000$ CPU hours.

Results

††\dagger††\dagger$\dagger$Reproducibility advisory. All the figures and tables in this paper, including those in the SI Appendix, are fully reproducible. Raw simulation results were stored in a mysql database from which all figures were generated. All our scripts, as well as a snapshot of the full result database, is available in the Data and Code Supplement [40].

Convergence rates.

Figure 3 shows typical convergence rates of the algorithms under study above their respective phase transitions for $\rho=0.2$ and $\beta=1$ contrasting AMP-SVST  with APG  (at $\delta=\delta^{nnm}(\rho,\beta)+0.1$, close to their common phase transition) and AMP-OPT  with NIHT  (at $\delta=\delta^{it}(\rho,\beta)+0.1$, close to the information-theoretical lower bound). When successful, matrix AMP algorithms exhibit roughly exponential convergence rates; for example, after $T=1000$ iterations, the relative error for AMP-OPT  is $10^{-14}$, compared with $10^{-5}$ for NIHT. Results for additional values of $\rho$ and $\beta$ are deferred to the SI Appendix (see section 7.3). Our results support the hypothesis that Matrix AMP algorithms converge exponentially fast when recovery is possible.

Phase transitions.

Figure 2 shows the estimated phase transition of the algorithms AMP-SVST, AMP-OPT, NIHT, as well as the asymptotic phase transition of APG  for $\beta=1$ and $N=100$. Similar results for $\beta=0.5$, $N=100$ appear in the SI Appendix, section 7.2 (for all three algorithms above), and for $N\in\left\{50,100,150\right\}$ and $\beta\in\left\{0.2,0.4,0.6,0.8,1.0\right\}$ (for AMP-OPT  and AMP-SVST). We found that AMP-SVST  exhibits a sharp phase transition, matching the phase transition of NNM and APG. Our results conclusively affirm the relation (9) for AMP-SVST. We further found that AMP-OPT  exhibits a phase transition that is less pronounced that previously studied AMP algorithms, see Discussion below. These results support the hypothesis that the phase transition of AMP-OPT  lies very close to the information-theoretic lower bound (12).

Discussion

Our contributions.

We present overwhelming evidence that the AMP framework, previously studied for sparse vector recovery and related problems, naturally extends to matrix recovery, with the notion of matrix rank replacing that of vector sparsity. We validate the correspondence between the phase transition of AMP-SVST  and the minimax MSE of the underlying shrinker $\eta^{soft}_{\lambda}$. Our results regarding both AMP-SVST  and AMP-OPT  suggests Matrix AMP as an appealing framework for designing matrix recovery algorithms. 2. 2.

We present two algorithms, AMP-SVST  and AMP-OPT, which converge exponentially fast when they succeed. To the best of our knowledge, these algorithms converge at a rate which meets or exceeds the state-of-the-art in matrix recovery from random linear measurements, sometimes by a significant margin. As discussed below, in contrast with previously suggested iterative algorithms for matrix recovery, both AMP-SVST  and AMP-OPT  offer clear “diagnostic” for whether recovery was successful. 3. 3.

We show that while algorithm AMP-OPT  improves on the state of the art in terms convergence rate, it also offers a near-optimal phase transition for $N,M\gtrsim 100$ (see important comments below).

Universality.

To test whether our results depend on the particular distribution used to generate the entries of the measurement matrix $\mathcal{A}$ we performed local evaluation of the phase transition of AMP-SVST  and AMP-OPT  using several distributions: Gaussian (light tail), Rademacher ($\pm 1$ with equal probability – no tail) and Student-t with $\nu=6$ degrees of freedom (heavy tail). All distributions are symmetric about [math], and were normalized to have variance $1/n$. For details of this experiment, see SI Appendix, section 7.4. We observe (see Figure 4 and additional results in the SI Appendix, section 7.4) that the phase transitions of AMP-SVST  and AMP-OPT  are the same, to high precision, under these three different distributions of the measurements matrix. Other authors have previously observed universality in other signal recovery settings: [12] observed universality in sparse vector recovery using $l_{1}$ norm minimization; [13] observed in passing universality for matrix recovery with NNM; [14] observed universality for AMP in various vector recovery-related problems.

State evolution.

We observe that Matrix AMP with our two denoisers behaves as predicted by the AMP formalism of State Evolution (SE) [15, 14]. Background and empirical evidence are presented in the SI Appendix, see section 6.

Other approaches to recovery by Approximate Message Passing.

The term Approximate Message Passing (AMP) is somewhat overloaded in the literature. Previous works in Approximate Message Passing algorithms for low-rank matrix recovery propose a generative model for low-rank matrices and derive involved AMP algorithms from first principles [43, 35, 36, 37, 24, 29, 38]. In contrast, this paper follows the framework of [15],

which delineates simple, highly structured iterative thresholding algorithms, whose design basically boils down to the design of the denoising nonlinearity. In SI Appendix Section 7.5 we demonstrate that AMP-OPT  compares favorably, for the case of random linear measurements, with P-BiG-AMP, a notable AMP algorithm for matrix recovery, derived from first principles.

Stopping conditions and sign of convergence.

State Evolution suggests that the estimated noise level $\hat{\sigma}_{t}$ in every iteration is in fact a good proxy for the current MSE: $\hat{\sigma}_{t}^{2}\approx\frac{1}{\delta NM}\left|\left|X_{t}-X\right|\right|^{2}$. In contrast with previously suggested iterative methods for matrix recovery, Matrix AMP algorithms therefore have the important advantage of offering continuous diagnostic measurements for the quality of the reconstruction $X_{t}$, and can choose to stop or continue accordingly. Background and empirical evidence are presented in the SI Appendix, section 6.4.

Finite-$N$ width of the AMP-OPT  phase transition.

Previously studied AMP algorithms for vector recovery problems demonstrated that, in finite problems, the transition from failed to successful recovery occurs in a transition region whose width shrinks with problem size. We observed that while the transition region of AMP-OPT  does shrink with $N$, it is wider and shrinks more slowly than the transition region of AMP-SVST  and previously studied AMP algorithms for vector recovery, see for instance Figure 4. Furthermore, evidence included in the SI Appendix suggests that below its phase transition, AMP-OPT  fails always, while above its phase transition it succeeds with probability that approaches $1$ rather slowly as $N$ grows, see SI Appendix Section 7.1. We believe this happens because the shrinker $\eta^{opt}$ is not Lipschitz continuous (see Figure 1) [3, 41, 42].

Optimality.

We conjecture that the information-theoretical lower bound $n/(NM)\geq\frac{r}{M}(1+M/N-r/M)$ and the denoising lower bound (11) are both asymptotically tight. Also, one naturally wonders whether the denoiser $\eta^{opt}$ in (15) is in fact optimal in the (asymptotic) minimax sense. Numerical evidence suggests that this is not the case, at least when $\rho$ is sufficiently large - see SI Appendix, section 4.3. Nonetheless, the worst-case asymptotic MSE of $\eta^{opt}$ seems to be very close to the minimax MSE.

Implementation with fast SVD.

It is natural to use fast randomized SVD methods such as [25] to accelerate the SVD step in the proposed algorithms. Note however that these methods introduce non-trivial complications and noise sensitivities which are likely to affect the phase transition.

Conclusion

This paper presents a near-optimal matrix recovery algorithm based on the Approximate Message Passing framework (AMP). The same ideas underlying AMP in previously studied vector recovery problems have been shown to hold intact in the matrix recovery problem. In particular, our results show that the AMP framework bridges the seemingly unrelated problems of matrix denoising on the one hand and matrix recovery from partial measurements on the other hand. Indeed, design of a near-optimal matrix denoiser $\eta^{opt}$ has lead to AMP-OPT, a near-optimal matrix recovery algorithm. For example, when $r=rank(X)$ is low, AMP-OPT  requires $n\approx r(N+M)$ measurements for recovery, whereas NNM requires $n\approx 3r(N+M)$. Our results point to a new approach for designing matrix recovery algorithms, via an interesting connection with statistical estimation and decision theory.

Acknowledgement

We thank the anonymous reviewers for many helpful comments. This work has been partially supported by Israeli Science Foundation grant no. 1523/16 and German-Israeli Foundation for Scientific Research and Development no. I-1100-407.1-2015. ER was partially supported by the HUJI Leibniz center.