1-Bit Matrix Completion

Introduction

The problem of recovering a matrix from an incomplete sampling of its entries—also known as matrix completion—arises in a wide variety of practical situations. In many of these settings, however, the observations are not only incomplete, but also highly quantized, often even to a single bit. In this paper we consider a statistical model for such data where instead of observing a real-valued entry as in the original matrix completion problem, we are now only able to see a positive or negative rating. This binary output is generated according to a probability distribution which is parameterized by the corresponding entry of the unknown low-rank matrix $\boldsymbol{M}$. The central question we ask in this paper is: “Given observations of this form, can we recover the underlying matrix?”

Questions of this form are often asked in the context of binary PCA or logistic PCA. There are a number of compelling algorithmic papers on these subjects, including , which suggest positive answers on simulated and real-world data. In this paper, we give the first theoretical accuracy guarantees under a generalized linear model. We show that $O(rd)$ binary observations are sufficient to accurately recover a $d\times d$, rank-$r$ matrix by convex programming. Our theory is inspired by the unquantized matrix completion problem and the closely related problem of 1-bit compressed sensing, described below.

Although these results are quite impressive, there is an important gap between the statement of the problem as considered in the matrix completion literature and many of the most common applications discussed therein. As an example, consider collaborative filtering and the now-famous “Netflix problem.” In this setting, we assume that there is some unknown matrix whose entries each represent a rating for a particular user on a particular movie. Since any user will rate only a small subset of possible movies, we are only able to observe a small fraction of the total entries in the matrix, and our goal is to infer the unseen ratings from the observed ones. If the rating matrix has low rank, then this would seem to be the exact problem studied in the matrix completion literature. However, there is a subtle difference: the theory developed in this literature generally assumes that observations consist of (possibly noisy) continuous-valued entries of the matrix, whereas in the Netflix problem the observations are “quantized” to the set of integers between 1 and 5. If we believe that it is possible for a user’s true rating for a particular movie 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, but this is somewhat unsatisfying because the ratings aren’t just quantized — there are also hard limits placed on the minimum and maximum allowable ratings. (Why should we suppose that a movie given a rating of 5 could not have a true underlying rating of 6 or 7 or 10?) The inadequacy of standard matrix completion techniques in dealing with this effect is particularly pronounced when we consider recommender systems where each rating consists of a single bit representing a positive or negative rating (consider for example rating music on Pandora, the relevance of advertisements on Hulu, or posts on sites such as MathOverflow). In such a case, the assumptions made in the existing theory of matrix completion do not apply, standard algorithms are ill-posed, and alternative theory is required.

2 Statistical learning theory and matrix completion

While the theory of matrix completion gained quite a bit of momentum following advances in compressed sensing, earlier results from Srebro et al. also addressed this problem from a slightly different perspective rooted in the framework of statistical learning theory. These results also deal with the binary setting that we consider here. They take a model-free approach and prove generalization error bounds; that is, they give conditions under which good agreement on the observed data implies good agreement on the entire matrix. For example, in , agreement is roughly measured via the fraction of predicted signs that are correct, but this can also be extended to other notions of agreement . From the perspective of statistical learning theory, this corresponds to bounding the generalization error under various classes of loss functions.

One important difference between the statistical learning approach and the path taken in this paper is that here we focus on parameter estimation — that is, we seek to recover the matrix $\boldsymbol{M}$ itself (or the distribution parameterized by $\boldsymbol{M}$ that governs our observations) — although we also prove generalization error bounds en route to our main results on parameter and distribution recovery. We discuss the relationship between our approach and the statistical learning approach in more detail in Section A.1 (see Remark 3). Briefly, our generalization error bounds correspond to the case where the loss function is the log likelihood, and do not seem to fit directly within the framework of the existing literature.

3 1-bit compressed sensing and sparse logistic regression

4 Challenges

In this paper, we extend the theory of matrix completion to the case of 1-bit observations. We consider a general observation model but focus mainly on two particular possibilities: the models of logistic and probit regression. We discuss these models in greater detail in Section 2.1, but first we note that several new challenges arise when trying to leverage results in 1-bit compressed sensing and sparse logistic regression to develop a theory for 1-bit matrix completion. First, matrix completion is in some sense a more challenging problem than compressed sensing. Specifically, some additional difficulty arises because the set of low-rank matrices is “coherent” with single entry measurements (see ). In particular, the sampling operator does not act as a near-isometry on all matrices of interest, and thus the natural analogue to the restricted isometry property from compressed sensing cannot hold in general—there will always be certain low-rank matrices that we cannot hope to recover without essentially sampling every entry of the matrix. For example, consider a matrix that consists of a single nonzero entry (which we might never observe). The typical way to deal with this possibility is to consider a reduced set of low-rank matrices by placing restrictions on the entry-wise maximum of the matrix or its singular vectors—informally, we require that the matrix is not too “spiky”.

We introduce an entirely new dimension of ill-posedness by restricting ourselves to 1-bit observations. To illustrate this, we describe one version of 1-bit matrix completion in more detail (the general problem definition is given in Section 2.1 below). Consider a $d\times d$ matrix $\boldsymbol{M}$ with rank $r$. Suppose we observe a subset $\Omega$ of entries of a matrix $\boldsymbol{Y}$. The entries of $\boldsymbol{Y}$ depend on $\boldsymbol{M}$ in the following way:

After considering this example, the problem might seem hopeless. However, an interesting surprise is that when we add noise to the problem (that is, when $\boldsymbol{Z}\neq\boldsymbol{0}$ is an appropriate stochastic matrix) the picture completely changes—this noise has a “dithering” effect and the problem becomes well-posed. In fact, we will show that in this setting we can sometimes recover $\boldsymbol{M}$ to the same degree of accuracy that is possible when given access to completely unquantized measurements! In particular, under appropriate conditions, $O(rd)$ measurements are sufficient to accurately recover $\boldsymbol{M}$.

5 Applications

The problem of 1-bit matrix completion arises in nearly every application that has been proposed for “unquantized” matrix completion. To name a few:

Recommender systems: As mentioned above, collaborative filtering systems often involve discretized recommendations . In many cases, each observation will consist simply of a “thumbs up” or “thumbs down” thus delivering only 1 bit of information (consider for example rating music on Pandora, the relevance of advertisements on Hulu, or posts on sites such as MathOverflow). Such cases are a natural application for 1-bit matrix completion.

Analysis of survey data: Another potential application for matrix completion is to analyze incomplete survey data. Such data is almost always heavily quantized since people are generally not able to distinguish between more than $7\pm 2$ categories . 1-bit matrix completion provides a method for analyzing incomplete (or potentially even complete) survey designs containing simple yes/no or agree/disagree questions.

Distance matrix recovery and multidimensional scaling: Yet another common motivation for matrix completion is to localize nodes in a sensor network from the observation of just a few inter-node distances . This is essentially a special case of multidimensional scaling (MDS) from incomplete data . In general, work in the area assumes real-valued measurements. However, in the sensor network example (as well as many other MDS scenarios), the measurements may be very coarse and might only indicate whether the nodes are within or outside of some communication range. While there is some existing work on MDS using binary data and MDS using incomplete observations with other kinds of non-metric data , 1-bit matrix completion promises to provide a principled and unifying approach to such problems.

Quantum state tomography: Low-rank matrix recovery from incomplete observations also has applications to quantum state tomography . In this scenario, mixed quantum states are represented as Hermitian matrices with nuclear norm equal to 1. When the state is nearly pure, the matrix can be well approximated by a low-rank matrix and, in particular, fits the model given in Section 2.2 up to a rescaling. Furthermore, Pauli-operator-based measurements give probabilistic binary outputs. However, these are based on the inner products with the Pauli matrices, and thus of a slightly different flavor than the measurements considered in this paper. Nevertheless, while we do not address this scenario directly, our theory of 1-bit matrix completion could easily be adapted to quantum state tomography.

6 Notation

We now provide a brief summary of some of the key notation used in this paper. We use $[d]$ to denote the set of integers $\{1,\ldots,d\}$. We use capital boldface to denote a matrix (e.g., $\boldsymbol{M}$) and standard text to denote its entries (e.g., $M_{i,j}$). Similarly, we let $\boldsymbol{0}$ denote the matrix of all-zeros and $\boldsymbol{1}$ the matrix of all-ones. We let $\|\boldsymbol{M}\|$ denote the operator norm of $\boldsymbol{M}$, $\left\|\boldsymbol{M}\right\|_{F}=\sqrt{\sum_{i,j}M_{i,j}^{2}}$ denote the Frobenius norm of $\boldsymbol{M}$, $\left\|\boldsymbol{M}\right\|_{*}$ denote the nuclear or Schatten-1 norm of $\boldsymbol{M}$ (the sum of the singular values), and $\left\|\boldsymbol{M}\right\|_{\infty}=\max_{i,j}\left|M_{i,j}\right|$ denote the entry-wise infinity-norm of $\boldsymbol{M}$. We will use the Hellinger distance, which, for two scalars $p,q\in$, is given by

This gives a standard notion of distance between two binary probability distributions. We also allow the Hellinger distance to act on matrices via the average Hellinger distance over their entries: for matrices $\boldsymbol{P},\boldsymbol{Q}\in^{d_{1}\times d_{2}}$, we define

Finally, for an event $\mathcal{E}$ ,$\mathds{1}_{[\mathcal{E}]}$ is the indicator function for that event, i.e., $\mathds{1}_{[\mathcal{E}]}$ is $1$ if $\mathcal{E}$ occurs and otherwise.

7 Organization of the paper

We proceed in Section 2 by describing the 1-bit matrix completion problem in greater detail. In Section 3 we state our main results. Specifically, we propose a pair of convex programs for the 1-bit matrix completion problem and establish upper bounds on the accuracy with which these can recover the matrix $\boldsymbol{M}$ and the distribution of the observations $\boldsymbol{Y}$. We also establish lower bounds, showing that our upper bounds are nearly optimal. In Section 4 we describe numerical implementations of our proposed convex programs and demonstrate their performance on a number of synthetic and real-world examples. Section 5 concludes with a brief discussion of future directions. The proofs of our main results are provided in the appendix.

The 1-bit matrix completion problem

We now consider two natural choices for $f$ (or equivalently, for $\boldsymbol{Z}$):

The logistic regression model, which is common in statistics, is captured by (2) with $f(x)=\frac{e^{x}}{1+e^{x}}$ and by (1) with $Z_{i,j}$ i.i.d. according to the standard logistic distribution.

The probit regression model is captured by (2) by setting $f(x)=1-\Phi(-x/\sigma)=\Phi(x/\sigma)$ where $\Phi$ is the cumulative distribution function of a standard Gaussian and by (1) with $Z_{i,j}$ i.i.d. according to a mean-zero Gaussian distribution with variance $\sigma^{2}$.

2 Approximately low-rank matrices

The particular choice of scaling, $\alpha\sqrt{rd_{1}d_{2}}$, arises from the following considerations. Suppose that each entry of $\boldsymbol{M}$ is bounded in magnitude by $\alpha$ and that $\operatorname*{rank}(\boldsymbol{M})\leq r$. Then

Thus, the assumption that $\left\|\boldsymbol{M}\right\|_{*}\leq\alpha\sqrt{rd_{1}d_{2}}$ is a relaxation of the conditions that $\operatorname*{rank}(\boldsymbol{M})\leq r$ and $\left\|\boldsymbol{M}\right\|_{\infty}\leq\alpha$. The condition that $\left\|\boldsymbol{M}\right\|_{\infty}\leq\alpha$ essentially means that the probability of seeing a $+1$ or $-1$ does not depend on the dimension. It is also a way of enforcing that $\boldsymbol{M}$ should not be too “spiky”; as discussed above this is an important assumption in order to make the recovery of $\boldsymbol{M}$ well-posed (e.g., see ).

Main results

We now state our main results. We will have two goals—the first is to accurately recover $\boldsymbol{M}$ itself, and the second is to accurately recover the distribution of $\boldsymbol{Y}$ given by $f(\boldsymbol{M})$.Strictly speaking, $f(\boldsymbol{M})\in^{d_{1}\times d_{2}}$ is simply a matrix of scalars, but these scalars implicitly define the distribution of $\boldsymbol{Y}$, so we will sometimes abuse notation slightly and refer to $f(\boldsymbol{M})$ as the distribution of $\boldsymbol{Y}$. All proofs are contained in the appendix.

In order to approximate either $\boldsymbol{M}$ or $f(\boldsymbol{M})$, we will maximize the log-likelihood function of the optimization variable $\boldsymbol{X}$ given our observations subject to a set of convex constraints. In our case, the log-likelihood function is given by

To recover $\boldsymbol{M}$, we will use the solution to the following program:

To recover the distribution $f(\boldsymbol{M})$, we need not enforce the infinity-norm constraint, and will use the following simpler program:

In many cases, $\mathcal{L}_{\Omega,\boldsymbol{Y}}(\boldsymbol{X})$ is a concave function and thus the above programs are convex. This can be easily checked in the case of the logistic model and can also be verified in the case of the probit model (e.g., see ).

2 Recovery of the matrix

We now state our main result concerning the recovery of the matrix $\boldsymbol{M}$. As discussed in Section 1.4 we place a “non-spikiness” condition on $\boldsymbol{M}$ to make recovery possible; we enforce this with an infinity-norm constraint. Further, some assumptions must be made on $f$ for recovery of $\boldsymbol{M}$ to be feasible. We define two quantities $L_{\alpha}$ and $\beta_{\alpha}$ which control the “steepness” and “flatness” of $f$, respectively:

In this paper we will restrict our attention to $f$ such that $L_{\alpha}$ and $\beta_{\alpha}$ are well-defined. In particular, we assume that $f$ and $f^{\prime}$ are non-zero in $[-\alpha,\alpha]$. This assumption is fairly mild—for example, it includes the logistic and probit models (as we will see below in Remark 1). The quantity $L_{\alpha}$ appears only in our upper bounds, but it is generally well behaved. The quantity $\beta_{\alpha}$ appears both in our upper and lower bounds. Intuitively, it controls the “flatness” of $f$ in the interval $[-\alpha,\alpha]$—the flatter $f$ is, the larger $\beta_{\alpha}$ is. It is clear that some dependence on $\beta_{\alpha}$ is necessary. Indeed, if $f$ is perfectly flat, then the magnitudes of the entries of $\boldsymbol{M}$ cannot be recovered, as seen in the noiseless case discussed in Section 1.4. Of course, when $\alpha$ is a fixed constant and $f$ is a fixed function, both $L_{\alpha}$ and $\beta_{\alpha}$ are bounded by fixed constants independent of the dimension.

with $C_{\alpha}:=C_{2}\alpha L_{\alpha}\beta_{\alpha}$. If $n\geq(d_{1}+d_{2})\log(d_{1}d_{2})$ then this simplifies to

Above, $C_{1}$ and $C_{2}$ are absolute constants.

Note that the theorem also holds if $\Omega=[d_{1}]\times[d_{2}]$, i.e., if we sample each entry exactly once or observe a complete realization of $\boldsymbol{Y}$. Even in this context, the ability to accurately recover $\boldsymbol{M}$ is somewhat surprising.

The logistic model satisfies the hypotheses of Theorem 1 with $\beta_{\alpha}=\frac{(1+e^{\alpha})^{2}}{e^{\alpha}}\approx e^{\alpha}$ and $L_{\alpha}=1$. The probit model has

where we can take $c_{1}=\pi$ and $c_{2}=8$. In particular, in the probit model the bound in (6) reduces to

Hence, when $\sigma<\alpha$, increasing the size of the noise leads to significantly improved error bounds—this is not an artifact of the proof. We will see in Section 3.4 that the exponential dependence on $\alpha$ in the logistic model (and on $\alpha/\sigma$ in the probit model) is intrinsic to the problem. Intuitively, we should expect this since for such models, as $\|\boldsymbol{M}\|_{\infty}$ grows large, we can essentially revert to the noiseless setting where estimation of $\boldsymbol{M}$ is impossible. Furthermore, in Section 3.4 we will also see that when $\alpha$ (or $\alpha/\sigma$) is bounded by a constant, the error bound (6) is optimal up to a constant factor. Fortunately, in many applications, one would expect $\alpha$ to be small, and in particular to have little, if any, dependence on the dimension. This ensures that each measurement will always have a non-vanishing probability of returning $1$ as well as a non-vanishing probability of returning $-1$.

The assumption made in Theorem 1 (as well as Theorem 2 below) that $\left\|\boldsymbol{M}\right\|_{*}\leq\alpha\sqrt{d_{1}d_{2}r}$ does not mean that we are requiring the matrix $\boldsymbol{M}$ to be low rank. We express this constraint in terms of $r$ to aid the intuition of researchers well-versed in the existing literature on low-rank matrix recovery. (If $\boldsymbol{M}$ is exactly rank $r$ and satisfies $\|\boldsymbol{M}\|_{\infty}\leq\alpha$, then as discussed in Section 2.2, $\boldsymbol{M}$ will automatically satisfy this constraint.) If one desires, one may simplify the presentation by replacing $\alpha\sqrt{d_{1}d_{2}r}$ with a parameter $\lambda$ and simply requiring $\left\|\boldsymbol{M}\right\|_{*}\leq\lambda$, in which case (6) reduces to

3 Recovery of the distribution

In many situations, we might not be interested in the underlying matrix $\boldsymbol{M}$, but rather in determining the distribution of the entries of $\boldsymbol{Y}$. For example, in recommender systems, a natural question would be to determine the likelihood that a user would enjoy a particular unrated item.

Surprisingly, this distribution may be accurately recovered without any restriction on the infinity-norm of $\boldsymbol{M}$. This may be unexpected to those familiar with the matrix completion literature in which “non-spikiness” constraints seem to be unavoidable. In fact, we will show in Section 3.4 that the bound in Theorem 2 is near-optimal; further, we will show that even under the added constraint that $\left\|\boldsymbol{M}\right\|_{\infty}\leq\alpha$, it would be impossible to estimate $f(\boldsymbol{M})$ significantly more accurately.

Furthermore, as long as $n\geq(d_{1}+d_{2})\log(d_{1}d_{2})$, we have

Above, $C_{1}$ and $C_{2}$ are absolute constants.

4 Room for improvement?

We now discuss the extent to which Theorems 1 and 2 are optimal. We give three theorems, all proved using information theoretic methods, which show that these results are nearly tight, even when some of our assumptions are relaxed. Theorem 3 gives a lower bound to nearly match the upper bound on the error in recovering $\boldsymbol{M}$ derived in Theorem 1. Theorem 4 compares our upper bounds to those available without discretization and shows that very little is lost when discretizing to a single bit. Finally, Theorem 5 gives a lower bound matching, up to a constant factor, the upper bound on the error in recovering the distribution $f(\boldsymbol{M})$ given in Theorem 2. Theorem 5 also shows that Theorem 2 does not suffer by dropping the canonical “spikiness” constraint.

Our lower bounds require a few assumptions, so before we delve into the bounds themselves, we briefly argue that these assumptions are rather innocuous. First, without loss of generality (since we can always adjust $f$ to account for rescaling $\boldsymbol{M}$), we assume that $\alpha\geq 1$. Next, we require that the parameters be sufficiently large so that

for an absolute constant $C_{0}$. Note that we could replace this with a simpler, but still mild, condition that $d_{1}>C_{0}$. Finally, we also require that $r\geq c$ where $c$ is either 1 or 4 and that $r\leq O(\min\{d_{1},d_{2}\}/\alpha^{2})$, where $O(\cdot)$ hides parameters (which may differ in each Theorem) that we make explicit below. This last assumption simply means that we are in the situation where $r$ is significantly smaller than $d_{1}$ and $d_{2}$, i.e., the matrix is of approximately low rank.

denote the set of matrices whose recovery is guaranteed by Theorem 1.

Fix $\alpha,r,d_{1},$ and $d_{2}$ to be such that $r\geq 4$ and (10) holds. Let $\beta_{\alpha}$ be defined as in (5), and suppose that $f^{\prime}(x)$ is decreasing for $x>0$. Let $\Omega$ be any subset of $[d_{1}]\times[d_{2}]$ with cardinality $n$, and let $\boldsymbol{Y}$ be as in (2). Consider any algorithm which, for any $\boldsymbol{M}\in K$, takes as input $Y_{i,j}$ for $(i,j)\in\Omega$ and returns $\widehat{\boldsymbol{M}}$. Then there exists $\boldsymbol{M}\in K$ such that with probability at least $3/4$,

as long as the right-hand side of (12) exceeds $r\alpha^{2}/\min(d_{1},d_{2})$. Above, $C_{1}$ and $C_{2}$ are absolute constants.Here and in the theorems below, the choice of $3/4$ in the probability bound is arbitrary, and can be adjusted at the cost of changing $C_{0}$ in (10) and $C_{1}$ and $C_{2}$. Similarly, $\beta_{\frac{3}{4}\alpha}$ can be replaced by $\beta_{(1-\epsilon)\alpha}$ for any $\epsilon>0$.

The requirement that the right-hand side of (12) be larger than $r\alpha^{2}/\min(d_{1},d_{2})$ is satisfied as long as $r\leq O(\min\{d_{1},d_{2}\}/\alpha^{2})$. In particular, it is satisfied whenever

for a fixed constant $C_{3}$. Note also that in the latent variable model in (1), $f^{\prime}(x)$ is simply the probability density of $Z_{i,j}$. Thus, the requirement that $f^{\prime}(x)$ be decreasing is simply asking the probability density to have decreasing tails. One can easily check that this is satisfied for the logistic and probit models.

Note that if $\alpha$ is bounded by a constant and $f$ is fixed (in which case $\beta_{\alpha}$ and $\beta_{\alpha^{\prime}}$ are bounded by a constant), then the lower bound of Theorem 3 matches the upper bound given in (6) up to a constant. When $\alpha$ is not treated as a constant, the bounds differ by a factor of $\sqrt{\beta_{\alpha}}$. In the logistic model $\beta_{\alpha}\approx e^{\alpha}$ and so this amounts to the difference between $e^{\alpha/2}$ and $e^{\alpha}$. The probit model has a similar change in the constant of the exponent.

4.2 Recovery from unquantized measurements

Next we show that, surprisingly, very little is lost by discretizing to a single bit. In Theorem 4, we consider an “unquantized” version of the latent variable model in (1) with Gaussian noise. That is, let $\boldsymbol{Z}$ be a matrix of i.i.d. Gaussian random variables, and suppose the noisy entries $M_{i,j}+Z_{i,j}$ are observed directly, without discretization. In this setting, we give a lower bound that still nearly matches the upper bound given in Theorem 1, up to the $\beta_{\alpha}$ term.

Fix $\alpha,r,d_{1},$ and $d_{2}$ to be such that $r\geq 1$ and (10) holds. Let $\Omega$ be any subset of $[d_{1}]\times[d_{2}]$ with cardinality $n$, and let $\boldsymbol{Z}$ be a $d_{1}\times d_{2}$ matrix with i.i.d. Gaussian entries with variance $\sigma^{2}$. Consider any algorithm which, for any $\boldsymbol{M}\in K$, takes as input $Y_{i,j}=M_{i,j}+Z_{i,j}$ for $(i,j)\in\Omega$ and returns $\widehat{\boldsymbol{M}}$. Then there exists $\boldsymbol{M}\in K$ such that with probability at least $3/4$,

as long as the right-hand side of (13) exceeds $r\alpha^{2}/\min(d_{1},d_{2})$. Above, $C_{1}$ and $C_{2}$ are absolute constants.

The requirement that the right-hand side of (13) be larger than $r\alpha^{2}/\min(d_{1},d_{2})$ is satisfied whenever

Following Remark 1, the lower bound given in (13) matches the upper bound proven in Theorem 1 for the solution to (4) up to a constant, as long as $\alpha/\sigma$ is bounded by a constant. In other words:

When the signal-to-noise ratio is constant, almost nothing is lost by quantizing to a single bit.

Perhaps it is not particularly surprising that 1-bit quantization induces little loss of information in the regime where the noise is comparable to the underlying quantity we wish to estimate—however, what is somewhat of a surprise is that the simple convex program in (4) can successfully recover all of the information contained in these 1-bit measurements.

Before proceeding, we also briefly note that our Theorem 4 is somewhat similar to Theorem 3 in . The authors in consider slightly different sets $K$: these sets are more restrictive in the sense that it is required that $\alpha\geq\sqrt{32\log n}$ and less restrictive because the nuclear-norm constraint may be replaced by a general Schatten-p norm constraint. It was important for us to allow $\alpha=O(1)$ in order to compare with our upper bounds due to the exponential dependence of $\beta_{\alpha}$ on $\alpha$ in Theorem 1 for the probit model. This led to some new challenges in the proof. Finally, it is also noteworthy that our statements hold for arbitrary sets $\Omega$, while the argument in is only valid for a random choice of $\Omega$.

4.3 Recovery of the distribution from 1-bit measurements

To conclude we address the optimality of Theorem 2. We show that under mild conditions on $f$, any algorithm that recovers the distribution $f(\boldsymbol{M})$ must yield an estimate whose Hellinger distance deviates from the true distribution by an amount proportional to $\alpha\sqrt{rd_{1}d_{2}/n}$, matching the upper bound of (9) up to a constant. Notice that the lower bound holds even if the algorithm is promised that $\|\boldsymbol{M}\|_{\infty}\leq\alpha$, which the upper bound did not require.

Fix $\alpha,r,d_{1},$ and $d_{2}$ to be such that $r\geq 4$ and (10) holds. Let $L_{1}$ be defined as in (5), and suppose that $f^{\prime}(x)\geq c$ and $c^{\prime}\leq f(x)\leq 1-c^{\prime}$ for $x\in$, for some constants $c,c^{\prime}>0$. Let $\Omega$ be any subset of $[d_{1}]\times[d_{2}]$ with cardinality $n$, and let $\boldsymbol{Y}$ be as in (2). Consider any algorithm which, for any $\boldsymbol{M}\in K$, takes as input $Y_{i,j}$ for $(i,j)\in\Omega$ and returns $\widehat{\boldsymbol{M}}$. Then there exists $\boldsymbol{M}\in K$ such that with probability at least $3/4$,

as long as the right-hand side of (14) exceeds $r\alpha^{2}/\min(d_{1},d_{2})$. Above, $C_{1}$ and $C_{2}$ are constants that depend on $c,c^{\prime}$.

The requirement that the right-hand side of (14) be larger than $r\alpha^{2}/\min(d_{1},d_{2})$ is satisfied whenever

for a constant $C_{3}$ that depends only on $c,c^{\prime}$. Note also that the condition that $f$ and $f^{\prime}$ be well-behaved in the interval $$is satisfied for the logistic model with$c=1/4$and$c^{\prime}=\frac{1}{1+e}\leq 0.269$. Similarly, we may take$c=0.242$and$c^{\prime}=0.159$ in the probit model.

Simulations

Before presenting the proofs of our main results, we provide algorithms and a suite of numerical experiments to demonstrate their usefulness in practice.The code for these algorithms, as well as for the subsequent experiments, is available online at http://users.ece.gatech.edu/${\sim}$mdavenport/. We present algorithms to solve the convex programs (3) and (4), and using these we can recover $\boldsymbol{M}$ (or $f(\boldsymbol{M})$) via 1-bit matrix completion.

We begin with the observation that both (3) and (4) are instances of the more general formulation

One possible solver for problems of the form (15) is the nonmonotone spectral projected-gradient (SPG) algorithm proposed by Birgin et al. . Another possibility is to use an accelerated proximal-gradient methods for the minimization of composite functions , which are useful for solving optimization problems of the form

where $f(x)$ and $g(x)$ are convex functions with $g(x)$ possibly non-smooth. This formulation reduces to (15) when choosing $g(x)$ to be the extended-real indicator function corresponding to $\mathcal{C}$:

Both algorithms are iterative and require at each iteration the evaluation of $f(x)$, its gradient $\nabla f(x)$, and an orthogonal projection onto $\mathcal{C}$ (i.e., the prox-function of $g(x)$)

For our experiments we use the SPG algorithm, which we describe in more detail below. The implementation of the algorithm is based on the SPGL1 code .

1.2 Spectral projected-gradient method

In basic gradient-descent algorithms for unconstrained minimization of a convex function $f(x)$, iterates are of the form $x_{k+1}=x_{k}-\alpha_{k}\nabla f(x_{k})$, where the step length $\alpha_{k}\in(0,1]$ is chosen such that sufficient descent in the objective function $f(x)$ is achieved. When the constraint $x\in\mathcal{C}$ is added, the basic scheme can no longer guarantee feasibility of the iterates. Projected gradient methods resolve this problem by including on orthogonal projections back onto the feasible set (17) at each iteration.

The nonmonotone SPG algorithm described in modifies the basic projected gradient method in two major ways. First, it scales the initial search direction using the spectral step-length $\gamma_{k}$ as proposed by Barzilai and Borwein . Second, it relaxes monotonicity of the objective values by requiring sufficient descent relative to the maximum objective over the last $t$ iterates (or $k$ when $k<t$). Two types of line search are considered. The first type is curvilinear and traces the following path:

The second type first determines a projected gradient step, and uses this to obtain the search direction $d_{k}$:

Next, a line search is done along the linear trajectory

In either case, once the step length $\alpha$ is chosen, we set $x_{k+1}=x(\alpha)$, and proceed with the next iteration.

In the 1-bit matrix completion formulation proposed in this paper, the projection onto $\mathcal{C}$ forms the main computational bottleneck. As a result, it is crucial to keep the number of projections to a minimum, and our implementation therefore relies primarily on the line search along the linear trajectory given by (19). The more expensive curvilinear line search is used only when the linear one fails. We have observed that this situation tends to arise only when $x_{k}$ is near optimal.

The optimality condition for (15) is that

Our implementation checks if (20) is approximately satisfied. In addition it imposes bounds on the total number of iterations and the run time.

1.3 Orthogonal projections onto the feasible sets

where the maximum is taken entrywise, and $\lambda\geq 0$ is the smallest value for which $\sum_{i=1}^{d}\max\{\sigma_{i}-\lambda,0\}\leq\tau$.

Unfortunately, no closed form solution is known for the orthogonal projection onto $\mathcal{C}_{2}$. However, the underlying problem

can be solved using iterative methods. In particular, we can rewrite (21) as

and apply the alternating-direction method of multipliers (ADMM) . The augmented Lagrangian for (22) with respect to the constraint $\boldsymbol{W}=\boldsymbol{Z}$ is given by

The ADMM iterates the following steps to solve (22):

Initialize $k=0$, and select $\mu_{k}$, $\boldsymbol{Y}_{k}$, $\boldsymbol{W}_{k}$, $\boldsymbol{Z}_{k}$ such that $\left\|\boldsymbol{W}_{k}\right\|_{\infty}\leq\kappa$ and $\left\|\boldsymbol{Z}_{k}\right\|_{*}\leq\tau$.

Minimize $\mathcal{L}_{\mu}(\boldsymbol{Y}_{k},W,\boldsymbol{Z}_{k})$ with respect to $\boldsymbol{W}$, which can be rewritten as

This is exactly the orthogonal projection of $\boldsymbol{B}=(\boldsymbol{X}+\boldsymbol{Y}_{k}+\mu\boldsymbol{Z}_{k})/(1+\mu)$ onto $\{\boldsymbol{W}\mid\left\|\boldsymbol{W}\right\|_{\infty}\leq\kappa\}$, and gives $\boldsymbol{W}_{k+1}(i,j)=\min\{\kappa,\max\{-\kappa,\boldsymbol{B}(i,j)\}\}$.

Minimize $\mathcal{L}_{\mu}(\boldsymbol{Y}_{k},\boldsymbol{W}_{k+1},\boldsymbol{Z})$ with respect to $\boldsymbol{Z}$. This gives

and simplifies to $\boldsymbol{Z}_{k+1}=\mathcal{P}_{\mathcal{C}_{1}}(\boldsymbol{W}_{k+1}-1/\mu_{k}\boldsymbol{Y}_{k})$.

Update $\boldsymbol{Y}_{k+1}=\boldsymbol{Y}_{k}-\mu(\boldsymbol{W}_{k+1}-\boldsymbol{Z}_{k+1})$, set $\mu_{k+1}=1.05\,\mu_{k}$, and increment $k$.

Return $\boldsymbol{Z}=\boldsymbol{Z}_{k}$ when $\left\|\boldsymbol{W}_{k}-\boldsymbol{Z}_{k}\right\|_{F}\leq\varepsilon$ and $\left\|\boldsymbol{Z}_{k}\right\|_{\infty}-\kappa\leq\varepsilon$ for some sufficiently small $\varepsilon>0$. Otherwise, repeat steps 1–4.

2 Synthetic experiments

To evaluate the performance of this algorithm in practice and to confirm the theoretical results described above, we first performed a number of synthetic experiments. In particular, we constructed a random $d\times d$ matrix $\boldsymbol{M}$ with rank $r$ by forming $\boldsymbol{M}=\boldsymbol{M}_{1}\boldsymbol{M}_{2}^{*}$ where $\boldsymbol{M}_{1}$ and $\boldsymbol{M}_{2}$ are $d\times r$ matrices with entries drawn i.i.d. from a uniform distribution on $[-\frac{1}{2},\frac{1}{2}]$. The matrix is then scaled so that $\|\boldsymbol{M}\|_{\infty}=1$. We then obtained 1-bit observations by adding Gaussian noise of variance $\sigma^{2}$ and recording the sign of the resulting value.

We begin by comparing the performance of the algorithms in (3) and (4) over a range of different values of $\sigma$. In this experiment we set $d=500$, $r=1$, and $n=0.15d^{2}$, and we measured performance of each approach using the squared Frobenius norm of the error (normalized by the norm of the original matrix $\boldsymbol{M}$) and averaged the results over 15 draws of $\boldsymbol{M}$.In evaluating these algorithms we found that it was beneficial in practice to follow the recovery by a “debiasing” step where the recovered matrix $\widehat{\boldsymbol{M}}$ is forced to be rank $r$ by computing the SVD of $\widehat{\boldsymbol{M}}$ and hard thresholding the singular values. In cases where we report the Frobenius norm of the error, we performed this debiasing step, although it does not dramatically impact the performance. The results are shown in Figure 1. We observe that for both approaches, the performance is poor when there is too little noise (when $\sigma$ is small) and when there is too much noise (when $\sigma$ is large). These two regimes correspond to the cases where the noise is either so small that the observations are essentially noise-free or when the noise is so large that each observation is essentially a coin toss. In the regime where the noise is of moderate power, we observe better performance for both approaches. Perhaps somewhat surprisingly, we find that for much of this range, the approach in (4) appears to perform almost as well as (3), even though we do not have any theoretical guarantees for (4). This suggests that adding the infinity-norm constraint as in (3) may have only limited practical benefit, despite the key role this constraint played in our analysis. By using the simpler program in (4) one can greatly simplify the projection step in the algorithm, so in practice this approach may be preferable.

We also conducted experiments evaluating the performance of both (3) and (4) as a function of $n$ for a particular choice of $\sigma$. The results showing the impact of $n$ on (4) are shown in Figure 2 (the results for (3) at this noise level are almost indistinguishable). In this experiment we set $d=200$, and chose $\sigma\approx 0.18$ such that $\log_{10}(\sigma)=0.75$, which lies in the regime where the noise is neither negligible nor overwhelming. We considered matrices with rank $r=3,5,10$ and evaluated the performance over a range of $n$. Figure 2(a) shows the performance in terms of the relative Frobenius norm of the error, and Figure 2(b) shows the performance in terms of the Hellinger distance between the recovered distributions. Consistent with our theoretical results, we observe a decay in the error (under both performance metrics) that appears to behave roughly on the order of $n^{-1/2}$.

3 Collaborative filtering

To evaluate the performance of our algorithm in a practical setting, we consider the MovieLens (100k) data set, which is available for download at http://www.grouplens.org/node/73. This data set consists of 100,000 movie ratings from 1000 users on 1700 movies, with each rating occurring on a scale from 1 to 5. For testing purposes, we converted these ratings to binary observations by comparing each rating to the average rating for the entire dataset (which is approximately 3.5). We then apply the algorithm in (4) (using the logistic model of $f(x)=\frac{e^{x}}{1+e^{x}}$) on a subset of 95,000 ratings to recover an estimate of $\boldsymbol{M}$. However, since our only source of data is the quantized ratings, there is no “ground truth” against which to measure the accuracy of our recovery. Thus, we instead evaluate our performance by checking to see if the estimate of $\boldsymbol{M}$ accurately predicts the sign of the remaining 5000 unobserved ratings in our dataset. The result of this simulation is shown in the first line of Table 1, which gives the accuracy in predicting whether the unobserved ratings are above or below the average rating of 3.5.

By comparison, the second line in the table shows the results obtained using a “standard” method that uses the raw ratings (on a scale from 1 to 5) and tries to minimize the nuclear norm of the recovered matrix subject to a constraint that requires the Frobenius norm of the difference between the recovered and observed entries to be sufficiently small. We implement this traditional approach using the TFOCS software package , and evaluate the results using the same error criterion as the 1-bit matrix completion approach—namely, we compare the recovered ratings to the average (recovered) rating. This approach depends on a number of input parameters: $\alpha$ in (3), the constraint on the Frobenius norm in the traditional case, as well as the internal parameter $\mu$ in TFOCS. We determine the parameter values by simply performing a grid search and selecting those values that lead to the best performance.

Perhaps somewhat surprisingly, the 1-bit approach performs significantly better than the traditional one, even though the traditional approach is given more information in the form of the raw ratings, instead of the binary observations.While not reported in the table, we also observed that the 1-bit approach is relatively insensitive to the choice of $\alpha$, so that this improvement in performance does not rely on a careful parameter setting. The intuition as to how this might be possible is that the standard approach is likely paying a significant penalty for actually requiring that the recovered matrix yields numerical ratings close to “1” or “5” when a user’s true preference could extend beyond this scale.

Discussion

Many of the applications of matrix completion consider discrete observations, often in the form of binary measurements. However, matrix completion from noiseless binary measurements is extremely ill-posed, even if one collects a binary measurement for each of the matrix entries. Fortunately, when there are some stochastic variations (noise) in the observations, recovery becomes well-posed. In this paper we have shown that the unknown matrix can be accurately and efficiently recovered from binary measurements in this setting. When the infinity norm of the unknown matrix is bounded by a constant, our error bounds are tight to within a constant and even match what is possible for undiscretized data. We have also shown that the binary probability distribution can be reconstructed over the entire matrix without any assumption on the infinity-norm, and we have provided a matching lower bound (up to a constant).

Our theory considers approximately low-rank matrices—in particular, we assume that the singular values belong to a scaled Schatten-1 ball. It would be interesting to see whether more accurate reconstruction could be achieved under the assumption that the unknown matrix has precisely $r$ nonzero singular values. We conjecture that the Lagrangian formulation of the problem could be fruitful for this endeavor. It would also be interesting to study whether our ideas can be extended to deal with measurements that are quantized to more than 2 (but still a small number) of different values, but we leave such investigations for future work.

Acknowledgements

We would like to thank Roman Vershynin for helpful discussions and Wenxin Zhou for pointing out an error in an earlier version of the proof of Theorem 6.

Appendix A Proofs of the main results

We now provide the proofs of the main theorems presented in Section 3. To begin, we first define some additional notation that we will need for the proofs. For two probability distributions $\mathcal{P}$ and $\mathcal{Q}$ on a finite set $A$, $D(\mathcal{P}\|\mathcal{Q})$ will denote the Kullback-Leibler (KL) divergence,

where $\mathcal{P}(x)$ denotes the probability of the outcome $x$ under the distribution $\mathcal{P}$. We will abuse this notation slightly by overloading it in two ways. First, for scalar inputs $p,q\in$, we will set

Second, for two matrices $\boldsymbol{P},\boldsymbol{Q}\in^{d_{1}\times d_{2}}$, we define

We first prove Theorem 2. Theorem 1 will then follow from an approximation argument. Finally, our lower bounds will be proved in Section A.3 using information theoretic arguments.

We will actually prove a slightly more general statement, which will be helpful in the proof of Theorem 1. We will assume that $\left\|\boldsymbol{M}\right\|_{\infty}\leq\gamma$, and we will modify the program (4) to enforce $\left\|\boldsymbol{X}\right\|_{\infty}\leq\gamma$. That is, we will consider the program

We will then send $\gamma\to\infty$ to recover the statement of Theorem 2. Formally, we prove the following theorem.

Above, $C_{1}$ and $C_{2}$ are absolute constants.

For the proof of Theorem 6, it will be convenient to work with the function

rather than with $\mathcal{L}_{\Omega,\boldsymbol{Y}}$ itself. The key will be to establish the following concentration inequality.

for some $r\leq\min\{d_{1},d_{2}\}$ and $\alpha\geq 0$. Then

where $C_{0}$ and $C_{1}$ are absolute constants and the probability and the expectation are both over the choice of $\Omega$ and the draw of $Y$.

We will prove this lemma below, but first we show how it implies Theorem 6. To begin, notice that for any choice of $\boldsymbol{X}$,

where the expectation is over both $\Omega$ and $\boldsymbol{Y}$. Next, note that by assumption $\boldsymbol{M}\in G$. Then, we have for any $\boldsymbol{X}\in G$

Moreover, from the definition of $\widehat{\boldsymbol{M}}$ we also have that $\widehat{\boldsymbol{M}}\in G$ and $\mathcal{L}_{\Omega,\boldsymbol{Y}}(\widehat{\boldsymbol{M}})\geq\mathcal{L}_{\Omega,\boldsymbol{Y}}(\boldsymbol{M})$. Thus

Applying Lemma 1, we obtain that with probability at least $1-C_{1}/(d_{1}+d_{2})$, we have

In this case, by rearranging and applying the fact that $\sqrt{d_{1}d_{2}}\leq d_{1}+d_{2}$, we obtain

Finally, we note that the KL divergence can easily be bounded below by the Hellinger distance:

This is a simple consequence of Jensen’s inequality combined with the fact that $1-x\leq-\log x$. Thus, from (26) we obtain

which establishes Theorem 6. Theorem 2 then follows by taking the limit as $\gamma\to\infty$.

The crux of the proof of the theorem involves bounding the quantity

The supremum is a function of the independent random variables $E_{i,j}\cdot\Delta_{i,j}$. To check the effectiveness of the method of bounded differences, suppose $E_{i,j}\Delta_{i,j}=1$ (for some $i,j$), and is replaced by $-1$. Then the empirical process can change by as much as $2\max_{X\in G}\left\|X\right\|_{\infty}=\alpha\sqrt{rd_{1}d_{2}}$, which is too large to yield an effective answer.

The fact that $\bar{\mathcal{L}}$ is not well-bounded in this way is also why the generalization error bounds of do not immediately generalize to provide results analogous to Theorem 2. To obtain this result, we must choose a loss function such that (27) reduces to the KL divergence, and unfortunately, this loss function is not well-behaved.

Finally, we also note that symmetrization and contraction principles, applied to bound empirical processes, are key tools in the theory of unquantized matrix completion. The empirical process that we must bound to prove Lemma 1 is a discrete analog of a similar process considered in Section 5 of .

We begin by noting that for any $h>0$, by using Markov’s inequality we have that

By a symmetrization argument (Lemma 6.3 in ),

where the $\varepsilon_{i,j}$ are i.i.d. Rademacher random variables and the expectation in the upper bound is with respect to both $\Omega$ and $\boldsymbol{Y}$ as well as with respect to the $\varepsilon_{i,j}$. To bound the latter term, we apply a contraction principle (Theorem 4.12 in ). By the definition of $L_{\gamma}$ and the assumption that $\|\widehat{\boldsymbol{M}}\|_{\infty}\leq\gamma$, both

are contractions that vanish at . Thus, up to a factor of $2$, the expected value of the supremum can only decrease when these are replaced by $X_{i,j}$ and $-X_{i,j}$ respectively. We obtain

where $\boldsymbol{E}$ denotes the matrix with entries given by $\varepsilon_{i,j}$, $\Delta_{\Omega}$ denotes the indicator matrix for $\Omega$ (so that $[\Delta_{\Omega}]_{i,j}=1$ if $(i,j)\in\Omega$ and 0 otherwise), and $\circ$ denotes the Hadamard product. Using the facts that the distribution of $\boldsymbol{E}\circ\boldsymbol{Y}$ is the same as the distribution of $\boldsymbol{E}$ and that $|\langle\boldsymbol{A},\boldsymbol{B}\rangle|\leq\|\boldsymbol{A}\|\|\boldsymbol{B}\|_{*}$, we have that

for some constant $C$. This in turn implies that

We first focus on the row sum $\sum_{j=1}^{d_{2}}\Delta_{i,j}$ for a particular choice of $i$. Using Bernstein’s inequality, for all $t>0$ we have

In particular, if we set $t\geq 6n/d_{1}$, then for each $i$ we have

where $W_{1},\ldots,W_{d_{1}}$ are i.i.d. exponential random variables.

Above, we have used the triangle inequality in the first line, followed by Jensen’s inequality in the second line. In the fifth line, (33), along with independence, allows us to introduce $\max_{i}W_{i}$. By standard computations for exponential random variables,

using the choice $h=\log(d_{1}+d_{2})\geq 1$ in the final line.

A similar argument bounds the column sums, and thus from (32) we conclude that

where the second and third inequalities both follow from Jensen’s inequality. Combining this with (29) and (31), we obtain

Plugging this into (28) we obtain that the probability in (28) is upper bounded by

provided that $C_{0}\geq 8(1+\sqrt{6})/e$, which establishes the lemma. ∎

A.2 Proof of Theorem 1

The proof of Theorem 1 follows immediately from Theorem 6 (with $\gamma=\alpha$) combined with the following lemma.

Let $f$ be a differentiable function and let $\left\|\boldsymbol{M}\right\|_{\infty},\|\widehat{\boldsymbol{M}}\|_{\infty}\leq\alpha$. Then

For any pair of entries $x=M_{i,j}$ and $y=\widehat{M}_{i,j}$, write

Using Taylor’s theorem to expand the quantity inside the square, for some $\xi$ between $x$ and $y$,

The lemma follows by summing across all entries and dividing by $d_{1}d_{2}$. ∎

A.3 Lower bounds

The proofs of our lower bounds each follow a similar outline, using classical information theoretic techniques that have also proven useful in the context of compressed sensing . At a high level, our argument involves first showing the existence of a set of matrices $\mathcal{X}$, so that for each $\boldsymbol{X}^{(i)}\neq\boldsymbol{X}^{(j)}\in\mathcal{X}$, $\|\boldsymbol{X}^{(i)}-\boldsymbol{X}^{(j)}\|_{F}$ is large. We will imagine obtaining measurements of a randomly chosen matrix in $\mathcal{X}$ and then running an arbitrary recovery procedure. If the recovered matrix is sufficiently close to the original matrix, then we could determine which element of $\mathcal{X}$ was chosen. However, Fano’s inequality will imply that the probability of correctly identifying the chosen matrix is small, which will induce a lower bound on how close the recovered matrix can be to the original matrix.

In the proofs of Theorems 3, 4, and 5, we will assume without loss of generality that $d_{2}\geq d_{1}$. Before providing these proofs, however, we first consider the construction of the set $\mathcal{X}$.

Let $K$ be defined as in (11), let $\gamma\leq 1$ be such that $\frac{r}{\gamma^{2}}$ is an integer, and suppose that $\frac{r}{\gamma^{2}}\leq d_{1}$. There is a set $\mathcal{X}\subset K$ with

For all $\boldsymbol{X}\in\mathcal{X}$, each entry has $|X_{i,j}|=\alpha\gamma$.

For all $\boldsymbol{X}^{(i)},\boldsymbol{X}^{(j)}\in\mathcal{X}$, $i\neq j$,

We use a probabilistic argument. The set $\mathcal{X}$ will by constructed by drawing

matrices $\boldsymbol{X}$ independently from the following distribution. Set $B=\frac{r}{\gamma^{2}}$. The matrix will consist of blocks of dimensions $B\times d_{2}$, stacked on top of each other. The entries of the first block (that is, $X_{i,j}$ for $(i,j)\in[B]\times[d_{2}]$) will be i.i.d. symmetric random variables with values $\pm\alpha\gamma$. Then $\boldsymbol{X}$ will be filled out by copying this block as many times as will fit. That is,

Now we argue that with nonzero probability, this set will have all the desired properties. For $\boldsymbol{X}\in\mathcal{X}$,

Further, because $\operatorname*{rank}{\boldsymbol{X}}\leq B$,

Thus $\mathcal{X}\subset K$, and all that remains is to show that $\mathcal{X}$ satisfies requirement 2.

For $\boldsymbol{X},\boldsymbol{W}$ drawn from the above distribution,

where the $\delta_{i,j}$ are independent $0/1$ Bernoulli random variables with mean $1/2$. By Hoeffding’s inequality and a union bound,

One can check that for $\boldsymbol{X}$ of the size given in (34), the right-hand side of the above tail bound is less than 1, and thus the event that $Z(\boldsymbol{X},\boldsymbol{W})>d_{2}B/4$ for all $\boldsymbol{X}\neq\boldsymbol{W}\in\mathcal{X}$ has non-zero probability. In this event,

where the second inequality uses the assumption that $d_{1}\geq B$ and the fact that $\lfloor x\rfloor\geq x/2$ for all $x\geq 1$. Hence, requirement (2) holds with nonzero probability and thus the desired set exists. ∎

A.3.2 Proof of Theorem 3

Before we prove Theorem 3, we will need the following lemma about the KL divergence.

Without loss of generality, we may assume that $x\leq y$. Indeed, $D(1-x\|1-y)=D(x\|y)$, and either $x\leq y$ or $1-x\leq 1-y$. Let $z=y-x$. A simple computation shows that

Thus, by Taylor’s theorem, there is some $\xi\in[0,z]$ so that

Since the right hand side is increasing in $\xi$, we may replace $\xi$ with $z$ and conclude

Now, for the proof of Theorem 3, we choose $\epsilon$ so that

where $C_{2}$ is an absolute constant to be specified later. We will next use Lemma 3 to construct a set $\mathcal{X}$, choosing $\gamma$ so that $\frac{r}{\gamma^{2}}$ is an integer and

We can make such a choice because $\epsilon\leq 1/32$ and $r\geq 4$. We verify that such a choice for $\gamma$ satisfies the requirements of Lemma 3. Indeed, since $\epsilon\leq 1/32$ and $\alpha\geq 1$ we have $\gamma\leq 1/4<1$. Further, we assume in the theorem that the right-hand side of (35) is larger than $Cr\alpha^{2}/d_{1}$ which implies that $r/\gamma^{2}\leq d_{1}$ for an appropriate choice of $C$.

Let $\mathcal{X}^{\prime}_{\alpha/2,\gamma}$ be the set whose existence is guaranteed by Lemma 3 with this choice of $\gamma$, and with $\alpha/2$ instead of $\alpha$. We will construct $\mathcal{X}$ by setting

Note that $\mathcal{X}$ has the same size as $\mathcal{X}^{\prime}_{\alpha/2,\gamma}$, i.e., $|\mathcal{X}|$ satisfies (34). $\mathcal{X}$ also has the same bound on pairwise distances

and every entry of $\boldsymbol{X}\in\mathcal{X}$ has

where $\alpha^{\prime}=(1-\gamma)\alpha$. Further, since for $\boldsymbol{X}^{\prime}\in\mathcal{X}^{\prime}_{\alpha/2,\gamma}$,

for $r\geq 4$ as in the theorem statement.

Now suppose for the sake of a contradiction that there exists an algorithm such that for any $\boldsymbol{X}\in K$, when given access to the measurements $\boldsymbol{Y}_{\Omega}$, returns an $\widehat{\boldsymbol{X}}$ such that

with probability at least $1/4$. We will imagine running this algorithm on a matrix $\boldsymbol{X}$ chosen uniformly at random from $\mathcal{X}$. Let

It is easy to check that if (37) holds, then $\boldsymbol{X}^{*}=\boldsymbol{X}$. Indeed, for any $\boldsymbol{X}^{\prime}\in\mathcal{X}$ with $\boldsymbol{X}^{\prime}\neq\boldsymbol{X}$, from (37) and (36) we have that

At the same time, since $\boldsymbol{X}\in\mathcal{X}$ is a candidate for $\boldsymbol{X}^{*}$, we have that

Thus, if (37) holds, then $\|\boldsymbol{X}^{*}-\widehat{\boldsymbol{X}}\|_{F}<\|\boldsymbol{X}^{\prime}-\widehat{\boldsymbol{X}}\|_{F}$ for any $\boldsymbol{X}^{\prime}\in\mathcal{X}$ with $\boldsymbol{X}^{\prime}\neq\boldsymbol{X}$, and hence we must have $\boldsymbol{X}^{*}=\boldsymbol{X}$. By assumption, (37) holds with probability at least $1/4$, and thus

We will show that this probability must in fact be large, generating our contradiction.

Because each entry of $\boldsymbol{Y}$ is independent,Note that here, to be consistent with the literature we are referencing regarding Fano’s inequality, $D$ is defined slightly differently than elsewhere in the paper where we would weight $D$ by $1/d_{1}d_{2}$.

Each term in the sum is either , $D(\alpha\|\alpha^{\prime})$, or $D(\alpha^{\prime}\|\alpha)$. By Lemma 4, all of these are bounded above by

and so, from the intermediate value theorem, for some $\xi\in[\alpha^{\prime},\alpha]$,

Using the assumption that $f^{\prime}(x)$ is decreasing for $x>0$ and the definition of $\alpha^{\prime}=(1-\gamma)\alpha$, we have

We now show that for appropriate values of $C_{0}$ and $C_{2}$, this leads to a contradiction. First suppose that $64n\epsilon^{2}\leq\beta_{\alpha^{\prime}}$. In this case we have

which together with (35) implies that $\alpha^{2}rd_{2}\leq 8$. If we set $C_{0}>8$ in (10), then this would lead to a contradiction. Thus, suppose now that $64n\epsilon^{2}>\beta_{\alpha^{\prime}}$. Then (40) simplifies to

Note $\beta$ is increasing as a function of $\alpha$ and $\alpha^{\prime}\geq 3\alpha/4$ (since $\gamma\leq 1/4$). Thus, $\beta_{\alpha^{\prime}}\geq\beta_{3\alpha/4}$. Setting $C_{2}\leq 1/512\sqrt{2}$ in (35) now leads to a contradiction, and hence (37) must fail to hold with probability at least $3/4$, which proves the theorem.

A.3.3 Proof of Theorem 4

for an absolute constant $C_{2}$ to be determined later. As in the proof of Theorem 3, we will consider running such an algorithm on a random element in a set $\mathcal{X}\subset K$. For our set $\mathcal{X}$, we will use the set whose existence is guaranteed by Lemma 3. We will set $\gamma$ so that $\frac{r}{\gamma^{2}}$ is an integer and

This is possible since $\epsilon\leq 1/4$ and $r,\alpha\geq 1$. One can check that $\gamma$ satisfies the assumptions of Lemma 3.

Now suppose that $\boldsymbol{X}\in\mathcal{X}$ is chosen uniformly at random, and let $\boldsymbol{Y}=\left.(\boldsymbol{X}+\boldsymbol{Z})\right|_{\Omega}$ as in the statement of the theorem. Let $\widehat{\boldsymbol{X}}$ be any estimate of $\boldsymbol{X}$ obtained from $\boldsymbol{Y}_{\Omega}$. We begin by bounding the mutual information $I(\boldsymbol{X};\widehat{\boldsymbol{X}})$ in the following lemma (which is analogous to [20, Equation 9.16]).

where $h$ denotes the differential entropy. Let $\mathbf{\xi}$ denote a matrix of i.i.d. $\pm 1$ entries. Then

and so, letting $\widetilde{\boldsymbol{X}}=\boldsymbol{X}\circ\mathbf{\xi}$,

Treating $\widetilde{\boldsymbol{X}}_{\Omega}+\boldsymbol{Z}_{\Omega}$ as a random vector of length $n$, we compute the covariance matrix as

We have that $h(\boldsymbol{Z}_{\Omega})=\frac{1}{2}\log\left((2\pi e)^{n}\sigma^{2n}\right)$, and so

Then the data processing inequality implies

We now proceed by using essentially the same argument as in the proof of Theorem 3. Specifically, we suppose for the sake of a contradiction that there exists an algorithm such that for any $\boldsymbol{X}\in K$, when given access to the measurements $\boldsymbol{Y}_{\Omega}$, returns an $\widehat{\boldsymbol{X}}$ such that

with probability at least $1/4$. As before, if we set

then we can show that if (42) holds, then $\boldsymbol{X}^{*}=\boldsymbol{X}$. Thus, if (42) holds with probability at least $1/4$ then

However, by Fano’s inequality, the probability that $\boldsymbol{X}\neq\widehat{\boldsymbol{X}}$ is at least

Plugging in $|\mathcal{X}|$ from Lemma 3 and $I(\boldsymbol{X};\widehat{\boldsymbol{X}})$ from Lemma 5, and using the inequality $\log(1+z)\leq z$, we obtain

Combining this with (43) and using the fact that $\gamma\leq 4\epsilon/\alpha$, we obtain

We now argue, as before, that this leads to a contradiction. Specifically, if $8n\epsilon^{2}/\sigma^{2}\leq 1$, then together with (41) this implies that $\alpha^{2}rd_{2}\leq 128$. If we set $C_{0}>128$ in (10), then this would lead to a contradiction. Thus, suppose now that $8n\epsilon^{2}/\sigma^{2}>1$, in which case we have

Thus, setting $C_{2}\leq 1/128$ in (41) leads to a contradiction, and hence (42) must fail to hold with probability at least $3/4$, which proves the theorem.

A.3.4 Proof of Theorem 5

The proof of Theorem 5 also mirrors the proof of Theorem 3. The main difference is the observation that the set constructed in Lemma 3 also works with the Hellinger distance. We begin as before by choosing $\epsilon$ so that

where $C_{2}$ is an absolute constant to be determined. Set $\gamma$ to be an integer so that

This is possible since by assumption $\alpha\geq 1$ and $\epsilon\leq\frac{c}{4}$. One can check that $\gamma$ satisfies the assumptions of Lemma 3.

As in the proof of Theorem 3, we will consider running such an algorithm on a random element in a set $\mathcal{X}\subset K$. For our set $\mathcal{X}$, we will use the set whose existence is guaranteed by Lemma 3. Note that since the Hellinger distance is bounded below by the Frobenius norm, we have that for all $\boldsymbol{X}^{(i)}\neq\boldsymbol{X}^{(j)}\in\mathcal{X}$,

Now suppose for the sake of a contradiction that there exists an algorithm such that for any $\boldsymbol{X}\in K$, when given access to the measurements $\boldsymbol{Y}_{\Omega}$, returns an $\widehat{\boldsymbol{X}}$ such that

with probability at least $1/4$. If we set

then we can show that if (45) holds, then $\boldsymbol{X}^{*}=\boldsymbol{X}$. Thus, if (45) holds with probability at least $1/4$ then

However, we may again apply Fano’s inequality as in (39). Using Lemma 4 we have

for some $|\xi|\leq\alpha\gamma$, where $L_{\alpha\gamma}$ is as in (5). By the assumption that $c^{\prime}<|f(x)|<1-c^{\prime}$ for $|x|<1$, and that

where $C^{\prime}=64c^{\prime}/(c^{2}(1-c^{\prime}))$. Thus, from (39), we have

We now argue once again that this leads to a contradiction. Specifically, if $C^{\prime}nL_{1}^{2}\epsilon^{2}\leq 1$, then together with (44) this implies that $\alpha^{2}rd_{2}\leq 128/c$. If we set $C_{0}>128/c$ in (10), then this would lead to a contradiction. Thus, suppose now that $C^{\prime}nL_{1}^{2}\epsilon^{2}>1$, in which case we have

Thus setting $C_{2}\leq c/32\sqrt{2C^{\prime}}$ in (44) leads to a contradiction, and hence (45) must fail to hold with probability at least $3/4$, which proves the theorem.

References