An overview of low-rank matrix recovery from incomplete observations

Introduction

Low-rank matrices arise in an incredibly wide range of settings throughout science and applied mathematics. To name just a few examples, we commonly encounter low-rank matrices in contexts as varied as:

ensembles of signals: the output of a sensor array or network, a collection of video frames, or a sequence of segments of a longer signal can often be highly correlated and represented using a low-rank matrix ;

system identification: low-rank (Hankel) matrices correspond to low-order linear, time-invariant systems ;

adjacency matrices: the connectivity structure of many graphs, such as those that arise in manifold learning and social networks, is often low rank ;

distance matrices: in many data embedding problems — such as those that arise in the context of multidimensional scaling , sensor localization , nuclear magnetic resonance spectroscopy , and others — the matrix of pairwise distances will typically have a rank dependent on the (low) dimension of the space in which the data lies;

item response data: low-rank models are frequently used in analyzing data sets containing the responses of various individuals to a range of items, such as survey data , educational data , the data generated by recommendation systems , and others;

machine learning: low-rank models are ubiquitous in machine learning, laying the foundation for both classical techniques such as principal component analysis as well as modern approaches to multi-task learning and natural language processing ;

quantum state tomography: a pure quantum state of $N$ ions can be described by a $2^{N}\times 2^{N}$ matrix with rank one .

In all of these settings, the matrices we are ultimately interested in can be extremely large. Moreover, as data becomes increasingly cheap to acquire, the potential size will continue to grow. This raises a number of challenges, but often a key obstacle is that fully observing the matrix of interest can prove to be an impossible task: it can be prohibitively expensive to fully sample the entire output of a sensor array; we might only be able to measure the strength of a few connections in a graph; and any particular user of a recommendation system will provide only a few ratings. In such settings we are left with a highly incomplete set of observations, and unfortunately, many of the most popular approaches to processing the data in the applications where low-rank matrices arise assume that we have a fully-sampled data set and are generally not robust to missing/incomplete data. In these situations we are confronted with the inverse problem of recovering the full matrix from our incomplete observations.

While such recovery is not always possible in general, when the matrix is low rank, it is possible to exploit this structure and to perform this kind of recovery in a surprisingly efficient manner. In fact, in recent years there has been tremendous progress in our understanding of how to solve such problems. While many of these applications have a relatively long history in which various existing approaches to dealing with incomplete observations have been independently developed, recent advances in the closely related field of compressive sensing have enabled a burst of progress in the last few years. We now have a unified framework which provides a strong base of theoretical results concerning when it is possible to recover a low-rank matrix from incomplete observations using efficient, practical algorithms .

In this paper we provide a survey of this developing field. We begin with a more formal mathematical statement of the problem of low-rank matrix recovery in Section 2, followed in Section 3 by an overview of some of the algorithms most commonly used in practice to solve these kinds of problems. We provide a brief overview of the existing theoretical guarantees for these algorithms in Sections 4, 5, and 6 for several concrete observation models with an emphasis on how many observations are required to reliably recover a low-rank matrix and what additional assumptions are potentially required. Finally, in Section 7 we describe an important application of these techniques to an important class of problems where we can solve quadratic and bilinear systems of equations by re-casting them as a simple problem of low-rank matrix recovery.

The Matrix Recovery Problem

We begin by carefully stating what we mean by low-rank matrix recovery. We observe a matrix ${\mathbf{X}}_{0}$, which we will assume to have size $M\times N$ and which we can express either exactly or approximately as having rank $R\ll M,N$. This means that we can write

Our discussion below will at times involve the adjoint of this operator, which is defined as

Our survey focuses on three basic variations of this measurement model. The first is taking ${\mathcal{A}}$ to be a random projection, where each of the $\boldsymbol{A}_{i}$ consisting of independent and identically distributed random variables. Although this model arises in only a limited number of practical situations, the theory is so streamlined that it can be understood almost from first principles (see Section 4). For our second model, ${\mathcal{A}}$ returns a subset of the entries of the target. Recovering from these samples is known as the matrix completion problem. In this case, each of the $\boldsymbol{A}_{i}$ has exactly one non-zero entry. The analysis for this problem, which we overview in Section 5, can also be extended to observing a subset of the expansion coefficients of ${\mathbf{X}}_{0}$ in a fixed (and known) orthobasis. The third model, which we will discuss in Section 7, is that the $\boldsymbol{A}_{i}$ are rank-1 matrices. These are encountered when each observation can be written as a quadratic or bilinear form in ${\mathbf{X}}_{0}$.

While these are the measurement models that have received the most attention in the literature, they are by no means the only interesting models. Other models inspired by applications in imaging and signal processing have also appeared recently in the literature (see for example ).

Algorithms for Matrix Recovery

We start by reviewing the classical problem of finding the best low-rank approximation to a given $M\times N$ matrix $\boldsymbol{X}_{0}$. By “best”, we mean closest in the sum-of-squares sense, and we formulate the problem as

where $\|\boldsymbol{X}-\boldsymbol{X}_{0}\|_{F}^{2}=\sum_{m,n}(X[m,n]-X_{0}[m,n])^{2}$ is the square of the standard Frobenius norm, and $R$ is the desired rank of the approximation. This problem is nonconvex, but is actually easy to solve explicitly using the singular value decomposition (SVD). In particular, if we decompose $\boldsymbol{X}_{0}$ as

where $K=\min(M,N)$, $\boldsymbol{U},\boldsymbol{V}$ are $M\times K$ and $N\times K$ matrices with orthonormal columns $\boldsymbol{u}_{1},\ldots,\boldsymbol{u}_{K}$ and $\boldsymbol{v}_{1},\ldots,\boldsymbol{v}_{K}$, and $\boldsymbol{\Sigma}$ is a diagonal $K\times K$ matrix with sorted entries $\sigma_{1}\geq\sigma_{2}\geq\cdots\geq\sigma_{K}\geq 0$, then the solution to (2) is found simply by truncating this expansion:

This is known as the Eckart-Young theorem; see [69, Chapter 7] for a detailed proof and discussion. For medium scale problems, computing the SVD to high precision is tractable, with computational complexity scaling as $O(K^{2}\max(M,N))$.

Computing the best low-rank approximation, then, is akin to thresholding the singular values: we take the matrix, compute its SVD, keep the large singular values while killing off the small ones, and then reconstruct. A variation of the program above makes this connection clearer. The Lagrangian of (2) is

As we vary the parameter $\lambda$ above, the solution to the program changes — in fact, the set of solutions produced for different $0\leq\lambda<\infty$ is exactly the same as the set of solutions for (2) produced for all $1\leq R\leq K$. Given $\lambda$, we solve (3) by computing the SVD, hard thresholding the singular values via

A common variation of the algorithm above involves replacing the hard threshold in (4) with a soft threshold. In this case we still set the singular values that are small to zero, but now the large values are shrunk:

This amounts to a more gradual phasing out of the terms that just cross the threshold. It turns out that this soft thresholding process can also be put in variational form; when $\gamma=\lambda/2$ the result of the procedure above is the solution to

where $\|\boldsymbol{X}\|_{*}$ is the nuclear norm, and is equal to the sum of the singular values of $\boldsymbol{X}$. $\|\boldsymbol{X}\|_{*}$ is also known as the trace norm, as it is equal to the trace when $\boldsymbol{X}$ is symmetric positive semidefinite. Unlike the rank, $\|\boldsymbol{X}\|_{*}$ is a convex function, and often appears as a convex proxy for rank in optimization problems . While in the approximation problem we are considering here, the solutions to (3) and (6) involve very similar computations, this will not be the case at all when we consider recovery from partial observations in the next section.

In the procedures above, the computational cost is dominated by computing the SVD. For matrices with $M,N$ on the order of $100$–$1000$, there are a number of exact methods with similar computational complexity — see [68, Chap. 45] for an overview. When $\boldsymbol{X}_{0}$ is large but is very well approximated by a matrix with modest rank, randomized algorithms can be used to compute an approximate SVD .

2 Low-rank recovery and nuclear norm minimization

The low-rank approximation problem described above readily admits a straightforward solution. However, in this survey we are concerned instead with the low-rank recovery problem where we are working from (possibly noisy) indirect observations, $\boldsymbol{y}\approx{\mathcal{A}}(\boldsymbol{X}_{0})$. In this case we would ideally like to solve the analog of (2),

Unfortunately, whereas (2) could be solved via a simple SVD, (3.2) is in general NP-hard.

In contrast, the nuclear norm minimization program remains tractable. In particular, with indirect observations, we replace (6) with

This is an unconstrained convex optimization program, and can be solved in a systematic way using a proximal algorithm . The solution(s) to (7) will obey, for any $\gamma>0$, the fixed point condition

One class of methods for solving (7) are based on the iteration

for some appropriately chosen sequence $\{\gamma_{k}\}$ . As discussed in the previous section, the subproblem (8) is solved by singular value soft-thresholding. Computing the SVD of the $\boldsymbol{X}_{k}-\gamma_{k}{\mathcal{A}}^{*}({\mathcal{A}}(\boldsymbol{X}_{k})-\boldsymbol{y})$ at each iteration is almost always the dominant cost, as it typically requires significantly more computation than applying ${\mathcal{A}}$ and ${\mathcal{A}}^{*}$.

State-of-the-art methods for solving (7) based on singular value thresholding are not too much more complicated. For example, the FISTA algorithm modifies the basic iteration in (9) through intelligent choices of the scaling coefficient $\gamma_{k}$ and by replacing $\boldsymbol{X}_{k}$ with a carefully chosen combination of $\boldsymbol{X}_{k}$ and $\boldsymbol{X}_{k-1}$. These small changes have almost no effect on the amount of computation done at each iteration, but they converge in significantly fewer iterations.

The simplicity of these proximal-type algorithms makes them very attractive for small to medium sized problems. However, as the number of rows and columns in the matrix gets to be several thousand, direct computation of the SVD becomes problematic. For specially structured ${\mathcal{A}}$, including the important case where ${\mathcal{A}}(\boldsymbol{X})$ returns a subset of the entries of $\boldsymbol{X}$, fast algorithms that take advantage of this structure to compute the SVD have been developed to solve (7) or closely related problems . In more general settings, techniques from randomized linear algebra have been applied to compute approximate SVDs .

The parameter $\lambda$ in (7) determines the trade-off between the fidelity of the solution to the measurements $\boldsymbol{y}$ and its conformance to the low-rank model. When we are very confident in the measurements, it might make sense to use them to define a set of linear equality constraints, solving

The output of this program will match (7) as $\lambda\rightarrow 0$. Some of the analytical results we review in Sections 4 and 5 reveal conditions under which (11) recovers a low-rank matrix ${\mathbf{X}}_{0}$ exactly given measurements $\boldsymbol{y}={\mathcal{A}}(\boldsymbol{X}_{0})$ as constraints.

3 Iterative hard thresholding

From the algorithmic point of view, iterative hard thresholding (IHT) algorithms are very similar to the proximal algorithms used to solve nuclear norm minimization in the previous section. However, when the target is very low rank, they tend to converge extremely quickly.

The basic iteration is as follows. From the current estimate $\boldsymbol{X}_{k}$, we first take a step in the direction of the gradient of $\|{\mathcal{A}}(\boldsymbol{X})-\boldsymbol{y}\|_{2}^{2}$, then project onto the set of rank $R$ matrices:

The $\operatorname{ProjectRank}_{R}$ operator computes the top $R$ left and right singular vectors and singular values — when $R$ is small compared to $M$ and $N$, this can be done in significantly less time than computing a full SVD , especially if the operator ${\mathcal{A}}$ and its adjoint ${\mathcal{A}}^{*}$ are structured in a such a way that there is a fast method for applying the matrices $\boldsymbol{Y}_{k+1}$ to a series of vectors. In these cases, the intermediate matrix $\boldsymbol{Y}_{k+1}$ is not computed explicitly, but each term in the first equation above can be handled efficiently in the SVD computation. The implementations of IHT in rely on existing software packages to do this.

In contrast to the nuclear norm minimization algorithm above, each of the iterates $\boldsymbol{X}_{k}$ IHT produces has a prescribed rank. The storage required for $\boldsymbol{X}_{k}$ is roughly $R(M+N)$, as opposed to the $MN$ required for a general $M\times N$ matrix. This difference is critical for large-scale applications.

4 Alternating projections

The alternating projections algorithm is another space efficient technique which stores the iterates in factored form. The algorithm is extraordinarily simple, and easy to interpret: looking for a $M\times N$ rank $R$ matrix that is consistent with $\boldsymbol{y}$:

is the same as looking for a $M\times R$ matrix $\boldsymbol{L}$ and a $N\times R$ matrix $\boldsymbol{R}$ whose product is consistent with $\boldsymbol{y}$:

This optimization problem is still non-convex, but with one of $\boldsymbol{L}$ or $\boldsymbol{R}$ fixed, it is a simple least-squares problem. This motivates the following iteration. Given current estimates $\boldsymbol{L}_{k}$,$\boldsymbol{R}_{k}$, we update using

Each step involves solving a linear system of equations with $RM$ or $RN$ variables for which we can draw on well-established algorithms in numerical linear algebra. Its simplicity and efficiency make it one of the most popular methods for large-scale matrix factorization , and it tends to outperform nuclear norm minimization , especially in cases where the rank $R$ is very small compared to $N,M$.

There are few general convergence guarantees for alternating projections, and the final solution tends to depend heavily on the initialization of $\boldsymbol{L}$ and $\boldsymbol{R}$. However, guarantees for the rate of convergence can be found in and recent work has provided some first theoretical results for conditions under which the iterations above converge to the true low-rank matrix (and a method for supplying a reliable starting point). These are discussed further in Sections 4 and 5 below.

Another advantage alternating projections is that the framework can be extended to handle structure on one or both of the factors $\boldsymbol{L}$ and $\boldsymbol{R}$. Typically, this means that the least-squares problems above are either regularized or constrained in a manner which encourages or enforces the desired structure. Extensions of the iterations above have been used successfully for problems including non-negative matrix factorization , sparse PCA, where we restrict the number of nonzero terms in $\boldsymbol{L},\boldsymbol{R}$ , and dictionary learning , where $\boldsymbol{L}$ is a well-conditioned matrix and $\boldsymbol{R}$ is sparse. Again, there are few strong theoretical guarantees for these algorithms, with notable exceptions in the recent works .

The optimization problems for alternating projections (12) and the Burer-Monteiro heurtistic (10) are similar (nonlinear) least-squares problems on the matrix factors $\boldsymbol{L},\boldsymbol{R}$. The algorithms used to solve them, though, have a distinct difference: instead of fixing for one factor and optimizing the other as in (13) above, solvers for (10) (e.g., ) typically take descent steps on $\boldsymbol{L}$ and $\boldsymbol{R}$ simultaneously. Convergence analysis for closely related local descent methods can be found in the recent works .

5 Other algorithms for matrix recovery

The methods above are by no means the only algorithms which have been proposed for low-rank matrix recovery. We close this section by briefly mentioning some other techniques.

Recent years have seen a renewed interest in Frank-Wolfe-type algorithms for minimizing norms defined by the convex hull of a set of atoms . These algorithms are of particular interest for minimizing the nuclear norm (where the atoms are rank-1 matrices), as they only require computing the leading singular vector at every iteration, rather than a full SVD . The nuclear norm problems in (7) and (11) can also be minimized by solving a series of weighted least-squares problems, each of which can be solved using standard techniques for linear systems of equations .

Beyond the nuclear norm, other proxies for rank exist. The max-norm is a convex function that results in a similar optimization program to (7), and is subject to similar heuristics as (10) for storage reduction . Alternatively, the logarithm of the determinant is a nonconvex, but smooth, proxy for rank when $\boldsymbol{X}$ is positive semidefinite (and can be applied to general matrices by embedding them in a PSD matrix). In practice, locally minimizing this function subject to convex constraints tends to produce low-rank solutions .

A greedy algorithm for low-rank recovery was presented in . This algorithm alternates between selecting an estimate for the $R$-dimensional subspace in which ${\mathbf{X}}_{0}$ lives, and projecting onto these subspaces; it is equipped with strong theoretical guarantees. Another class of algorithms evolves the left and right singular vectors along the Grassman manifold . These algorithms, which are specialized for the matrix completion problem, are computationally efficient, and have equally strong theoretical performance guarantees.

Matrix Recovery from Gaussian Observations

The theory of low-rank recovery is clean and elegant when the measurement operator ${\mathcal{A}}(\cdot)$ is a random projection. Applications where this is a good model for the observations are limited, but looking at it as an abstract problem gives us real mathematical insight about why low-rank recovery works.

The discussion in this section will center on “Gaussian” ${\mathcal{A}}(\cdot)$, where the entries of each of the $\boldsymbol{A}_{m}$ are independent and identically distributed normal random variables with zero mean and variance $L^{-1}$ — this variance is chosen so that $\operatorname{E}[\|{\mathcal{A}}(\boldsymbol{X})\|_{2}^{2}]=\|\boldsymbol{X}\|_{F}^{2}$ for any fixed $M\times N$ matrix $\boldsymbol{X}$.

We first examine the fundamental question of whether we can distinguish different rank-$R$ matrices viewed through the lens of the operator ${\mathcal{A}}$. One way to formalize this is asking whether ${\mathcal{A}}$ preserves the distances between all such matrices; this is certainly true if there exists a $0\leq\delta<1$ such that

for all $\boldsymbol{X}_{1},\boldsymbol{X}_{2}$ of rank $R$ or smaller. This condition is known as the matrix restricted isometry property (matrix-RIP), and is the matrix analog to the restricted isometry property from compressive sensing . The first immediate consequence is that all matrices of rank $R$ or less have unique images, since if ${\mathcal{A}}(\boldsymbol{X}_{1})={\mathcal{A}}(\boldsymbol{X}_{2})$ for $\boldsymbol{X}_{1}\not=\boldsymbol{X}_{2}$, then the lower bound would be violated. Qualitatively, the upper and lower bounds tell us that two rank $R$ matrices are as distinguishable from their measurements as they would be if they were observed directly. With (14) established, we can use any number of techniques to recover a low-rank matrix from measurements through ${\mathcal{A}}$; we discuss some specific guarantees below.

This can be established through relatively simple probabilistic methods. We sketch the argument below; the detailed proof in smartly combines ideas from and .

To begin, notice that it is enough to show that $(1-\delta)\leq\|{\mathcal{A}}(\boldsymbol{X})\|_{2}^{2}\leq(1+\delta)$ for all $\boldsymbol{X}$ of rank $2R$ and unit Frobenius norm. Then there are three basic steps for establishing (15).

For any arbitrary fixed $M\times N$ matrix $\boldsymbol{X}$ with $\|\boldsymbol{X}\|_{F}^{2}=1$,

for $t\leq\frac{1}{2}$, where $C$ and $c$ are reasonable constants that can be calculated explicitly (standard calculations yield $C\leq 2$ and $c\geq 1/8$). Since the entries of ${\mathcal{A}}(\boldsymbol{X})$ are independent Gaussian random variables, $\|{\mathcal{A}}(\boldsymbol{X})\|_{2}^{2}$ is a chi-squared random variable with $L$ degrees of freedom, and the inequality above follows from standard tail bounds .

We want (16) to hold not just for a single matrix, but uniformly over the set of all unit-norm rank $2R$ matrices. It is straightforward to get a uniform result over any finite set $\mathcal{Q}$ of such matrices by using a union bound:

Even though the set of all unit-norm rank $2R$ matrices is infinite, the next step shows us that it is enough to consider a finite subset.

Let $\mathcal{R}$ denote the set of $M\times N$ matrices $\boldsymbol{X}$ with $\text{rank}(\boldsymbol{X})\leq R$ and $\|\boldsymbol{X}\|_{F}=1$. Now let $\mathcal{R}_{\epsilon}$ be any finite $\epsilon$-approximation to $\mathcal{R}$ — this means that $\mathcal{R}_{\epsilon}\subset\mathcal{R}$ and for every $\boldsymbol{X}\in\mathcal{R}$ there is an $\bar{\boldsymbol{X}}$ in $\mathcal{R}_{\epsilon}$ that is within $\epsilon$: $\|\boldsymbol{X}-\bar{\boldsymbol{X}}\|_{F}\leq\epsilon$. Then a short argument shows that for $\epsilon=\delta/(4\sqrt{2})$,

Thus the supremum over the infinite set on the right can be replaced with the maximum over the finite set on the left. Using the result from step 1, we now have

The size of $\mathcal{R}_{\epsilon}$ can be estimated as a function of $\epsilon$. The bound in reads

The uniformity of the matrix-RIP, that it holds for all pairs of rank $R$ matrices simultaneously, results in stability guarantees for each of these algorithms when the measurements are made in the presence of noise, or the target matrix is only approximately low rank.

2 Convex geometry and Gaussian widths

Recent works analyze the performance of this program using very intuitive geometrical principles. The target matrix $\boldsymbol{X}_{0}$ is a member of two different convex sets; it is in the nuclear norm ball of radius $\|\boldsymbol{X}_{0}\|_{*}$,

and it is in the affine space consisting of all $M\times N$ matrices that have the same measurements,

By definition, $\boldsymbol{X}_{0}$ is the unique solution to (11) if and only if it is the only matrix in the intersection of $\mathcal{B}$ and $\mathcal{S}$; this is illustrated in Figure 1(a). Generating ${\mathcal{A}}$ with a Gaussian distribution is the same as choosing the orientation of $\mathcal{S}$ uniformly at random. The local geometry of the tip of $\mathcal{B}$ is determined by the rank of $\boldsymbol{X}_{0}$, smaller ranks make this point more singular, and hence decrease the probability of an intersection. The dimension of the set $\mathcal{S}$ is $MN-L$, the same as the null space of ${\mathcal{A}}$; as $L$ increases, $\mathcal{S}$ gets smaller, and the probability of an intersection decreases.

We can make both of these statements precise. The collection of directions that lead into the ball (i.e., decrease the nuclear norm) $\mathcal{B}$ from the point $\boldsymbol{X}_{0}$ is called the tangent cone of $\mathcal{B}$ at $\boldsymbol{X}_{0}$:

Asking if there is a better feasible point in (11) than $\boldsymbol{X}_{0}$ is exactly the same as asking if the subspace $\mathcal{S}-\boldsymbol{X}_{0}=\operatorname{Null}({\mathcal{A}})$ (the affine set $\mathcal{S}$ shifted to the origin) intersects the cone $\mathcal{T}_{\mathcal{B}}(\boldsymbol{X}_{0})$ at any place other than the origin.

When does a randomly chosen subspace intersect a cone only at the origin? A clean answer is given in , which builds on the classic work . This answer depends on the notion of the Gaussian width of a set $\mathcal{C}$, which is defined as

where $\boldsymbol{G}$ is a matrix whose entries are independent zero-mean Gaussian random variables with unit variance, and the expectation is taken with respect to $\boldsymbol{G}$. The quantity $w(\mathcal{C})$ can be interpreted as the amount we expect a set $\mathcal{C}$ to align with a randomly drawn vector. It is shown in that a randomly chosen subspace will intersect a cone $\mathcal{C}$ only at the origin with high probability if the codimension of the subspace is at least as large as the square of the Gaussian width:

It is clear that for subsets of $M\times N$ matrices, $\omega(\mathcal{C})\leq\sqrt{MN}$, and if $\mathcal{C}$ is a $D$ dimensional subspace, a quick calculation shows that $w(\mathcal{C})=\sqrt{D}$. The codimension of $\operatorname{Null}({\mathcal{A}})$ is always $L$, so measuring the Gaussian width of the tangent cone gives an immediate sufficient condition on the number of measurements needed for accurate recovery. The bound on the Gaussian width $w(\mathcal{T}_{\mathcal{B}}(\boldsymbol{X}_{0}))^{2}~{}\leq~{}3R(M+N-R)$ gives us the sharp sufficient condition of

The relation in (20) is very close to being necessary as well. In , it is shown that if the number of observations is not too far below $w(\mathcal{T}_{\mathcal{B}}(\boldsymbol{X}_{0}))^{2}$, then the probability that (11) fails to recover $\boldsymbol{X}_{0}$ is very close to $1$. These papers contain an impressive suite of numerical experiments showing that the success or failure of equality-constrained nuclear norm minimization (11) can be predicted accurately from the parameters $M,N,R$.

It should be mentioned that the analysis in and applies to many different types of structured recovery problems based on convex optimization; the low-rank recovery results discussed above are an important special case.

Matrix Completion

As we have just seen, the theory of low-rank matrix recovery is particularly elegant in the case where the measurement operator ${\mathcal{A}}(\cdot)$ is a (Gaussian) random projection operator. While this provides some insight into the kind of behavior we can hope for in many applications, it is also far from representative of the type of observations we often encounter in practice. In particular, in many settings of interest the $\boldsymbol{A}_{i}$ are highly structured — in the case of matrix completion, the $\boldsymbol{A}_{i}$ will have only a single nonzero value of 1 corresponding to the row and column of the observed element. An equivalent way to think about this type of measurement is that we only observe the entries of $\boldsymbol{X}_{0}$ on a subset $\Omega$ of the complete set of entries. This kind of measurement model arises in a variety of practical settings. For example, we might have a large number of questions we would like to potentially ask a number of users (as in a large-scale survey or recommendation system) but in practice we might only expect to receive responses to a few of these questions from any given user. Similarly, in many large-scale graphs (such as graphs representing the strength of social connections or the distances between sensors or other items) we might only be able to measure/observe the strength of a few connections in a graph.

Even if the elements of $\Omega$ are chosen at random, this scenario has some significant differences from the case where the $\boldsymbol{A}_{i}$ are Gaussian. It is clear that we cannot expect the theory developed using Gaussian widths to be of much use, since our observations are not Gaussian, but there is a more fundamental problem that arises in the case of matrix completion.

In particular, the most immediate challenge we encounter when developing a theory of matrix completion is that it is no longer possible to obtain the kind of uniform guarantees that apply to all low-rank matrices described in Section 4. To see why, consider a matrix with rank one, but where one (or both) singular vectors are sparse, meaning that their energy is concentrated on just a few entries. When this occurs, as illustrated in Figure 2, the resulting matrix will also have its energy mostly concentrated on just a few entries, in which case most entries give us very little information and the recovery problem is highly ill-posed unless nearly all the entries are observed. (Another way to see this is to realize that if only a few entries of such a matrix are observed, the matrix is very likely to live in the nullspace of the measurement operator ${\mathcal{A}}$.) More generally, if any particular column (or row) is approximately orthogonal to the span of the remaining columns (or rows), then it will be impossible to estimate without essentially observing the entire column (or row).

2 Recovery guarantees

There is now a rich literature providing a range of guarantees under which it is possible to recover a matrix $\boldsymbol{X}_{0}$ from randomly chosen entries under the assumption that $\boldsymbol{X}_{0}$ is incoherent and/or satisfies certain similar conditions. As a representative example, we will describe the guarantees that are possible when using nuclear norm minimization as a recovery technique, as first developed in and further refined in .

To state the main conclusion of this literature, we will assume for the moment that $\boldsymbol{X}_{0}$ is an $M\times N$ matrix of rank $R$ with singular value decomposition $\boldsymbol{U}\boldsymbol{\Sigma}\boldsymbol{V}^{*}$. We will also assume that $\mu(\boldsymbol{U}),\mu(\boldsymbol{V})\leq\mu_{0}$ and further that the matrix $\boldsymbol{E}=\boldsymbol{U}\boldsymbol{V}^{*}$ has a maximum entry bounded by $\mu_{1}\sqrt{R/MN}$ .As argued in , as a simple consequence of the Cauchy-Schwarz inequality, the assumption on $\boldsymbol{E}$ will always hold with $\mu_{1}=\mu_{0}\sqrt{R}.$ In addition, a more refined analysis using alternative concentration bounds can eliminate the need for any assumption on $\mu_{1}$ (see, for example, ). Moreover, given a limited amount of a priori information about the underlying matrix, it is also possible to obtain results that omit any dependence on the incoherence by sampling certain rows/columns more heavily . We will also assume, without loss of generality, that $M>N$. Then if $L$ entries of $\boldsymbol{X}_{0}$ are observed with locations sampled uniformly at random with

$\boldsymbol{X}_{0}$ will be the solution to (11) with high probability. We note that the required number of observations in the matrix completion setting exceeds that required in the Gaussian case (i.e., $O(R(M+N))$) in two natural ways. First, $L$ scales with the level of incoherence as quantified by $\mu_{0}$ and $\mu_{1}$. For example, in the case where $\mu(\boldsymbol{U})$ or $\mu(\boldsymbol{V})$ approach their maximal value then the bound in (22) reduces to the requirement that we observe nearly every entry. The second difference is that we have an additional $\log$ factor. While the power of $2$ on the $\log$ may not be strictly necessary, some logarithmic dependence on the dimension of the matrix is a necessary consequence of the random observation model. In particular, as a consequence of the classic coupon collector problem, we need at least $O(M\log M)$ observations simply to ensure that we observe each row at least once. (If $N>M$, the same argument applies with columns instead of rows.)

Finally, we also note that the above result applies to the specific case of noise-free observations of a matrix with rank at most $R$. Similar results can also be established that guarantee approximate recovery in the case of noisy observations and approximately low-rank matrices with the amount of recovery error being naturally determined by the amount of noise and degree of approximation error .

The full proof of the exact recovery result described above is somewhat involved, but has been significantly simplified in the work of and subsequently in to the point where the key ingredients are relatively straightforward. In light of the discussion of Section 4.1 one might expect that a possible avenue of attack would be to show that by selecting elements of our matrix at random, we obtain a measurement operator ${\mathcal{A}}$ that satisfies something like the matrix-RIP in (14) but which holds only for matrices satisfying our incoherence assumptions. In fact, it is indeed possible to pursue this route (and the incoherence assumption is vital to ensure that $\operatorname{E}[\|{\mathcal{A}}(\boldsymbol{X})\|_{2}^{2}]=\|\boldsymbol{X}\|_{F}^{2}$ as is required at the outset in Section 4). This is essentially the approach taken in . However, the difference between two incoherent matrices is not necessarily itself incoherent, which leads to some significant challenges in an RIP-based analysis. Moreover, this approach fails to yield the kind of exact recovery guarantee described above in the exactly low rank case and is primarily of interest in the noisy setting.

The more challenging part of the argument is the second step, in which one must show that given random samples of a matrix $\boldsymbol{X}_{0}$ satisfying the required coherence properties, such a dual certificate must exist (with high probability). The original proof of this in and the subsequent improvement in involved rather intricate analysis, made especially difficult by the fact that when we observe $L$ distinct entries of a matrix, our observations are not fully independent (since we cannot sample the same entry twice). In this analysis is dramatically simplified by two observations: (i) it is actually sufficient to merely obtain an approximate dual certificate, which can be easier to construct, and (ii) one can alternatively consider an observation model of sampling with replacement. Under this model, one can analyze the adjoint operator ${\mathcal{A}}^{*}(\cdot)$ by treating its output as a sum of independent random matrices, which enables the use of the powerful concentration inequalities recently developed in that bound the deviation of this sum from its expected value. This allows one to analyze the relatively simple “golfing scheme” of , which consists of an iterative construction of an approximate dual certificate. See for a condensed description of this approach and for an overview (and other applications) of the matrix concentration inequalities used in this analysis. It is also worth noting that this line of analysis also provides theory for alternative low-rank matrix recovery scenarios beyond simply the matrix completion case. Indeed, the arguments in were originally developed to address the problem of quantum state tomography, where the goal is to efficiently determine the quantum state of a system via a small number of observations . This problem can be posed as a matrix recovery problem, but where the $\boldsymbol{A}_{i}$ are constructed from Pauli matrices. Perhaps somewhat surprisingly, the theory described above can also be readily adapted to this scenario.

Finally, we also note that similar guarantees hold for a range of alternative algorithmic approaches. For example, both the approaches based on alternating minimization and spectral methods with iterations similar to the proximal method described in Section 3 have been shown to provide exact reconstruction under similar coherence assumptions, but at the cost of a slight increase in the required number of observations. In particular, the results of all involve a dependence on the condition number of ${\mathbf{X}}_{0}$ in which the number of required observations grows when the $R^{\text{th}}$ singular value gets too small. By a clever iterative scheme, shows that it is possible to eliminate this dependence in the context of alternating minimization. However, just as in the Gaussian measurement case, the best-known guarantees for alternating minimization also involve a polynomial dependence on the rank $R$ as opposed to the linear dependence which can be obtained using other approaches. Nevertheless, the substantial computational advantages of these approaches means that they provide an attractive alternative in practice.

Nonlinear Observation Models

Although the theoretical results described in Sections 4 and 5 are quite impressive, there is an important gap between the observation model described by (1) and many common applications of low-rank matrix recovery. As an example, consider the matrix completion problem in the context of a recommendation system where $\boldsymbol{X}_{0}$ represents a matrix whose entries each represent a rating for a particular user on a particular item. In most practical recommendation systems (or indeed, any system soliciting any kind of feedback from people), the observations are “quantized”, for example, to the set of integers between 1 and 5. If we believe that it is possible for a user’s true rating to be, for example, 4.5, then we must account for the impact of this “quantization noise” on our recovery. Of course, one could potentially treat quantization simply as a form of bounded noise and rely on the existing stability guarantees mentioned in Sections 4 and 5, but this is somewhat unsatisfying because the ratings are not simply quantized — there are also hard limits placed on the minimum and maximum allowable ratings. (Why should we suppose that an item given a rating of 5 could not have a true underlying rating of 6 or 7 or 10? And note that if this is indeed the case, then the “quantization error” can be potentially extremely large.) In such a situation, it can be much more advantageous to directly consider a nonlinear observation model of the form

where $Q(\cdot)$ is a scalar function that captures the impact of quantization, or any other potential nonlinearity of interest. We describe some concrete examples below.

The inadequacy of standard low-rank matrix recovery techniques in dealing with this effect is particularly pronounced when we consider problems where each observation is quantized to a single-bit. In particular, suppose that our observations are given by

In such a case, the assumptions made in the standard theory of matrix recovery do not apply, standard algorithms are ill-posed, and an alternative theory is required. To see why, simply observe that in the noise-free setting (where $z_{i}=0$), we could rescale ${\mathbf{X}}_{0}$ arbitrarily without changing any observations.In fact, in the noise-free setting the situation is even worse than one might suspect. Even if the normalization is fixed/known a priori, the problem remains highly ill-posed. See for further discussion.

What is perhaps somewhat surprising is that when considering the noisy setting the situation completely changes — the noise has a “dithering” effect and the problem becomes well-posed. In fact, it is possible to show that one can sometimes recover ${\mathbf{X}}_{0}$ to the same degree of accuracy that is possible when given access to completely unquantized measurements, and that given sufficiently many measurements it is possible to recover ${\mathbf{X}}_{0}$ to an arbitrary level of precision. For the specific case of matrix completion, this observation model is analyzed in with the main conclusion being that it is possible to recover any ${\mathbf{X}}_{0}$ belonging to a certain class of approximately low-rank matrices up to an error proportional to $\sqrt{R(M+N)/L}$, so that by taking $L=CR(M+N)$ one can drive the recovery error to be arbitrarily small. The recovery algorithm in is a simple modification of (7), but where the fidelity constraint $\|{\mathcal{A}}(\boldsymbol{X})-\boldsymbol{y}\|_{F}^{2}$ is replaced by the negative log-likelihood of $\boldsymbol{X}$ given observations from the model in (24). A limitation of this approach is that it requires knowledge of the distribution of the noise $z_{i}$, but note that for many common noise distributions this still results in a convex optimization problem which can be solved with variants of the same algorithms described in Section 3. Similar results are described in , which again uses a penalized maximum-likelihood estimator but with a different regularizer than the nuclear norm. See also and which suggest potential improvements for the case of exactly low-rank matrices. Note that many of these results are of interest even in the case where every entry of the matrix is observed, as this provides a theoretically justified way to reveal the “underlying” low-rank matrix given only quantized data.

While most of the existing research into the one-bit observation model has focused on the model in (24), and specifically in the matrix completion setting, it is important to note that it would not be difficult to extend this literature to handle the Gaussian observation model described in Section 4 and/or related one-bit observation models. As an example, an alternative model might involve a setting where we can only tell if the magnitude of our observations is “small” or “large”, as captured by a model along the lines of

A similar theory could be developed for this setting. See for example applications and one possible approach.

2 Comparisons

Another application of the one-bit observation models in (24) or (25) arises when one considers scenarios involving comparisons between different entries in the matrix ${\mathbf{X}}_{0}$. This might occur in the context of paired comparisons, where a person evaluates a pair of items and indicates whether they are similar or dissimilar, or whether one is preferred to the other. This type of data frequently arises when dealing with judgements made by human subjects, since people are typically more accurate and find it easier to make such judgements than to assign numerical scores . These settings can be readily accommodated by the models in (24) or (25) by slightly modifying the $\boldsymbol{A}_{i}$. In particular, when comparing $X_{0}[m,n]$ to $X_{0}[m^{\prime},n^{\prime}]$ one can set $\boldsymbol{A}_{i}=\boldsymbol{e}_{m}\boldsymbol{e}_{n}^{T}-\boldsymbol{e}_{m^{\prime}}\boldsymbol{e}_{n^{\prime}}^{T}$ so that, for example, (24) becomes

A similar theory can be developed for low-rank recovery under this model. For example, see for discussion of paired comparisons and for a generalization to ordinal comparisons among groups of entries in ${\mathbf{X}}_{0}$.

3 Categorical observations and other noise models

While the discussion above has focused on the one-bit case where $Q(\cdot)=\text{sign}(\cdot)$, the algorithms developed for this case can often be readily extended to handle more general nonlinearities. For example, in the case of more general quantization schemes to arbitrary finite alphabets, one must simply be able to compute the log-likelihood function in order to be able to compute the regularized maximum-likelihood estimate as analyzed in . This has applications to multi-bit quantization schemes (e.g., ratings being quantized to the integers from 1 to 5) and also to handle general forms of categorical data (e.g., group membership, race/ethnicity, zip code, etc.). Moreover, in many cases the analysis can also be extended to handle this kind of observation . Finally, it is also worth noting that this body of work has also established a variety of both algorithmic and analytical techniques for handling a variety of complex probabilistic observation models. As a result, this has laid the foundation for a number of works which consider more general noise models beyond simple bounded perturbations or Gaussian measurement noise. For example, explores the impact of the signal-dependent (Poisson) noise that arises in dealing with applications involving count data. These techniques are further generalized in to general exponential noise families. See also for a treatment of matrix factorization problems allowing for a quite general family of noise models.

Lifting

The progress in recovering a low-rank matrix from an incomplete set of linear measurements described above has also affected the way we think about solving quadratic and bilinear systems of equations. There is a simple, but perhaps until recently under-appreciated, way to re-cast a system of quadratic equations as a system of linear equations whose solution obeys a rank constraint. This method, known as lifting, is best illustrated with a small, concrete example. Consider the following system of $6$ quadratic equations in three unknowns $v_{1},v_{2},v_{3}$:

These equations are quadratic in the entries of the vector $\boldsymbol{v}$, but they are linear in the entries of the matrix

For example, the first two equations above can be written as

The reader can verify that $v_{1}=-1,v_{2}=3,v_{3}=5$ is indeed a solution to the equations above. This solution is also unique up to a sign change.

Quadratic and bilinear problems that arise in applications typically have special structure in the equation coefficients $\boldsymbol{A}_{m}$. One type of problem that has received considerable attention in the literature recently is recovering a vector $\boldsymbol{v}_{0}\in R^{N}$ from observations of the magnitude of a series of linear measurements:

By squaring the observations in (27), we see that

Finally, there are mathematical guarantees for other algorithms for solving the nonlinear inverse problem in (27). Algorithms very similar to the alternating minimization algorithm from Section 3.4 have been used for phase retrieval since the 1970s . In , a theoretical analysis of alternating projections for phase retrieval was given (along with an intelligent method for initializing the algorithm) that gives a guarantee of effectiveness when $L\gtrsim N\log^{3}N$. Recent work in shows that a local descent algorithm again coupled with a smart initialization is as effective for phase retrieval as the convex relaxations above while being far more computationally efficient. Algorithms of this type also have performance guarantees for recovering rank-$R$ matrices from measurements of the form (28) .

An extended survey of the recent work on this problem, including a more in depth discussion of most of the results above, can be found in .

2 Blind deconvolution

Perhaps an even more important and prevalent problem, especially in the signal processing and communications communities, is blind deconvolution. Here we observe the (discrete-time) convolution of two signals,

The DTFT of the observations $\boldsymbol{y}$, after we zero-pad it outside its support, is the point-by-point multiplication of the DTFTs of $\boldsymbol{s}$ and $\boldsymbol{h}$ (also after zero-padding). Since all three signals are zero outside of $\{0,\ldots,L-1\}$, we can evaluate the DTFT of the observations at $L$ equally spaced frequencies between $-\pi$ and $\pi$ by computing the vector inner products

This says that the dimension of our linear models (which determines their expressive power) can be within a logarithmic factor of the ambient dimension. General identifiablility results for deterministic subspace models are discussed in .

The lifting technique described above for blind deconvolution is just one of many methods that have been proposed for this important problem (see the books or the survey for an overview of algorithms used in digital communications, imaging, and other areas). The lifting technique, however, puts the problem into the realm of optimization. This makes it very natural to add (convex) constraints for modeling prior knowledge about the signal, or integrate indirect or partial measurements. In , for example, it is shown that if the image is modulated before being blurred by an unknown kernel, the recovery problem is actually very well-posed. Recovery in this scenario is possible without any prior knowledge of the image, and the restrictions on the blur kernel are very mild. In , the lifting framework is applied to the closely related problem of auto-calibrating sensor arrays. Techniques for encouraging $\boldsymbol{s}$ and/or $\boldsymbol{h}$ to be sparse have recently been studied in .

References