Matrix Completion via Max-Norm Constrained Optimization

Introduction

The max-norm was recently proposed as an alternative convex surrogate to the rank of the matrix. For collaborative filtering problems, the max-norm has been shown to be empirically superior to the trace-norm Srebro, Rennie and Jaakkola (2004). Foygel and Srebro (2011) 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 (Srebro, Sridharan and Tewari, 2010). Further, a max-norm constrained maximum likelihood method was considered by Cai and Zhou (2013) for one-bit matrix completion, where instead of observing real-valued entries of an unknown matrix one is only able to see binary outputs, i.e. yes/no, true/false, agree/disagree (Davenport et al., 2014). Theoretical guarantees are obtained in general non-uniform sampling models, and numerical studies show that the max-norm based approach is comparable to and sometimes slightly outperform the corresponding trace-norm method.

Matrix completion has been well analyzed in the uniform sampling model, where observed entries are assumed to be sampled randomly and uniformly. In such a setting, the trace-norm regularized approach has been shown to have good theoretical and numerical performance. However, in some applications such as collaborative filtering, the uniform sampling model is unrealistic. For example, in the Netflix problem, the uniform sampling model is equivalent to assuming all users are equally likely to rate each movie and all movies are equally likely to be rated by any user. From a practical point of view, invariably some users are more active than others and some movies are more popular and thus rated more frequently. Hence, the sampling distribution is in fact non-uniform in the real world. In such a setting, Salakhutdinov and Srebro (2010) showed that the standard trace-norm relaxation can sometimes behave poorly, and suggested a weighted trace-norm penalty, which incorporates the knowledge of true sampling distribution in its construction. 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, a practically available method relies on the empirically-weighted trace-norm (Foygel et al., 2011). It is also worth noticing that, when the sampling probabilities are bounded from below and above, the trace-norm penalized estimator is minimax optimal up to a logarithmic factor (Klopp, 2014). We refer to Fang et al. (2015b) for further numerical evaluations of the trace-norm regularized method under various non-uniform sampling schemes.

In this paper, we employ the max-norm as a convex relaxation for the rank to study matrix completion based on noisy observations in a general, unspecified sampling model. The rate of convergence for the max-norm constrained least squares estimator is obtained. Information-theoretical methods are used to establish a matching minimax lower bound in the general non-uniform sampling model. Together, the minimax upper and lower bounds yield the optimal rate of convergence for the Frobenius norm loss. It is shown that the max-norm regularized approach indeed provides a unified and robust approximate recovery guarantee with respect to sampling schemes. In the uniform sampling model as a special case, our results also show that the extra logarithmic factors appeared in the error rates obtained by Srebro, Sridharan and Tewari (2010) and Foygel and Srebro (2011) could be avoided after a careful analysis to match the minimax lower bound with the upper bound (see Theorems 3.1 and 3.3 and the discussions in Section 3).

The max-norm constrained minimization problem is a convex program. To solve general convex programs that involve either a max-norm constraint or a max-norm penalization, a first-order algorithm was proposed by Lee et al. (2010), which is computationally effective and outperforms the semi-definite programming (SDP) method of Srebro, Rennie and Jaakkola (2004). In principle, the method of Lee et al. (2010) is based on nonconvex relaxations. Therefore, their algorithm is only guaranteed to find a stationary point, and statistical properties of such solutions are difficult to analyze. Recently, Fang et al. (2015b) proposed a scalable algorithm based on the alternating direction of multipliers method to efficiently solve the max-norm constrained optimization problem with guaranteed rate of convergence to the global optimum. In summary, the max-norm constrained empirical risk minimization problem can indeed be implemented in polynomial time as a function of the sample size and matrix dimensions.

The remainder of the paper is organized as follows. After introducing basic notation and definitions, Section 2 collects a few useful results on the max-norm, trace-norm and Rademacher complexity that will be needed in the rest of the paper. Section 3 introduces the model and the estimation procedure and then investigates the theoretical properties of the estimator. Both minimax upper and lower bounds are given. The results show that the max-norm constrained minimization method achieves the optimal rate of convergence over the parameter space. Comparison with past work is also given. Computation and implementation issues are discussed in Section 4. A brief discussion is given in Section 5, and the proofs of the main results and key technical lemmas are placed in Section 6.

Notations and Preliminaries

In this section, we begin with some notation that will be used throughout the paper, and then collect some known results on the max-norm, trace-norm and Rademacher complexity that will be applied repeatedly later.

The factor of equivalence is the Grothendieck’s constant $K_{G}\in(1.67,1.79)$. Based on these properties, the max-norm regularization is expected to be more effective when dealing with uniformly bounded data (Lee et al., 2010).

Of the same spirit as the definition of the max-norm in (1.1), the trace-norm has the following equivalent characterization in terms of matrix factorizations:

See, for example, Srebro and Shraibman (2005). It is easy to see that

with $K_{G}\in(1.67,1.79)$ denoting the Grothendieck’s constant. Moreover, $\mathcal{M}_{\pm}$ is a finite class with cardinality $|\mathcal{M}_{\pm}|=2^{d-1}$, where $d=d_{1}+d_{2}$.

2 Rademacher complexity

A technical tool used in our analysis involves data-dependent estimates of the Rademacher and Gaussian complexities of a function class. We refer to Bartlett and Mendelson (2002) and Srebro and Shraibman (2005) for a detailed introduction of these concepts.

where $\boldsymbol{\varepsilon}=(\varepsilon_{1},\varepsilon_{2},\ldots,\varepsilon_{n})^{\intercal}$ is a Rademacher sequence. The Rademacher complexity with respect to a distribution $\mathcal{P}$ is the expectation, over an independent and identically distributed (i.i.d.) sample of $|S|$ points drawn from $\mathcal{P}$, denoted by

Replacing $\varepsilon_{1},\ldots,\varepsilon_{n}$ with independent standard normal variables $g_{1},\ldots,g_{n}$ leads to the definition of (empirical) Gaussian complexity.

Considering a matrix as a function from the index pairs to the entry values, Srebro and Shraibman (2005) obtained upper bounds on the Rademacher complexity of the unit balls under both the trace-norm and the max-norm. Specifically, for any $d_{1},d_{2}>2$ and any sample of size $2<|S|<d_{1}d_{2}$, the empirical Rademacher complexity of the max-norm unit ball is bounded by

Max-Norm Constrained Empirical Risk Minimization

for some $\sigma>0$. The noise variables $\xi_{t}$ are independent with zero mean and unit variance. By expressing the model as in (3.1), it is implicitly assumed that the noise on the entry is drawn independently each time.

When $\Pi$ corresponds to the uniform distribution, $\|A\|_{\Pi}=(d_{1}d_{2})^{-1/2}\|A\|_{F}$.

2 Max-norm constrained least squares estimator

Given a collection of observations $Y_{S}=\{Y_{i_{t}j_{t}}\}_{t=1}^{n}$ from the observation model $(\ref{mc-md})$, we estimate the unknown $M^{*}\in\mathcal{K}(\alpha,R)$ for some $\alpha,R>0$ by the minimizer of the empirical risk with respect the quadratic loss function

The minimization procedure requires that all the entries of $M^{*}$ are bounded in magnitude by a prespecified constant $\alpha$. This condition enforces that $M^{*}$ should not be too “spiky”, and a too large bound may jeopardize exactness of the estimation. See, for example, Koltchinskii, Lounici and Tsybakov (2011), Negahban and Wainwright (2012) and Klopp (2014). On the other hand, as argued in Lee et al. (2010), the max-norm regularization is expected to be more effective particularly for uniformly bounded data, which is our main motivation for using the max-norm constrained estimator.

Although the max-norm constrained minimization problem $(\ref{max-est})$ is a convex program, fast and efficient algorithms for solving large-scale optimization problems that incorporate the max-norm have only been developed recently in Lee et al. (2010) and Fang et al. (2015b). We will show in Section 4 that the convex optimization problem $(\ref{max-est})$ can be implemented in polynomial time as a function of the sample size $n$ and dimensions $d_{1}$ and $d_{2}$.

3 Upper bounds

In this section, we state our main results regarding the recovery of an approximately low-rank (low-max-norm) matrix $M^{*}$ using max-norm constrained empirical risk minimization.

Suppose that the noise sequence $\{\xi_{t}\}_{t=1}^{n}$ are independent sub-exponential random variables; that is, there is a constant $K>0$ such that

The parameters $\alpha,R>0$ are such that $M^{*}\in\mathcal{K}(\alpha,R)$. Then, for a sample size $n$ satisfying $d\leq n\leq d_{1}d_{2}$,

with probability greater than $1-2e^{-d}$, where $C>0$ is an absolute constant. If, in addition, assumption $(\ref{ass1})$ is satisfied, then for a sample size $n$ with $d\leq n\leq d_{1}d_{2}$,

holds with probability at least $1-2e^{-d}$.

The upper bounds given in Theorem 3.1 hold with high probability. The rate of convergence under expectation can be obtained as a direct consequence. More specifically, for a sample size $n$ with $d\leq n\leq d_{1}d_{2}$, we have

In view of the upper bound in (6.1), when the noise level $\sigma$ is comparable to or dominated by $\alpha$, the rate is of order $\alpha R\,(\frac{d}{n})^{1/2}$. To fully understand how the random noise affects the estimation accuracy particularly when $\sigma$ is much smaller than $\alpha$, we provide a complementary result in Theorem 3.2 which generalizes Theorem 9 in Foygel and Srebro (2011) to the general non-uniform sampling model.

Assume that the conditions of Theorem 3.1 are satisfied and $\sigma\leq\alpha$. Then,

holds with probability at least $1-2n^{-1}$ over a random sample of size $n$ satisfying $d\leq n\leq d_{1}d_{2}$, where $C>1$ is a constant.

An interesting consequence of Theorem 3.2 is that, in the noiseless case where $\sigma=0$ and a random subset of the entries of $M^{*}$ are perfectly observed, then for any prespecified tolerance level $\epsilon>0$, the target matrix $M^{*}$ can be approximately recovered in the sense that $\|\widehat{M}_{\max}-M^{*}\|_{\Pi}^{2}\leq\epsilon$ whenever the sample size n\gtrsim\max\big{\{}\frac{R^{2}d}{\epsilon}(\log n)^{3},\frac{\alpha^{2}}{\epsilon}(\log n)^{3/2}\big{\}}.

4 Information-theoretic lower bounds

To derive the lower bound, we assume that the sampling distribution $\Pi$ satisfies

for a positive constant $\mu\geq 1$. Clearly, when $\mu=1$, it amounts to say that the sampling distribution is uniform.

Suppose that the noise sequence $\{\xi_{t}\}_{t=1}^{n}$ are i.i.d. standard normal random variables, the sampling distribution $\Pi$ satisfies the condition $(\ref{ass2})$ and the quintuple $(n,d_{1},d_{2},\alpha,R)$ satisfies

Then the minimax $\|\cdot\|_{F}$-risk is lower bounded as

In particular, for a sample size $n\geq\frac{1}{\alpha^{2}\mu}R^{2}d$,

Assume that both $\nu$ and $\mu$, respectively appeared in $(\ref{ass1})$ and $(\ref{ass2})$, are bounded above by universal constants, then comparing the lower bound (3.14) with the upper bound $(\ref{mc-ubd3})$ shows that if the sample size $n>(\frac{R}{\alpha})^{2}d$, the optimal rate of convergence is $R\sqrt{d/n}$; that is,

and the max-norm constrained least-squares estimator (3.5) is rate-optimal. The requirement here on the sample size $n>(\frac{R}{\alpha})^{2}(d_{1}+d_{2})$ is weak. If, in addition, $d_{1}=d_{2}$, condition (3.12) is reduced to $\alpha^{2}d^{-1}\lesssim R^{2}\lesssim\sigma\alpha d$, which is a mild constraint since $R^{2}$ is of order $\alpha^{2}r_{0}$ in the exact low-rank case where $r_{0}=\mbox{rank}(M^{*})$.

The proof of Theorem 3.3 uses information-theoretic methods. A key technical tool for the proof is the following lemma which guarantees the existence of a suitably large packing set for $\mathcal{K}(\alpha,R)$ in the Frobenius norm.

Let $r=(\frac{R}{\alpha})^{2}$ and let $\gamma\leq 1$ be such that $r\leq\gamma^{2}(d_{1}\wedge d_{2})$ is an integer. Then, there exists a subset $\mathcal{M}\subseteq\mathcal{K}(\alpha,R)$ with cardinality

For any two distinct $M^{i},M^{j}\in\mathcal{M}$,

The proof of Lemma 3.1 is based on an adaptation of the arguments used to prove Lemma 3 in Davenport et al. (2014), which for self-containment, is given in Section 6.4.

5 Comparison to past work

We now compare the results established in this section with those known in the literature for matrix completion under uniform or general sampling schemes.

It is now well-known that the exact recovery of a low-rank matrix in the noiseless case requires the “incoherence conditions” on the target matrix $M^{*}$ (Candés and Recht, 2009; Candés and Tao, 2010; Recht, 2011; Gross, 2011). In this paper, we consider instead a general setting of approximately low-rank matrices, and prove that approximate recovery is still possible without enforcing exact structural assumptions.

Accordingly, define the weighted norms as

where $W_{{\rm r}}=d_{1}\cdot\mbox{diag}(\pi_{1\cdot},\ldots,\pi_{d_{1}\cdot})$ and $W_{{\rm c}}=d_{2}\cdot\mbox{diag}(\pi_{\cdot 1},\ldots,\pi_{\cdot d_{2}})$

Negahban and Wainwright (2012) proposed the following estimator of $M^{*}$ based on the trace-norm penalized minimization:

In the context of low-trace-norm (approximately low-rank) matrix recovery where the true matrix $M^{*}$ satisfies (3.17), they proved that for properly chosen $\lambda_{n}$ depending on $\sigma$ (see, e.g. Corollary 2 therein), there exist absolute positive constants $c_{1}$–$c_{3}$ such that

holds with probability at least $1-c_{2}\exp(-c_{3}\log d)$.

Since only a relatively small sample of the entries of $M^{*}$ is observed, these estimates may not be accurate enough. The max-norm constrained minimization approach, on the other hand, is proved (Theorem 3.1) to be effective in the presence of non-uniform sampling distributions. The method does not require either a product distribution or the knowledge of the exact true sampling distribution. From this point of view, the max-norm constrained method indeed yields a more robust approximate recovery guarantee, with respect to the sampling distributions.

Suppose that the noise sequence $\{\xi_{t}\}_{t=1}^{n}$ are i.i.d. $N(0,1)$ random variables and the sampling distribution $\Pi$ is uniform on $[d_{1}]\times[d_{2}]$. Then the following inequalities hold with probability at least $1-3d^{-1}$:

The optimum $\widehat{M}_{\max}$ to the convex program $(\ref{max-est})$ satisfies

The minimum $\widehat{M}_{{\rm tr}}$ to the SDP $(\ref{NW-est})$ with all weighted norms replaced by the standard ones and with a properly chosen $\lambda_{n}$ satisfies

The upper bound $(\ref{compare1})$ follows immediately from $(\ref{mc-ubd})$ in Theorem 3.1, and $(\ref{compare2})$ is a straightforward extension of Theorem 7 in Klopp (2014) on exact low-rank matrix recovery to the case of low-trace-norm matrix reconstruction. The proof is essentially the same and thus is omitted.

Foygel and Srebro (2011) analyzed the recovery guarantee for $\widehat{M}_{\max}$ based on an excess risk bound for empirical risk minimization with a smooth loss function recently developed in Srebro, Sridharan and Tewari (2010). Specifically, assuming a uniform sampling model with sub-exponential noise and that the target matrix $M^{*}\in\mathcal{K}(\alpha,R)$, they proved that with high probability,

After a more delicate analysis, our result shows that the additional logarithmic factors in $(\ref{compare3})$ purely arise from an artifact of the proof technique and thus can be avoided. Moreover, in view of the lower bounds given in Theorem 3.3, we see that the max-norm constrained least square estimator $\widehat{M}_{\max}$ achieves the optimal rate of convergence for recovering approximately low-rank matrices over the parameter space $\mathcal{K}(\alpha,R)$ under the Frobenius norm loss. To our knowledge, the best known rate for the trace-norm regularized estimator given in $(\ref{compare2})$ is near-optimal up to logarithmic factors in a minimax sense, over a larger parameter space $\mathcal{K}_{{\rm tr}}(\alpha,R)$.

5.2 Uniform/non-uniform sampling distributions

See, for example, Theorem 9 in Lee, Shraibman and Spalek (2008). As a result, by considering a max-norm penalized estimator that solves

all the possible marginal probabilities are taken into account, and therefore the solution is expected to be more robust with respect to the unknown sampling distributions.

Foygel et al. (2011) proved the error bound for the excess risk of the empirically-weighted trace-norm constrained estimator when the loss function is Lipschitz. It is interesting to investigate whether the results similar to those in Negahban and Wainwright (2012) hold for the empirically-weighted trace-norm constrained and penalized estimators when the quadratic loss function is used.

It is also worth noting that, under condition (3.6) and when the sampling distribution is nearly uniform in the sense that

for some constants $\nu,L\geq 1$, Klopp (2014) showed that the trace-norm penalized estimator

Computational Algorithms

Due to Srebro, Rennie and Jaakkola (2004), the max-norm of a $d_{1}\times d_{2}$ matrix $M$ can be computed via a semi-definite program:

Correspondingly, we can reformulate $(\ref{max-est})$ as the following SDP problem

where the objective function $f$ is given by

This SDP can be solved using standard interior-point methods, though are fairly slow and do not scale to matrices with large dimensions. For large-scale problems, an alternative factorization method based on $(\ref{eq1.1})$, as described below, is preferred (Lee et al., 2010).

where $\{(i_{1},j_{1}),\ldots,(i_{n},j_{n})\}\subseteq([d_{1}]\times[d_{2}])^{n}$ is a training set of row-column indices, $U_{i}$ and $V_{j}$ denote the $i$th row of $U$ and the $j$th row of $V$, respectively. This problem, however, is non-convex since it involves a constraint on all product factorizations $UV^{\intercal}$. When the size of the problem $k$ is large enough, Burer and Choi (2006) proved that this reformulated problem has no local minima. To solve this problem fast and efficiently, Lee et al. (2010) suggested the following first-order method.

where $\tau>0$ is a stepsize parameter. If $\|\widetilde{U}^{t+1}(\widetilde{V}^{t+1})^{\intercal}\|_{\infty}>\alpha$, we replace

otherwise we keep it still. Next, compute updates according to

2 An alternating direction method of multipliers based approach

The first-order algorithm described in Section 4.1 is computationally efficient and fast. However, (4.1) is in principle a non-convex optimization problem and thus the algorithm is only guaranteed to find a stationary point. Recently, an alternating direction method of multipliers (ADMM) based approached was proposed by Fang et al. (2015b) to solve the convex program (3.5) efficiently with strong theoretical guarantee. Furthermore, it was shown in Fang et al. (2015a) that the worst-case rate of convergence of the ADMM method is of order $1/t$, where $t$ denotes the iteration counter. We briefly summarize this ADMM approach here for the sake of readability.

In this notation, the problem (4.1) can be equivalently formulated as

More specifically, consider the augmented Lagrangian function of (4.5) that is given by

for $W\in\mathcal{P}$ and $X\in\mathcal{S}^{d}_{+}=\{S\in\mathcal{S}^{d}:S\succeq 0\}$, where $Z$ denotes the dual variable and $\rho>0$ is prespecified. The ADMM is used to solve (4.5) iteratively as follows: Initialize $(W^{0},X^{0},Z^{0})$ and $\rho>0$; at the $(t+1)$-th iteration, update $(W,X,Z)$ according to

3 Implementation

In order to estimate $R_{0}$ directly from a missing data matrix, it can be seen from $(\ref{eq2.2})$ that $\alpha_{0}\sqrt{r_{0}}$ is a sharp upper bound on $R_{0}$ and is more amenable to estimation. Fortunately, it is possible to convincingly specify $\alpha_{0}$ beforehand in many real-life applications. When dealing with the Netflix data, for instance, $\alpha_{0}$ can be chosen as the highest rating index; in the structure-from-motion problem, $\alpha_{0}$ depends on the range of the camera field of view, which in most cases is sufficiently large to capture the feature point trajectories. In case where the percentage of missing entries is low, the largest magnitude of the observed entries can be used as an alternative for $\alpha_{0}$.

As for $r_{0}$, we recommend the rank estimation approach recently developed in Juliá et al. (2011), which was shown to be effective in computer vision problems. Recall that in the structure-from-motion problem, each column of the data matrix corresponds a trajectory along the frames of a given feature point, and can be regarded as a signal vector with missing coordinates. Due to the rigidity of the moving objects, it was noted in Juliá et al. (2011) that the behavior of observed and missing data is the same and thus they both generate an analogous (frequency) spectral representation. Motivated by this observation, the proposed approach is based on the study of changes in frequency spectra on the initial matrix after missing entries are recovered.

In general, choosing the tuning parameter $R>0$ in (3.5) adaptively is a difficult problem. In the regression case, it can be done by the Scaled LASSO method (Sun and Zhang, 2012). It is unclear whether a similar approach would work for matrix completion problems. By convexity and strong duality, the optimization program in (3.5) is equivalent to

for a properly chosen $\lambda$. In fact, for any $R>0$ specified in (3.5), there exists a $\lambda>0$ such that the solutions to the two problems (3.5) and (4.2) coincide. In practice, we suggest to solve (4.2) using the ADMM method described in Section 4.2 with $\lambda$ obtained via cross-validation, in a way similarly to that for LASSO or the trace-norm penalized $M$-estimator studied in Negahban and Wainwright (2011).

Next we describe an implementation of the max-norm constrained matrix completion procedure, which incorporates the rank estimation approach in Juliá et al. (2011). Assume without loss of generality that $\alpha_{0}$ is known.

Given the observed partial matrix $M_{S}$, the initial matrix $M_{{\rm ini}}$ is obtained by adding the average of the corresponding column to the missing entries of $M_{S}$. Applying the Fast Fourier Transform (FFT) to the columns of $M_{{\rm ini}}$ and taking its modulus, i.e. $F:=|{\rm FFT}(M_{{\rm ini}})|$.

Set an initial rank $r=2$ and an upper bound $r_{\max}$. Clearly, $r_{\max}\leq\min(d_{1},d_{2})$ and it can be computed automatically by adding a criteria for stopping the iteration.

For the current value of $r$, using the computational algorithms given in Section 4 with $R=\alpha_{0}\sqrt{r}$ to solve the max-norm constraint optimization $(\ref{max-est})$. The resulting estimated full matrix is denoted by $\widehat{M}_{r}$.

Apply the FFT to $\widehat{M}_{r}$ as in step 1. Write $F_{r}=|{\rm FFT}(\widehat{M}_{r})|$ and compute the error $e(r)=\|F-F_{r}\|_{F}$.

If $r<r_{\max}$, set $r=r+1$ and go to step 3.

and the corresponding $\widehat{M}_{r^{*}}$ is the final estimate of $M^{*}$. Clearly, the above procedure can be modified by replacing the rank $r$ with the max-norm $R$. A suitable initial value for the max-norm is $R=\alpha_{0}\sqrt{2}$ and at each iteration, increase $R=R+\delta$ with a fixed step size $\delta>0$. An upper-bound $R_{\max}$ could be automatically computed by adding some criteria for stopping the iteration.

Discussions

This paper considers the approximate recovery of approximately low-rank matrices, in particular low-max-norm matrices in contrary to low-trace-norm matrices. The max-norm ball with radius 1 is nearly equivalent to the convex hull of rank-1 matrices, and therefore is an alternative convex surrogate for the rank. A max-norm constrained empirical risk minimization method is proposed and its theoretical properties are studied along with computational algorithms. Allowing for unknown non-uniform sampling which is an important relaxation of the uniform assumption in practice, it is shown that the method is rate-optimal and can be solved efficiently in polynomial time.

When the underlying matrix has exactly rank $r$, it is known that using the trace-norm based approach leads to a mean square error of order $O\{rd(\log d)/n\}$ (Keshavan and Montanari, 2010; Koltchinskii, Lounici and Tsybakov, 2011; Negahban and Wainwright, 2012; Klopp, 2014), where $d=d_{1}+d_{2}$. In the ideal uniform sampling model, the trace-norm regularized method is arguably the mostly preferable one as it achieves optimal rate of convergence (up to a logarithmic factor) and is computationally feasible. The sampling scheme considered in this paper is unspecified and is allowed to be highly non-uniform, which brings additional randomness and uncertainty to the recovery problem. Therefore, we are essentially dealing with a much more complex model, and the max-norm constraint is not only introduced as a convex relaxation for low-rankness according to (2.2) but also takes into account the effect of non-uniform sampling.

Proofs

We prove the main results, Theorems 3.1 and 3.3, in this section. The proofs of a few key technical lemmas including Lemma 3.1 are also given.

For ease of exposition, we write $\widehat{M}=\widehat{M}_{\max}$ as long as there is no ambiguity. To illustrate the main idea, we first consider the case where $\xi_{1},\ldots,\xi_{n}$ are i.i.d. normal random variables and prove that there exists an absolute constant $C$ such that for any $t\in(0,1)$ and a sample size $n$ satisfying $2<n\leq d_{1}d_{2}$,

holds with probability greater than $1-t-e^{-d}$. The case of sub-exponential noise can be obtained via a straightforward adaptation of the arguments for Gaussian noise.

To begin with, noting that $\widehat{M}$ is optimal and $M^{*}$ is feasible for the convex optimization problem $(\ref{max-est})$, we thus have the basic inequality that

This, combined with our model assumption $Y_{i_{t}j_{t}}=M^{*}_{i_{t}j_{t}}+\sigma\xi_{t}$ yields that

where $\widehat{\Delta}=\widehat{M}-M^{*}\in\mathcal{K}(2\alpha,2R)$ is the error matrix. By (6.2), the major challenges in proving Theorem 3.1 consist of two parts, bounding the left-hand side of $(\ref{ineq1})$ from below in a uniform sense and the right-hand side of $(\ref{ineq1})$ from above.

Step 1. (Upper bound). Recalling that $\{\xi_{t}\}_{t=1}^{n}$ is a sequence of $N(0,1)$ random variables and that $S=\{(i_{1},j_{1}),\ldots,(i_{n},j_{n})\}$ is drawn i.i.d. according to $\Pi$ on $[d_{1}]\times[d_{2}]$, we define

Due to Pisier (1989), we obtain that for any realization of the training set $S$ and for any $\delta>0$, with probability at least $1-\delta$ over $\boldsymbol{\xi}=\{\xi_{t}\}_{t=1}^{n}$,

Thus it remains to estimate the following expectation over the class of matrices $\mathcal{K}(\alpha,R)$:

As a direct consequence of $(\ref{eq2.3})$, we have

where $\mathcal{M}_{\pm}$ contains rank-one sign matrices with cardinality $|\mathcal{M}_{\pm}|=2^{d-1}$. For each $M\in\mathcal{M}_{\pm}$, $\sum_{t=1}^{n}\xi_{t}M_{i_{t}j_{t}}$ is a Gaussian random variable with mean zero and variance $n$. Then, the expectation of the Gaussian maximum in (6.5) can be bounded by

Since this upper bound holds uniformly over all realizations of $S$, we conclude that with probability at least $1-\delta$ over both the random samples $S$ and the noise $\boldsymbol{\xi}=\{\xi_{t}\}_{t=1}^{n}$,

In the case of sub-exponential noise, i.e. $\{\xi_{t}\}_{t=1}^{n}$ satisfies the assumption $(\ref{sub-exp})$, it follows from $(\ref{eq2.3})$ that

For any realization of the training set $S=\{(i_{1},j_{1}),\ldots,(i_{n},j_{n})\}$ and for any $M\in\mathcal{M}_{\pm}$ fixed, it follows from a Bernstein-type inequality for sub-exponential random variables (Vershynin, 2012) that

where $c>0$ is an absolute constant. By the union bound, it can be easily verified that for a sample size $n\geq d$,

holds with probability at least $1-e^{-d}$ for some absolute constant $C>0$.

Step 2. (Lower bound). For the given sampling distribution $\Pi$, note that

Here, $\delta$ can be regarded as a tolerance parameter. The goal is to show that there exists some function $f_{\beta}$ such that with high probability, the following inequality

holds uniformly over $M\in\mathcal{C}(\beta,\delta)$.

Proof of $(\ref{RSC})$. Instead, we will prove a stronger result that with exponentially high probability,

holds for all $M\in\mathcal{C}(\beta,\delta)$, based on a straightforward adaptation of the peeling argument used in Negahban and Wainwright (2012). Taking $\varrho=\frac{3}{2}$, define a sequence of subsets

In fact, if there exists some $M\in\mathcal{C}(\beta,\delta)$ satisfying

Therefore, the main task is to show that the latter event occurs with small probability. To this end, define the maximum deviation for each $S\subseteq([d_{1}]\times[d_{2}])^{n}$ that

The following lemma shows that $n^{-1}\|M_{S}\|_{2}^{2}$ does not deviate far from its expectation uniformly for all $M\in\mathcal{B}(D)$.

There exists a universal positive constant $C_{1}$ such that, for any $D>0$,

In view of the above lemma, we take $f_{\beta}(n,d_{1},d_{2})=C_{1}\beta\sqrt{d/n}$ and consider the following sequence of events

with $c_{0}=\log(3/2)/26$, where we used the elementary inequality that

Consequently, for a sample size $n\leq d_{1}d_{2}$ satisfying $\exp(-c_{0}n\delta)\leq\frac{1}{2}$, or equivalently, $n>(c_{0}\delta)^{-1}\log 2$, we obtain that with probability greater than $1-2\exp(-c_{0}n\delta)$,

holds for all $M\in\mathcal{C}(\beta,\delta)$.

Step 3. Now we combine the results in Step 1 and Step 2 to finish the proof. On one hand, it follows from $(\ref{step1})$ that for a sample size $2<n\leq d_{1}d_{2}$,

holds with probability at least $1-e^{-d}$. On the other hand, set $\widetilde{\Delta}=\widehat{\Delta}/(2\alpha)$ such that $\|\widetilde{\Delta}\|_{\infty}\leq 1$ and $\|\widetilde{\Delta}\|_{\max}\leq R/\alpha:=\beta$, or equivalently, $\widetilde{\Delta}\in\mathcal{K}(1,\beta)$. For any $0<t<1$, applying $(\ref{step2})$ with $\delta=\frac{\log(2/t)}{c_{0}n}$ implies that for a sample size $n$ with $2<n\leq d_{1}d_{2}$,

holds with probability at least $1-t$. The last two displays, joint with the basic inequality $(\ref{ineq1})$ lead to the final conclusion $(\ref{mc-ubd})$ after a simple rescaling. Similarly, using the upper bound $(\ref{step1-gen})$, instead of $(\ref{step1})$, together with the lower bound $(\ref{step2})$ proves $(\ref{mc-ubd-gen})$ in the case of sub-exponential noise. ∎

Here, we prove the concentration inequality given in Lemma 6.1. The argument is based on some basic techniques of probability in Banach spaces, including symmetrization, contraction inequality and Bousquet’s version of Talagrand concentration inequality as well as the upper bound $(\ref{eq2.4})$ on the empirical Rademacher complexity of the max-norm ball.

where $\{\varepsilon_{i}\}_{i=1}^{n}$ is an i.i.d. Rademacher sequence, independent of $S$. Given an index set $S=\{(i_{1},j_{1}),\ldots,(i_{n},j_{n})\}$, since $|M_{i_{t}j_{t}}|\leq 1$, using Ledoux-Talagrand contraction inequality (Ledoux and Talagrand, 1991) implies that for $d=d_{1}+d_{2}$,

where we used inequality $(\ref{eq2.4})$ in the last step. Since the “worst-case” Rademacher complexity is uniformly bounded, we have

Next, applying Bousquet’s version of Talagrand’s concentration inequality for empirical processes indexed by bounded functions (Bousquet, 2003) yields that for every $t>0$,

with probability at least $1-e^{-t}$. The conclusion $(\ref{prob-be})$ thus follows by taking $t=nD/26$. ∎

2 Proof of Theorem 3.2

The proof is based on a general result in Srebro, Sridharan and Tewari (2010) on excess risk bounds for learning with a smooth loss. Recall that the noisy response is of the form $Y_{i_{t}j_{t}}=M^{*}_{i_{t}j_{t}}+\xi_{t}$ for $t=1,2,\ldots$, where the location $(i_{t},j_{t})$ of the entry is drawn from $[d_{1}]\times[d_{2}]$ according to $\Pi$ and the noise $\xi_{t}$ on the entry is drawn independently each time. For every $d_{1}\times d_{2}$ matrix $M$, define the quadratic loss function

and its empirical counterpart $\widehat{\mathcal{L}}(M)=\frac{1}{n}\sum_{t=1}^{n}(M_{i_{t}j_{t}}-Y_{i_{t}j_{t}})^{2}+\sigma^{2}$ for a given i.i.d. sample $\{(i_{t},j_{t}),Y_{i_{t}j_{t}}=M^{*}_{i_{t}j_{t}}+\xi_{t}\}_{t=1}^{n}$. In this notation, our estimator $\widehat{M}_{\max}$ can be written as $\widehat{M}_{\max}=\mathop{\rm arg\min}_{M\in\mathcal{K}(\alpha,R)}\widehat{\mathcal{L}}(M)$.

In view of Definition 2.1, define the worst-case Rademacher complexity as

where $\mathcal{K}=\mathcal{K}(\alpha,R)$.

For any $B>0$, let $\mathcal{E}_{B}$ be the event that $\max_{1\leq t\leq n}|\xi_{t}|\leq B$ holds. On $\mathcal{E}_{B}$, applying Theorem 1 in Srebro, Sridharan and Tewari (2010) by taking $H=2$ and $b=5\alpha^{2}+4\alpha\sigma B$ that, for any $0<\delta<1$,

holds with probability at least $1-\delta$ over a random sample $\{(i_{t},j_{t})\}_{t=1}^{n}$ of size $n$, where $C_{1}>0$ is an absolute constant. By (2.4), the worst-case Rademacher complexity $R_{n}(\mathcal{K})$ is bounded by $6R\sqrt{d/n}$. Moreover, note that $\min_{M\in\mathcal{K}(\alpha,R)}{\mathcal{L}}(M)=\mathcal{L}(M^{*})=\sigma^{2}$ and

Putting the above calculations together, we obtain that on the event $\mathcal{E}_{B}$,

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

Finally, it follows from Borell’s inequality that for every $t>0$,

In particular, taking $\delta=n^{-1}$ in both (6.15) and (6.16) proves (3.10). ∎

3 Proof of Theorem 3.3

By construction in Lemma 3.1, setting $\delta=\gamma\alpha\sqrt{d_{1}d_{2}/2}$ we see that $\mathcal{M}$ is a $\delta$-packing set of $\mathcal{K}(\alpha,R)$ in the Frobenius norm. Next, a standard argument (Yang and Barro, 1999; Yu, 1997) yields a lower bound on the $\|\cdot\|_{F}$-risk in terms of the error in a multi-way hypothesis testing problem. More specifically,

where $K(M^{i}\|M^{i})$ denotes the Kullback-Leibler divergence between distributions $(Y_{S}|M^{i})$ and $(Y_{S}|M^{j})$. For the observation model $(\ref{mc-md})$ with i.i.d. Gaussian noise, we have

provided that $r(d_{1}\vee d_{2})\geq 48$ and $\gamma^{4}\leq\frac{\sigma^{2}}{128\alpha^{2}}\frac{r(d_{1}\vee d_{2})}{\mu n}$. If $\frac{\sigma^{2}}{128\alpha^{2}}\frac{r(d_{1}\vee d_{2})}{\mu n}>1$, we choose $\gamma=1$ so that

Otherwise, as long as the parameters $(n,d_{1},d_{2},\alpha,R)$ satisfy $(\ref{quater})$, taking

4 Proof of Lemma 3.1

Next, we show that above random procedure succeeds in generating a set having all desired properties, with non-zero probability. For $1\leq i\leq N$, it is easy to see that

Consequently, $M^{i}\in\mathcal{K}(\alpha,R)$ and it remains to show that the set $\{M^{i}\}_{i=1}^{N}$ satisfies property (ii). In fact, for any $1\leq i\neq j\leq N$,

where $\delta_{kl}$ are independent $0/1$ Bernoulli random variables with mean $1/2$. Using Hoeffding’s inequality gives

Because there are less than $N^{2}/2$ such index pairs in total, the above inequality, together with the union bound implies that with probability at least $1-\frac{N^{2}}{2}\exp(-Bd_{2}/8)\geq 1/2$,

holds for all $i\neq j$. This completes the proof of Lemma 3.1. ∎

Acknowledgements

We thank the editors and an anonymous referee for their careful reviews and constructive comments.

References