A Max-Norm Constrained Minimization Approach to 1-Bit Matrix Completion

Introduction

Matrix completion, which aims to recover a low-rank matrix from a subset of its entries, has been an active area of research in the last few years. It has a range of successful applications. In some real-life situations, however, the observations are highly quantized, sometimes even to a single bit and thus the standard matrix completion techniques do not apply. Take the Netflix problem as an example, the observations are the ratings of movies, which are quantized to the set of integers from $1$ to $5$. In the more extreme case such as recommender systems, only a single bit of rating standing for a “thumbs up” or “thumbs down” is recorded at each occurrence. Another example of applications is targeted advertising, such as the relevance of advertisements on Hulu. Each user who is watching TV shows on Hulu is required to answer yes/no to the question“Is this ad relevant to you?”. Noise effect should be considered since there are users who just click no to all the advertisements. In general, people would prefer to have advertisement catered to them, rather than to endure random advertisement. Targeted marketing that utilizes customer needs tends to serve better than random, scattershot advertisements. Similar idea has already been employed in mail system . Other examples from recommender systems include rating music on Pandora and posts on Reddit or MathOverflow, in which each observation consists of a single bit representing a positive or negative rating. Similar problem also arises in analyzing incomplete survey designs containing simple agree/disagree questions in the analysis of survey data, and distance matrix recovery in multidimensional scaling using binary and incomplete data . See for more detailed discussions.

To be more specific, consider an arbitrary unknown $d_{1}\times d_{2}$ target matrix $M^{*}$ with rank at most $r$. Suppose a subset $S=\{(i_{1},j_{1}),...,(i_{n},j_{n})\}$ of entries of a binary matrix $Y$ is observed, where the entries of $Y$ depend on $M^{*}$ in the following way:

Foygel and Srebro (2011) first used the max-norm for matrix completion under the uniform sampling distribution. Their results are direct consequences of a recent bound on the excess risk for a smooth loss function, such as the quadratic loss, with a bounded second derivative . Matrix completion under a non-degenerate random sampling model was considered by the present authors in an earlier paper . It was shown that the max-norm constrained minimization method is rate-optimal and it yields a more stable approximate recovery guarantee, with respect to the sampling distributions, than trace-norm based approaches.

Davenport, et al. (2012) analyzed 1-bit matrix completion under the uniform sampling model, where observed entries are assumed to be sampled randomly and uniformly. In such a setting, the trace-norm constrained approach has been shown to achieve minimax rate of convergence. However, in certain application such as collaborative filtering, the uniform sampling model is over idealized. In the Netflix problem, for instance, the uniform sampling model is equivalent to assuming all users are equally likely to rate every movie and all movies are equally likely to be rated by any user. In practice, inevitably some users are more active than others and some movies are more popular and thus rated more frequently. Therefore, the sampling distribution is in fact non-uniform. In such a scenario, Salakhutdinov and Srebro (2010) showed that the standard trace-norm relaxation can behave very poorly, and suggested to use a weighted variant of the trace-norm, which takes the sampling distribution into account. Since the true sampling distribution is most likely unknown and can only be estimated based on the locations of those entries that are revealed in the sample, what commonly used in practice is the empirically-weighted trace norm. Foygel, et al. (2011) provided rigorous recovery guarantees for learning with the standard weighted, smoothed weighted and smoothed empirically-weighted trace-norms. In particular, they gave upper bounds on excess error, which show that there is no theoretical disadvantage of learning with smoothed empirical marginals as compared to learning with smoothed true marginals.

In this paper we study matrix completion based on noisy 1-bit observations under a general (non-degenerate) sampling model using the max-norm as a convex relaxation for the rank. The rate of convergence for the max-norm constrained maximum likelihood estimate is obtained. A matching minimax lower bound is established under the general non-uniform sampling model using information-theoretical methods. The minimax upper and lower bounds together yield the optimal rate of convergence for the Frobenius norm loss. As a comparison with the max-norm constrained optimization approach, we also analyze the recovery guarantee of the weighted trace-norm constrained method in the setting of non-uniform sampling distributions. Our result includes an additional logarithmic factor, which might be an artifact of the proof technique. To sum up, the max-norm regularized approach indeed provides a unified and stable approximate recovery guarantee with respect to the sampling distributions, while previously used approaches are based on different variants of the trace-norm which may sometimes seem artificial to practitioners.

When the noise distribution is Gaussian or more generally log-concave, the negative log-likelihood function for $M$, given the measurements, is convex, hence computing the max-norm constrained maximum likelihood estimate is a convex optimization problem. The computational effectiveness of this method is also studied, based on a first-order algorithm developed in for solving convex programs involving a max-norm constraint, which outperforms the semi-definite programming method of Srebro, et al. (2004). It will be shown in Section 4 that the convex optimization problem can be implemented in polynomial time as a function of the sample size and the matrix dimensions.

The rest of the paper is organized as follows. Section 2 begins with the basic notation and definitions, and then states a collection of useful results on the matrix norms, Rademacher complexity and distances between matrices that will be needed throughout the paper. Section 3 introduces the 1-bit matrix completion model and the estimation procedure and investigates the theoretical properties of the estimator. Both minimax upper and lower bounds are established. The results show that the max-norm constraint maximum likelihood estimator is rate-optimal over the parameter space. Section 3 also gives a comparison of our results with previous work. Computational algorithms are discussed in Section 4, and numerical performance of the proposed algorithm is considered in Section 5. The proofs of the main results are given in Section 7. The paper is concluded with a brief discussion in Section 6.

Notations and Preliminaries

In this section, we introduce basic notation and definitions that will be used throughout the paper, and state some known results on the max-norm, trace-norm and Rademacher complexity that will be used repeatedly later.

Recall the definition $(\ref{maxnorm})$ of the max-norm, the trace-norm can be analogously defined in terms of matrix factorization as

By the elementary inequality $\|M_{m\times n}\|_{F}\leq\sqrt{m}\|M_{m\times n}\|_{2,\infty}$, we see that

Furthermore, as was noticed in Lee, et al. (2010), the max-norm, which is defined in $(\ref{maxnorm})$, is comparable with a trace-norm more precisely in the following sense :

2 Rademacher complexity

Considering matrices as functions from index pairs to entry values, a technical tool used in our proof involves data-dependent estimates of the Rademacher complexity of the classes that consist of low trace-norm and low max-norm matrices. We refer to Bartlett and Mendelson (2002) for a detailed introduction of this concept.

where $\varepsilon=(\varepsilon_{1},...,\varepsilon_{n})$ is a Rademacher sequence. The Rademacher complexity with respect to the distribution $\mathcal{P}$ is the expectation, over a sample $S$ of $|S|$ points drawn i.i.d. according to $\mathcal{P}$, denoted by

The following properties regarding $\hat{R}_{S}(\mathcal{F})$ are useful.

If $\mathcal{F}\subseteq\mathcal{G}$, $\hat{R}_{S}(\mathcal{F})\leq\hat{R}_{S}(\mathcal{G})$.

$\hat{R}_{S}(\mathcal{F})=\hat{R}_{S}(\mbox{\rm conv}(\mathcal{F}))=\hat{R}_{S}(\mbox{\rm absconv}(\mathcal{F}))$, where $\mbox{\rm conv}(\mathcal{F})$ is the class of convex combinations of functions from $\mathcal{F}$, and $\mbox{\rm absconv}(\mathcal{F})$ denotes the absolutely convex hull of $\mathcal{F}$, that is, the class of convex combinations of functions from $\mathcal{F}$ and $-\mathcal{F}$.

In particular, we are interested in calculating the Rademacher complexities of the trace-norm and max-norm balls. To this end, define for any radius $R>0$ that

First, recall that any matrix with unit trace-norm is a convex combination of unit-norm rank-one matrices, and thus

where $K>0$ denotes a universal constant.

The unit max-norm ball, on the other hand, can be approximately characterized as a convex hull. Due to the Grothendieck’s inequality, it was shown in that

where $\mathcal{M}_{\pm}:=\{M\in\{\pm 1\}^{d_{1}\times d_{2}}:\mbox{rank}(M)=1\}$ is the class of rank-one sign matrices, and $K_{G}\in(1.67,1.79)$ is the Grothendieck’s constant. It is easy to see that $\mathcal{M}_{\pm}$ is a finite class with cardinality $|\mathcal{M}_{\pm}|=2^{d-1}$, $d=d_{1}+d_{2}$. For any $d_{1},d_{2}>2$ and any sample of size $2<|S|\leq d_{1}d_{2}$, the empirical Rademacher complexity of the unit max-norm ball is bounded by

3 Discrepancy

In order to get both upper and lower prediction error bounds on the weighted squared Frobenius norm between the proposed estimator, given by $(\ref{max-est})$ below, and the target matrix described via model $(\ref{1b})$, we will need the following two concepts of discrepancies between matrices as well as their connections. In particular, we will focus on element-wise notion of discrepancy between two $d_{1}\times d_{2}$ matrices $P$ and $Q$.

First, for two matrices $P$, $Q:[d_{1}]\times[d_{2}]\rightarrow^{d_{1}\times d_{2}}$, their Hellinger distance is given by

where $d_{H}^{2}(p;q)=(\sqrt{p}-\sqrt{q})^{2}+(\sqrt{1-p}-\sqrt{1-q})^{2}$ for $p,q\in$. Next, the Kullback-Leibler divergence between two matrices $P$, $Q:[d_{1}]\times[d_{2}]\rightarrow^{d_{1}\times d_{2}}$ is defined by

The relationship between the two “distances” is as follows. For any two scalars $p,q\in$, we have

which in turn implies that, for any two matrices $P$, $Q:[d_{1}]\times[d_{2}]\rightarrow^{d_{1}\times d_{2}}$,

The proof of $(\ref{dis0})$ is based on the Jensen’s inequality and an elementary inequality that $1-x\leq-\log x$ for any $x>0$.

Max-Norm Constrained Maximum Likelihood Estimate

In this section, we introduce the max-norm constrained maximum likelihood estimation procedure for 1-bit matrix completion and investigates the theoretical properties of the estimator. The results are also compared with other results in the literature.

Instead of assuming the uniform sampling distribution , here we allow a general sampling distribution $\Pi=\{\pi_{kl}\}$, satisfying $\sum_{(k,l)\in[d_{1}]\times[d_{2}]}\pi_{kl}=1$, according to which we make $n$ independent random choices of entries. The drawback of the setting is that, with fairly high probability, some entries will be sampled multiple times. Intuitively it would be more practical to assume that entries are sampled without replacement, or equivalently, to sample $n$ of the $d_{1}d_{2}$ binary entries observed with noise without replacing. Due to the requirement that the drawn entries be distinct, the $n$ samples are not independent. This dependence structure turns out to impede the technical analysis of the learning guarantees. To avoid this complication, we will use the i.i.d. approach as a proxy for sampling without replacement throughout this paper. As has been noted in , between sampling with and without replacement both in a uniform sense, that is, making $n$ independent uniform choices of entries versus choosing a set $S$ of entries uniformly at random over all subsets that consist of exactly $n$ entries, the latter is indeed as good as the former. See Sect. 7.4 below for more details.

Next we list three natural choices for $F$, or equivalently, for the distribution of $\{Z_{i,j}\}$.

(Logistic regression/Logistic noise): The logistic regression model is described by $(\ref{1b-md})$ with

and equivalently by $(\ref{1b})$ with $Z_{i,j}$ i.i.d. following the standard logistic distribution.

(Probit regression/Gaussian noise): The probit regression model is described by $(\ref{1b-md})$ with

where $\Phi$ denotes the cumulative distribution function of $N(0,1)$, and equivalently by $(\ref{1b})$ with $Z_{i,j}$ i.i.d. following $N(0,\sigma^{2})$.

(Laplace noise): Another interesting case is that $Z_{i,j}$’s are i.i.d. Laplace noise (Laplace$(0,b)$), with

Davenport, et al. (2012) have focused on approximately low-rank matrices recovery by considering the following class of matrices

2 Max-norm constrained maximum likelihood estimate

Now, given a collection of observations $Y_{S}=\{Y_{i_{t},j_{t}}\}_{t=1}^{n}$ from the observation model $(\ref{1b-md})$, the negative log-likelihood function can be written as

Then we consider estimating the unknown $M^{*}\in K_{\max}(\alpha,R)$ by maximizing the empirical likelihood function subject to a max-norm constraint, i.e.,

The optimization procedure requires that all the entries of $M_{0}$ are bounded in absolute value by a pre-defined constant $\alpha$. This condition is reasonable while also critical in approximate low-rank matrix recovery problems by controlling the spikiness of the solution. Indeed, the measure of the “spikiness” of matrices is much less restrictive than the incoherence conditions imposed in exact low-rank matrix recovery. See, e.g. .

When the noise distribution is log-concave so that the log-likelihood is a concave function, the max-norm constrained minimization problem $(\ref{max-est})$ is a convex program and we recommend a fast and efficient algorithm developed in for solving large-scale optimization problems that incorporate the max-norm. We will show in Section 4 that the convex optimization problem $(\ref{max-est})$ can indeed be implemented in polynomial time as a function of the sample size $n$ and the matrix dimensions $d_{1}$ and $d_{2}$.

3 Upper bounds

To establish an upper bound on the prediction error of estimator $\hat{M}_{\max}$ given by $(\ref{max-est})$, we need the following assumption on the unknown matrix $M^{*}$ as well as the regularity conditions on the function $F$ in $(\ref{1b-md})$.

Condition U: Assume that there exist positive constants $R$ and $\alpha$ such that

$F$ and $F^{\prime}$ are non-zero in $[-\alpha,\alpha]$, and

In particular under condition $(U2)$, the quantity

is well-defined. As prototypical examples, we specify below the quantities $L_{\alpha}$, $\beta_{\alpha}$ and $U_{\alpha}$ in the cases of Logistic, Gaussian and Laplace noise:

(Logistic regression/Logistic noise): For $F(x)=e^{x}/(1+e^{x})$, we have

(Probit regression/Gaussian noise): For $F(x)=\Phi(x/\sigma)$, straightforward calculations show that

(Laplace noise): For a Laplace$(0,b)$ distribution function, we have

Now we are ready to state our main results concerning the recovery of an approximately low-rank matrix $M^{*}$ using the max-norm constrained maximum likelihood estimate. We write hereafter $d=d_{1}+d_{2}$ for brevity.

Suppose that Condition U holds and assume that the training set $S$ follows a general weighted sampling model according to the distribution $\Pi$. Then there exists an absolute constant $C$ such that, for a sample size $2<n\leq d_{1}d_{2}$ and for any $\delta>0$, the minimizer $\hat{M}_{\max}$ of the optimization program $(\ref{max-est})$ satisfies

with probability at least $1-\delta$. Here $\|\cdot\|_{\Pi}$ denotes the weighted Frobenius norm with respect to $\Pi$, i.e.,

While using the trace-norm to study this general weighted sampling model, it is common to assume that each row and column is sampled with positive probability (Nagahban and Wainwright, 2012; Klopp, 2012), though in some applications this assumption does not seem realistic. More precisely, assume that there exists a positive constant $\mu\geq 1$ such that

Then, under condition $(\ref{ass1})$ and the conditions of Theorem 3.1,

holds with probability at least $1-4/d$, where $C>0$ denotes an absolute constant.

Klopp (2012) studied the problem of standard matrix completion with noise, also in the case of general sampling distribution, using the trace-norm penalized approach. However, the Assumption 1 therein requires that the distribution $\pi_{kl}$ over entries is bounded from above, which is quite restrictive especially in the Netflix problem. It is worth noticing that this upper bound condition on sampling distribution is not required in both results $(\ref{mn-up})$ and $(\ref{mn-up2})$.

it was shown that for any $\delta\in(0,1)$ and a sample size $2<n\leq d_{1}d_{2}$,

holds with probability greater than $1-\exp(-d)-\delta$, where $C^{\prime}>0$ is a universal constant.

In 1-bit observations case when $Z_{i,j}\stackrel{{\scriptstyle i.i.d.}}{{\sim}}N(0,\sigma^{2})$, it is equivalent that the function $F$ in model $(\ref{1b-md})$ is given by $F(\cdot)=\Phi(\cdot/\sigma)$. According to $(\ref{prob-ex})$, we have

holds with probability at least $1-\delta$.

Comparing the upper bounds in $(\ref{TZ12.ubd})$ and $(\ref{1-bit.ubd})$ and note that $\alpha\vee\sigma\leq\alpha+\sigma\leq 2(\alpha\vee\sigma)$, we see that there is no essential loss of recovery accuracy by discretizing to binary measurements as long as $\frac{\alpha}{\sigma}$ is bounded by a constant . On the other hand, as the signal-to-noise ratio $\frac{\alpha}{\sigma}\geq 1$ increases, the error bounds deteriorate significantly. In fact, the case $\alpha\gg\sigma$ essentially amounts to the noiseless setting, in which it is impossible to recover $M^{*}$ based on any subset of the signs of its entries.

4 Information-theoretic lower bounds

We now establish minimax lower bounds by using information-theoretic techniques. The lower bounds given in Theorem 3.2 below show that the rate attained by the max-norm constrained maximum likelihood estimator is optimal up to constant factors.

Assume that $F^{\prime}(x)$ is decreasing and $\frac{F(x)(1-F(x))}{(F^{\prime}(x))^{2}}$ is increasing for $x>0$, and let $S$ be any subset of $[d_{1}]\times[d_{2}]$ with cardinality $n$. Then, as long as the parameters $(R,\alpha)$ satisfy

the minimax risk for estimating $M$ over the parameter space $K_{\max}(\alpha,R)$ satisfies

In fact, the lower bound $(\ref{minimax-lbd0})$ is a special case of the following general result, which will be proved in Sect. 7.2. Let $\gamma^{*}>0$ be the solution of the following equation

Then the minimax risk for estimating $M$ over the parameter space $K_{\max}(\alpha,R)$ satisfies

To see the existence of $\gamma^{*}$ defined above, setting

then it is easy to see that $h(\gamma)$ is strictly increasing and $g(\gamma)$ is decreasing for $\gamma\in(0,1)$ with $h(0)=0$ and $g(0)>0$. Therefore, equation $(\ref{ga*})$ has a unique solution $\gamma^{*}\in(0,\frac{1}{2}]$, i.e. $h(\gamma^{*})=g(\gamma^{*})$.

Assume that $\mu$ and $\alpha$ are bounded above by universal constants and let the function $F$ be fixed, so that both $L_{\alpha}$ and $\beta_{\alpha}$ are bounded. Also notice that $\beta_{(1-\gamma^{*})\alpha}\geq\beta_{\alpha/2}$ since $\gamma^{*}\leq 1/2$. Then comparing the lower bound $(\ref{minimax-lbd})$ with the upper bound $(\ref{mn-up2})$ shows that if the sample size $n\geq\frac{R^{2}\beta_{\alpha/2}}{4\alpha^{4}}(d_{1}+d_{2})$, the optimal rate of convergence is $R\sqrt{\frac{d_{1}+d_{2}}{n}}$, i.e.

and the max-norm constrained maximum likelihood estimate (3.10) is rate-optimal. If the target matrix $M^{*}$ is known to have rank at most $r$, we can take $R=\alpha\sqrt{r}$, such that the requirement here on the sample size $n\geq\frac{\beta_{\alpha/2}}{4\alpha^{2}}r(d_{1}+d_{2})$ is weak and the optimal rate of convergence becomes $\alpha\sqrt{\frac{r(d_{1}+d_{2})}{n}}$.

5 Comparison to prior work

In this paper, we study a matrix completion model proposed in , in which it is assumed that a binary matrix is observed at random from a distribution parameterized by an unknown matrix which is (approximately) low-rank. It is noteworthy that some earlier papers on collaborative filtering or matrix completion, including Srebro, et al. (2004) and references therein, also dealt with binary observations that are assumed to be noisy versions of the underlying matrix, in Logistic or Bernoulli conditional model. The goal there is to predict directly the quantized values, or equivalently, to reconstruct the sign matrix, instead of the underlying real values, therefore the non-identifiability issue could be avoided.

We next turn to a detailed comparison of our results for 1-bit matrix completion to those obtained in , also for approximately low-rank matrices. Using the trace-norm as a proxy to rank, Davenport, et al. (2012) have studied 1-bit matrix completion under the uniform sampling distribution over the parameter space

for some $\alpha>0$ and $r\leq\min\{d_{1},d_{2}\}$ is a positive integer. To recover the unknown $M^{*}\in K_{*}(\alpha,r)$, given a collection of observations $Y_{S}$ where $S$ follows a Bernoulli model, i.e. every entry $(k,l)\in[d_{1}]\times[d_{2}]$ is observed independently with equal probability $\frac{n}{d_{1}d_{2}}$, they propose the following trace-norm constrained MLE

and prove that for a sample size $n\geq d\log(d)$, $d=d_{1}+d_{2}$, with high probability,

Comparing to $(\ref{mn-up2})$ with $R=\alpha\sqrt{r}$, it is easy to see that under the uniform sampling model, the error bounds in (rescaled) Frobenius norm for the two estimates $\hat{M}_{\max}$ and $\hat{M}_{\rm tr}$ are of the same order. Moreover, Theorem 3 in and Theorem 3.2, respectively, provide lower bounds showing that both $\hat{M}_{\rm tr}$ and $\hat{M}_{\max}$ achieve the minimax rate of convergence for recovering approximately low-rank matrices over the parameter spaces $K_{*}(\alpha,r)$ and $K_{\max}(\alpha,R)$ respectively.

As mentioned in the introduction, the uniform sampling distribution assumption is restrictive and not valid in many applications including the well-known Netflix problem. When the sampling distribution is non-uniform, it was shown in Salakhutdinov and Srebro (2010) that the standard trace-norm regularized method might fail, specifically in the setting where the row and column marginal distributions are such that certain rows or columns are sampled with very high probabilities. Moreover, it was proposed to use a weighted variant of the trace-norm, which incorporates the knowledge of the true sampling distribution in its construction, and showed experimentally that this variant indeed leads to superior performance. Using this weighted trace-norm, Negahban and Wainright (2012) provided theoretical guarantees on approximate low-rank matrix completion in general sampling case while assuming that each row and column is sampled with positive probability (See condition $(\ref{ass1})$). In addition, requiring that the probabilities to observe an element from any row or column are of order $O((d_{1}\wedge d_{2})^{-1})$, Klopp (2012) analyzed the performance of the trace-norm penalized estimators, and provided near-optimal (up to a logarithmic factor) bounds which are similar to the bounds in this paper.

Next we provide an analysis of the performance of the weighted trace-norm in 1-bit matrix completion. Given the knowledge of the true sampling distribution, we establish an upper bound on the error in recovering $M^{*}$, which comparing to $(\ref{DP-ubd})$, includes an additional $\log^{1/2}(d)$ factor. We do not rule out the possibility that this logarithmic factor might be an artifact of the technical tools used in proof described below. The proof in for the trace-norm regularization in uniform sampling case may also be extended to the weighted trace-norm method under the general sampling model, by using the matrix Bernstein inequality instead of Seginer’s theorem. The extra logarithmic factor, however, is still inevitable based on this argument. We will not pursue the details in this paper.

Given a sampling distribution $\Pi=\{\pi_{kl}\}$ on $[d_{1}]\times[d_{2}]$, define its row- and column-marginals as

respectively. Under the condition $(\ref{ass1})$, we have

As in , consider the following weighted trace-norm with respect to the distribution $\Pi$:

where $(M_{w})_{k,l}:=\sqrt{\pi_{k\cdot}\pi_{\cdot l}}M_{k,l}$. Notice that if $M$ has rank at most $r$ and $\|M\|_{\infty}\leq\alpha$, then

Analogous to the previous studied class $K_{*}(\alpha,r)$ containing the low trace-norm matrices, define

and consider estimating the unknown $M^{*}\in K_{\Pi,*}$ by solving the following optimization problem:

The following theorem states that the weighted trace-norm regularized approach can be nearly as good as the max-norm regularized estimator (up to logarithmic and constant factors), under a general weighted sampling distribution. The theoretical performance of the weighted trace-norm is first studied by Foygel, et al. (2011) in the standard matrix completion problems under arbitrary sampling distributions.

Suppose that Condition U holds but with $M^{*}\in K_{\Pi,*}$, assume that the training set $S$ follows a general weighted sampling model according to the distribution $\Pi$ satisfying $(\ref{ass1})$. Then there exists an absolute constant $C>0$ such that, for a sample size $n\geq\mu\min\{d_{1},d_{2}\}\log(d)$ and any $\delta>0$, the minimizer $\hat{M}_{w,tr}$ of the optimization program $(\ref{weighted})$ satisfies

Since the construction of weighted trace-norm $\|\cdot\|_{w,*}$ highly depends on the underlying sampling distribution which is typically unknown in practice, the constraint $M^{*}\in K_{\Pi,*}$ seems to be artificial. The max-norm constrained approach, on the contrary, does not require the knowledge of the exact sampling distribution and the error bound in weighted Frobenius norm, as shown in $(\ref{mn-up})$, holds even without prior assumption on $\Pi$, e.g., $(\ref{ass1})$.

To clarify the major difference between the principles behind $(\ref{DP-ubd})$ and $(\ref{wei-up})$, we remark that one of the key technical tools used in is a bound of Seginer (2000) on the spectral norm of a random matrix with i.i.d. zero mean entries (corresponding to the uniform sampling distribution), i.e. for any $h\leq 2\log(\max\{d_{1},d_{2}\})$,

where $a_{k\cdot}$ (resp. $a_{\cdot l}$) denote the rows (resp. columns) of $A$ and $K$ is a universal constant. Under the non-uniform sampling model, we will deal with a matrix with independent entries that are not necessarily identically distributed, to which case an alternative result of Latala (2005) can be applied, i.e.

or instead, resorting to the matrix Bernstein inequality. Using either inequality would thus bring an additional logarithmic factor, appeared in $(\ref{wei-up})$.

It is also worth noticing that though the sampling distribution is not known exactly in practice, its empirical analogues are expected to be stable enough as an alternative. According to Forgel, et al. (2011), given a random sample $S=\{(i_{t},j_{t})\}_{t=1}^{n}$, consider the empirical marginals

as well as the smoothed empirical marginals

The smoothed empirically-weighted trace-norm $\|\cdot\|_{\check{w},*}$ can be defined in the same spirit as in the definition $(\ref{wtn})$ of weighted trace-norm, only with $\{\pi_{ij}\}$ replaced by $\{\check{\pi}_{ij}\}$. Then the unknown matrix can be estimated via regularization on the $\check{\pi}$-weighted trace-norm, that is,

Adopting [9, Theorem 4] to the current 1-bit problem will lead to a learning guarantee similar to $(\ref{wei-up})$.

Computational Algorithm

Problems of the form $(\ref{max-est})$ can now be solved using a variety of algorithms, including interior point method , Frank-Wolfe-type algorithm and projected gradient method . The first two are convex methods with guaranteed convergence rates to the global optimum, though can be slow in practice and might not scale to matrices with hundreds of rows or columns. We describe in this section a simple first order method due to Lee, et al. (2010), which is a special case of a projected gradient algorithm for solving large-scale convex programs involving the max-norm. This method is non-convex, but as long as the size of the problem is large enough, it is guaranteed that each local minimum is also a global optimum, due to Burer and Monteiro (2003).

Then the global optimum of $(\ref{max-est})$ is equal to that of

where $\tau>0$ is a stepsize parameter and $t=0,1,2,...$. Next, we project $(U(\tau),V(\tau))$ onto $\mathcal{M}_{k}(R)$ according to

otherwise we keep it still. The resulting update is then denoted by $(U^{t+1},V^{t+1})$.

It is important to note that the choice of $k$ must be large enough, at least as big as the rank of $M^{*}$. Suppose that, before solving $(\ref{max-est})$, we know that the target matrix $M^{*}$ has rank at most $r^{*}$. Then it is best to solve $(\ref{est2'})$ for $k=r^{*}+1$ in the sense that, if we choose $k\leq r^{*}$, then $(\ref{est2'})$ is not equivalent to $(\ref{max-est})$, and if we take $k>r^{*}+1$, then we would be solving a larger program than necessary. In practice, we do not know the exact value of $r^{*}$ in advance. Nevertheless, motivated by Burer and Monteiro (2003), we suggest the following scheme to solve the problem which avoids solving $(\ref{est2'})$ for $r\gg r^{*}$:

Choose an initial small $k$ and compute a local minimum $(U,V)$ of $(\ref{est2'})$, using above projected gradient method.

where $S\subset[d_{1}]\times[d_{2}]$ is a training set of row-column indices, $u_{i}$ and $v_{j}$ denote the $i$-th row of U and $j$-th row of V, respectively. The stochastic gradient method says that at $t$-th iteration, we only need to pick one training pair $(i_{t},j_{t})$ at random from $S$, then update $g(u_{i_{t}}^{T}v_{j_{t}};Y_{i_{t},j_{t}})$ via the previous procedure. More precisely, if $\|u_{i_{t}}\|_{2}^{2}>R$, we project it back so that $\|u_{i_{t}}\|_{2}^{2}=R$, otherwise we do not make any change (do the same for $v_{j_{t}}$). Next, if $|u_{i_{t}}^{T}v_{j_{t}}|>\alpha$, replace $u_{i_{t}}$ and $v_{i_{t}}$ with $\sqrt{\alpha}u_{i_{t}}/|u_{i_{t}}^{T}v_{j_{t}}|^{1/2}$ and $\sqrt{\alpha}v_{i_{t}}/|u_{i_{t}}^{T}v_{j_{t}}|^{1/2}$ respectively, otherwise we keep everything still. At the $t$-th iteration, we do not need to consider any other rows of $U$ and $V$. This simple algorithm could be computationally as efficient as optimization with the trace-norm.

Numerical results

In this section, we report the simulation results for low-rank matrix recovery based on 1-bit observations. In all cases presented below, we solved the convex program $(\ref{est2'})$ by using our implementation in MATLAB of the projected gradient algorithm proposed in Sect. 4 for a wide range of values of the step-size parameter $\tau$.

We first consider a rank-$2$, $d\times d$ target matrix $M^{*}$ with eigenvalues $\{d/\sqrt{2},d/\sqrt{2},0,...,0\}$, so that $\|M^{*}\|_{F}/d=1$. We choose to work with the Gaussian conditional model under uniform sampling. Let $Y_{S}$ be the noisy binary observations with $S=\{(i_{1},j_{1}),...,(i_{t},j_{t})\}$, that is, for $(i,j)\in S$,

where $\Omega^{+}=\{(i,j)\in S:Y_{i,j}=1\}$ and $\Omega^{-}=\{(i,j)\in S:Y_{i,j}=-1\}$. In Figure 1, averaging the results over 20 repetitions, we plot the squared Frobenius norm of the error (normalized by the dimension) $\|\hat{M}-M^{*}\|_{F}^{2}/d^{2}$ versus a range of sample sizes $s=|S|$, with the noise level $\sigma$ taken to be $\alpha/2$, for three different matrix sizes, $d\in\{80,120,160\}$. Naturally, in each case, the Frobenius error decays as $s$ increases, although larger matrices require larger sample sizes, as reflected by the upward shift of the curves as $d$ is increased.

Next, we compare the performance of the max-norm based regularization with that of the trace-norm using the same criterion as in . More specifically, the target matrix $M^{*}$ is constructed at random by generating $M=LR^{T}$, where $L$ and $R$ are $d\times r$ matrices with i.i.d. entries drawn from Uniform $[-1/2,1/2]$, so that rank$(M^{*})=r$. It is then scaled such that $\|M^{*}\|_{\infty}=1$, while in the last case, $M^{*}$ is formed such that $\|M^{*}\|_{F}/d=1$. As before, we focus on the Gaussian conditional model but with noise level $\sigma$ varies from $10^{-3}$ to $10$, and set $d=500$, $r=1$ and $s=0.15d^{2}$, which is exactly the same case studied in . We plot in Figure 2 the squared Frobenius norm of the error (normalized by the norm of the underlying matrix $M^{*}$) over a range of different values of noise level $\sigma$ on a logarithmic scale. As evident in Figure 2, the max-norm based regularization performs slightly but consistently better than the trace-norm, except on the one point where $\sigma=\log_{10}(0.25)$. Also, we see that for both methods, the performance is poor when the noise is either too little or too much.

In the third experiment, we consider matrices with dimension $d=200$ and choose a moderate level of noise, that is, $\sigma=\log_{10}(-0.75)$, according to previous experiences. Figure 3 plots the relative Frobenius norm of the error versus the sample size $s$ for three different matrix ranks, $r\in\{3,5,10\}$. Indeed, larger rank means larger intrinsic dimension of the problem, and thus increases the difficulty of any reconstruction procedure.

Discussion

This paper studies the problem of recovering a low-rank matrix based on highly quantized (to a single bit) noisy observation of a subset of entries. The problem was first formulated and studied by Davenport, et al. (2012), where the authors consider approximately low-rank matrices in terms that the singular values belong to a scaled Schatten-1 ball. When the infinity norm of the unknown matrix $M^{*}$ is bounded by a constant and its entries are observed uniformly in random, they show that $M^{*}$ can be recovered from binary measurements accurately and efficiently.

Our theory, on the other hand, focuses on approximately low-rank matrices in the sense that unknown matrix belongs to certain max-norm ball. The unit max-norm ball is nearly the convex hull of rank-1 matrices whose entries are bounded in magnitude by 1, thus is a natural convex relaxation of low-rank matrices, particularly with bounded infinity norm. Allowing for non-uniform sampling, we show that the max-norm constrained maximum likelihood estimation is rate-optimal up to a constant factor, and that the corresponding convex program may be solved efficiently in polynomial time. An interesting question naturally arises that whether it is possible to push the theory further to cover exact low-rank matrix completion from noisy binary measurements.

In our previous work , we suggest to use max-norm constrained least square estimation to study standard matrix completion (based on noisy observations) under a general sampling model. Similar errors bounds are obtained, which are tight to within a constant. Comparing both results in the case of Gaussian noise demonstrates that as long as the signal-to-noise ratio remains constant, almost nothing is lost by quantizing to a single bit.

Proofs

Since $\hat{M}_{\max}$ is optimal and $M^{*}$ is feasible to the optimization problem $(\ref{max-est})$, we have

Since $M^{*}$ has a fixed value which does not depend on $S$, the empirical likelihood term $\mathcal{D}_{S}(M^{*};Y)$ is an unbiased estimator of $\mathcal{D}_{\Pi}(M^{*};Y)$, i.e.

which in turn implies that that with probability at least $1-\delta$ over choosing a subset $S$ according to $\Pi$,

Moreover, observe that the left-hand side of $(\ref{up-1})$ is equal to

This, combined with $(\ref{riskbd1})$, $(\ref{up-2})$ and $(\ref{up-1})$ implies that for any $\delta>0$, the following inequality holds with probability at least $1-\delta$ over $S$:

This, together with $(\ref{dis})$ and Lemma 7.1 below gives $(\ref{mn-up})$.

Let $F$ be an arbitrary differentiable function, and $s,t$ are two real numbers satisfying $|s|,|t|\leq\alpha$. Then

The proof of Theorem 3.1 is now completed.

2 Proof of Theorem 3.2

The proof for the lower bound follows an information-theoretic method based on Fano’s inequality , as used in the proof of Theorem 3 in . To begin with, we have the following lemma which guarantees the existence of a suitably large packing set for $K_{\max}(\alpha,R)$ in the Frobenius norm. The proof follows from Lemma 3 of with a simple modification, see, e.g., the proof of Lemma 3.1 in .

Let $r=(R/\alpha)^{2}$ and $\gamma\leq 1$ be such that $\frac{r}{\gamma^{2}}\leq\min(d_{1},d_{2})$ is an integer. There exists a subset $\mathcal{S}(\alpha,\gamma)\subset K_{\max}(\alpha,R)$ with cardinality

For any $N\in\mathcal{S}(\alpha,\gamma)$, rank$(N)\leq\frac{r}{\gamma^{2}}$ and $N_{k,l}\in\{\pm\gamma\alpha/2\}$, such that

For any two distinct $N^{k},N^{l}\in\mathcal{S}(\alpha,\gamma)$,

Then we construct the packing set $\mathcal{M}$ by setting

provided that $r\geq 4$. Therefore, $\mathcal{M}$ is indeed a $\delta$-packing of $K_{\max}(\alpha,R)$ in the Frobenius metric with

i.e. for any two distinct $M,M^{\prime}\in\mathcal{M}$, we have $\|M-M^{\prime}\|_{F}\geq\delta$.

Next, a standard argument (e.g. ) yields a lower bound on the $\|\cdot\|_{F}$-risk in terms of the error in a multi-way hypothesis testing problem. More concretely,

where $I(M^{\star};Y_{S})$ denotes the mutual information between the random parameter $M^{\star}$ in $\mathcal{M}$ and the observation matrix $Y_{S}$. Following the proof of Theorem 3 in , we could bound $I(M^{\star};Y_{S})$ as follows:

where the last inequality holds provided that $F^{\prime}(x)$ is decreasing on $(0,\infty)$. Substituting this into the Fano’s inequality $(\ref{fi})$ yields

Recall that $\gamma^{*}>0$ solves the equation $(\ref{ga*})$, i.e.

Requiring $\frac{64\log(2)(\gamma^{*})^{2}}{d_{1}\vee d_{2}}\leq r\leq(d_{1}\wedge d_{2})(\gamma^{*})^{2}$, which is guaranteed by $(\ref{r})$, to ensure that this probability is least $1/4$. Consequently, we have

which in turn implies $(\ref{minimax-lbd})$.

3 Proof of Theorem 3.3

The proof of Theorem 3.3 modifies the proof of Theorem 3.1, therefore we only outline the key steps in the following. Let $\{A_{1},...,A_{n}\}=\{(i_{1},j_{1}),...,(i_{n},j_{n})\}$ be independent random variables taking values in $[d_{1}]\times[d_{2}]$ according to $\Pi$, and recall that

According to and the proof of Theorem 3.1, it suffices to derive an upper bound on

where $\varepsilon_{t}$ are i.i.d. Rademacher random variables. Then it follows from $(\ref{tn.cvhull})$ that

Putting pieces together, we conclude that

which in turn yields that for any $\delta\in(0,1)$, inequality

holds with probability at least $1-\delta$, provided that $n\geq\mu\min\{d_{1},d_{2}\}\log(d)$.

4 An extension to sampling without replacement

In this paper, we have focused on sampling with replacement. We shall show here that in the uniform sampling setting, the results obtained in this paper continue to hold if the (binary) entries are sampled without replacement. Recall that in the proof of Theorem 3.1, we let $A_{1},...,A_{n}$ be random variables taking values in $[d_{1}]\times[d_{2}]$, $S=\{A_{1},...,A_{n}\}$ and assume the $A_{t}$’s are distributed uniformly and independently, i.e. $S\sim\Pi=\{\pi_{kl}\}$ with $\pi_{kl}\equiv\frac{1}{d_{1}d_{2}}$. The purpose now is to prove that the arguments remain valid when the $A_{t}$’s are selected without replacement, denoted by $S\sim\Pi_{0}$. In this notation, we have

By Lemma 3 in and $(\ref{riskbd1})$, for any $\delta>0$,

Acknowledgements

We would like to thank Yaniv Plan for helpful discussions and for pointing out the importance of allowing non-uniform sampling. A part of this work was done when the second author were visiting the Wharton Statistics Department of the University of Pennsylvania. He wishes to thank the institution and particularly the first author for their hospitality.

References