Near-optimal matrix recovery from random linear
measurements

arXiv:1705.09958v1 [cs.IT] 28 May 2017

Elad Romanov1 and Matan Gavish
1

∗1

School of Computer Science and Engineering, The Hebrew University,
Jerusalem, Israel

Abstract
In matrix recovery from random linear measurements, one is interested in recovering
an unknown M -by-N matrix X0 from n < M N measurements yi = T r(A>
i X0 ) where
each Ai is an M -by-N measurement matrix with i.i.d random entries, i = 1, . . . , 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 often 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 recovery tradeoff between the information content of the object
X0 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 the fastest algorithm known, but also near-optimal
in terms of the matrices it successfully recovers.

∗

To whom correspondence should be addressed. E-mail: gavish@cs.huji.ac.il

1

Contents
1 Introduction
1.1 Phase transitions in matrix recovery . . . . . . . . . . .
1.2 Phase transitions and minimax MSE . . . . . . . . . .
1.3 Matrix denoising and optimal singular value shrinkage
1.4 Unification of lower bounds . . . . . . . . . . . . . . .
1.5 Detailed description of the proposed algorithms . . . .
1.6 Main hypotheses . . . . . . . . . . . . . . . . . . . . .

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

.
.
.
.
.
.

3
4
6
7
8
8
9

2 Methods

10

3 Results

11

4 Discussion

14

5 Conclusion

15

2

1

Introduction

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, yi = T r(A>
i X) (i = 1, . . . , n), where A1 , . . . , An 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 ≥ M N . 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 [4, 8, 18, 22, 30, 6] and efficient
algorithms [3, 24, 37, 20, 36, 21]. Applications have been developed in fields ranging widely,
from video and image processing [5, 40] and system identification [23, 30] to quantum state
tomography [19] and collaborative filtering [22, 4]; see also the recent survey [9].
It is convenient to denote the linear measurement scheme by a linear operator A :
RM ×N → Rn , so that y = A(X) are the n linear measurements available to the scientist. As a simple model for random linear measurements, one takes Ai to be white Gaussian
vectors, or equivalently, assumes that the entries of A are independent and identically distributed N (0, 1/n) random variables. These are known as Gaussian measurements; more
generally we will assume that the entries of 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 M N 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
Xt+1 = ηt (Xt + µt A∗ zt )
zt = y − AXt ,

(1)

where Xt is the current estimate for X, µt is a step size and zt ∈ Rn is the current residual
vector. Here, ηt : RM ×N → RM ×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 η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 [14]. These algorithms offer significant improvement
over variations of (1). Starting at X0 = 0, we propose the Matrix AMP iteration
Xt+1 = ηt (Xt + A∗ zt )
zt = y − AXt + bt zt−1

3

(2)

where ηt are a sequence of singular value denoisers and
1
∇ · ηt−1 (Xt−1 + A∗ zt−1 ) .
(3)
n
Here, ∇ · ηt is the divergence of the matrix function η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 η can yield a full-blown AMP
matrix recovery algorithm by setting
bt =

ηt (W ) = σ̂t η(W/σ̂t )

(4)

in (2), where σ̂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
[14, 15, 13, 28, 26].
The first algorithm we present, AMP-SVST, is a Matrix AMP iteration based on the
famous soft thresholding nonlinearity
ηλsof t (y) = (y − λ)+

(5)

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 [17]. The
explicit form of shrinker, which we denote by η 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 ηλsof t and the optimal singular value shrinker η opt into the Matrix AMP iteration
(2), one obtains matrix recovery algorithms which are superior, to the best of our knowledge,
to any known matrix recovery algorithm in the case of random linear measurements, in the
two following 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 [14, 13]. In contrast, state of the art first order algorithms
based on variations of (1) exhibit significantly slower convergence rates.
(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 ≈ r(N +M ) measurements, which
compares favorably with the ≈ 3r(N + M ) measurements required by the popular Nuclear
Norm Minimization algorithm, and agrees with the information-theoretic lower bound.

1.1

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
4

10

Soft tresholding
Hard tresholding
Optimal shrinkage

8
6
4
2
0

0

2

4

6

8

10

Figure 1: Three common singular value shrinkers. AMP-SVST is based on iterative soft
thresholding (5), whereas AMP-OPT uses a variation of the optimal shrinker given in
(10). Hard thresholding is used in the iterative hard thresholding family of algorithms, see
e.g [36]. All the plotted shrinkers are tuned for β = 1.
of the object to be recovered grows), n must grow as well. Many variations of (1) solve the
Nuclear Norm Minimization (NNM) convex program
minW ||W ||∗
subject to A(W ) = y ,
where ||·||∗ 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 [30, 12, 29, 1]. Following these authors, we
consider a sequence of recovery problems obeying a proportional growth model as N → ∞,
r
n
M
→ β ∈ (0, 1],
→ ρ ∈ (0, 1],
→ δ ∈ (0, 1] ,
(6)
N
M
NM
where β, ρ and δ are known as the aspect ratio, rank fraction and undersampling ratio,
respectively. Letting δ ∗ (r, M, N ) = n(r, M, N )/(M N ) be the location of the phase transition
we consider the formal limit
δ ∗ (ρ, β) = lim δ ∗ (N, M, r) .
N →∞

(7)

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
5

1.0

0.8

δ

0.6

0.4

0.2

0.0
0.0

Information-theoretic lower bound
NNM phase transition
PˆT : AMP-OPT
PˆT : AMP-SVST
PˆT : NIHT (Tanner and Wei 2013)

0.1

0.2

0.3

0.4

0.5

0.6

ρ

Figure 2: Evaluated phase transition of the proposed algorithms AMP-SVST and AMPOPT; evaluated phase transition of state-of-the-art algorithm NIHT; information theoretical lower bound (12); and theoretical asymptotic phase transition of NNM (8), explicitly
provided in the Supporting Information [32]. Gaussian observations, N = M = 100 (β = 1).
See Methods section below for details. Similar figures for other values of β appear in the
Supporting Information.
transition curves, an approach that has proven very successful in compressed sensing [16, 11,
15, 25].

1.2

Phase transitions and minimax MSE

It has been empirically [12] and theoretically [1, 29] demonstrated that NNM (and any
iterative singular value shrinkage algorithm converging to the NNM solution) exhibits a
phase transition, which we will denote by δ nnm (ρ, β). These authors show that δ nnm (ρ, β) is
given exactly by the minimax Mean Square Error (MSE) of matrix denoising by singular value
soft thresholding over matrices of asymptotic rank fraction ρ and aspect ratio β. Formally,
δ nnm (ρ, β) = lim min max
N →∞

λ

r(X)≤ρM

1
E ηλsof t (Y ) − X
MN

2
F

.

(8)

Here, Y = X + Z with Z a matrix with i.i.d N (0, 1) entries, and ||·||2F 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 β and ρ, see [10, 12] and the
Supporting Information [32]. This striking identity is the only known instance in matrix
6

recovery of a deep connection between decision theory (specifically, minimax MSE) and
compressed sensing phase transitions, a phenomenon that has been first observed in sparse
vector recovery problems [14]. 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 η, and under the assumption known as AMP State Evolution, [13, 2]
found that the phase transition δ ∗ (ρ) of the AMP iteration, corresponding to η, matches
M(ρ|η), the asymptotic minimax MSE of denoising a signal with rank-fraction ρ using η.
Formally,
δ ∗ (ρ) = M(ρ|η) .

(9)

This fundamental connection has been observed for multiple vector recovery problems, with
the rank fraction ρ replaced by some other measure of sparsity or information content.
In this paper we observe that this fundamental connection holds for the Matrix AMP
iteration (2) with two very different choices of shrinker η. When η = ηλsof t (with λ 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
ηλsof t , namely, will have the same phase transition curve as NNM, given by (8).

1.3

Matrix denoising and optimal singular value shrinkage

To the best of our knowledge, the identity (9) has not been observed in AMP-based matrix
recovery algorithms so far. Still, it 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 [17, 35]. The
optimal singular value shrinker can be shown to be
η

spiked

(y) =

 q
 1 (y 2
y

− β − 1)2 − 4β

0

√
y ≥1+ β
√ .
y ≤1+ β

(10)

Adapting the shrinker η spiked to the proportional growth model (6) requires technical detail
and is deferred to the Supporting Information [32]. We denote the adapted version of (10),
formally defined in (15) below, by η opt . The shrinkers η opt and ηλsof t 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.

7

1.4

Unification of lower bounds

There is reason to believe that the minimax MSE of the η opt is not only appealing, but in fact
near-optimal. Indeed, [10] show that for any measurable matrix denoiser X̂ : RM ×N → RM ×N
(not necessarily based on singular value shrinkage) the minimax MSE of η is lower-bounded
by
2
r
1
r
r2 + r
max
E X − X̂(X + Z) ≥
+
−
.
(11)
F
rank(X)≤r N M
M
N
MN
Moreover, [10] show that this lower bound, if achieved, is achieved by a singular value
shrinker. Empirical evidence, included in the Supporting Information [32], suggest that the
minimax MSE of η 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
r
r2 + r
r
+
−
= ρ(1 + β − βρ) ,
N →∞ M
N
MN
lim

and that the Matrix AMP iteration based on η 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 δ it (ρ, β). Indeed, as the set of rank-r M -by-N matrices is a smooth
manifold of dimension r(M + N − r) embedded in RM ×N (see, e.g, [38]), faithful recovery
of any matrix in the manifold requires at least n ≥ r(N + M − r) linear measurements. It
follows that
r(N + M − r)
= ρ(1 + β − βρ) .
(12)
δ it (ρ, β) = lim
N →∞
NM
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 η opt is shown to have a
phase transition close to the curve δ(ρ) = ρ(1 + β − βρ), 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.

1.5

Detailed description of the proposed algorithms

Complete specification of the proposed algorithms requires full specification of ∇ · ηt from
(2), as well as σ̂t and η from (4). To estimate the noise level σ̂t we can use [14, 13]
σ̂t =

1
Φ−1 (0.75)

median(|zt |) ,

(13)

where zt is defined in (2). The divergence term in (2) for a general singular value shrinker is
given by [7, 10]
8

∇ · η(X) =

M
X

X η(σi )σi − η(σj )σj
dη
(σi ) + 2
σi2 − σj2
i=1 dσ
i<j

(14)

M
X

η(σi )
+ (N − M )
i=1 σi
where σ1 > . . . > σ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η/dσ is shown next for each
choice of shrinker.
√
Algorithm 1: AMP-SVST. We use ηλsof t ((5)) with a tuning parameter λ = N λ∗ (ρ, β)
achieving the asymptotic minimax MSE in the right hand side of (8). Details of calculating
λ∗ (ρ, β) are deferred to the Supporting Information [32]. The (weak) derivative of ηλsof t is
simply dηλsof t /dσ = 1{σ>λ} .
Algorithm 2: AMP-OPT. We use a properly calibrated version of η spiked from (10),
defined by
!
√
σ
√
,
(15)
η opt (σ) = α N · η spiked
α N
√
√
√
where α = ( 1 − βρ + β − βρ)/(1 + β). The full derivation of (15) and α is deferred to
the Supporting Information [32]. Now, the derivative of η spiked from (10) is given by




2σ(σ 2 − β − 1) 
1q
dη spiked
(σ) =  − 2 (σ 2 − β − 1)2 − 4β + q
dσ
σ
(σ 2 − β − 1)2 − 4β
when σ > 1 +

1.6

√


√ 
β, so that the derivative of η opt is given by (dη spiked /dσ) σ/(α N ) .

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 AMPOPT is close to the lower bound (12).

9

Convergence rates. One of the appealing properties of AMP algorithms for sparse vector recovery is their exponential rate of convergence [14, 13] – see Supporting Information
for more details. We hypothesize that when N M is large both AMP-SVST and AMPOPT similarly demonstrate an exponential rate of convergence whenever recovery is possible.

2

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. We compared the proposed algorithms AMP-SVST and
AMP-OPT with two state-of-the-art algorithms based on variations of iterative singular
value shrinkage:
1. Normalized iterated hard thresholding (NIHT) [36]: An alternating projection algorithm, with an adaptive stepsize, based on hard thresholding of singular values.
2. Accelerated proximal gradient singular value thresholding (APG) [37]: An accelerated
proximal gradient descent algorithm for NNM, based on soft thresholding of singular
values.
Both algorithms are described in detail in the Supporting Information. 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 ∆t :=
||Xt − X||F / ||X||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 Supporting Information.
Evaluating phase transitions. Our analysis follows the methods of [11, 27, 12]. For each
of the proposed algorithms we assigned a suspected phase transition δ ∗ . The suspected phase
transition for AMP-SVST was δ ∗ = δ nnm from (8), and for AMP-OPT and NIHT we
used δ ∗ = δ it from (12). For several problem dimensions N and parameter combinations M =
βN, r = ρM we performed many random recovery experiments, using undersampling ratios δ
above and below δ ∗ . 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 Supporting Information). At each point in parameter space, we chronicled
the empirical probability of success π(r|n, N, M ). We then fitted a logit curve log(π̂/(1 −
π̂)) = a + b(δ − δ ∗ ). Setting π = 1/2, we obtain an estimate δ̂ = δ ∗ − a/b of the real phase
10

transition at ρ, β. Due to space consideration, we defer further details of our methodology
to the Supporting Information.
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 [12, 27]. Our framework is written in Python and is
made available on the Data and Code Supplement [31]. We used the Spark parallelization
framework [39] to orchestrate the parallel computations required, and executed the code on
large Amazon Web Services machines using roughly 15, 000 CPU hours.

3

Results

Reproducibility advisory. All the figures and tables in this paper, including those in
the Supporting Information, 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 [31].
0

log10 (kXt − Xk/kXk)

−2
−4
−6
−8

−10
−12

NIHT
AMP-OPT
NNM using APG
AMP-SVST

−14
−16

0

200

400

600

800

1000

# iteration

Figure 3: Comparison of convergence rates: relative error (log scale) over iteration number.
Both matrix AMP algorithms exhibit exponential convergence in iteration number until machine precision is reached (NIHT appears to converge exponentially as well, yet significantly
slower). See methods section for further details. N = M = 50, ρ = 0.2, at undersampling
ratio δ ∗ + 0.1.

11

Convergence rates. Figure 3 shows typical convergence rates of the algorithms under
study above their respective phase transitions for ρ = 0.2 and β = 1 contrasting AMPSVST with APG and AMP-OPT with NIHT. 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 ρ and β are deferred to the Supporting Information. 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
β = 1 and N = 100. Similar results for β = 0.5, N = 100 appear in the Supporting Information (for all three algorithms above), and for N ∈ {50, 100, 150} and β ∈
{0.2, 0.4, 0.6, 0.8, 1.0} (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).

12

Empirical success probability

1.0
gaussian
rademacher
t(6)

0.8

0.6

0.4

0.2

0.0
−0.05

0.00

PˆT − P T = 0.0145
PˆT − P T = 0.0112
PˆT − P T = 0.0144
0.05
∗

0.10

0.15

δ−δ

Empirical success probability

1.0
gaussian
rademacher
t(6)

0.8

0.6

0.4

0.2

0.0
−0.06

−0.04

−0.02

PˆT − P T = 0.00413
PˆT − P T = 0.00244
PˆT − P T = 0.00404
0.00
0.02
0.04
∗

0.06

0.08

0.10

δ−δ

Figure 4: Evidence of universality. Empirical success probability (with logit fit and estimated
phase transition) for various probability distributions of the entries of the measurement
matrix A. Top: AMP-OPT. Bottom: AMP-SVST. M = N = 100, ρ = 0.3.

13

4

Discussion

Our contributions.
1. 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 ηλsof t . Our results regarding both AMP-SVST and AMPOPT suggests Matrix AMP as an appealing framework for designing matrix recovery
algorithms.
2. We present two algorithms, AMP-SVST and AMP-OPT, which converge exponentially fast when they succeed. This makes them, to the best of our knowledge,
the fastest matrix recovery algorithms known 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. We show that algorithm AMP-SVST can safely replace NNM methods for N, M &
100 as it offers an identical phase transition, with significantly faster convergence.
4. 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 & 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 A we performed local evaluation of the phase
transition of AMP-SVST and AMP-OPT using several distributions: Gaussian (light
tail), Rademacher (±1 with equal probability – no tail) and Student-t with ν = 6 degrees
of freedom (heavy tail). All distributions are symmetric about 0, and were normalized
to have variance 1/n. For details of this experiment, see Supporting Information. We
observe (see Figure 4 and additional results in the Supporting Information) 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: [11] observed universality in sparse
vector recovery using l1 norm minimization; [12] observed in passing universality for matrix
recovery with NNM; [13] 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) [14, 13]. Background and empirical
evidence are presented in the Supporting Information.

14

Stopping conditions and sign of convergence. State Evolution suggests that the
estimated noise level σ̂t in every iteration is in fact a good proxy for the current MSE:
σ̂t2 ≈ δN1M ||Xt − X||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 Xt , and can choose
to stop or continue accordingly. Background and empirical evidence are presented in the
Supporting Information.
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. Furthermore, evidence included in
the Supporting Information suggests that the below its phase transition, AMP-OPT fails
always, while above its phase transition it succeeds with probability that depends on N .
We believe this happens because the shrinker η opt is not Lipschitz continuous (see Figure 1).
Indeed, current theoretical results on AMP for vector recovery problems require the denoiser
to be Lipschitz continuous [2, 33, 34].
Optimality. We conjecture that the information-theoretical lower bound n/(N M ) ≥ Mr (1+
M/N − r/M ) and the denoising lower bound (11) are both asymptotically tight. Also, one
naturally wonders whether the denoiser η opt in (15) is in fact optimal in the (asymptotic)
minimax sense. Numerical evidence suggests that this is not the case, at least when ρ is
sufficiently large - see Supporting Information. Nonetheless, the worst-case asymptotic MSE
of η opt seems to be very close to the minimax MSE.

5

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 η opt has lead to AMPOPT, a near-optimal matrix recovery algorithm. For example, when r = rank(X) is low,
AMP-OPT requires n ≈ r(N +M ) measurements for recovery, whereas NNM requires n ≈
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.

15

Acknowledgement
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.

References
[1] D. Amelunxen, M. Lotz, M. B. McCoy, and J. A. Tropp. Living on the edge: Phase
transitions in convex programs with random data. Information and Inference, page
iau005, 2014.
[2] M. Bayati and A. Montanari. The dynamics of message passing on dense graphs, with
applications to compressed sensing. IEEE Trans. Info. Theory, 57(2):764–785, 2011.
[3] J.-F. Cai, E. J. Candès, and Z. Shen. A singular value thresholding algorithm for matrix
completion. SIAM Journal on Optimization, 20(4):1956–1982, 2010.
[4] E. Candes and B. Recht. Exact matrix completion via convex optimization. Communications of the ACM, 55(6):111–119, 2012.
[5] E. J. Candès, X. Li, Y. Ma, and J. Wright. Robust principal component analysis?
Journal of the ACM (JACM), 58(3):11, 2011.
[6] E. J. Candes and Y. Plan. Matrix completion with noise. Proc. IEEE, 98(6):925–936,
2010.
[7] E. J. Candes, C. A. Sing-Long, and J. D. Trzasko. Unbiased risk estimates for singular
value thresholding and spectral estimators. IEEE Trans. Sig. Proc., 61(19):4643–4657,
2013.
[8] E. J. Candès and T. Tao. The power of convex relaxation: Near-optimal matrix completion. IEEE Transactions on Information Theory, 56(5):2053–2080, 2010.
[9] M. A. Davenport and J. Romberg. An overview of low-rank matrix recovery from incomplete observations. IEEE Journal of Selected Topics in Signal Processing, 10(4):608–622,
2016.
[10] D. Donoho and M. Gavish. Minimax risk of matrix denoising by singular value thresholding. The Annals of Statistics, 42(6):2413–2440, 2014.
[11] D. Donoho and J. Tanner. Observed universality of phase transitions in high-dimensional
geometry, with implications for modern data analysis and signal processing. Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and
Engineering Sciences, 367(1906):4273–4293, 2009.
16

[12] D. L. Donoho, M. Gavish, and A. Montanari. The phase transition of matrix recovery
from gaussian measurements matches the minimax mse of matrix denoising. Proceedings
of the National Academy of Sciences, 110(21):8405–8410, 2013.
[13] D. L. Donoho, I. Johnstone, and A. Montanari. Accurate prediction of phase transitions
in compressed sensing. IEEE Trans. Info. Theory, 59(6):3396–3433, 2013.
[14] D. L. Donoho, A. Maleki, and A. Montanari. Message-passing algorithms for compressed
sensing. Proceedings of the National Academy of Sciences, 106(45):18914–18919, 2009.
[15] D. L. Donoho, A. Maleki, and A. Montanari. The noise-sensitivity phase transition in
compressed sensing. IEEE Trans. on Info. Theory, 57(10):6920–6941, 2011.
[16] D. L. Donoho and J. Tanner. Sparse nonnegative solution of underdetermined linear
equations by linear programming. Proceedings of the National Academy of Sciences of
the United States of America, 102(27):9446–9451, 2005.
[17] M. Gavish and D. L. Donoho. Optimal shrinkage of singular values. IEEE Trans. Info.
Theory, 63(4):2137–2152, 2017.
[18] D. Gross. Recovering low-rank matrices from few coefficients in any basis. IEEE Transactions on Information Theory, 57(3):1548–1566, 2011.
[19] D. Gross, Y.-K. Liu, S. T. Flammia, S. Becker, and J. Eisert. Quantum state tomography
via compressed sensing. Physical review letters, 105(15):150401, 2010.
[20] P. Jain, R. Meka, and I. S. Dhillon. Guaranteed rank minimization via singular value
projection. In Advances in Neural Information Processing Systems, pages 937–945, 2010.
[21] P. Jain, P. Netrapalli, and S. Sanghavi. Low-rank matrix completion using alternating
minimization. In Proc. 45th annual ACM symposium Theory computing, pages 665–674.
ACM, 2013.
[22] R. H. Keshavan, A. Montanari, and S. Oh. Matrix completion from noisy entries.
Journal of Machine Learning Research, 11(Jul):2057–2078, 2010.
[23] Z. Liu and L. Vandenberghe. Interior-point method for nuclear norm approximation
with application to system identification. SIAM Mat. Anal. and App., 31(3):1235–1256,
2009.
[24] S. Ma, D. Goldfarb, and L. Chen. Fixed point and bregman iterative methods for matrix
rank minimization. Mathematical Programming, 128(1):321–353, 2011.
[25] A. Maleki and D. L. Donoho. Optimally tuned iterative reconstruction algorithms for
compressed sensing. IEEE Journal of Selected Topics in Signal Processing, 4(2):330–341,
2010.

17

[26] C. Metzler, A. Maleki, and R. Baraniuk. From denoising to compressed sensing. IEEE
Trans. Info. Theory, 62(9):5117–5144, 2016.
[27] H. Monajemi, S. Jafarpour, M. Gavish, D. L. Donoho, S. Ambikasaran, S. Bacallado,
D. Bharadia, Y. Chen, Y. Choi, M. Chowdhury, et al. Deterministic matrices matching
the compressed sensing phase transitions of gaussian random matrices. Proc. of the
Natl. Acad. Sci., 110(4):1181–1186, 2013.
[28] A. Montanari. Graphical models concepts in compressed sensing. Compressed Sensing:
Theory and Applications, pages 394–438, 2012.
[29] S. Oymak and B. Hassibi. Sharp mse bounds for proximal denoising. Foundations of
Computational Mathematics, 16(4):965–1029, 2016.
[30] B. Recht, M. Fazel, and P. A. Parrilo. Guaranteed minimum-rank solutions of linear
matrix equations via nuclear norm minimization. SIAM review, 52(3):471–501, 2010.
[31] E. Romanov. Code and data supplement for “near-optimal matrix recovery from random
linear measurements”, https://purl.stanford.edu/rt605yk2478, 2017.
[32] E. Romanov and M. Gavish. Supporting information for “near-optimal matrix recovery
from random linear measurements”, https://purl.stanford.edu/rt605yk2478, 2017.
[33] C. Rush, A. Greig, and R. Venkataramanan. Capacity-achieving sparse regression codes
via approximate message passing decoding. In Information Theory (ISIT), 2015 IEEE
International Symposium on, pages 2016–2020. IEEE, 2015.
[34] C. Rush and R. Venkataramanan. Finite-sample analysis of approximate message passing. In Information Theory (ISIT), 2016 IEEE International Symposium on, pages
755–759. IEEE, 2016.
[35] A. A. Shabalin and A. B. Nobel. Reconstruction of a low-rank matrix in the presence
of gaussian noise. Journal of Multivariate Analysis, 118:67–76, 2013.
[36] J. Tanner and K. Wei. Normalized iterative hard thresholding for matrix completion.
SIAM Journal on Scientific Computing, 35(5):S104–S125, 2013.
[37] K.-C. Toh and S. Yun. An accelerated proximal gradient algorithm for nuclear norm
regularized linear least squares problems. Pacific Journal of optimization, 6(615-640):15,
2010.
[38] B. Vandereycken. Low-rank matrix completion by riemannian optimization. SIAM
Journal on Optimization, 23(2):1214–1236, 2013.
[39] M. Zaharia, M. Chowdhury, M. J. Franklin, S. Shenker, and I. Stoica. Spark: Cluster
computing with working sets. HotCloud, 10(10-10):95, 2010.

18

[40] X. Zhou, C. Yang, H. Zhao, and W. Yu. Low-rank modeling and its applications in
image analysis. ACM Computing Surveys (CSUR), 47(2):36, 2015.

19

