# The Phase Transition of Matrix Recovery from Gaussian Measurements Matches
the Minimax MSE of Matrix Denoising

David L. Donoho 111Department of Statistics, Stanford University Matan Gavish
111Department of Statistics, Stanford University Andrea Montanari
111Department of Statistics, Stanford University 222Department of Electrical
Engineering, Stanford University

###### Abstract

Let $X_{0}$ be an unknown $M$ by $N$ matrix. In matrix recovery, one takes
$n<MN$ linear measurements $y_{1},\dots,y_{n}$ of $X_{0}$, where $y_{i}={\rm
Tr}(a_{i}^{T}X_{0})$ and each $a_{i}$ is a $M$ by $N$ matrix. For measurement
matrices with Gaussian i.i.d entries, it known that if $X_{0}$ is of low rank,
it is recoverable from just a few measurements. A popular approach for matrix
recovery is Nuclear Norm Minimization (NNM): solving the convex optimization
problem $\text{min }\|X\|_{*}\text{ subject to }y_{i}={\rm Tr}(a_{i}^{T}X)$
for all $1\leq i\leq n$, where $\|\cdot\|_{*}$ denotes the nuclear norm,
namely, the sum of singular values. Empirical work reveals a _phase
transition_ curve, stated in terms of the undersampling fraction
$\delta(n,M,N)=n/(MN)$, rank fraction $\rho=r/N$ and aspect ratio $\beta=M/N$.
Specifically, a curve $\delta^{*}=\delta^{*}(\rho;\beta)$ exists such that, if
$\delta>\delta^{*}(\rho;\beta)$, NNM typically succeeds, while if
$\delta<\delta^{*}(\rho;\beta)$, it typically fails.

An apparently quite different problem is matrix denoising in Gaussian noise,
where an unknown $M$ by $N$ matrix $X_{0}$ is to be estimated based on direct
noisy measurements $Y=X_{0}+Z$, where the matrix $Z$ has iid Gaussian entries.
It has been empirically observed that, if $X_{0}$ has low rank, it may be
recovered quite accurately from the noisy measurement $Y$. A popular matrix
denoising scheme solves the unconstrained optimization problem $\text{min
}\|Y-X\|_{F}^{2}/2+\lambda\|X\|_{*}$. When optimally tuned, this scheme
achieves the asymptotic minimax MSE ${\cal
M}(\rho)=\lim_{N\rightarrow\infty}\inf_{\lambda}\sup_{{\rm
rank}(X)\leq\rho\cdot N}MSE(X,\hat{X}_{\lambda})$.

We report extensive experiments showing that the phase transition
$\delta^{*}(\rho)$ in the first problem (Matrix Recovery from Gaussian
Measurements) coincides with the minimax risk curve ${\cal M}(\rho)$ in the
second problem (Matrix Denoising in Gaussian Noise): $\delta^{*}(\rho)={\cal
M}(\rho)$, for any rank fraction $0<\rho<1$.

Our experiments considered matrices belonging to two constraint classes: real
$M$ by $N$ matrices, of various ranks and aspect ratios, and real symmetric
positive semidefinite $N$ by $N$ matrices, of various ranks. Different
predictions ${\cal M}(\rho)$ of the phase transition location were used in the
two different cases, and were validated by the experimental data.

###### Contents

  1. 1 Introduction
  2. 2 Methods
    1. 2.1 Generation of Problem Instances
    2. 2.2 Convex Optimization
    3. 2.3 Probability of Exact Recovery
    4. 2.4 Estimating the Probability of Exact Recovery
    5. 2.5 Asymptotic Phase Transition Hypothesis
    6. 2.6 Empirical Phase Transitions
    7. 2.7 Experimental Design
    8. 2.8 Computing
  3. 3 Results
  4. 4 Discussion: Existing Literature and Our Contributions
  5. 5 Conclusions
  6. A Asymptotic Minimax MSE Formula
  7. B Summary of Empirical Results
  8. C Data Deposition
  9. D Code Deposition

##  1 Introduction

Let $X_{0}$ be an unknown $M$ by $N$ matrix. How many measurements must we
obtain in order to ‘completely know’ $X_{0}$? While it seems that $MN$
measurements must be necessary, in recent years intense research in applied
mathematics, optimization and information theory, has shown that, when $X_{0}$
is of low rank, we may efficiently recover it from a relatively small number
of linear measurements by convex optimization [3, 1, 2]. Applications have
been developed in fields ranging widely, for example from video and image
processing [4], to quantum state tomography [5], to collaborative filtering
[1, 6].

Specifically, let ${\cal A}:{\mathbb{R}}^{M\times N}\to{\mathbb{R}}^{n}$ be a
linear operator and consider measurements $y={\cal A}(X_{0})$. If $n<MN$, the
problem of inferring $X_{0}$ from $y$ may be viewed as attempting to solve an
underdetermined system of equations. Under certain circumstances, it has been
observed that this (seemingly hopeless) task can be accomplished by solving
the so-called nuclear norm minimization problem

$(P_{nuc})\qquad\text{min }\|X\|_{*}\;\;\text{ subject to }y={\cal A}(X)\,.$
(1)

Here the _nuclear norm_ $\|X\|_{*}$ is the sum of singular values of $X$. For
example, it was found that if $X_{0}$ is sufficiently low rank, with a
principal subspace in a certain sense incoherent to the measurement operator
${\cal A}$, then the solution $X_{1}=X_{1}(y)$ to $(P_{nuc})$ is precisely
$X_{0}$. Such incoherence can be obtained by letting ${\cal A}$ be random, for
instance if ${\cal A}(X_{0})_{i}={\rm Tr}(a_{i}^{T}X_{0})$ with
$a_{i}\in{\mathbb{R}}^{m\times n}$ having i.i.d. Gaussian entries. In this
case we speak of _“matrix recovery from Gaussian measurements”_ [3].

A key phrase from the previous paragraph: ‘if $X_{0}$ is _sufficiently_ low
rank’. Clearly there must be a quantitative trade-off between the rank of
$X_{0}$ and the number of measurements required, such that higher rank
matrices require more measurements. In the Gaussian measurements model, with
$N$ sufficiently large, empirical work by Recht, Xu and Hassibi [8, 7], Fazel,
Parillo and Recht [3], Tanner and Wei [9] and Oymak and Hassibi [10],
documents a _phase transition_ phenomenon. For matrices of a given rank, there
is a fairly precise number of required samples, in the sense that a transition
from non recovery to complete recovery takes place sharply as the number of
samples varies across this value. For example, in Figure 1 below we report
results obtained in our own experiments, showing that, for reconstructing
matrices of size 60 by 60 which are of rank 20, 2600 Gaussian measurements are
sufficient with very high probability, but 2400 Gaussian measurements are
insufficient with very high probability.

$\delta$ |  0.67 |  0.68 |  0.68 |  0.68 |  0.68 |  0.69 |  0.69 |  0.69 |  0.69 |  0.70 |  0.71 |  0.72 |  0.73
---|---|---|---|---|---|---|---|---|---|---|---|---|---
$n$ |  2400 |  2437 |  2446 |  2455 |  2465 |  2474 |  2483 |  2492 |  2502 |  2511 |  2538 |  2575 |  2612
$\hat{\pi}(r|n,N)$ |  0.00 |  0.00 |  0.30 |  0.20 |  0.30 |  0.45 |  0.60 |  0.80 |  0.60 |  0.90 |  1.00 |  1.00 |  1.00
Figure 1: Data from typical Phase Transition experiment. Here $N=60$, $r=20$,
and the number $n$ of Gaussian measurements varies. Note: our formula predicts
an asymptotic phase transition at $\delta^{*}=0.6937$, corresponding to
$n=2497$. And, indeed, the success probability is close to $1/2$ at that $n$.
All runs involved $T=20$ Monte Carlo trials.

In this paper, we present a simple and explicit formula for the phase
transition curve in matrix recovery from Gaussian measurements. The formula
arises in an apparently unrelated problem: matrix de-noising in Gaussian
noise. In this problem, we again let $X_{0}$ denote an $M$ by $N$ matrix, and
we observe $Y=X_{0}+Z$, where $Z$ is Gaussian iid noise $Z_{ij}\sim{\cal
N}(0,1)$. We consider the following nuclear norm de-noising scheme:

$\displaystyle(P_{nuc,\lambda})\qquad\min\Big{\\{}\frac{1}{2}\|Y-X\|_{F}^{2}+\lambda\|X\|_{*}\Big{\\}}\,.$
(2)

In this problem the measurements $Y$ are direct, so in some sense complete,
but noisy. The solution $\hat{X}_{\lambda}(Y)$ can be viewed as a shrinkage
estimator. In the basis of the singular vectors $U_{Y}$ and $V_{Y}$ of $Y$,
the solution $\hat{X}_{\lambda}(Y)$ is diagonal, and the diagonal entries are
produced by soft thresholding of the singular values of $Y$.

Because the measurements $y$ in the matrix recovery problem are noiseless but
incomplete, while the measurements $Y$ in the matrix denoising problem are
complete but noisy, the problems seem quite different. Nevertheless, we show
here that there is a deep connection between the two problems.

Let us quantify performance in the denoising problem by the minimax MSE,
namely

$\mathcal{M}(\rho;M,N)=\min_{\lambda}\max_{rank(X)\leq\rho
N}MSE(X_{0}\,,\,\hat{X}_{\lambda}(Y)),$

where MSE refers to the dimension-normalized mean-squared error

$\frac{1}{MN}E\|X-\hat{X}_{\lambda}\|_{F}^{2}$

and subscript $F$ denotes Frobenius norm. The asymptotic minimax MSE
$\mathcal{M}(\rho;\beta)=\lim_{N\rightarrow\infty}\mathcal{M}(\rho;\beta N,N)$
has been derived in [11]. Explicit formulas for the curve
$\rho\mapsto\mathcal{M}(\rho;\beta)$ appear in the Appendix. A parametric form
is given for the case of asymptotically square matrices, $\beta=1$. Figures 2
and 3 depict the various minimax MSE curves.

Figure 2:  The two asymptotic minimax MSE’s $M(\rho|{\bf X})$: ${\bf X}=M_{N}$
(black), ${\bf X}=Symm_{N}$ (red), in the case of square matrices. For non
square matrices, see curves in Figure 3 Figure 3:  Curves: Asymptotic Minimax
MSE’s for nonsquare matrix settings $Mat_{M,N}$; varying shape factor
$\beta=M/N$. Points: Empirical phase transition locations for matrix recovery
from incomplete measurements, see Table 8.

We can now state our main hypothesis for matrix recovery from Gaussian
measurements.

Main Hypothesis: Asymptotic Phase Transition Formula. _Consider a sequence of
matrix recovery problems with parameters $\\{(r,n,M,N)\\}_{N\geq 1}$ having
limiting fractional rank $\rho=\lim_{N\to\infty}r/N$, limiting aspect ratio
$\beta=\lim_{N\to\infty}M/N$, and limiting incompleteness fraction
$\delta=\lim_{N\to\infty}n/(MN)$. In the limit of large problem size $N$, the
solution $X_{1}(y)$ to the nuclear norm minimization problem $(P_{nuc})$ is
correct with probability converging to one if $\delta>{\cal M}(\rho;\beta)$
and incorrect with probability converging to one if $\delta<{\cal
M}(\rho;\beta)$. _ In short: _The asymptotic phase transition
$\delta^{*}(\rho,\beta)$ in Gaussian matrix recovery is equal to the
asymptotic minimax MSE ${\cal M}(\rho;\beta)$._

In particular, for the case of small rank $r$, by studying the small $\rho$
asymptotics of Eq. (14), we obtain that reconstruction is possible from $n\geq
2r(M+N+\sqrt{MN})(1+O(r/N))$ measurements, but not from substantially less.

This brief announcement tests this hypothesis by conducting a substantial
computer experiment generating large numbers of random problem instances. We
use statistical methods to check for disagreement between the hypothesis and
the predicted phase transition. To bolster the solidity of our results, we
conduct the experiment in two different settings: $(1)$ the matrix $X_{0}$ is
a general $M$ by $N$ matrix, for various rank fractions $\rho$ and aspect
ratios $\beta$; $(2)$ $X_{0}$ is a symmetric positive definite matrix, for
various rank fractions $\rho$. In the latter case the positive semidefinite
constraint is added to the convex program $(P_{nuc})$. As described below,
there are different asymptotic MSE curves for the two settings. We demonstrate
an empirically accurate match in each of the cases, showing the depth and
significance of the connection we expose here.

In the discussion and conclusion we connect our result with related work in
the field of sparsity-promoting reconstructions, where the same formal
identity between a minimax MSE and a phase transition boundary has been
observed, and in some cases even proved. We also discuss recent rigorous
evidence towards establishing the above matrix recovery phase transition
formula.

##  2 Methods

We investigated the hypothesis that the asymptotic phase transition boundary
agrees with the proposed phase transition formula to within experimental
error.

For notational simplicity we will focus here on the case $M=N$, and defer the
case of non-square matrices to the SI. Hence, we will drop throughout the main
text the argument $\beta=1$. The asymptotic phase plane at point
$(\rho,\delta)$ is associated to triples $(r,n,N)$, where $\rho=r/N\in[0,1]$
is the rank fraction, and $\delta=\delta(n,N|{\bf X})=n/dim({\bf X})$ is the
under sampling ratio, where $dim({\bf X})$ is the dimension of the underlying
collection of matrices ${\bf X}$. We performed a sequence of experiments, one
for each tuple, in which we generated random rank-$r$ $N$ by $N$ matrices
$X_{0}\in{\bf X}$, random measurement matrices ${\cal A}=A$ of size $n\times
N^{2}$, and obtained random problem instances $(y,A)$. We then applied a
convex optimization procedure, obtaining a putative reconstruction
$\hat{X}=\hat{X}(y,A)$. We declared a reconstruction successful when the
Frobenius norm was smaller than a threshold. Our raw empirical observations
consist of a count of empirical successes and sample sizes at a selected set
of positions $(\rho,\delta)$ and a selected set of problem sizes $N$. From
these raw counts we produce fitted success probabilities $\hat{\pi}(r|n,N,{\bf
X})$, The finite-$N$ phase transition is the place where the true underlying
probability of successful reconstruction take the value $50$%. We tested the
hypothesis that the finite-$N$ transition was consistent with the proposed
asymptotic phase transition formula.

This section discusses details of data generation and data analysis.

###  2.1 Generation of Problem Instances

Each problem instance $(y,A)$ was generated by, first, generating a random
rank $r$ matrix $X_{0}$, then, generating a random measurement matrix
$A=A_{n,N^{2}}$ and then applying $y=A\cdot vec(X_{0})$.

We considered problem instances of two specific types, corresponding to
matrices $X_{0}\in{\bf X}$, with ${\bf X}$ one of two classes of matrices

  * •

${Mat}_{N}$: all $N\times N$ matrices with real-valued entries

  * •

${Sym}_{N}$: all $N\times N$ real symmetric matrices which are nonnegative-
semidefinite

In the case ${\bf X}={Mat}_{N}$, we consider low-rank matrices
$X_{0}=UV^{\prime}$ where $U$ and $V$ are each $N$ by $r$ partial orthogonal
matrices in the Stiefel manifold $St(N,r)$. The matrices are uniformly
distributed on $St(N,r)$. In the case ${\bf X}={Sym}_{N}$, we consider low-
rank matrices $X_{0}=UU^{\prime}$ where $U$ is an $N$ by $r$ partial
orthogonal matrix in $St(N,r)$, and again the matrix is uniformly distributed.

For measurement matrices $A$, we use Gaussian random matrices satisfying
$A_{i,j}\sim N(0,1/n)$.

###  2.2 Convex Optimization

For a given problem instance $(y,A)$, we attempt to recover the underlying
low-rank object $X_{0}$ from the measurements $y$ by convex optimization. Each
of our choices ${\bf X}$ gives rise to an associated optimization problem:

$(P^{\bf X}_{nuc})\qquad\text{min }||X||_{*}\text{ subject to }y=A\cdot
vec(X),\qquad X\in{\bf X}.$

Here ${\bf X}$ is one of these two classes of matrices ${Mat}_{N}$ or
${Sym}_{N}$. The two resulting optimization problems can each be reduced to a
so-called semidefinite programming problem; see [13, 12].

###  2.3 Probability of Exact Recovery

Since both the measurement matrix $A$, and the underlying low-rank object
$X_{0}$ are random, $(y,A)$ is a random instance for $(P^{\bf X}_{nuc})$. The
probability of exact recovery is defined by

$\pi(r|n,N,{\bf X})=\text{Prob}\\{X_{0}\mbox{ is the unique solution of
}(P^{\bf X}_{nuc})\\}.$

Clearly $0\leq\pi\leq 1$; for fixed $N$, $\pi$ is monotone decreasing in $r$
and monotone increasing in $n$. Also
$\pi(r|n,N,{Mat}_{N})<\pi(r|n,N,{Sym}_{N})$.

###  2.4 Estimating the Probability of Exact Recovery

Our procedure follows [14, 15]. For a given matrix type ${\bf X}$ and rank $r$
we conduct an experiment whose purpose is to estimate $\pi(r|n,N,{\bf X})$
using $T$ Monte Carlo trials. In each trial we generate a random instance
$(y,A)$ which we supply to a solver for $(P^{\bf X}_{(nuc)})$, obtaining the
result $X_{1}$. We compare the result $X_{1}$ to $X_{0}$. If the relative
error $\|X_{0}-X_{1}\|_{F}/\|X_{0}\|_{F}$ is smaller than a numerical
tolerance, we declare the recovery a success; if not, we declare it a failure.
(In this paper, we used an error tolerance of $0.001$.) We thus obtain $T$
binary measurements $Y_{i}$ indicating success or failure in reconstruction.
The empirical success fraction is then calculated as

$\hat{\pi}(r|n,N,T,{\bf
X})=\frac{\\#\\{\text{successes}\\}}{\\#\\{\text{trials}\\}}=\frac{1}{T}\sum_{i=1}^{T}Y_{i}\,.$

These are the raw observations generated by our experiments.

###  2.5 Asymptotic Phase Transition Hypothesis

Consider a sequence of tuples $(r,n,N)$ with $r/N\rightarrow\rho$ and
$n/N\rightarrow\delta$. We assume that there is an asymptotic phase transition
curve $\delta^{*}(\rho|{\bf X})$, i.e. a curve obeying

$\pi(r|{n,N},{\bf
X})\rightarrow\left\\{\begin{array}[]{ll}1&\delta<\delta^{*}(\rho|{\bf X})\\\
0&\delta>\delta^{*}(\rho|{\bf X})\end{array}\right.$ (3)

For many convex optimization problems the existence of such an asymptotic
phase transition is rigorously proven; see the Discussion below.

The hypothesis we investigate concerns the value of $\delta^{*}(\rho|{\bf
X})$; specifically, whether

$\delta^{*}(\rho|{\bf X})={\cal M}(\rho|{\bf X}).$ (4)

Here ${\cal M}(\rho|{Mat})$ (respectively ${\cal M}(\rho|{Sym})$ ) is the
minimax MSE for SVT for general matrices (respectively, positive definite
ones). Formulas for ${\cal M}$ were derived by the Authors in [11];
computational details are provided in the Appendix.

###  2.6 Empirical Phase Transitions

The empirical phase transition point is estimated by fitting a smooth function
$\hat{\pi}(n/N)$ (in fact a logistic function) to the empirical data
$\hat{\pi}(r|n,N,{\bf X})$ using the glm() command in the R statistics
language. In fact we fit the logistic model that
$\mbox{logit}(\pi)\equiv\log(\frac{\pi}{1-\pi})=a+b\Delta$, where
$\Delta(\delta|\rho)=\delta-M(\rho)$ is the offset between $\delta$ and the
predicted phase transition. The coefficients $a$ and $b$ are called the
intercept and slope, and will be tabulated below. The intercept gives the
predicted logit exactly at $\Delta=0$, i.e. $\delta={\cal M}(\rho)$. The
empirical phase transition is located at $\hat{\delta}(r,N,M,{\bf X})={\cal
M}(\rho)-a/b$. This is the value of $\delta=\delta(n,N|{\bf X})$ solving

$\hat{\pi}(\delta)=1/2.$

Under the hypothesis (4) we have

$\lim_{N\rightarrow\infty,r/N\rightarrow\rho}\lim_{T\rightarrow\infty}\hat{\delta}(r,N,T,{\bf
X})={\cal M}(\rho|{\bf X}).$

Consequently, in data analysis we will compare the fitted values
$\hat{\delta}(r,N,T,{\bf X})$ with ${\cal M}(r/N|{\bf X})$.

###  2.7 Experimental Design

To address our main hypothesis regarding the agreement of phase transition
boundaries, we measure $\hat{\pi}$ at points $\delta=n/N$ and $\rho=r/N$ in
the phase plane $(\delta,\rho)$ which we expect to be maximally informative
about the location of the phase transition. In fact the informative locations
in binomial response models correspond to points where the probability of
response is in the middle range $(1/10,9/10)$ [16]. As a rough approximation
to such an optimal design, we sample at equispaced $\delta\in[{\cal
M}(\rho|{\bf X})-0.05,{\cal M}(\rho|{\bf X})+0.05]$.

###  2.8 Computing

We primarily used the MATLAB computing environment, and the popular CVX convex
optimization package [17]. A modeling system for disciplined convex
programming by Boyd, Grant and others, supporting two open source interior-
point solvers: SeDuMi and SDPT3 [19, 18].

We also studied the robustness of our results across solvers. Zulfikar Ahmed
translated our code into Python and used the general purpose solver package
CVXOPT by Anderson and Vandeberghe [20].

##  3 Results

Our experimental data have been deposited at [21] they are contained in a text
file with more than 100,000 lines, each line reporting one batch of Monte
Carlo experiments at a given $r,n,M,N$ and ${\bf X}$. Each line also documents
the number of Monte Carlo trials $T$, and the observed success fraction
$\hat{\pi}$. The file also contains metadata identifying the solver and the
researcher responsible for the run.

In all cases, we observed a transition from no observed successes at
$\delta={\cal M}(\rho)-0.05$ to no observed failures at $\delta={\cal
M}(\rho)+0.05$. Figure 3 shows results we obtained at non square matrix
ensembles, with varying $\beta=M/N$. The minimax MSE curves $M(\rho|\beta)$
vary widely, but the observed PT’s track the curves closely.

$\rho$ | $N$ | $T$ | ${\cal M}(\rho)$ | $\hat{\delta}(\rho)$ | $a$ | $b$ | $Z$
---|---|---|---|---|---|---|---
1/10 | 40 | 400 | 0.351 | 0.352 | -0.128 | 182.978 | -0.581
50 | 400 | 0.351 | 0.350 | 0.282 | 200.131 | 1.213
60 | 400 | 0.351 | 0.352 | -0.096 | 221.212 | -0.398
80 | 400 | 0.351 | 0.350 | 0.415 | 295.049 | 1.452
100 | 400 | 0.351 | 0.349 | 0.641 | 383.493 | 1.900
Figure 4: Results with ${\bf X}={Mat}_{N}$, with $\rho=1/10$ and Varying $N$.
$T$ total number of reconstructions; $a$, $b$: fitted logistic regression
parameters; $Z$: traditional Z-score of logistic regression intercept. See
Table 6 for the complete table.

Figure 4 shows a small subset of our results in the square case ${\bf
X}={Mat}_{N}$, to explain our empirical results; the full tables are given in
SI for the square, non square and symmetric positive definite cases. In the
square case, the empirical phase transition agrees in all cases with the
formula ${\cal M}(\rho)$ to two digits accuracy. Table 7 shows that, in the
symmetric nonnegative-definite case ${\bf X}={Sym}_{N}$, the empirical phase
transition falls within $[{\cal M}(\rho)-0.01,{\cal M}(\rho)+0.01]$.

Previous empirical studies of phase transition behavior in sparse recovery
show that, even in cases where the asymptotic phase transition curve is
rigorously proven and analytically available, such large-$N$ theory cannot be
expected to match empirical finite-$N$ data to within the usual naive standard
errors [14, 15]. Instead, one observes a finite transition zone of width
$\approx c_{1}/N^{1/2}$ and a small displacement of the empirical phase
transition away from the asymptotic Gaussian phase transition, of size
$\approx c_{2}/N$. Hence, the strict literal device of testing the hypothesis
that $E\hat{\delta}={\cal M}(\rho)$ is not appropriate in this setting 111As
shown in Figure 4, and in Tables 6,7,8,9,and 10, our results in many cases do
generate $Z$ scores for the intercept term in the logistic regression which
are consistent with traditional acceptance of this hypothesis. However,
traditional acceptance is not needed, in order for the main hypothesis to be
valid, and because of finite-$N$ scaling effects indicated above, would not
ordinarily be expected to hold.

A precise statement of our hypothesis uses the fitted logistic parameters
$\hat{a}=\hat{a}(r,N,M,{\bf X})$ and $\hat{b}=\hat{b}(r,N,T,{\bf X})$, defined
above. The asymptotic relation222 $=_{P}$ denotes convergence in probability.

$\lim_{N\rightarrow\infty}\lim_{T\rightarrow\infty}\frac{\hat{a}(r,N,T,{\bf
X})}{\hat{b}(r,N,T,{\bf X})}=_{P}0,\qquad r=\lfloor\rho N\rfloor,$ (5)

implies $\delta(\rho|{\bf X})={\cal M}(\rho|{\bf X})$ . Now note that in
Figure 4 the coefficient $b$ scales directly with $N$ and takes the value
several hundred for large $N$. This means that, in these experiments, the
transition from complete failure to complete success happens in a zone of
width $<1/100$. Notice also that $a$ stays bounded, for the most part between
$0$ and $2$. This means that the response probability evaluated exactly at
${\cal M}(\rho)$ obeys typically:

$\mbox{logit}(\hat{\pi}({\cal M}(\rho))\in[0,2],$

Figure 5:  Finite-$N$ scaling effects. Fitted slopes $b$ from Tables 6 and 7
versus matrix sizes $N$. Left Panel: asymmetric square matrices ${\bf
X}={Mat}_{N}$. Right Panel: symmetric nonnegative-definite matrices ${\bf
X}={Sym}_{N}$. Superimposed lines have formulas $b=10+3N$. Note: in Left
Panel, data for $\rho\in\\{3/4,9/10\\}$ were excluded because the slopes $b$
were very large (700 and 2700, respectively).

Figure 5, Panel A presents the fitted slopes $b$ and problem sizes $N$,as well
as an empirical fit. All the data came from Table 6, but we omitted results
for $\rho\in\\{3/4,9/10\\}$ because those slopes were large multiples of all
other slopes. Similarly Figure 5, Panel B presents the slopes from Table 7. In
each plot, the line $b=6+3N$ is overlaid for comparison. Both plots show a
clear tendency of the slopes to increase with $N$.

The combination of linear growth of $b$ with $N$ and non-growth of $a$
333Actually, even sub linear growth of $a$ implies the result. implies Eq.
(5), and acceptance of our main hypothesis.

#### Nongaussian Measurements.

This paper studies matrix recovery from Gaussian measurements of the unknown
matrix, and specifically does not study recovery from partial entry wise
measurements typically called ’matrix completion’. Entrywise measurements is
yield a phase transition at a different location [9]. Our conclusions do
extend to certain nonGaussian measurements based on ${\cal A}$ with
independent and identically distributed entries that are equiprobable $\pm 1$
(a.k.a. Rademacher matrices). See Table 10.

##  4 Discussion: Existing Literature and Our Contributions

Phase transitions in the success of convex optimization at solving non-convex
problems have been observed previously. Donoho and Tanner considered linear
programs useful in the recovery of sparse vectors, rigorously established the
existence of asymptotic phase transitions, and derived exact formulas for the
asymptotic phase transition [25, 22, 14, 23, 24] as well as finite $N$
exponential bounds. Their work considered the recovery of $k$-sparse vectors
with $N$ entries from $n$ Gaussian measurements. The phase diagram in those
papers could be stated in terms of variables $(\kappa,\delta)$, where, as in
this paper, $\delta=n/N$ is the under sampling fraction, and $\kappa=k/N$ is
the sparsity fraction, which plays a role analogous to the role played here by
the rank fraction $\rho$. The proofs in those papers were obtained by
techniques from high-dimensional combinatorial geometry, and the formulas for
the asymptotic phase transition were implicit, written in terms of special
functions.

Donoho, Maleki, and Montanari [26] later developed a new so-called Approximate
Message Passing (AMP) approach to the sparse recovery problem, which gave new
formulas for phase boundaries, confirmed rigorously by Bayati and Montanari in
[27, 28]. While the previous formulas involved combinatorial geometry, the new
(necessarily equivalent) ones involved instead minimax decision theory. An
extensive literature on AMP algorithms has since developed see, e.g. [30, 29],
implying, among other things, universality in the sparse recovery phase
transition [31].

Donoho, Johnstone and Montanari [32] generalized previous AMP-based phase
transition results to block-sparse, monotone, and piecewise constant signals.
They provided evidence that the observed phase transition
$\delta^{*}(\,\cdot\,)$ of the associated convex optimization reconstruction
methods obeyed formulas of the form

$\delta^{*}(\kappa)={\cal M}(\kappa),\;\;0<\kappa<1.$ (6)

Here $\kappa$ is a variable measuring generalized sparsity and ${\cal
M}(\kappa)$ the minimax MSE of an appropriate denoiser based on direct noisy
measurements. The main result in this brief report fits in this general
framework whereby the sparsity $\kappa$ is identified with the fractional rank
$\rho$, and the minimax MSE symbol ${\cal M}$ applies to the singular value
thresholding denoiser. Our main result is therefore an extension of DJM-style
formulas from the sparse recovery setting problem to the low rank recovery
setting.

Much mathematical study of matrix recovery [1, 8, 3, 33] has focused on
providing rigorous bounds which show the existence of a region of success,
without however establishing a phase transition phenomenon, or determining its
exact boundary. A relatively accurate result by Candès and Recht (CR) for the
case of Gaussian measurements [34], implies that $n\geq 3r(M+N)[1+O(r/N)]$
measurements are sufficient. Our formulas for the square case $M=N$ show that

${\cal M}(\rho|Mat_{N})\sim 6\rho,\qquad\rho\rightarrow 0.$

This agrees with the CR bound in the very low-rank case. However, our relation
$\delta^{*}(\rho)={\cal M}(\rho)$ is apparently noticeably more accurate than
the CR formula at finite $N$. Table 9 presents experiments where the rank is
fixed at $r=1$,$2$,$3$, or $4$, and $N$ varies between $40$ and $90$. Even
though in such cases the corresponding $\rho=r/N$ is rather small, for example
$1/90$ in the case $r=1$ and $N=90$, the empirical PT in our experiments
agrees much more closely with the nonlinear formula ${\cal M}(\rho)$ than it
does with the linear formula $6\rho$. Also, in the non square case $\beta\neq
1$, ${\cal M}\sim 2\rho(1+\beta+\sqrt{\beta})$, which is strictly smaller than
$6\rho$ for $\beta<1$, so the CR formula is noticeably less accurate than the
$\delta^{*}={\cal M}$ formula in the non square case.

Of the mathematical methods developed to identify phase transitions but which
are not based on combinatorial geometry or approximate message passing, the
most precise are based on the ‘Escape Through the Mesh’ (ETM) technique of
Yoram Gordon [35, 10]. ETM was used to prove upper bounds on the number of
Gaussian measurements needed for reconstruction of sparse signals by Stojnic
[35] and for low-rank matices by Oymak and Hassibi [10]. In particular, [10]
studies one of the two cases studied here and observes in passing that in the
square case, ETM gives bounds that seemingly agree with actual phase
transition measurements. Very recently, building on the same approach, a DJM-
style inequality $\delta^{*}(\kappa)\leq{\cal M}(\kappa)$ has been announced
by Oymak and Hassibi for a wide range of convex optimization problems in [36],
including nuclear norm minimization; [36] also presented empirical evidence
for $\delta^{*}(\rho)={\cal M}(\rho)$ in the square case ${\bf X}={Mat}_{N}$.

Our Contributions. This paper presents explicit formulas for the minimax MSE
of singular value thresholding in various cases, and shows that the new
formulas match the appropriate empirical phase transitions in a formal
comparison. Compared to earlier work, we make here the following
contributions:

  * •

A Broad Range of Phase Transition Studies for non square, square, and
symmetric nonnegative-definite matrix recovery from Gaussian measurements. We
also made certain nonGaussian measurements and observed similar results.

  * •

A Broad Range of Prediction formulas. We make available explicit formulas for
Minimax MSE in the square symmetric nonnegative-definite, or asymmetric case,
as well as non square cases.

  * •

Careful Empirical Technique, including the following:

    * Reproducibility. Publication of the code and data underlying our conclusions.

    * Validation. Matlab/CVX results were re-implemented in Python / CVXOPT, with similar results. Code was executed on 3 different computing clusters, with similar results.

    * Study of Finite-$N$-scaling. We studied tendencies of $a$,$b$ as $N$ varies, and observed behavior consistent with our asymptotic phase transition hypothesis. Studies at a single $N$ could only have shown that an empirical phase transition was near to a given theoretical curve, at a given $N$, but not shown the scaling behavior with $N$ that the main hypothesis properly involves.

##  5 Conclusions

For the problem of matrix recovery from Gaussian measurements, our
experiments, as well as those of others, document the existence of a
finite-$N$ phase transition. We compared our measured empirical phase
transition curve with a formula from the theory of matrix denoising and
observed a compelling match. Although the matrix recovery and matrix denoising
problems are superficially different, this match evidences a deeper
connection, such that mean squared error properties of a denoiser in a noise
removal problem give precisely the exact recovery properties of a matrix
recovery rule in a noiseless, but incomplete data problem.

This connection suggests both new limits on what is possible in the matrix
recovery problem, but also new ways of trying to reach those limits.

## Acknowledgments

Thanks to Iain Johnstone for advice at several crucial points. This work was
partially supported by NSF DMS 0906812 (ARRA). MG was supported by a William
R. and Sara Hart Kimball Stanford Graduate Fellowship and a Technion EE Sohnis
Promising Scientist Award. AM was partially supported by the NSF CAREER award
CCF-0743978 and the grant AFOSR/DARPA FA9550-12-1-0411.

## References

  * Candès and Recht [2008] Emmanuel J Candès and Benjamin Recht.  Exact Matrix Completion via Convex Optimization.  _Foundations of Computational Mathematics_ , 9(6):717–772, 2008.  URL http://arxiv.org/abs/0805.4471.
  * Gross [2011] David Gross.  Recovering Low-Rank Matrices From Few Coefficients in Any Basis.  _IEEE Transactions on Information Theory_ , 57(3):1548–1566, March 2011.  ISSN 0018-9448.  doi: 10.1109/TIT.2011.2104999.  URL http://ieeexplore.ieee.org/lpdocs/epic03/wrapper.htm?arnumber=5714248.
  * Recht et al. [2010a] Benjamin Recht, Maryam Fazel, and PA Parrilo.  Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization.  _SIAM review_ , 52(3):471–501, 2010a.  URL http://epubs.siam.org/doi/abs/10.1137/070697835.
  * Candès et al. [2009] Emmanuel J Candès, Xiaodong Li, Yi Ma, and John Wright.  Robust Principal Component Analysis ?  2009\.  URL http://arxiv.org/pdf/0912.3599.
  * Gross et al. [2010] David Gross, Yi-Kai Liu, Steven Flammia, Stephen Becker, and Jens Eisert.  Quantum State Tomography via Compressed Sensing.  _Physical Review Letters_ , 105(15):1–4, October 2010.  ISSN 0031-9007.  doi: 10.1103/PhysRevLett.105.150401.  URL http://link.aps.org/doi/10.1103/PhysRevLett.105.150401.
  * Keshavan et al. [2010] RH Keshavan, Andrea Montanari, and S Oh.  Matrix completion from noisy entries.  _The Journal of Machine Learning_ , 11:2057–2078, 2010\.  URL http://dl.acm.org/citation.cfm?id=1859890.1859920.
  * Recht et al. [2008] Benjamin Recht, W Xu, and Babak Hassibi.  Necessary and sufficient conditions for success of the nuclear norm heuristic for rank minimization.  In _Proceedings of the 47th IEEE conference on Decision and Control Cancun, Mexico_ , 2008.  URL http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=4739332.
  * Recht et al. [2010b] Benjamin Recht, Weiyu Xu, and Babak Hassibi.  Null space conditions and thresholds for rank minimization.  _Mathematical Programming_ , 127(1):175–202, October 2010b.  ISSN 0025-5610.  doi: 10.1007/s10107-010-0422-2.  URL http://www.springerlink.com/index/10.1007/s10107-010-0422-2.
  * Tanner and Wei [2012] Jared Tanner and KE Wei.  Normalized iterative hard thresholding for matrix completion.  2012\.  URL http://people.maths.ox.ac.uk/tanner/papers/TaWei_NIHT.pdf.
  * Oymak and Hassibi [2010] Samet Oymak and Babak Hassibi.  New Null Space Results and Recovery Thresholds for Matrix Rank Minimization.  page 28, November 2010.  URL http://arxiv.org/abs/1011.6326.
  * Donoho and Gavish [2013] David L. Donoho and Matan Gavish.  Minimax Risk of Matrix Denoising by Singular Value Thresholding.  Technical report, Stanford University Department of Statistics, 2013.
  * Fazel [2002] Maryam Fazel.  _Matrix rank minimization with applications_.  PhD thesis, Stanford University, 2002.  URL http://www.inatel.br/docentes/dayan/TP524/Artigos/MatrixRankMinimizationwithApplications.pdf.
  * Fazel et al. [2001] Maryam Fazel, H Hindi, and Stephen P. Boyd.  A rank minimization heuristic with application to minimum order system approximation.  In _American Control Conference, 2001_ , pages 4734–4739, 2001.  URL http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=945730.
  * Donoho and Tanner [2009a] David L. Donoho and Jared 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, Series A_ , 367(1906):4273–93, November 2009a.  ISSN 1364-503X.  doi: 10.1098/rsta.2009.0152.  URL http://www.ncbi.nlm.nih.gov/pubmed/19805445.
  * Monajemi et al. [2013] Hatef Monajemi, Sina Jafarpour, Matan Gavish, and David L. Donoho.  Deterministic matrices matching the compressed sensing phase transitions of Gaussian random matrices.  _Proceedings of the National Academy of Sciences of the United States of America_ , 110(4):1181–6, January 2013.  ISSN 1091-6490.  doi: 10.1073/pnas.1219540110.  URL http://www.pnas.org/cgi/content/long/110/4/1181.
  * Kalish [1990] Leslie A. Kalish.  Efficient design for estimation of median lethal dose and quantal dose-response curves.  _Biometrics_ , 46(3):737–48, September 1990.  ISSN 0006-341X.  URL http://www.ncbi.nlm.nih.gov/pubmed/2242412.
  * Grant and Boyd [2010] M. Grant and Stephen P. Boyd.  CVX: Matlab Software for Disciplined Convex Programming, 2010.  URL cxvr.com/cvx.
  * Sturm [1999] JF Sturm.  Using SeDuMi 1.02, a MATLAB toolbox for optimization over symmetric cones.  _Optimization methods and software_ , 1999.  URL http://www.tandfonline.com/doi/abs/10.1080/10556789908805766.
  * Toh et al. [1999] K.C. Toh, M.J. Todd, and Tutuncu R.H.  SDPT3 - a Matlab software package for semidefinite programming, Optimization Methods and Software, 1999.
  * Andersen et al. [2012] D. Andersen, J. Dahl, and L. Vandenberghe.  CVXOPT: Python Software for Convex Optimization, 2012.  URL http://abel.ee.ucla.edu/cvxopt.
  * Donoho et al. [2013] David L. Donoho, Matan Gavish, and Andrea Montanari.  Data for the article The Phase Transition of Matrix Recovery from Gaussian Measurements Matches the Minimax MSE of Matrix Denoising.  http://purl.stanford.edu/tz124hw0000, 2013.  URL purl.stanford.edu/tz124hw0000.  accessed 15 February 2013.
  * Donoho [2005] David L. Donoho.  High-Dimensional Centrally Symmetric Polytopes with Neighborliness Proportional to Dimension.  _Discrete & Computational Geometry_, 35(4):617–652, December 2005.  ISSN 0179-5376.  doi: 10.1007/s00454-005-1220-0.  URL http://www.springerlink.com/index/10.1007/s00454-005-1220-0.
  * Donoho and Tanner [2009b] David L. Donoho and Jared Tanner.  Counting faces of randomly projected polytopes when the projection radically lowers dimension.  _Journal of the American Mathematical Society_ , 22(1):1–53, 2009b.
  * Donoho and Tanner [2010] David L. Donoho and Jared Tanner.  Precise Undersampling Theorems.  _Proceedings of the IEEE_ , 98(6):913–924, June 2010.  ISSN 0018-9219.  doi: 10.1109/JPROC.2010.2045630.  URL http://ieeexplore.ieee.org/lpdocs/epic03/wrapper.htm?arnumber=5458001.
  * Donoho and Tanner [2005] David L. Donoho and Jared Tanner.  Neighborliness of randomly projected simplices in high dimensions.  _Proceedings of the National Academy of Sciences of the United States of America_ , 102(27):9452–7, July 2005.  ISSN 0027-8424.  doi: 10.1073/pnas.0502258102.  URL http://www.pnas.org/cgi/content/long/102/27/9452.
  * Donoho et al. [2009] David L. Donoho, Arian Maleki, and Andrea Montanari.  Message-passing algorithms for compressed sensing.  _Proceedings of the National Academy of Sciences of the United States of America_ , 106(45):18914–9, November 2009.  ISSN 1091-6490.  doi: 10.1073/pnas.0909892106.  URL http://www.pnas.org/cgi/content/long/106/45/18914.
  * Bayati and Montanari [2011] Mohsen Bayati and Andrea Montanari.  The Dynamics of Message Passing on Dense Graphs, with Applications to Compressed Sensing.  _IEEE Transactions on Information Theory_ , 57(2):764–785, February 2011.  ISSN 0018-9448.  doi: 10.1109/TIT.2010.2094817.  URL http://ieeexplore.ieee.org/lpdocs/epic03/wrapper.htm?arnumber=5695122.
  * Bayati and Montanari [2012] Mohsen Bayati and Andrea Montanari.  The LASSO Risk for Gaussian Matrices.  _IEEE Transactions on Information Theory_ , 58(4):1997–2017, April 2012.  ISSN 0018-9448.  doi: 10.1109/TIT.2011.2174612.  URL http://ieeexplore.ieee.org/lpdocs/epic03/wrapper.htm?arnumber=6069859.
  * Rangan [2011] Sundeep Rangan.  Generalized approximate message passing for estimation with random linear mixing.  In _2011 IEEE International Symposium on Information Theory Proceedings_ , pages 2168–2172. Ieee, July 2011.  ISBN 978-1-4577-0596-0.  doi: 10.1109/ISIT.2011.6033942.  URL http://ieeexplore.ieee.org/lpdocs/epic03/wrapper.htm?arnumber=6033942.
  * Vila and Schniter [2011] Jeremy Vila and Philip Schniter.  Expectation-maximization Bernoulli-Gaussian approximate message passing.  In _Record of the Forty Fifth Asilomar Conference on Signals, Systems and Computers (ASILOMAR), 2011_ , number July, pages 799–803, 2011.  ISBN 9781467303231.  URL http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=6190117.
  * Bayati et al. [2012] Mohsen Bayati, Marc Lelarge, and Andrea Montanari.  Universality in Polytope Phase Transitions and Message Passing Algorithms.  July 2012.  URL http://arxiv.org/abs/1207.7321.
  * Donoho et al. [2011] David L. Donoho, Iain M Johnstone, and Andrea Montanari.  Accurate Prediction of Phase Transitions in Compressed Sensing via a Connection to Minimax Denoising.  2011\.  URL http://http//arxiv.org/abs/1111.1041v2.
  * Oymak et al. [2011] Samet Oymak, Karthik Mohan, Maryam Fazel, and Babak Hassibi.  A simplified approach to recovery conditions for low rank matrices.  _2011 IEEE International Symposium on Information Theory Proceedings_ , pages 2318–2322, July 2011.  doi: 10.1109/ISIT.2011.6033976.  URL http://ieeexplore.ieee.org/lpdocs/epic03/wrapper.htm?arnumber=6033976.
  * Candès and Recht [2011] Emmanuel J Candès and Benjamin Recht.  Simple bounds for recovering low-complexity models.  _Mathematical Programming_ , (June 2011):1–11, 2011.  URL http://www.springerlink.com/index/655751h7568x3621.pdf.
  * Stojnic [2009] Mihailo Stojnic.  Various thresholds for l1-optimization in compressed sensing.  July 2009.  URL http://arxiv.org/abs/0907.3666.
  * Oymak and Hassibi [2012] Samet Oymak and Babak Hassibi.  On a Relation between the Minimax Denoising and the Phase Transitions of Convex Functions, 2012.  URL http://www.its.caltech.edu/~soymak/Relation.pdf.
  * Marcenko and Pastur [1967] VA Marcenko and LA Pastur.  Distribution of eigenvalues for some sets of random matrices.  _Mathematics USSR Sbornik_ , 1(4):457–483, 1967\.  URL http://iopscience.iop.org/0025-5734/1/4/A01.

##  Appendix A Asymptotic Minimax MSE Formula

The following provides explicit formulas for the matrix denoising minimax
curves $\mathcal{M}(\rho,\beta|Mat)$ and $\mathcal{M}(\rho|Sym)$ used above.
Please see [11] for the derivations. Computer programs that efficiently
calculate these quantities are provided in [21]. Let

$\displaystyle
P_{\gamma}(x;k)=\frac{1}{2\pi\gamma}\intop_{x}^{\gamma_{+}}t^{k-1}\sqrt{(\gamma_{+}-t)(t-\gamma_{-}})\,dt\,,$
(7)

where $\gamma_{\pm}=\left(1\pm\sqrt{\gamma}\right)^{2}$, denote the
complementary incomplete moments of the Marc̆enko-Pastur distribution [37].
Define

$\displaystyle\mathbf{M}_{\alpha}(\Lambda;\rho,\tilde{\rho})=\rho+\tilde{\rho}-\rho\tilde{\rho}+(1-\tilde{\rho})\Big{[}\rho\Lambda^{2}+$
$\displaystyle+\alpha(1-\rho)\Big{(}P_{\gamma}(\Lambda^{2};1)-2\Lambda
P_{\gamma}(\Lambda^{2};\tfrac{1}{2})+\Lambda^{2}P_{\gamma}(\Lambda^{2};0)\Big{)}\Big{]}\,,$

with
$\gamma=\gamma(\rho,\tilde{\rho})=(\tilde{\rho}-\rho\tilde{\rho})/(\rho-\rho\tilde{\rho})$.

#### Case ${\bf X}={Mat}_{M,N}$.

The minimax curve is given by
$\mathcal{M}(\rho,\beta|{Mat})=\inf_{\Lambda}\mathbf{M}_{1}(\Lambda;\rho,\beta\rho).$
The following minimaxity interpretation is proved in [11]

$\displaystyle\lim_{N\to\infty}\inf_{\lambda}\sup_{\stackrel{{\scriptstyle
X_{0}\in{\mathbb{R}}^{\lfloor\beta N\rfloor\times
N}}}{{rank(X_{0})\leq\rho\beta
N}}}MSE(\hat{X}_{\lambda},X_{0})=\mathcal{M}(\rho,\beta|{Mat})\,.$ (9)

#### Case ${\bf X}={Sym}_{N}$.

The minimax curve is given by
$\mathcal{M}(\beta|{Sym})=\inf_{\Lambda}\mathbf{M}_{1/2}(\Lambda;\rho,\rho).$
The following minimaxity interpretation is proved in [11]

$\displaystyle\lim_{N\to\infty}\inf_{\lambda}\sup_{\stackrel{{\scriptstyle
X_{0}\in S_{N}}}{{rank(X_{0})\leq\rho
N}}}MSE(\hat{X}_{\lambda},X_{0})=\mathcal{M}(\rho|{Sym})\,.$ (10)

#### Computing $\mathcal{M}(\rho,\beta|{\bf X})$.

The map $\Lambda\mapsto\mathbf{M}_{\alpha}(\Lambda;\rho,\tilde{\rho})$ is
convex. Solving $d\mathbf{M}_{\alpha}/d\Lambda=0$ we get that
$\text{argmin}_{\Lambda}\,\mathbf{M}_{\alpha}(\Lambda;\rho,\tilde{\rho})$ is
the unique root of the equation

$\displaystyle\Lambda^{-1}P_{\gamma}(\Lambda^{2};\tfrac{1}{2})-P_{\gamma}(\Lambda^{2};0)=\frac{\rho}{\alpha(1-\rho)}\,.$
(11)

The right hand side of (11) is decreasing in $\Lambda$ and the solution is
determined numerically by binary search.

For square matrices ($\rho=\tilde{\rho}$) (A) can be expressed using
elementary trigonometric functions. In [11] it is shown that

$\displaystyle\mathbf{M}_{\alpha}(\Lambda;\rho,\rho)=\rho(2-\rho)+(1-\rho)\big{[}\rho\Lambda^{2}+$
(12) $\displaystyle+\alpha(1-\rho)\left(Q_{2}\left(\Lambda\right)-2\lambda
Q_{1}\left(\Lambda\right)+\Lambda^{2}Q_{0}\left(\Lambda\right)\right)\big{]}\,.$

where

$\displaystyle Q_{0}(x)$ $\displaystyle=$
$\displaystyle\frac{1}{\pi}\intop_{x}^{2}\sqrt{4-x^{2}}=1-\frac{x}{2\pi}\sqrt{4-x^{2}}-\frac{2}{\pi}atan(\frac{x}{\sqrt{4-x^{2}}})$
$\displaystyle Q_{1}(x)$ $\displaystyle=$
$\displaystyle\frac{1}{\pi}\intop_{x}^{2}x\sqrt{4-x^{2}}=\frac{1}{3\pi}(4-x^{2})^{3/2}$
$\displaystyle Q_{2}(x)$ $\displaystyle=$
$\displaystyle\frac{1}{\pi}\intop_{x}^{2}x^{2}\sqrt{4-x^{2}}=1-\frac{1}{4\pi}x\sqrt{4-x^{2}}(x^{2}-2)-\frac{2}{\pi}asin(\frac{x}{2})$

are the complementary incomplete moments of the Quarter Circle law. Moreover

$\displaystyle\text{argmin}_{\Lambda}\mathbf{M}_{\alpha}(\Lambda;\rho,\rho)=2\cdot\sin\left(\theta_{\alpha}(\rho)\right)\,,$
(13)

where $\theta_{\alpha}(\rho)\in[0,\pi/2]$ is the unique solution to the
transcendental equation

$\displaystyle\theta+cot(\theta)\cdot\left(1-\frac{1}{3}cos^{2}(\theta)\right)=\frac{\pi(1+\alpha^{-1}\rho-\rho)}{2(1-\rho)}\,,$
(14)

which is a simplified version of (11).

#### Parametric representation of the minimax curves.

For square matrices ($\rho=\tilde{\rho}$) the minimax curves
$\mathcal{M}(\rho,1|{Mat})$ and $\mathcal{M}(\rho|{Sym})$ admit a parametric
representation in the $(\rho,\mathcal{M})$ plane using elementary
trigonometric functions, see [11]. As $\theta$ ranges over $[0,\pi/2]$,

$\displaystyle\rho(\theta)$ $\displaystyle=$ $\displaystyle
1-\frac{\pi/2}{\theta+(\cot(\theta)\cdot(1-\frac{1}{3}cos^{2}(\theta)))}$
$\displaystyle\mathcal{M}(\theta)$ $\displaystyle=$ $\displaystyle
2\rho(\theta)-\rho^{2}(\theta)+4\rho(\theta)(1-\rho(\theta))sin^{2}(\theta)$
$\displaystyle+$
$\displaystyle\frac{4}{\pi}(1-\rho)^{2}\left[(\pi-2\theta)(\frac{5}{4}-cos(\theta)^{2})+\frac{sin(2\theta)}{12}(cos(2\theta)-14)\right]$

is a parametric representation of $\mathcal{M}(\rho,1|{Mat})$, and similarly

$\displaystyle\rho(\theta)$ $\displaystyle=$ $\displaystyle
1-\frac{\theta+(\cot(\theta)\cdot(1-\frac{1}{3}cos^{2}(\theta)))-\pi/2}{\theta+(\cot(\theta)\cdot(1-\frac{1}{3}cos^{2}(\theta)))+\pi/2}$
$\displaystyle\mathcal{M}(\theta)$ $\displaystyle=$ $\displaystyle
2\rho(\theta)-\rho^{2}(\theta)+4\rho(\theta)(1-\rho(\theta))sin^{2}(\theta)$
$\displaystyle+$
$\displaystyle\frac{2}{\pi}(1-\rho)^{2}\left[(\pi-2\theta)(\frac{5}{4}-cos(\theta)^{2})+\frac{sin(2\theta)}{12}(cos(2\theta)-14)\right]$

is a parametric representation of $\mathcal{M}(\rho|{Sym})$.

##  Appendix B Summary of Empirical Results

$\rho$ | N | T | ${\cal M}(\rho)$ | $\hat{\delta}(\rho)$ | $a$ | $b$ | $Z$
---|---|---|---|---|---|---|---
1/10 | 40 | 400 | 0.351 | 0.352 | -0.128 | 182.978 | -0.581
50 | 400 | 0.351 | 0.350 | 0.282 | 200.131 | 1.213
60 | 400 | 0.351 | 0.352 | -0.096 | 221.212 | -0.398
80 | 400 | 0.351 | 0.350 | 0.415 | 295.049 | 1.452
100 | 400 | 0.351 | 0.349 | 0.641 | 383.493 | 1.900
1/8 | 32 | 400 | 0.414 | 0.412 | 0.262 | 119.424 | 1.457
48 | 400 | 0.414 | 0.413 | 0.262 | 204.025 | 1.120
64 | 400 | 0.414 | 0.412 | 0.560 | 270.147 | 2.006
80 | 400 | 0.414 | 0.413 | 0.296 | 290.536 | 1.057
1/6 | 36 | 400 | 0.507 | 0.505 | 0.356 | 150.509 | 1.749
60 | 400 | 0.507 | 0.506 | 0.137 | 216.067 | 0.571
1/4 | 32 | 400 | 0.655 | 0.651 | 0.398 | 97.107 | 2.375
48 | 400 | 0.655 | 0.653 | 0.356 | 191.373 | 1.554
64 | 400 | 0.655 | 0.653 | 0.507 | 242.289 | 1.935
80 | 400 | 0.655 | 0.651 | 1.312 | 333.842 | 3.584
1/3 | 36 | 400 | 0.765 | 0.765 | -0.006 | 148.439 | -0.029
60 | 400 | 0.765 | 0.762 | 1.145 | 317.633 | 3.348
90 | 400 | 0.765 | 0.762 | 1.487 | 384.289 | 3.610
$1/2$ | 50 | 400 | 0.905 | 0.903 | 0.658 | 260.939 | 2.361
$3/4$ | 40 | 400 | 0.989 | 0.986 | 1.535 | 709.400 | 3.837
$9/10$ | 50 | 400 | 0.999 | 0.998 | 2.478 | 2980.216 | 3.435
Figure 6: Results with ${\bf X}={Mat}_{N}$.  $\rho$ | N | T | ${\cal M}(\rho)$ | $\hat{\delta}(\rho)$ | $a$ | $b$ | $Z$
---|---|---|---|---|---|---|---
1/10 | 40 | 800 | 0.315 | 0.310 | 0.787 | 148.605 | 4.268
60 | 800 | 0.315 | 0.312 | 0.519 | 186.020 | 2.640
| 80 | 800 | 0.315 | 0.311 | 0.995 | 265.114 | 3.864
1/8 | 40 | 800 | 0.371 | 0.369 | 0.270 | 141.743 | 1.585
56 | 800 | 0.371 | 0.369 | 0.230 | 176.063 | 1.212
| 80 | 800 | 0.371 | 0.369 | 0.368 | 244.721 | 1.628
1/7 | 35 | 800 | 0.407 | 0.403 | 0.623 | 126.659 | 3.730
49 | 800 | 0.407 | 0.404 | 0.496 | 140.034 | 2.878
| 70 | 800 | 0.407 | 0.404 | 0.581 | 184.556 | 2.898
1/6 | 36 | 800 | 0.453 | 0.447 | 0.702 | 116.847 | 4.607
60 | 800 | 0.453 | 0.447 | 1.097 | 193.706 | 5.133
| 90 | 800 | 0.453 | 0.449 | 1.170 | 311.016 | 4.246
1/5 | 40 | 800 | 0.511 | 0.505 | 0.786 | 134.698 | 5.448
50 | 800 | 0.511 | 0.507 | 0.707 | 164.398 | 4.528
| 80 | 800 | 0.511 | 0.507 | 0.859 | 205.621 | 4.778
1/4 | 36 | 800 | 0.588 | 0.579 | 0.967 | 110.989 | 7.120
60 | 800 | 0.588 | 0.582 | 1.155 | 181.174 | 6.379
| 80 | 800 | 0.588 | 0.583 | 1.450 | 271.651 | 6.066
1/3 | 36 | 800 | 0.694 | 0.685 | 1.085 | 124.413 | 3.926
60 | 800 | 0.694 | 0.688 | 1.034 | 192.690 | 5.710
| 90 | 800 | 0.694 | 0.689 | 1.212 | 263.543 | 5.481
1/2 | 36 | 800 | 0.844 | 0.837 | 1.306 | 180.027 | 6.965
60 | 800 | 0.844 | 0.838 | 1.872 | 320.651 | 6.425
| 90 | 800 | 0.844 | 0.840 | 1.884 | 438.138 | 5.513
Figure 7: Results with ${\bf X}={Sym}_{N}$.

$\beta$ | $\rho$ | ${\cal M}(\rho)$ | $\hat{\delta}(\rho)$ | $a$ | $b$ | $Z$ | $\sqrt{M\cdot N}$
---|---|---|---|---|---|---|---
$1/4$ | 0.100 | 0.241 | 0.241 | 0.020 | 215.506 | 0.085 | 60.000
0.125 | 0.290 | 0.290 | -0.037 | 322.468 | -0.127 | 64.000
0.143 | 0.323 | 0.321 | 0.468 | 238.477 | 1.813 | 70.000
0.167 | 0.365 | 0.362 | 0.884 | 247.690 | 3.109 | 60.000
0.200 | 0.421 | 0.422 | -0.153 | 246.541 | -0.597 | 60.000
0.250 | 0.498 | 0.496 | 0.339 | 214.134 | 1.405 | 64.000
0.333 | 0.610 | 0.607 | 0.856 | 284.629 | 2.834 | 60.000
| 0.500 | 0.788 | 0.786 | 0.635 | 297.803 | 2.124 | 64.000
$1/3$ | 0.100 | 0.255 | 0.255 | -0.105 | 269.362 | -0.395 | 51.962
0.143 | 0.340 | 0.340 | 0.031 | 210.274 | 0.133 | 48.497
0.167 | 0.383 | 0.381 | 0.373 | 169.227 | 1.733 | 51.962
0.200 | 0.440 | 0.437 | 0.589 | 197.335 | 2.457 | 51.962
0.250 | 0.517 | 0.517 | 0.113 | 218.479 | 0.470 | 55.426
0.333 | 0.629 | 0.626 | 0.661 | 219.930 | 2.584 | 51.962
0.500 | 0.801 | 0.799 | 0.572 | 274.122 | 2.036 | 55.426
$1/2$ | 0.200 | 0.475 | 0.474 | 0.127 | 159.697 | 0.616 | 42.426
0.250 | 0.554 | 0.553 | 0.258 | 201.035 | 1.114 | 45.255
0.286 | 0.604 | 0.603 | 0.307 | 167.955 | 1.434 | 49.497
0.333 | 0.665 | 0.662 | 0.486 | 155.704 | 2.327 | 42.426
0.400 | 0.738 | 0.737 | 0.292 | 184.091 | 1.308 | 42.426
0.500 | 0.827 | 0.825 | 0.460 | 228.287 | 1.826 | 45.255
0.667 | 0.930 | 0.927 | 0.729 | 248.019 | 3.073 | 42.426
$3/5$ | 0.100 | 0.296 | 0.295 | 0.146 | 191.950 | 0.646 | 38.730
0.125 | 0.352 | 0.351 | 0.242 | 219.767 | 1.000 | 61.968
0.143 | 0.389 | 0.388 | 0.105 | 179.344 | 0.484 | 54.222
0.167 | 0.436 | 0.433 | 0.425 | 198.788 | 1.814 | 46.476
0.200 | 0.495 | 0.494 | 0.237 | 146.409 | 1.196 | 38.730
0.250 | 0.575 | 0.573 | 0.361 | 164.339 | 1.704 | 46.476
0.333 | 0.686 | 0.684 | 0.509 | 258.045 | 1.886 | 46.476
0.500 | 0.843 | 0.842 | 0.153 | 156.137 | 0.751 | 46.476
$2/3$ | 0.100 | 0.305 | 0.306 | -0.168 | 221.047 | -0.697 | 36.742
0.125 | 0.363 | 0.362 | 0.064 | 151.402 | 0.321 | 39.192
0.143 | 0.401 | 0.399 | 0.263 | 142.052 | 1.346 | 34.293
0.167 | 0.448 | 0.446 | 0.271 | 157.783 | 1.323 | 36.742
0.200 | 0.509 | 0.510 | -0.209 | 155.201 | -1.034 | 36.742
0.250 | 0.589 | 0.587 | 0.315 | 152.493 | 1.555 | 39.192
0.333 | 0.700 | 0.697 | 0.408 | 190.613 | 1.768 | 36.742
0.500 | 0.854 | 0.853 | 0.257 | 232.865 | 1.034 | 39.192
$3/4$ | 0.100 | 0.317 | 0.317 | -0.043 | 163.597 | -0.207 | 34.641
0.125 | 0.376 | 0.375 | 0.237 | 192.933 | 1.047 | 55.426
0.143 | 0.415 | 0.412 | 0.562 | 178.363 | 2.476 | 48.497
0.167 | 0.463 | 0.461 | 0.304 | 174.599 | 1.405 | 41.569
0.200 | 0.525 | 0.524 | 0.159 | 120.459 | 0.886 | 34.641
0.250 | 0.606 | 0.604 | 0.276 | 154.259 | 1.353 | 41.569
0.333 | 0.716 | 0.714 | 0.378 | 180.963 | 1.682 | 41.569
0.500 | 0.867 | 0.866 | 0.341 | 283.911 | 1.227 | 41.569
$4/5$ | 0.100 | 0.324 | 0.323 | 0.220 | 174.269 | 1.019 | 44.721
0.125 | 0.384 | 0.381 | 0.423 | 160.001 | 1.999 | 44.721
0.143 | 0.423 | 0.422 | 0.273 | 248.331 | 1.055 | 62.610
0.167 | 0.472 | 0.472 | -0.025 | 190.094 | -0.112 | 40.249
0.200 | 0.534 | 0.533 | 0.204 | 195.806 | 0.889 | 44.721
0.250 | 0.616 | 0.616 | -0.076 | 197.138 | -0.333 | 44.721
0.333 | 0.726 | 0.725 | 0.213 | 203.942 | 0.910 | 40.249
0.500 | 0.875 | 0.873 | 0.412 | 215.351 | 1.689 | 44.721
Figure 8: Results with non square matrices ${\bf X}={Mat}_{M,N}$. Each row
based on $T=400$ Monte Carlo trials.

N | r | $\rho$ | $6\rho$ | ${\cal M}(\rho)$ | $\hat{\delta}(\rho)$ | $a$ | $b$ | $Z$
---|---|---|---|---|---|---|---|---
90 | 1 | 0.011 | 0.067 | 0.059 | 0.054 | 3.595 | 828.597 | 4.627
80 | 1 | 0.013 | 0.075 | 0.065 | 0.059 | 2.766 | 534.431 | 5.117
70 | 1 | 0.014 | 0.086 | 0.072 | 0.069 | 1.444 | 432.528 | 3.962
60 | 1 | 0.017 | 0.100 | 0.081 | 0.078 | 1.374 | 434.911 | 3.542
50 | 1 | 0.020 | 0.120 | 0.094 | 0.091 | 1.159 | 348.687 | 3.226
90 | 2 | 0.022 | 0.133 | 0.103 | 0.101 | 1.225 | 611.400 | 2.499
40 | 1 | 0.025 | 0.150 | 0.114 | 0.114 | -0.169 | 328.547 | -0.513
40 | 1 | 0.025 | 0.150 | 0.114 | 0.114 | -0.169 | 328.547 | -0.513
70 | 2 | 0.029 | 0.171 | 0.127 | 0.125 | 6.383 | 2746.893 | 0.007
30 | 1 | 0.033 | 0.200 | 0.145 | 0.145 | 0.013 | 137.638 | 0.051
30 | 1 | 0.033 | 0.200 | 0.145 | 0.145 | 0.013 | 137.638 | 0.051
30 | 1 | 0.033 | 0.200 | 0.145 | 0.145 | 0.013 | 137.638 | 0.051
80 | 3 | 0.037 | 0.225 | 0.160 | 0.158 | 0.493 | 345.100 | 1.111
50 | 2 | 0.040 | 0.240 | 0.169 | 0.166 | 1.071 | 364.190 | 1.981
70 | 3 | 0.043 | 0.257 | 0.179 | 0.177 | 3.442 | 1883.417 | 0.006
90 | 4 | 0.044 | 0.267 | 0.184 | 0.183 | 2.044 | 1773.906 | 0.006
20 | 1 | 0.050 | 0.300 | 0.203 | 0.203 | -0.055 | 94.130 | -0.219
20 | 1 | 0.050 | 0.300 | 0.203 | 0.203 | -0.055 | 94.130 | -0.219
20 | 1 | 0.050 | 0.300 | 0.203 | 0.203 | -0.055 | 94.130 | -0.219
20 | 1 | 0.050 | 0.300 | 0.203 | 0.203 | -0.055 | 94.130 | -0.219
70 | 4 | 0.057 | 0.343 | 0.226 | 0.227 | -2.755 | 2732.891 | -0.001
50 | 3 | 0.060 | 0.360 | 0.235 | 0.230 | 6.462 | 1302.371 | 0.008
30 | 2 | 0.067 | 0.400 | 0.256 | 0.255 | 0.106 | 152.849 | 0.290
30 | 2 | 0.067 | 0.400 | 0.256 | 0.255 | 0.106 | 152.849 | 0.290
40 | 3 | 0.075 | 0.450 | 0.281 | 0.280 | 0.114 | 191.765 | 0.249
50 | 4 | 0.080 | 0.480 | 0.296 | 0.294 | 0.361 | 225.965 | 0.622
20 | 2 | 0.100 | 0.600 | 0.351 | 0.345 | 0.387 | 60.051 | 1.345
20 | 2 | 0.100 | 0.600 | 0.351 | 0.345 | 0.387 | 60.051 | 1.345
20 | 2 | 0.100 | 0.600 | 0.351 | 0.345 | 0.387 | 60.051 | 1.345
30 | 4 | 0.133 | 0.800 | 0.434 | 0.434 | -0.005 | 115.737 | -0.011
20 | 3 | 0.150 | 0.900 | 0.472 | 0.480 | -0.690 | 86.275 | -1.510
20 | 4 | 0.200 | 1.200 | 0.572 | 0.587 | -0.994 | 64.706 | -2.057
Figure 9: Results with low rank square matrices ${\bf X}={Mat}_{N,N}$. Each row based on $T=400$ Monte Carlo trials. $\rho$ | $N$ | ${\cal M}(\rho)$ | $\hat{\delta}(\rho)$ | $a$ | $b$ | $Z$
---|---|---|---|---|---|---
$1/10$ | 40.000 | 0.351 | 0.350 | 0.187 | 170.554 | 0.878
60.000 | 0.351 | 0.350 | 0.239 | 286.765 | 0.863
80.000 | 0.351 | 0.351 | 0.166 | 290.560 | 0.597
$1/8$ | 40.000 | 0.414 | 0.413 | 0.145 | 146.920 | 0.736
56.000 | 0.414 | 0.412 | 0.552 | 215.088 | 2.220
80.000 | 0.414 | 0.413 | 0.667 | 375.205 | 1.992
$1/7$ | 35.000 | 0.456 | 0.457 | -0.158 | 141.742 | -0.813
49.000 | 0.456 | 0.455 | 0.158 | 216.040 | 0.659
70.000 | 0.456 | 0.454 | 0.911 | 394.564 | 2.527
$1/6$ | 36.000 | 0.507 | 0.505 | 0.370 | 159.357 | 1.768
60.000 | 0.507 | 0.505 | 0.570 | 245.893 | 2.137
90.000 | 0.507 | 0.505 | 0.689 | 335.299 | 2.168
$1/5$ | 40.000 | 0.572 | 0.571 | 0.089 | 172.047 | 0.417
60.000 | 0.572 | 0.568 | 0.912 | 248.273 | 3.188
80.000 | 0.572 | 0.570 | 0.598 | 285.542 | 2.073
$1/3$ | 36.000 | 0.765 | 0.759 | 0.752 | 124.479 | 3.825
60.000 | 0.765 | 0.763 | 0.608 | 234.115 | 2.324
81.000 | 0.765 | 0.762 | 1.053 | 288.525 | 3.306
$1/2$ | 36.000 | 0.905 | 0.903 | 0.487 | 229.398 | 1.914
60.000 | 0.905 | 0.902 | 1.403 | 379.105 | 3.518
90.000 | 0.905 | 0.903 | 11.365 | 3971.012 | 0.008
Figure 10: Results with Rademacher measurements of square matrices ${\bf
X}={Mat}_{N,N}$. Each row based on $T=400$ Monte Carlo trials.

##  Appendix C Data Deposition

The data have been deposited in a text file at [21]. A typical fragment of the
file is given here:

    Line Project             Experiment      M  N  S Instance rank  rho              delta            Err0       Err1  Err2
    3381 Nuc_CVX_20121129dii Nuc_CVX_N12S01m 12 12 1 m 4 0.333333333333333 0.73395061728395 0.016630719182629 0 0.444444444444444
    3382 Nuc_CVX_20121129dii Nuc_CVX_N12S01m 12 12 1 m 4 0.333333333333333 0.739213775178687 0.0159010475527301 0 0.465277777777778
    3383 Nuc_CVX_20121129dii Nuc_CVX_N12S01m 12 12 1 m 4 0.333333333333333 0.744476933073424 0.0161172486497232 0 0.416666666666667
    3384 Nuc_CVX_20121129dii Nuc_CVX_N12S01m 12 12 1 m 4 0.333333333333333 0.749740090968161 0.00149080158529591 1 1
    3385 Nuc_CVX_20121129dii Nuc_CVX_N12S01m 12 12 1 m 4 0.333333333333333 0.755003248862898 0.0340839945130298 0 0.173611111111111
    3386 Nuc_CVX_20121129dii Nuc_CVX_N12S01m 12 12 1 m 4 0.333333333333333 0.760266406757635 0.0235186056093925 0 0.361111111111111
    3387 Nuc_CVX_20121129dii Nuc_CVX_N12S01m 12 12 1 m 4 0.333333333333333 0.765529564652372 0.0136400454215757 0 0.506944444444444
    3388 Nuc_CVX_20121129dii Nuc_CVX_N12S01m 12 12 1 m 4 0.333333333333333 0.770792722547109 1.11644848368808e-09 1 1
    3389 Nuc_CVX_20121129dii Nuc_CVX_N12S01m 12 12 1 m 4 0.333333333333333 0.776055880441845 0.00194118165102407 0 1
    3390 Nuc_CVX_20121129dii Nuc_CVX_N12S01m 12 12 1 m 4 0.333333333333333 0.781319038336582 0.0059943065062774 0 0.902777777777778
    3391 Nuc_CVX_20121129dii Nuc_CVX_N12S01m 12 12 1 m 4 0.333333333333333 0.786582196231319 1.83878516274726e-09 1 1
    3392 Nuc_CVX_20121129dii Nuc_CVX_N12S01m 12 12 1 m 4 0.333333333333333 0.791845354126056 7.85225405546251e-09 1 1
    3393 Nuc_CVX_20121129dii Nuc_CVX_N12S01m 12 12 1 m 4 0.333333333333333 0.797108512020793 4.47138355029567e-10 1 1
    3394 Nuc_CVX_20121129dii Nuc_CVX_N12S01m 12 12 1 m 4 0.333333333333333 0.80237166991553 1.03566257175607e-08 1 1
    3395 Nuc_CVX_20121129dii Nuc_CVX_N12S01m 12 12 1 m 4 0.333333333333333 0.807634827810266 2.62954112892325e-09 1 1
    3396 Nuc_CVX_20121129dii Nuc_CVX_N12S01m 12 12 1 m 4 0.333333333333333 0.812897985705003 0.00120410713874671 1 1
    3397 Nuc_CVX_20121129dii Nuc_CVX_N12S01m 12 12 1 m 4 0.333333333333333 0.81816114359974 1.03259685844506e-08 1 1
    3398 Nuc_CVX_20121129dii Nuc_CVX_N12S01m 12 12 1 m 4 0.333333333333333 0.823424301494477 2.52202641519081e-10 1 1
    3399 Nuc_CVX_20121129dii Nuc_CVX_N12S01m 12 12 1 m 4 0.333333333333333 0.828687459389214 1.94437637948599e-09 1 1
    3400 Nuc_CVX_20121129dii Nuc_CVX_N12S01m 12 12 1 m 4 0.333333333333333 0.83395061728395 1.69829115472996e-10 1 1

The fields have the following meaning

  * •

Line – Line number in file; in the above example, lines 3381-3400.

  * •

Project – File identifier – allows identification of code and logs that
generated these data; in the above fragment, ’Nuc_CVX_20121129dii’.

  * •

Experiment – File identified – allows identification of code and logs that
generated these data; in the above fragment ’Nuc_CVX_N12S01m’.

  * •

M,N – matrix size of $X_{0}$, i.e $M$ by $N$ matrix; in the above fragment
$M=N=12$.

  * •

S – number of matrices in a stack (see below)

  * •

Instance – alphabetic code a-t, identifying one of 20 identical runs which
generated this result; in the above fragment, ’m’.

  * •

rank – integer rank of matrix; in the above fragment ’m’.

  * •

rho – fraction in $[0,1]$, $\rho=\mbox{rank}/N$; in the above fragment, $1/3$.

  * •

delta – $\delta=n/(MN)$ in case of asymmetric matrix, or $\delta=2*n/(N(N+1))$
in case of symmetric matrix.

  * •

Err0 – $\|\hat{X}-X_{0}\|_{F}/(NM)^{1/2}$.

  * •

Err1 – 1 iff $\|\hat{X}-X_{0}\|_{F}/\|X_{0}\|_{F}<tol$, and $0$ otherwise

  * •

Err2 – fraction of entries with discrepancy $|\hat{X}(i,j)-X_{0}(i,j)|<tol$.

Additional concepts:

  * •

Numerical Tolerance. In our experiments, we used a numerical error tolerance
parameter $tol=0.001$.

  * •

Our experiments also extensively covered cases where $X_{0}$ is a ’stack of
matrices’, i.e. a 3-way array $M\times N\times S$, where $S$ is the number of
items in the stack. The only cases of interest for this paper are $S=1$.

##  Appendix D Code Deposition

There are two types of code deposition.

  * •

Reproduction from Data Deposition. The code that actually makes the figures
and tables we presented in this paper, starting from the data deposition. This
is deposited at [21]. The code we actually ran to create our figures and
tables is a set of R scripts, and was run on a Mac OS X. We believe the same
code runs with minimal changes on a LINUX environment.

  * •

RunMyCode Deposition. For readers who wish simply to compute the value of the
Minimax Mean-Squared Error over each of the matrix classes we considered, we
offer a Minimax MSE calculator at RunMyCode.org.

  * •

Full Code and Results Deposition. At [21], we also offer a literal dump of the
code we ran and all the logs and result files we obtained.

We believe the first two items are self-documenting. The third item can be
explained as follows. Our database of the experiment and all results is
contained in a unix directory tree rooted at exp.

Inside directory exp one finds further directories, as indicated below:

    Nuc_CVX_20120605a
    Nuc_CVX_20120607a
    ...
    Nuc_CVX_20121107a
    Nuc_CVX_20121113a
    ...
    Nuc_CVX_20121120a
    Nuc_CVX_20121121a
    ...
    SDP_CVX_20121230h
    SDP_CVX_20121230i
    ..
    Nuc_CVX_20130110a
    Nuc_CVX_20130110c
    ...
    Nuc_CVX_20130117Cc
    Nuc_CVX_20130117Cd
    ...
    Nuc_CVX_20130126Df
    Nuc_CVX_20130126Dg
    Nuc_CVX_20130126Dh

These directory names are precisely the Project names used in the data
deposition. Say we look inside one of these directories, for example
’Nuc_CVX_20121121b’. We will find a directory called bin containing software,
and a list of further directories. We excerpt from a 2-column listing of those
directories:

    Nuc_CVX_N05S1aΨΨΨNuc_CVX_N30S1c
    Nuc_CVX_N05S1bΨΨΨNuc_CVX_N30S1d
    Nuc_CVX_N05S1cΨΨΨNuc_CVX_N30S1e
    ...
    Nuc_CVX_N10S1jΨΨΨNuc_CVX_N35S1l
    Nuc_CVX_N10S1kΨΨΨNuc_CVX_N35S1m
    Nuc_CVX_N10S1lΨΨΨNuc_CVX_N35S1n
    ...
    Nuc_CVX_N15S1hΨΨΨNuc_CVX_N40S1j
    Nuc_CVX_N15S1iΨΨΨNuc_CVX_N40S1k
    Nuc_CVX_N15S1jΨΨΨNuc_CVX_N40S1l
    ...
    Nuc_CVX_N20S1oΨΨΨNuc_CVX_N45S1q
    Nuc_CVX_N20S1pΨΨΨNuc_CVX_N45S1r
    Nuc_CVX_N20S1qΨΨΨNuc_CVX_N45S1s
    ...
    Nuc_CVX_N25S1nΨΨΨNuc_CVX_N50S1p
    Nuc_CVX_N25S1oΨΨΨNuc_CVX_N50S1q
    Nuc_CVX_N25S1pΨΨΨNuc_CVX_N50S1r
    ...

These directory names are precisely the ’Experiment’ values seen in the data
deposition. Inside the bin directory we find the matlab code used in common by
all the above experiments.

    Nuc_CVX_20121121b$ ls bin
    aveRank.mΨΨΨsolveNuc_CVX_Stack_Arb.m
    predNucPT.mΨΨΨstackRankRMatrix.m
    rankRMatrix.m

Inside one of the experiment directories we find experiment-specific files (in
both Matlab and Bash script) as well as output files. The subdirectory logs
contains logs that were created while the jobs were running.

    Nuc_CVX_20121121b $ ls -R -C1 Nuc_CVX_N25S1o
    bashMain.sh
    PTExperiment.mΨ
    matlabMain.mΨ
    results.mat
    randState.matΨ

    Nuc_CVX_N25S1o/logs:
    Nuc_CVX_N25S1o.stderrΨΨΨ
    Nuc_CVX_N25S1o.stdout
    runMatlab.20121121b_Nuc_CVX_N25S1o.log

Here the .mat files were produced by matlab during the running of the
experiment.

  * •

randState.mat preserves the state of the random number generator at the
beginning of that experiment.

  * •

results.mat gives the results of the individual problem instances; typically
in the form $(r,n,N,M,Err0,Err1,Err2)$.


# A spinwave Ising machine

Artem Litvinenko1∗, Roman Khymyn1, Victor H. González1, Ahmad A. Awad1,
Vasyl Tyberkevych2, Andrei Slavin2 and Johan Åkerman1∗

1Department of Physics, University of Gothenburg, Fysikgränd 3, 412 96
Gothenburg, Sweden
2Department of Physics, Oakland University, 48309, Rochester, Michigan, USA,

∗To whom correspondence should be addressed; E-mails:
<EMAIL_ADDRESS><EMAIL_ADDRESS>
> We demonstrate a spin-wave-based time-multiplexed Ising Machine (SWIM),
> implemented using a 5 $\mu$m thick Yttrium Iron Garnet (YIG) film and off-
> the-shelf microwave components. The artificial Ising spins consist of 34–68
> ns long 3.125 GHz spinwave RF pulses with their phase binarized using a
> phase-sensitive microwave amplifier. Thanks to the very low spinwave group
> velocity, the 7 mm long YIG waveguide can host an 8-spin MAX-CUT problem and
> solve it in less than 4 $\mu$s while consuming only 7 $\mu$J. Using a real-
> time oscilloscope, we follow the temporal evolution of each spin as the SWIM
> minimizes its energy and find both uniform and domain-propagation-like
> switching of the spin state. The SWIM has the potential for substantial
> further miniaturization, scalability, speed, and reduced power consumption,
> and may become a versatile platform for commercially feasible optimization
> problem solvers with high performance.

## Introduction

As Moore’s law comes to an end due to physical limitations, while the amount
processing data is constantly growing, novel analog and digital computing
paradigm are being investigated. Large part of data processing tasks consist
of combinatorial optimization problems and require special purpose
accelerators for power efficient and fast computing. Combinatorial
optimization is essential to select the optimal path from an enormous range of
choices and appears in various social and industrial fields such as optimizing
financial transactions (?), logistics, including travel (?) and channel
assignment for wireless communications (?), genetic engineering and molecular
design for drug discovery (?). Conventional computers based on Von Neumann
architectures are inefficient in solving hard combinatorial optimization
problems due to the factorial growth of all possible combinations to be
evaluated using brute force. Fortunately, Ising machines have emerged as a
promising non-von-Neumann computing scheme that can accelerate computation of
NP-hard optimization problems.

An Ising machine maps an NP-problem onto the Ising Hamiltonian of an array of
$N$ binarized and interacting physical entities, commonly referred to as spins
$s_{i}=\pm 1$.

$H(s_{1},...,s_{N})=-\sum_{i<j}J_{ij}s_{i}s_{j}-\sum_{i=1}h_{i}s_{i}$ (1)

The coding of a specific problem is performed by adjusting the pairwise
coupling $J_{ij}$ terms such that the lowest energy spin glass configuration
represents the solution to the NP-problem. This is possible because spin glass
solutions belong to the NP-hard complexity class (?), and many NP-complete and
NP-hard problems can be reformulated as Ising Hamiltonians (?). Then, the
energy minimization of the entire array can be employed as a solution search
algorithm in a process called annealing (?).

Various types of Ising machines (?) have been proposed in the recent decades
exploiting novel physical building blocks such as spintronic devices (?, ?, ?,
?), memristor crossbars (?), quantum superconducting junctions (?, ?), metal-
insulator relaxation oscillators (?), and degenerate optical parametric
oscillators (?, ?, ?) as well as employing conventional CMOS technology with
analog electric oscillators (?, ?), FPGAs (?, ?). In particular, optical
Coherent Ising Machines have attracted great attention due to their high
computational speed, time-multiplexing method that provides all-to-all Ising
spin connections and the largest amount of supported Ising spins amongst all
other implementations.

Optical CIMs use time-multiplexed degenerate optical parametric oscillators
(DOPOs) that propagate in a ring cavity in the form of short light-wave
pulses. All-to-all coupling between DOPOs can be implemented either
physically, via external optical delay lines (?, ?, ?, ?), or electronically,
using FPGA blocks (?, ?, ?, ?). In contrast, spatially-resolved IMs based on
physical oscillators are inherently limited to low density connectivity as the
number of all-to-all interconnections grows as $O(n^{2})$ and hence represent
an unsolvable problem for 2.5-D mass electronics production technology. DOPOs
phase degeneracy allows one to exploit the bistable DOPOs phase as equivalent
two-state Ising spins. CIM belongs to a annealing scheme class known as
dynamical system solvers as their collective steady-state oscillation
minimizes the individual phases described by the Kuramoto model.

$\dot{\Phi}_{i}=\omega_{0}+K\sum_{j}J_{ij}\sin(\Phi_{i}-\Phi_{j})+Kh_{i}\sin(\Phi_{i}-\omega_{0}t)$
(2)

Where $K$ is a global coupling factor. See that eq.2 reduces to eq.1 for
$\sin(\Phi_{i}-\Phi_{j})=\pm 1$, promoting binarization which can be done by
the parametric pumping at double frequency which will result in an additional
term in eq.2. In the rotating frame, the steady-state solution
($\dot{\Phi}_{i}=0$) necessitates $\Phi_{i}-\Phi_{j}=0,\pi$ and
$\Phi_{i}-\omega_{0}t=0,\pi$ for $N$ mutually coupled oscillators. This state
will have the properties of a spin glass solution and thus the architecture
can be used as scalable NP-solvers, as shown in a progressive evolution from
100 to 100,000 optical spins in ring-based CIMs. Nevertheless, the commercial
feasibility of optical CIMs remains elusive as the technology requires optical
tables, kilowatts of power, and kilometers of optical fibers, which altogether
blocks its further development from a proof-of-principle demonstration to a
miniaturized commercially viable device.

In this work, we present a spinwave Ising machine (SWIM) with artificial spin
states implemented via the phase of spinwave RF pulses propagating in an
Yttrium iron garnet (YIG) thin film. Spinwaves are being actively used for
implementation of signal processing (?, ?, ?, ?, ?, ?, ?, ?), frequency
synthesis (?, ?, ?), logic (?, ?, ?, ?, ?, ?, ?) and computing (?, ?, ?, ?, ?)
applications due to their inherent non-linearities and exceptionally slow
propagation velocities. In comparison to optical pulses, the use of spinwaves
allows for 5–7 orders of magnitude reduction in delay line length. Another
advantage comes from the simplicity and efficiency of microwave current-to-
spinwave transducers and the off-the-shelf availability of power-efficient
microwave electronic components and circuits for the signal processing
required in a time-multiplexing architecture of Ising machines. In particular,
it allows the use of electronic phase-sensitive amplification blocks consuming
only milliwatts of power in contrast to optical kilowatt-power parametric
phase-sensitive amplifiers. Altogether, the use of propagating spinwaves for
the design of a time-multiplexed Ising machine brings the concept from a
laboratory demonstrator to a commercially feasible technology with potential
for CMOS integration.

## Design and characteristics of the spinwave Ising Machine

The SWIM employs a time-multiplexing approach similar to optical CIMs (?, ?,
?, ?, ?, ?, ?, ?) but in the microwave frequency domain, which leads to
certain signal processing modifications and substantial performance
improvements. As shown in Fig. 1, it consists of two blocks: a) an electronic
part for linear and parametric amplification of the RF pulses, RF pulse
interconnection, and RF pulse measurement, and b) a physical YIG spinwave
waveguide where all the spinwave RF pulses are excited by a first transducer
(antenna), allowed to propagate along the film length, and finally transformed
back to a microwave RF pulses by a second transducer. In contrast to optical
CIMs, where optical pulses propagate and are amplified in an optical ring
system, never leaving the ring, the SWIM is a multi-physical system where
amplification of the propagating RF pulses and their further signal processing
is performed outside of the YIG waveguide in a power-efficient electronic
system. This is possible because the carrier frequency of the propagating
spinwaves is around 3 GHz and, therefore, can be easily handled by inexpensive
and readily available commercial RF components.

Spin waves are a fundamental type of excitation in magnetic systems (?, ?, ?,
?, ?), representing collective precession of elementary magnetic moments
coupled via exchange and dipole-dipole interactions. The spin-wave excitation
frequency depends on the strength and orientation of the static magnetic field
applied to the film and on the magnetic properties of its material, such as
the saturation magnetization ($M_{S}$), the exchange stiffness ($A$), and
gyromagnetic ratio ($\gamma$). In this work, the static magnetic field is
applied at angle $\theta=$ 53∘ _w.r.t._ to the film plane, with its in-pane
component parallel to the antennas, _i.e._ perpendicular to the spinwave
propagation along the waveguide. In this configuration, the propagating waves
are of a mixed type between magnetostatic surface spin waves (MSSW) and
forward volume spin waves (FVMSW). The frequency in this configuration is
described by the dispersion relation (?, ?):

$\displaystyle\omega_{0}$
$\displaystyle=\gamma\mu_{0}\sqrt{(H_{int}+M_{s}l_{ex}^{2}k^{2})(H_{int}+M_{s}l_{ex}^{2}k^{2}+M_{s}F_{0})},$
(3) $\displaystyle F_{0}$
$\displaystyle=P_{0}+\cos^{2}\theta_{int}\left[1-P_{0}\right.+M_{s}\left.P_{0}(1-P_{0})/\left(H_{int}+M_{s}l_{ex}^{2}k^{2}\right)\right],$
(4)

where $k$ is the full wavevector of the spin wave, $P_{0}=1-(1-\exp(-kd))$ and
$d=$ 5 $\mu$m being the thickness of the film,
$l_{ex}=\sqrt{2A/(\mu_{0}M_{s}^{2})}\simeq 1.7\cdot 10^{-8}$ m for YIG, and
$H_{int}$ and $\theta_{int}$ define the value and out of plane angle of the
internal field; these can be found from the strength and angle of the applied
field ($H$ and $\theta$) using the solution of the magnetostatic problem (?,
?). In our case $\mathit{H}=$ 0.04 T and the choice of $\theta=53^{\circ}$ is
motivated by an empirically found high excitation efficiency and low total
losses in the YIG delay line when using simple thin copper wire antennas.
Moreover, spinwaves in such a configuration are strongly unidirectional and
will be excited most efficiently in only one direction, which prohibits multi-
transit spinwave echo signals and, therefore, improves the frequency stability
of the SWIM.

Our time-multiplexed SWIM can be considered as a ring oscillator circuit.
According to the circuit design theory, a circuit oscillates when the
Barkhausen stability criteria are satisfied:

$\displaystyle\beta\mathbf{A}$ $\displaystyle=1,$ (5)
$\displaystyle\angle\beta\mathbf{A}$ $\displaystyle=2\pi
n,\>n\in\\{0,1,2...\\}$ (6)

where $\mathbf{A}$ is the total amplification in the loop, $\beta$ is the
total loss, and $\angle\beta\mathbf{A}$ is the phase accumulated along the
loop. The Barkhausen criteria are further narrowed by adding a parametric
phase-sensitive amplifier (PSA), which limits stable oscillations only at
either phase 0 or $\pi$ relative to a pumping reference signal at twice the
oscillating frequency and consequently phase binarize the system. The PSA also
induces second-harmonic frequency locking (?, ?) to the external reference
signal,

$\displaystyle\omega_{0}t-\omega_{ref}t/2$ $\displaystyle={\pi
n},\>n\in\\{0,1,2...\\}$ (7) $\displaystyle\omega_{0}$
$\displaystyle=\omega_{ref}/2$ (8)

which further improves the frequency stability.

These phase-binarized oscillation conditions are valid for continuous
oscillations in the loop. In order to define separate time-multiplexed and
phase-binarized aritificial Ising spins, we introduce an RF switch ($\bf{1}$)
that is triggered by a square wave signal with variable frequency to control
the spinwave RF pulse length, $\tau_{p}$. The frequency of the switching is
determined by the total delay in the ring and the required number of supported
artificial spins:

$f_{sw}=N/{\tau}_{delay}$ (9)

The RF switch also prevents propagating spinwave RF pulses from spreading due
to spinwave dispersion and non-constant delay time over the occupied frequency
range.

Figure 1: The spinwave Ising machine (SWIM). (a) Schematic of the SWIM with
the in-plane magnetized YIG waveguide in the center together with peripheral
microwave components, described in the text, and a real-time oscilloscope for
direct studies of the time-evoluation of all spins; couplings between spins
are implemented using three delay cables of different lengths. (b) Measurement
of $S_{21}$ of the YIG delay line showing the spinwave spectrum. BW${}_{sw}=$
60 MHz is the spinwave spectral bandwidth measured at –3 dB. (c) The
frequency-dependent delay time of the YIG delay line. $\tau_{delay}=$ 271 ns
is the average delay time over BWsw with a total deviation of
$\Delta\tau_{delay}=$ 30 ns. (d) Amplification of the parametric phase-
sensitive amplifier as a function of the relative phase between the
propagating RF pulses and the reference signal $\omega_{ref}=$ 6.25 GHz, with
a phase sensitivity of $\mathit{A_{PSA}}=$ 6.1 dB.

One major advantage of using spinwaves is their 5–7 orders of magnitude slower
propagation than the vacuum speed of light, which allows us to reduce the
length of the YIG waveguide to millimeters compared to the kilometers-long
optical fibers of CIM. As our artificial Ising spins consist of spinwave RF
pulses of length $\tau_{p}$, the maximum number of supported spins is
$N=\tau_{delay}/\tau_{p}$, where the waveguide delay time
$\tau_{delay}=\mathit{l_{YIG}}/v_{g,sw}$ is the ratio of the waveguide length
($\mathit{l_{YIG}}$) and the spinwave group velocity ($v_{g,sw}$). An
approximate theoretical value of the minimal pulse length is given by
convolution of _i_) the 3-dB bandwidth $BW_{sw}$ of the spinwave spectrum,
which controls the slew rate of the spinwave RF pulses, and _ii_) the delay
time deviation $\Delta{\tau}_{delay}$ within $BW_{sw}$, which broadens the
pulses as they propagate between the two transducers.

Fig. 1(b) shows the spinwave spectrum of the YIG spinwave waveguide in the
form of its $S_{21}$-parameter. The bandwidth of the spinwave generation
spectrum is BW${}_{sw}=$ 60 MHz measured at the –3 dB level. Fig. 1(c) shows
the corresponding frequency-dependent delay time, with an average value across
BWsw of ${\tau}_{delay}=$ 271 ns and a deviation of 30.2 ns.The limit of
minimal pulse width derived from bandwidth is 16.8 ns. However, the delay time
deviation causes a significant broadening of propagating spinwave RF pulses
limiting the total minimal pulse width by 34 ns.

## Defining and solving MAX-CUT problems

In this section, we evaluate the computational performance of the SWIM on
simple 4-spin and 8-spin MAX-CUT optimization problems. While we have
successfully populated the waveguide with up to 10 spins, in agreement with
the theoretical estimate above, we here limit the problem size to eight spins
to avoid any undesired direct interaction between spinwave RF pulses as they
begin to partially overlap when more closely packed.

The first experiment demonstrates the principle of operation and routes for
obtaining a solution to the simple 4-spin MAX-CUT problem with
antiferromagnetic nearest-neighbor couplings described by the following Ising
matrix:

Figure 2: Experimental demonstration of 4-spin SWIM solving a MAX-CUT problem.
(a) Transformation of the initial non-optimal state to an optimal solution for
4-spin MAX-CUT optimization problem. The number of cuts in each case is
defined by the number of couplings which cross the line which separates spins
if they were regrouped according to their spin values. Time traces of SWIM
control and computational signals in the initial non-optimal state (b) and the
final solution (c). Top panel: a control signal of switch 1 with a repetition
frequency of 14.8 MHz and duty cycle of 50$\%$ that forms oscillations in the
SWIM ring oscillator into separate propagating RF pulses. Middle panel: an
amplified signal of RF propagating pulses $V_{ampl}$. Bottom panel: a
calculated instantaneous phase signal. The time traces in middle panels are
colored according to a value of the calculated instantaneous signal to
visually show the pulse instantaneous phase.

$J_{ij}=\left(\begin{array}[]{cccc}0&-1&0&0\\\ 0&0&-1&0\\\ 0&0&0&-1\\\
-1&0&0&0\end{array}\right)$ (10)

This matrix can be realized using a single coupling delay cable (3) with a
delay of 68 ns, _i.e._ of exactly one period of the pulse repetition time.
Each microwave RF pulse that passes through Coupler 2 generates both a
spinwave RF pulse in the YIG delay line and a coupling microwave RF pulse
propagating through the coupling delay cable 1. After 68 ns, the coupling
microwave RF pulse reaches Coupler 3 and combines with a different RF
microwave RF pulse converted from a propagating spinwave RF pulse, effectively
implementing the coupling between these two spins. The negative sign of the
coupling coefficients $J_{ij}$ is realized using an additional $180^{\circ}$
degree of phase shift implemented with a variable phase shifter.

For a demonstration of the computational dynamics, the SWIM is first placed in
a randomly chosen steady state $s_{j}^{1}=\\{+1-1-1-1\\}$ with the coupling
cable $\bf{3}$ disconnected. The value +1 corresponds to $180^{\circ}$ of
phase difference between the RF pulse and the reference signal
$\mathit{f_{ref}}$, while the value -1 corresponds to $0^{\circ}$ phase
difference. The time traces of the signal of propagating RF pulses and their
instantaneous phase are shown in Fig. 2(b), where the color of the time trace
corresponds to the instantaneous phase. As can be seen, the propagating
spinwave RF pulses are well-defined and separated from each other; their
individual phases are similarly well-defined and uniform in time. The initial
$s_{j}^{1}$ represents a non-optimal MAX-CUT solution with the number of cuts
equal to 2 (see Fig. 2(a)&(b)). At $t=0$ we connect the coupling cable
$\bf{3}$ and let the SWIM evolve for 3.76 $\mu$s or 12 circulation periods,
during which it reaches a stable state and finds the optimal number of cuts.
Fig. 2(c) shows the resulting spin state during the 12th circulation period,
where the third spin has now switched its phase to $\pi$ under the influence
of the coupling matrix. The SWIM has hence changed its state to
$s_{j}^{12}=\\{+1-1+1-1\\}$, which represents an optimal solution with the
number of cuts equal to 4 (see Fig. 2(a)).

Similarly, we confirmed the SWIM’s capability to solve 8-spin MAX-CUT
problems. The frequency $f_{sw}$ of RF switch $\bf{1}$ was hence increased to
29.6 MHz according to eq.9 to increase the spin capacity to eight.

As the corresponding repetition period decreased to 34 ns, a shorter coupling
delay cable with 34 ns delay was used for this experiment.

Figure 3: Experimental demonstration of 8-spin MAX-CUT problem computation.
(a) An optimal solution for 8-spin MAX-CUT optimization problem with 8 cuts.
(b) Time traces of SWIM computational signals in the final optimal solution
state. Top panel: an amplified signal of RF propagating pulses $V_{ampl}$. The
total duration of the time trace signals corresponds to the delay time in YIG
delay line. Bottom panel: a calculated instantaneous phase signal

As can be seen in Fig. 2(b) the spinwave RF pulses start to show some minor
overlap, which could nevertheless be neglected, as the SWIM had reached the
correct solution when interrogated after 12 cycles.

Figure 4: Spin switching scenarios for 4- and 8-spin MAX-CUT problem
computations. (a,b) Evolution of the 3rd artificial spin from an initial non-
optimal state to an optimal solution for 4-spin and 8-spin MAX-CUT
optimization problem. Top panels: an amplified signal $V_{ampl}$ within time
intervals which correspond to the 3rd propagating RF pulse. Bottom panels: a
calculated instantaneous phase signal. (c,d) Envelopes (left panels) and
instantaneous phase signals (right panels) of the 3rd propagating RF pulse
signal within 12 circulation periods plotted in the form of overlapping traces
in a relative scale. The color of the time traces corresponds to the
circulation number.

To further visualize the routes to a solution in greater detail, time traces
and the instantaneous phases for the 3rd spin for both 4- and 8-spin MAX-CUT
problems were captured at twelve consecutive circulation steps and compiled
into Fig. 4(a,b) with scale brakes. The signal envelopes and instantaneous
phases are also compiled in Fig. 4(c,d) in the form of overlapping traces with
different colors corresponding to the number of circulations. Interestingly,
the SWIM demonstrates two different scenarios of switching between
$c_{3}^{1}=-1$ and $c_{3}^{12}=+1$. The first scenario is observed for the
4-spin MAX-CUT problem where the pulses are longer (40 ns). During the first 4
circulations, the spinwave RF pulse gradually separates into two distinct
domains of approximately equal width, where only the trailing domain reverses
its instantaneous phase. Withing the next 4-5 circulations the amplitude of
the switching domain starts to decrease. The corresponding domain wall can
also be clearly seen in the instantaneous phase. From this point on, the
domain wall starts to propagate through the pulse until a fully switched
single domain pulse is established at about the 10th cycle. The exact origin
of this behavior is yet to be identified.

A different switching mechanism can be observed for the 8-spin MAX-CUT problem
with the 34 ns long coupling delay cable $\bf{3}$. In this case, the evolution
of the instantaneous phase occurs via a gradual and uniform shift from 0 to
$\pi$. The difference between the two mechanisms is particularly clear when
the evolutions of the envelope and instantaneous phase are plotted in
Fig.4c&d. It is possible that the shorter pulse length in the 8-spin case
promotes a more uniform response to the coupling. Regardless of its origin,
this gradual uniform switching mechanism might open up the possibility to use
the SWIM concept as a potential platform to implement XY-Ising machines (?).

A key parameter for benchmarking the performance of an Ising machine is its
time-to-solution. In order to achieve optimal solutions in the simple cases of
4- and 8-spin MAX-CUT problems, a single system run appears sufficient.
However, for larger combinatorial problems with dense coupling graphs, there
may exist a number of local minima in the Ising Hamiltonian. Therefore, a
larger SWIM will require classical annealing schemes with multiple system runs
and randomly chosen or idle initial conditions in order to perform statistical
analysis on the probability of different stable solutions. The shorter the
time-to-solution, the more such runs can be carried out in a given time.

In our experiments, we periodically switch off the amplification in the loop
for 2 $\mu$s to ensure that the SWIM returns to an idle state, after which we
again switch on the loop amplification. Fig. 4(a) shows the temporal evolution
of the SWIM prior to and after shutdown. It first takes 1.9 $\mu$s for the
SWIM to develop circulating RF pulses to saturation level. It then takes an
additional 1.5 $\mu$s to reach a stable state that represents the optimal
solution. The rapid development of the amplitude of the circulating pulses is
explained by a large small-signal over-amplification in the loop of the ring-
based SWIM circuit ($\sim$6 dB) and leads to system evolution into a non-
optimal solution before the signal reaches amplitude saturation. Reducing this
over-amplification will lead to slower development of oscillation and it might
improve the time-to-solution for large size optimization problems as at low
amplitude of the RF pulses the SNR is lower and the system would more easily
switch between its states. Another approach that could be attempted, similar
to CIM (?), is to gradually increase the amplification ratio until the system
starts oscillating in the lowest energy state that corresponds to the minimum
state of the Ising Hamiltonian.

Figure 5: Measurement of the time-to-solution parameter for 4-spin MAX-CUT
problem computation. (a) SWIM time traces for the measurement of time-to-
saturation and time-to-solution parameters. Top panel: a control signal
$V_{ctrl\>ampl}$ for the linear amplifier. Middle panel: an amplified signal
of RF propagating pulses $V_{ampl}$. Bottom panel: phase values in the center
of each propagating RF pulse sampled at each circulation period. The linear
amplifier switches down at the moment $1\>\mu s$ with blackout period of
$2\>\mu s$ to ensure signal suppression. (b) Circulation number and time-to-
solution parameters for 4-spin MAX-CUT problem as a function of the coupling
strength.

The SWIM design allows for the control of the overall coupling strength
between spins by attenuating the signal coming from the coupling delay lines
with a variable attenuator (see Fig. 1(a)). To benchmark the SWIM dynamic
range, a series of ten consecutive time-to-solution measurements was conducted
at different coupling strengths for the 4-spin MAX-CUT problem. A statistical
analysis is presented in Fig. 5(b). For this very simple problem, the SWIM
demonstrates almost constant time-to-solution values for coupling strengths
from 0 to –12 dB with about 20 circulations or 5 $\mu$s. The very weak trend
of a longer time-to-solution at weaker coupling strength is then more
pronounced at –15 dB. The large SWIM dynamical range of 15 dB should allow the
mapping of a wide range of optimization problems.

## Discussion and outlook

The SWIM concept has a high potential for further scaling in terms of spin
capacity and physical size and the improvement of power consumption and speed
performance parameters. Exploitation of non-linear spinwave solitons (?, ?, ?,
?) and proper choice of magnetization angle could compensate spinwave
dispersion and flatten delay time over a wide range of frequencies allowing to
use shorter spinwave RF pulses effectively increasing the number of supported
Ising spins with the same length of YIG waveguide. The use of exchange-
dominated spinwaves with exceptionally slow propagation velocity (?, ?, ?, ?)
in nanoscale YIG films could help to further miniaturize the YIG waveguide.

The use of miniaturized microwave signal processing components opens a
possibility of SWIM scaling by parallelizing the system using multiple ring
circuits comprising identical phase sensitive amplifiers, YIG delay lines,
etc. Parallelizing SWIM would improve the time-to-solution parameter by the
number of the parallel rings as it decreases circulation time. Such a solution
would require an FPGA (?) as a highly parallel control system for the
interconnection of Ising spins.

## Conclusion

We have demonstrated a spinwave time-multiplexed Ising machine (SWIM) and
characterized its computational performance and functionality. We have shown
how the SWIM can solve 4- and 8-spin NP-hard MAX-CUT problems in about 3.5
$\mu$s, which is comparable to D-wave’s performance and faster than CIM. The
multiphysics design allows the use of power-efficient and conventional low-
power microwave components consuming 2 W of power, which amounts to an energy
consumption of only 7 $\mu$J. This outperforms both D-wave’s Advantage and CIM
by 3–5 orders of magnitude. Thanks to the exceptionally slow spin wave
propagation speed, the 7 mm long YIG waveguide can host up to 11 spins and
further scaling is envisioned using both slower spinwave modes and patterned
YIG waveguides. The vast variety of nonlinear spinwave modes and phenomena,
and the direct availability of more optimized microwave signal processing
components, bodes well for the SWIM scaling potential in terms of number of
supported Ising spins, device size, and power consumption, making SWIM a
commercially feasible platform for solving a wide range of optimization
problems.

## References

  * 1. K. Tatsumura, R. Hidaka, M. Yamasaki, Y. Sakai, H. Goto, 2020 IEEE International Symposium on Circuits and Systems (ISCAS) (IEEE, 2020), pp. 1–5.
  * 2. T. Zhang, Q. Tao, J. Han, 2021 18th International SoC Design Conference (ISOCC) (IEEE, 2021), pp. 288–289.
  * 3. M. Hasegawa, H. Ito, H. Takesue, K. Aihara, IEICE Transactions on Communications 104, 210 (2021).
  * 4. R. Liu, X. Li, K. S. Lam, Current opinion in chemical biology 38, 117 (2017).
  * 5. F. Barahona, Journal of Physics A: Mathematical and General 15, 3241 (1982).
  * 6. A. Lucas, Frontiers in Physics 2 (2014).
  * 7. N. Mohseni, P. L. McMahon, T. Byrnes, Nature Reviews Physics 4, 363 (2022).
  * 8. N. Mohseni, P. L. McMahon, T. Byrnes, Nature Reviews Physics 4, 363 (2022).
  * 9. D. I. Albertsson, et al., Applied Physics Letters 118, 112404 (2021).
  * 10. A. Houshang, et al., Physical Review Applied 17, 014003 (2022).
  * 11. B. C. McGoldrick, J. Z. Sun, L. Liu, Physical Review Applied 17, 014006 (2022).
  * 12. B. Sutton, K. Y. Camsari, B. Behin-Aein, S. Datta, Scientific Reports 7, 44370 (2017).
  * 13. F. Cai, et al., Nature Electronics 3, 409 (2020).
  * 14. M. W. Johnson, et al., Nature 473, 194 (2011).
  * 15. P. I. Bunyk, et al., IEEE Transactions on Applied Superconductivity 24, 1 (2014).
  * 16. S. Dutta, et al., 2019 IEEE International Electron Devices Meeting (IEDM) (2019), pp. 37.8.1–37.8.4. ISSN: 2156-017X.
  * 17. T. Honjo, et al., Science Advances 7, eabh0952 (2021).
  * 18. T. Inagaki, et al., Science 354, 603 (2016).
  * 19. P. L. McMahon, et al., Science 354, 614 (2016).
  * 20. T. Wang, L. Wu, J. Roychowdhury, Proceedings of the 56th Annual Design Automation Conference 2019, DAC ’19 (Association for Computing Machinery, New York, NY, USA, 2019).
  * 21. S. Patel, L. Chen, P. Canoza, S. Salahuddin, arXiv preprint arXiv:2008.04436 (2020).
  * 22. K. Tatsumura, M. Yamasaki, H. Goto, Nature Electronics 4, 208 (2021).
  * 23. K. Takata, et al., Scientific reports 6, 1 (2016).
  * 24. Y. Yamamoto, et al., npj Quantum Information 3, 1 (2017).
  * 25. Y. Haribara, S. Utsunomiya, Y. Yamamoto, A Coherent Ising Machine for MAX-CUT Problems: Performance Evaluation against Semidefinite Programming and Simulated Annealing (Springer Japan, Tokyo, 2016), pp. 251–262.
  * 26. A. Marandi, Z. Wang, K. Takata, R. L. Byer, Y. Yamamoto, Nature Photonics 8, 937 (2014).
  * 27. S. Kako, et al., Advanced Quantum Technologies 3, 2000045 (2020).
  * 28. Á. Papp, W. Porod, Á. I. Csurgay, G. Csaba, Scientific reports 7, 1 (2017).
  * 29. A. B. Ustinov, A. V. Drozdovskii, B. A. Kalinikos, Applied physics letters 96, 142513 (2010).
  * 30. A. B. Ustinov, B. A. Kalinikos, Applied Physics Letters 93, 102504 (2008).
  * 31. A. B. Ustinov, B. A. Kalinikos, E. Lähderanta, Journal of Applied Physics 113, 113904 (2013).
  * 32. A. Sadovnikov, V. Gubanov, S. Sheshukova, Y. P. Sharaevskii, S. Nikitov, Physical Review Applied 9, 051002 (2018).
  * 33. A. A. Grachev, A. V. Sadovnikov, S. A. Nikitov, Nanomaterials 12, 1520 (2022).
  * 34. I. Ustinova, et al., Technical Physics Letters 42, 891 (2016).
  * 35. Y. L. Etko, A. Ustinov, Technical Physics Letters 37, 1015 (2011).
  * 36. A. Litvinenko, S. Grishin, Y. P. Sharaevskii, V. Tikhonov, S. Nikitov, Technical Physics Letters 44, 263 (2018).
  * 37. A. Litvinenko, et al., Physical Review Applied 15, 034057 (2021).
  * 38. V. V. Vitko, A. A. Nikitin, A. B. Ustinov, B. A. Kalinikos, IEEE Magnetics Letters 9, 1 (2018).
  * 39. A. B. Ustinov, E. Lähderanta, M. Inoue, B. A. Kalinikos, IEEE Magnetics Letters 10, 1 (2019).
  * 40. A. A. Nikitin, et al., Applied Physics Letters 106, 102405 (2015).
  * 41. Q. Wang, et al., Nature Electronics 3, 765 (2020).
  * 42. A. Khitun, M. Bao, K. L. Wang, Journal of Physics D: Applied Physics 43, 264005 (2010).
  * 43. A. Khitun, K. L. Wang, Journal of Applied Physics 110, 034306 (2011).
  * 44. A. Khitun, M. Bao, K. L. Wang, IEEE Transactions on Magnetics 44, 2141 (2008).
  * 45. C. Davies, et al., IEEE Transactions on Magnetics 51, 1 (2015).
  * 46. S. Watt, M. Kostylev, A. B. Ustinov, Journal of Applied Physics 129, 044902 (2021).
  * 47. S. Watt, M. Kostylev, A. B. Ustinov, B. A. Kalinikos, Physical Review Applied 15, 064060 (2021).
  * 48. M. Balinskiy, A. Khitun, AIP Advances 12, 045307 (2022).
  * 49. Y. Khivintsev, et al., Journal of Magnetism and Magnetic Materials 545, 168754 (2022).
  * 50. Z. Guo, et al., Proceedings of the IEEE (2021).
  * 51. A. Serga, A. Chumak, B. Hillebrands, Journal of Physics D: Applied Physics 43, 264002 (2010).
  * 52. V. Cherepanov, I. Kolokolov, V. L’vov, Physics Reports 229, 81 (1993).
  * 53. H. Glass, Proceedings of the IEEE 76, 151 (1988).
  * 54. M. Haidar, et al., Journal of Applied Physics 117, 17D119 (2015).
  * 55. A. Navabi, et al., Physical review applied 11, 034046 (2019).
  * 56. B. Kalinikos, A. Slavin, Journal of Physics C: Solid State Physics 19, 7013 (1986).
  * 57. S. Muralidhar, et al., Physical Review Letters 126, 037204 (2021).
  * 58. R. Lebrun, et al., Phys. Rev. Lett. 115, 017201 (2015).
  * 59. A. Litvinenko, et al., Physical Review Applied 16, 024048 (2021).
  * 60. K. P. Kalinin, A. Amo, J. Bloch, N. G. Berloff, Nanophotonics 9, 4127 (2020).
  * 61. A. B. Ustinov, N. Y. Grigor’eva, B. A. Kalinikos, JETP letters 88, 31 (2008).
  * 62. A. B. Ustinov, B. A. Kalinikos, V. E. Demidov, S. O. Demokritov, Physical Review B 80, 052405 (2009).
  * 63. P. A. Kolodin, et al., Phys. Rev. Lett. 80, 1976 (1998).
  * 64. M. Kostylev, B. Kalinikos, Technical Physics 45, 277 (2000).
  * 65. A. Chumak, A. Serga, B. Hillebrands, Journal of Physics D: Applied Physics 50, 244001 (2017).
  * 66. A. Mahmoud, et al., Journal of Applied Physics 128, 161101 (2020).
  * 67. V. Tikhonov, A. Litvinenko, Applied Physics Letters 115, 072410 (2019).

## Author Contributions

A.L. conceived the concept and designed the circuit; A.L., V.G., and A.A.
performed the measurements and analyzed the data; A.L, R.K. performed
theoretical calculations with help from A.S. and V.T.; J.Å. managed the
project; all co-authors contributed to the manuscript, the discussion, and
analysis of the results.