Robust Low-Rank Matrix Estimation

Introduction

Netflix, Spotify, Apple Music, Amazon and many other on-line services offer an almost infinite amount of songs or films to their users. Clearly, a single person will never be able to watch every film or to listen to every available song. For this reason, an elaborate recommendation system is necessary in order to allow the users to choose content that already match his or her preferences. Many models and estimation methods have been proposed to address this question. Matrices provide an appropriate way of modelling this problem. Imagine that the plethora of films/songs is identified with the rows of a matrix, call it $B^{*}$, and the users with its columns. One entry of the matrix corresponds to the rating given to film $``i"$ (row) by user $``j"$ (column). This matrix will have many missing entries. These entries are bounded and we can expect the rows of $B^{*}$ to be very similar to each other. It is therefore sensible to assume that $B^{*}$ has a low rank. The challenge is now to predict the missing ratings/fill in the empty entries of $B^{*}$. Define for this purpose the set of observed (possibly noisy) entries

where $\eta$ is for instance the mean highest rating, and $R_{n}$ is some convex empirical error measure that is defined by the data, e.g.

Since the rank of a matrix is not convex we use the nuclear norm as its convex surrogate. This leads us to a relaxed convex optimization problem. For $B\in\mathcal{B}$ we

for some $\tau>0$. The model described above can be considered as a special case of the trace regression model.

In the trace regression model (see e.g. Rohde and Tsybakov 2011) one considers the observations $(X_{i},Y_{i})$ satisfying

where $\varepsilon_{i}$ are i.i.d. random errors. The matrices $X_{i}$ are so-called masks. They are assumed to lie in

where $e_{k}(q)$ is the $q$-dimensional $k$-th unit vector and $e_{l}(p)$ is the $p$-dimensional $l$-th unit vector. We will assume that the $X_{i}$ are i.i.d. with

for all $i\in\left\{1,\dots,n\right\}$, $k\in\left\{1,\dots q\right\}$, and $j\in\left\{1,\dots,p\right\}$. However, we point out that it is not necessary for our estimators to know this distribution. This knowledge will only be used in the proofs of the theoretical results.

The trace regression model together with the space $\chi$ and the distribution on $\chi$ is equivalent to the matrix completion case. The entries of the vector $Y$ can be identified with the observed entries as those in the matrix $A$.

From this, it can be seen that we are in a high-dimensional setting since the number of observations $n$ must be smaller than or equal to the total number of entries of $A$. The setup described above is then called uniform sampling matrix completion. A very similar setup was first considered in Srebro, Rennie and Jaakkola 2005 and Srebro and Shraibman 2005.

As in the standard regression setting, parameter estimation in the trace regression model is also done via empirical risk minimization. Using the Lagrangian form for $B\in\mathcal{B}$ we

where $R_{n}(B)=1/n\sum_{i=1}^{n}\rho(Y_{i}-\trace(X_{i}B))$, $\rho$ is a convex loss function and $\lambda>0$ is the tuning parameter. The loss function is often chosen to be the quadratic loss (or one of its modifications) as in Koltchinskii, Lounici and Tsybakov 2011; Negahban and Wainwright 2011; Negahban and Wainwright 2012; Rohde and Tsybakov 2011 and many others. In Lafond 2015 the case of an error distribution belonging to an exponential family is considered. As long as the errors are assumed to be light tailed as it is the case for i.i.d. Gaussian errors the least squares estimator will perform very well. However, the ratings are heavily subject to frauds (e.g. by the producer of a film). It is necessary to take this fact into account also in the estimation procedure. One might also be interested in estimating the median or another quantile of the ratings. For this purpose, M-estimators based on different losses than the quadratic loss are usually chosen.

2 Proposed estimators

In this paper, we consider the absolute value loss and the Huber loss. The first robust estimator is then given by

defines the Huber loss function. The tuning parameter $\kappa>0$ is assumed to be given for our estimation problem. The possible values for the Huber parameter $\kappa$ depend on the distribution of the errors as shown in Lemma 2.1. In practice, one usually estimates $\kappa$ and $\lambda$ with methods such as cross-validation. Notice that it could happen that the estimators defined in Equations 1.6 and 1.7 are not unique since the objective functions are not strictly convex. As will be shown, the rates depend on the Lipschitz constants of the loss functions and on $\eta$. Typically, the Lipschitz constants of the absolute value loss as well as of the Huber loss induce smaller constants in the rates compared to the Lipschitz constant of the truncated quadratic loss.

Assuming $B=B^{0}$ in Corollary 3.1 the upper bound is typically of the form

where $\lesssim$ means that some multiplicative constants (depending on the tuning parameter $\kappa$) are omitted.

where $0<r<1$ and $\rho_{r}^{r}$ is some reasonably small constant.

3 Related Literature

A first study with robust matrix estimation was made in Chandrasekaran et al. 2011 in a setting with no missing entries. In order to avoid identifiability issues the authors introduce “incoherence” conditions on the low-rank component. These conditions make sure that the low-rank component itself is not too sparse. The locations of the corruptions are assumed to be fixed. In the context of Principal Component Analysis (PCA) which is a special case of the matrix regression model robustness was investigated in Candès et al. 2011. The authors assume that the matrix to be estimated is decomposed in a low-rank matrix and a sparse matrix. In contrast to Chandrasekaran et al. 2011 the non-zero entries of the sparse matrix are assumed to be drawn randomly following a uniform distribution. Following this line of research Li 2011 apply these conditions to the matrix completion problem with randomly observed entries. In a parallel work Chen et al. 2013 consider the case where the indices of the observed entries may be simultaneously both random and deterministic. In these papers only noiseless robust matrix completion is considered.

Cambier and Absil 2016 study computational aspects of robust matrix completion (in the previously mentioned setting). A method relying on Riemannian optimization is proposed. The authors assume the rank of the matrix to be estimated to be known.

In Foygel et al. 2011 weighted nuclear norm penalized estimators with (possibly nonconvex) Lipschitz continuous loss functions are studied from a learning point of view. The partially observed entries are assumed to follow a possibly non-uniform distribution on the set $\chi$. In contrast, our derivations rely among other properties on the convexity of the risk (i.e. the margin conditions).

In contrast to the previously mentioned papers on robust matrix completion, we consider (possibly heavy-tailed) random errors that affect the observations but not the truth.

4 Organization of the paper

Preliminaries

In this section the assumptions on the loss functions, the risk, and the distribution of the errors are presented. In particular, Assumptions 1-3 below are on the curvature of the (theoretical) risk. They are used to derive the deterministic sharp and non-sharp oracle inequalities. It is important to notice that the curvature of the risk mainly depends on the properties of the distribution of the errors. Assumptions 4 and 5 below will be shown to be sufficient for Assumptions 2 and 3 to hold, respectively.

The first assumption is about the loss function.

The next two assumptions ensure the identifiability of the parameters by requiring a sufficient convexity of the risk around the target.

One-point-margin condition. There is an increasing strictly convex function $G$ with $G(0)=0$ such that for all $B\in\mathcal{B}$

where $R$ is the theoretical risk function.

Two-point-margin condition. There is an increasing strictly convex function $G$ with $G(0)=0$ such that for all $B,B^{\prime}\in\mathcal{B}$ we have

where $R$ is the theoretical risk function and $[\dot{R}(B^{\prime})]_{kl}=\frac{\partial}{\partial B_{kl}}R(B)|_{B=B^{\prime}}$.

Assumption 1 is crucial when it comes to the application of the Contraction Theorem which in turn allows us to apply the dual norm inequality to find a bound for the random part of the oracle bounds. Assumptions 2 and 3 are essential in the proofs of the (deterministic) results. In particular, in addition to the differentiability of the empirical risk $R_{n}$, Assumption 3 is responsible for the sharpness of the first oracle bound that will be proved. The margin conditions are strongly related to the shape of the distribution function and the corresponding density of the errors.

For the specific application to the Huber loss and absolute value loss estimators we show that mild conditions on the distribution of the errors ensure a sufficient curvature of the risk for both loss functions under study.

Assumption 3 holds under a weak condition on the distribution function of the errors:

Assume that there exists a constant $C_{1}>0$ such that the distribution function $F$ with density with respect to Lebesgue measure $f$ of the errors fulfills

Assumption 4 implies Assumption 3 with $G(u)=u^{2}/(2C_{1}^{2}pq)$.

The following assumption guarantees that Assumption 2 holds.

Suppose $\varepsilon_{1},\dots,\varepsilon_{n}$ are i.i.d. with median zero and density $f$ with respect to Lebesgue measure. Assume that for $C_{2}>0$

Assumption 5 implies Assumption 2 with $G(u)=u^{2}/(2C_{2}^{2}pq)$.

Another important fact is that when the distribution of the errors is assumed to be symmetric around zero $B^{*}=B^{0}$. This phenomenon is discussed in Section 4 of the Supplement.

2 Properties of the nuclear norm

Consider the singular value decomposition of the oracle $B$ with rank $s^{\star}$ given by

where $P$ is a $p\times q$ matrix, $Q$ a $q\times p$ matrix and $\Lambda$ a $q\times q$ diagonal matrix containing the ordered singular values $\Lambda_{1}\geq\dots\geq\Lambda_{q}$. Then the nuclear norm is given by

The matrix $B^{+}$ is called “active” part of the oracle $B$, whereas the matrix $B^{-}$ is called the “non-active” part. The singular value decomposition of $B^{+}$ is given by

We observe that the integer $s$ is not necessarily the rank of the oracle $B$. The choice of $s$ is free. One may choose a value that trades off the roles of the “active” part $B^{+}$ and “non-active” part $B^{-}$, see Lemma 4.1. The following lemma is adapted from Lemma 7.2 and Lemma 12.5 in van de Geer 2016.

We then say that the triangle property holds at $B^{+}$.

From now on, we write $\Omega^{+}=\Omega_{B^{+}}^{+}$ and $\Omega^{-}=\Omega_{B^{+}}^{-}$. Equation 2.8 is proved in Appendix A.

Hence, the property that our estimators should mimic is not the rank of the oracle but rather the fact that the “non-active” part is zero under the semi-norm induced by the active part.

Moreover, we define the norm $\underline{\Omega}$ as

Notice that the semi-norms $\Omega^{+}$ and $\Omega^{-}$ form a complete pair, meaning that $\underline{\Omega}:=\Omega^{+}+\Omega^{-}$ is a norm.

The estimation error in several different norms can thus be “computed” in general (semi-)norms.

Let $X_{1},\dots,X_{n}$ be i.i.d. $q\times p$ matrices that satisfy for some $\alpha\geq 1$ (and all $i$)

Then for a constant $C$ and for all $t>0$

This theorem is used with the tail summation property of the expectation in the derivations of the tail bounds in Section 3 of the Supplementary Material.

Oracle inequalities

We first give two deterministic sharp and non-sharp oracle inequalities. The connection to the empirical process parts and to the specific loss functions follows in Subsection 3.3. Let $B^{0}=\underset{B^{\prime}\in\mathcal{B}}{\argmin}\ R(B^{\prime})$ be the target. It is assumed that $q\leq p$.

Here, we assume that the loss function is differentiable and Lipschitz continuous. The next lemma gives a connection between the empirical risk and the penalization term.

Suppose that $R_{n}$ is differentiable. Then for all $B\in\mathcal{B}$

The following theorem is inspired by Theorem 7.1 in van de Geer 2016. In contrast to this theorem, we need to bound the empirical process part differently. In view of the application to the matrix completion problem, we assume a specific bound on the empirical process.

Suppose that Assumptions 1 and 3 hold, that the loss function is differentiable, and let $H$ be the convex conjugate of $G$. Assume further for all $B^{\prime}\in\mathcal{B}$ that for $\lambda_{\varepsilon}>0$ and $\lambda_{*}>0$

Take $\lambda>\lambda_{\varepsilon}$. Let $0\leq\delta<1$ be arbitrary, and define

In the proof of this theorem the differentiability of the loss function and Assumption 3 are crucial. Without this property an additional term arising from the one-point-margin condition would appear in the upper bound. This term would then lead to a non-sharp bound.

2 Non-sharp oracle inequality

Instead of bounding an empirical process term depending on the derivative of the empirical and theoretical risks we need to consider differences of these functions.

Suppose that Assumptions 1 and 2 hold. Let $H$ be the convex conjugate of $G$. Suppose further that for $\lambda_{\varepsilon}>0$, $\lambda_{*}>0$, and all $B^{\prime}\in\mathcal{B}$

Let $0<\delta<1$, take $\lambda>\lambda_{\varepsilon}$ and define

It has to be noticed that the above bound is “good” only if $R(B)-R(B^{0})$ is already small. The main cause for the non-sharpness is Assumption 2 that leads to an additional term in the upper bound of the inequality.

3 Applications to specific loss functions

We now apply the deterministic sharp and non-sharp oracle inequalities to the case of the Huber loss and absolute value loss, respectively. We assume in both cases that the distribution of the errors is symmetric around $0$ so that $B^{0}=B^{*}$. This is discussed in detail in Section B.4.

We first consider the case that arises by choosing the Huber loss. Theorem 3.1 together with Lemma 2.1 and the first claim of Lemma 3.2 in the Supplement imply the following corollary. It is useful to notice that the Lipschitz constant of the Huber loss is $2\kappa$.

Let $B=B^{+}+B^{-}$ where $B^{+}$ and $B^{-}$ are defined in Equation 2.4. Let Assumption 4 be satisfied.

and $\lambda_{*}=8\eta(4\eta+2\kappa)p\log(p+q)/(3n)+\lambda_{\varepsilon}$.

Assume that $\lambda>\lambda_{\varepsilon}$. Take $0\leq\delta<1$,

Choose $j_{0}:=\lceil\log_{2}(7q\sqrt{pq}\eta)\rceil$ and define

Then we have with probability at least $1-\alpha$ that

Assumption 4 guarantees that the risk function is sufficiently convex. From this assumption we also obtain a bound for the possible values of the tuning parameter $\kappa$. We can also see that the results hold for errors with a heavier tail than the Gaussian. The choice of the noise level $\lambda_{\varepsilon}$ and consequently of the tuning parameter $\lambda$ results from the the probability inequalities for the empirical process in Section B.3. The quantity $\lambda_{*}$ is also a consequence of the bound on the empirical process part. However, it does not affect the asymptotic rates.

Absolute value loss - non-sharp oracle inequality

The next corollary is an application to the case of the absolute value loss. Theorem 3.2 combined with Lemma 2.2 and the second claim of Lemma 3.2 in the Supplement lead to the following corollary. The Lipschitz constant in this case is $1$.

Let the oracle $B$ be as in 2.4. Suppose that Assumption 5 is satisfied. For a constant $C_{0}>0$ let

and $\lambda_{*}=8\eta p\log(p+q)/(3n)+\lambda_{\varepsilon}$. Take $0<\delta<1$ and $\lambda>\lambda_{\varepsilon}$. Choose $j_{0}:=\lceil\log_{2}(7q\sqrt{pq}\eta)\rceil$ and define

Then we have with probability at least $1-\alpha$ that

Also in this case, the choices of $\lambda_{\varepsilon}$ and $\lambda_{*}$ are a consequence of the probability bounds.

Asymptotics and Weak sparsity

The results in Section 3 are valid for finite values of the dimension of the matrix $p,q$, the rank, and the number of observed entries $n$. A question that is answered in this section is how the estimation errors of the proposed estimators behave when $n,p,$ and $q$ are allowed to grow.

As mentioned in Negahban and Wainwright 2012, practical reasons motivate the assumption that the matrix $B^{0}$ is not exactly low-rank but only approximately. In relation to the matrix completion problem one observes that the ratings given by the users are unlikely to be exactly equal but rather very similar. This translates to a matrix that is not low-rank. However, it is sensible to assume that the matrix is almost low-rank. The notion of weak sparsity quantifies this assumption by assuming that for some $0<r<1$ and $\rho>0$

where $\Lambda_{1}^{0},\dots,\Lambda_{q}^{0}$ are the singular values of $B^{0}$. For $r=0$ we have under the convention that $0^{0}=0$ that

where $s_{0}$ is the rank of $B^{0}$. The following lemma gives a bound of the non-active part of the matrix $B$ that appears in the oracle bounds.

We first consider the asymptotic behavior of our estimators in the case of an exactly low-rank matrix and deduce from this the asymptotics for the case of an approximately low-rank matrix.

By Corollary 3.1, assuming that $q\log(1+q)=o\left(\frac{n}{\log(p+q)}\right)$, and therefore using the choice for the noise level

We choose for simplicity the oracle to be the matrix $B^{0}$ itself with $s_{0}=\mbox{rank}(B^{0})$. Then, we make use of the two point margin condition that is shown to hold in Lemma 2.1.The resulting rate is then given by

where $\kappa$ is the Huber parameter and $C_{1}$ is the constant from Lemma 2.1.

The rate (4.4) depends on $\eta$ as in Koltchinskii, Lounici and Tsybakov 2011 and on the Lipschitz constant of the loss function which is typically smaller than $\eta$. If $C_{1}^{2}=O(\eta)$, the constant in front of the rate is of order $O(\eta^{4})$. This is a “worst-case” scenario that shows the cost that is paid when allowing for very general error distributions as in our case. We emphasize that in this case the distribution of the errors is not required to have a density.

In addition to the rate obtained for the Frobenius norm, we are also able to derive rates for the estimation error measured in nuclear norm. From Corollary 3.1 and Equation 2.9 under the previous conditions it follows that

1.2 non-sharp

By Corollary 3.2 it is known that the assumption $q\log(1+q)=o\left(\frac{n}{\log(p+q)}\right)$ leads to the choice $\lambda_{\varepsilon}\asymp\sqrt{\log(p+q)/nq}$. Therefore, we have

What can be observed comparing the rates in Equations (4.3) and (4.6) is the presence of the additional term $R(B)-R(B^{0})$ in the non-sharp case in contrast to the sharp case. We choose again the oracle to be the matrix $B^{0}$ itself. By the one point margin condition derived in Lemma 2.2 we see that the rate of convergence in this case is given by

where the constant $C_{2}$ comes from Lemma 2.2.

If $C_{2}^{2}=O(\eta)$ a comparison with the rates obtained in Koltchinskii, Lounici and Tsybakov 2011 shows that the rates agree. In contrast to the rate obatined for the Huber loss (Equation 4.4), the distribution of the errors is assumed to have a density. This leads to a constant of order $O(\eta^{2})$ in a “worst-case” scenario. It is a natural consequence of the stronger assumption on the distribution of the errors. This is comparable to the constant obtained in Koltchinskii, Lounici and Tsybakov 2011.

In analogy to the previous case, we are able to derive a rate for the estimation error measured in nuclear norm:

The rates are indeed very slow but this is not surprising given that per entry the number of observations is about $n/(pq)$. The price to pay for the estimation of the reduced number of parameters $ps_{0}$ is given by the term $\log(p+q)$.

2 Weak sparsity

In what follows, the asymptotic behavior of the proposed estimators is discussed when applied to an estimation problem where one aims at estimating a matrix that is not exactly low-rank. With Lemma 4.1 and the rates given in the previous section we are able to derive an explicit rate also for the approximatley low-rank case. For this purpose, we assume that Equation 4.1 holds.

The following corollary gives rates for the estimation error of the Huber estimator when used for estimation of a not exactly low-rank matrix.

With $q\log(1+q)=o\left(\frac{n}{\log(p+q)}\right)$ we choose

2.2 Absolute value estimator

Using the oracle inequality under the weak sparsity assumption we obtain the following result.

With $q\log(1+q)=o\left(\frac{n}{\log(p+q)}\right)$ we choose

Then we have for the Frobenius norm of the estimation error

Simulations

In this section, the robustness of the Huber estimator 1.7 is empirically demonstrated. In Subsection 5.3, the Huber estimator is compared with the estimator proposed in Klopp, Lounici and Tsybakov 2016 under model 1.4 and 1.9 with each Student t and standard Gaussian noise. The sample size ranges in all simulations for all dimensions considered here from $3p\log(p)s_{0}$ to $pq$. Between minimal and maximal sample size there are in each case $10$ points. To illustrate the rate derived in Section 4 we compute the error $\|\hat{B}_{H}-B^{0}\|_{F}^{2}$ for different dimensions of the problem under increasing number of observations.

To compute the solution of the optimization problem 1.7 functions from the Matlab library cvx (CVX Research 2012) were used.

Throughout this section the error is assumed to have the following shape

To verify the robustness of the estimator 1.7 and the rate of convergence that was derived in Section 4 we use the Student t distribution with $3$ degrees of freedom. Every point in the plots corresponds to an average of $25$ simulations. The value of the tuning parameter is set to

A comparison with $\lambda_{\varepsilon}$ from Corollary 3.1 indicates that $\lambda$ is rather small. For the settings we consider in this section we found that this value for $\lambda$ is more appropriate. As done in Candès and Plan 2010, for a better comparison between the error curves of our estimator and the oracle rate in Equation (4.4) this rate was multiplied with $1.68$ in the case of Student t distributed errors and with $1.1$ in the case of Gaussian errors.

The variance of the Student t distribution with $\nu>2$ degrees of freedom is given by

Figure 1(a) shows a comparison between the Huber estimator 1.7 with the estimator that uses the quadratic loss in the case of Student t with $3$ degrees of freedom distributed errors. As expected, the estimator that uses the quadratic loss is not robust against the corrupted entries. On the other hand, we can see in Figure 1(b) that the Huber estimator performs almost as well as the quadratic loss estimator in the case of Gaussian errors with variance $1$. In agreement with the theory, the rate of the estimator is very close to the oracle rate for sufficiently large sample sizes. The value of $\kappa$ that we used in the simulations is $1.345$. The maximal rating $\eta$ is chosen to be $\eta=10$.

2 Changing the problem size

In order to confirm/verify the theoretical results, we proceed similarly to what was done in Negahban and Wainwright 2011 and Negahban and Wainwright 2012 in the corresponding cases. Here, we consider three different problem sizes: $p,q\in\left\{30,50,80\right\}$. In Figure 2(a) we observe that as the problem gets harder, i.e. as the dimension of the matrix increases, also the sample size needs to be larger. Figure 2(b) shows that by rescaling the sample size by $n/(3ps_{0}\log(p))$ the rate of convergence agrees very well with the theoretical one. It is assumed that the rank of the matrices is $s_{0}=2$ for all cases. Every point corresponds to an average of $25$ simulations. The maximal rating $\eta$ and the tuning parameter $\kappa$ are chosen as before.

3 Comparison with a low-rank + sparse estimator

In this subsection, we compare the performance of the Huber estimator 1.7 with the performance of the low-rank matrix estimator proposed by Klopp, Lounici and Tsybakov 2016 1.10. We first compare the estimators $\hat{B}_{H}$ and $\hat{L}$ with the observations $Y_{i}$ generated according to the model 1.4 with standard Gaussian and Student t with $3$ degrees of freedom distributed errors. Equation (21) in Klopp, Lounici and Tsybakov 2016 suggests that the tuning parameters are chosen as follows:

where $\lambda_{1}$ and $\lambda_{2}$ are the tuning parameters of the estimator 1.10. Also in this case it has to be noticed that the tuning parameters are smaller than the theoretical values given in their paper.

In Figure 3(a) the Huber estimator 1.7 is compared with the low-rank plus sparse estimator 1.10 under the model 1.4 with i.i.d. Student t noise with $3$ degrees of freedom. As expected, these estimators perform comparably well under the trace regression model 1.4. In Figure 3(b) the same estimators are compared under the model 1.4 with i.i.d. standard Gaussian noise. Also in this case, we see that both estimators achieve approximately the same error. These observations are not surprising since the theoretical analysis of Section 3 could be carried over by adapting the (semi-)norms to the different penalization.

We now consider the model proposed in Klopp, Lounici and Tsybakov 2016 where around $5\%$ of the observed entries are taken to be only one rating. This is the case of malicious users who systematically rate only one particular movie with the same rating. We refer to Section 2.3 of Klopp, Lounici and Tsybakov 2016 for more details on this setting. In Figure 4(a) we see that the Huber estimator outperforms the low-rank plus sparse estimator with Student t noise with $3$ degrees of freedom. This might be due to the quadratic loss function and to the choice of the tuning parameters. In Figure 4(b) where Gaussian noise is considered we observe that both estimators perform almost equally well.

Discussion

In this paper, we have derived sharp and non-sharp oracle inequalities for two robust nuclear norm penalized estimators of the noisy matrix completion problem. The robust estimators were defined using the well-known Huber loss for which the sharp oracle inequality has been derived and the absolute value loss for which we have shown a non-sharp oracle inequality. For both types of oracle inequalities we proved a general deterministic result first and added then the part arising from the empirical process. We have also shown how to apply the oracle inequalities to the case where we only assume weak sparsity, i.e. approximately low-rank matrices. It is worth pointing out that our estimators do not require the distribution on the set of matrices (1.5) to be known in contrast to e.g. Koltchinskii, Lounici and Tsybakov 2011. In our case, the distribution on the set of matrices (1.5) is only needed in the theoretical analysis. The proofs of the oracle inequalities rely on the properties of the nuclear norm, and for the empirical process part on the Concentration, Symmetrization, and Contraction Theorems. A main tool in this context was also the bound on the largest singular value of a matrix with finite Orlicz norm. Our simulations, in the case of the Huber loss, showed a very good agreement with the convergence rates obtained by our theoretical analysis. We saw that the oracle rate is attained up to constants in presence of non-Gaussian noise and that the robust estimation procedure outperforms the quadratic loss function.

It is left to future research to establish a sharp oracle inequality also for the case of a non-differentiable robust loss function. The Contraction inequality used in this paper for the Huber loss requires that also the derivative of the loss is Lipschitz continuous. This is not the case for the absolute value loss. Thanks to the convexity of the loss function it might be possible to derive a sharp result also for this case.

Appendix A Proofs of main results

Using the triangle property at $B^{+}$ with $B^{\prime}=B^{+}$ we obtain

By the triangle property at $B^{+}$ with $B^{\prime}=B=B^{+}+B^{-}$ we have that

By the triangle inequality it follows using $\Omega^{-}(B^{+})=0$ that

For an arbitrary $B$ we have again by the triangle inequality

Applying the triangle property at $B^{+}$ we find that

Apply now twice the triangle inequality (first inequality) to find that

where it was used that $\Omega(B^{-})=0$ and that $\Omega^{-}(B)\leq\Omega^{-}(B^{-})\leq\|B^{-}\|_{\mbox{nuclear}}$. ∎

Let $B\in\mathcal{B}$. Define for $0<t<1$

Letting $t\rightarrow 0$ the claim follows. ∎

The first order Taylor expansion of $R$ at $\hat{B}$ is given by

then by the two-point-margin condition 3 we find that

By the two-point inequality (Lemma 3.1) we have that

We then have by the convex conjugate inequality

We start the proof with the following inequality using the fact that $\hat{B}$ is the minimizer of the objective function.

Then, by adding and subtracting $R(\hat{B})$ on the left hand side and $R(B)$ on the right hand side we obtain

Applying Assumption 3.1, the definition of $\underline{\Omega}$, and Lemma 2.8 we obtain

Since later on we apply Assumption 2 we subtract on both sides of the above inequality $R(B^{0})$.

It is then useful to make the following case distinction that allows us to obtain an upper bound for the estimation error. Case 1 If $\overline{\lambda}\Omega^{+}(\hat{B}-B)\leq\frac{(1-\delta)}{\delta}\left(\lambda_{*}+R(B)-R(B^{0})+2\lambda\|B^{-}\|_{\mbox{nuclear}}\right)$, then

By multiplying Equaiton A on both sides with $\delta$ we arrive at

Case 2 If $\overline{\lambda}\Omega^{+}(\hat{B}-B)\geq\frac{(1-\delta)}{\delta}\left(\lambda_{*}+R(B)-R(B^{0})+2\lambda\|B^{-}\|_{\mbox{nuclear}}\right)$, then

We then obtain using the definition of $\Omega^{+}$ in Lemma 2.3

Invoking the convex conjugate inequality and Assumption 2 we get

Combining the two cases we have for the estimation error

and for the second claim we conclude that

The theoretical risk function arising from the Huber loss is given by

Suppose that $X_{i}$ has its only $1$ at entry $(k,j)$. Then $XB=(B)_{jk}$. Define

The second derivative of $r(x,B)$ with respect to $B_{jk}$ is then given by

Therefore, the Taylor expansion around $B^{\prime}$ is given by

We can see that Assumption 3 holds with $G(u)=u^{2}/(2C_{1}^{2}pq)$. ∎

For the (theoretical) risk function $R$ arising from the absolute value loss we have

Using the tower property of the conditional expectation we obtain

Suppose that $X_{i}$ has its only $1$ at entry $(k,j)$. Then $XB=(B)_{jk}$. Define

The Taylor expansion of $r(x,B)$ around $B^{0}$, assuming that $B^{0}$ minimizes $r$, is given by

which means that the one point margin Condition 2 is satisfied with $G(u)=u^{2}/(2C_{2}^{2}pq)$. ∎

Appendix B Supplemental Material

This supplemental material contains an application to real data sets, the proofs of the lemmas in Section 2 of the main text, and a section on the bound of the empirical process part of the estimation problem.

In Section 5 we have shown several synthetic data examples. The convex optimization problems there were solved using the semidefinite programming (SDP) toolbox CVX Research 2012. These algorithms work very well for comparably low-dimensional optimization problems. When real datasets are considered, due to the much larger problem sizes different algorithms are needed. To solve the optimization problem with real data a proximal gradient algorithm is used. The algorithm is given in pseudocode.

We define $F_{1}$ to be the empirical risk for the Huber loss

We define $F_{2}$ to be the empirical risk for the quadratic loss

where $\dot{\rho}_{H}(Y_{i}-\trace(X_{i}B))$ is given in the proof of Lemma 2.1 in the main paper.

For the nuclear norm the proximity operator has a closed form: let $W=U\text{diag}(\sigma_{1},\dots,\sigma_{\min(p,q)})V^{\prime}$ be the singular value decomposition of $W$, then

It is known that the solution of the optimization problem $\hat{B}_{H}$ satisfies the following fixed point equation

The same holds also for the quadratic loss function. To compute the proximal operator in Algorithm 1 the function proxanuclearnorm from the Matlab toolbox Unlocbox 2016 was used. The algorithm is a Nesterov-type Accelerated Proximal Gradient algorithm. We refer to Section 4.3 in Parikh and Boyd 2014 and references therein. In particular, the present form of Algorithm 1 goes back to Beck and Teboulle 2009. The Huber constant is chosen to be $\kappa=2$ in all the examples. Algorithm 1 is applied to the Huber loss ($i=1$) and to the quadratic loss ($i=2$). The tuning parameter for the quadratic loss case is smaller than for the Huber loss.

This method was applied to the MovieLens data set which consists of $p=943$ users, $q=1682$ movies, and $100^{\prime}000$ observed ratings, we call the set containing the indices of the observed ratings $\Gamma$. The minimal and maximal ratings are $1$ and $5$ respectively. Every user has rated at least $20$ movies. The data set is available at MovieLens 1998. Training and testing sets of different sizes are drawn randomly without replacement from the set of observed ratings. The testing error is computed as follows: with $100^{\prime}000=n_{\text{test}}+n_{\text{train}}$ and using the disjoint union $\Gamma=\Gamma_{\text{test}}\cup\Gamma_{\text{train}}$ (where $|\Gamma_{\text{test}}|=n_{\text{test}}$ and $|\Gamma_{\text{train}}|=n_{\text{train}}$) we have

It can be observed that the test error of the quadratic loss estimator is slightly smaller than the test error for the Huber loss estimator. This might be due to the fact that there are no heavily corrupted entries in the data. On the other hand, this indicates that the Huber estimator is able to “adapt” also to the usual setting without heavy corruptions.

We have applied Algorithm 1 also to the MovieLens data set with $1^{\prime}000^{\prime}209$ observed ratings from $6^{\prime}040$ users on $3^{\prime}952$ movies. The minimal and maximal ratings are $1$ and $5$ respectively. Every user has rated at least $20$ movies. The data set is available at MovieLens 2003.

B.2 Proofs of Lemmas in Section 2

The following lemma shows “equivalence” of the nuclear norm and Frobenius norm. This fact is useful since it is more common and meaningful to measure the estimation error in the Frobenius norm. The rest of the lemma contains technicalities that are used in the sequel to derive the triangle property.

where $\left\|A\right\|_{F}$ is the Frobenius norm of $A$ and $\Lambda_{\mbox{\emph{max}}}(A)$ its largest singular value.

Consider the singular value decomposition (SVD) of $A$ with $s:=\mbox{rank}(A)$

The first claim follows from Hölder’s inequality applied in two different manners for the lower and upper bounds.

For the second claim consider the $p-$dimensional $j-$th unit vector $e_{j}$. Define

Then, by the invariance of the trace under cyclic permutations we obtain the claimed result. ∎

In order to show that the triangle property holds we also need the dual norm of the nuclear norm so that to apply the dual norm inequality (see Lemma B.2). The subdifferential of the nuclear norm is then used to deduce the triangle property.

The following lemma gives the dual norm and the subdifferential of the nuclear norm. It cites results of Lange 2013 and Watson 1992.

where $\Lambda_{\max}(A)$ is the largest singular value of $A$ . Moreover, the subgradient of the nuclear norm is given by

The derivation of the dual norm of the nuclear norm can be found in Example 14.3.6. from Lange 2013. A justification for the second claim can be found in Watson 1992. ∎

Let $Z\in\partial\left\|B^{+}\right\|_{\mbox{nuclear}}$, i.e.

Recall the definition of the subdifferential of the nuclear norm

For $Z\in\partial\left\|B^{+}\right\|_{\mbox{nuclear}}$ we have

To prove the first assertion we need to bound the right hand side of the above inequality. For simplicity consider first $\trace\left(Z_{1}^{T}B^{\prime}\right)$. Using the invariance of the trace under cyclic permutations we have

since $\Lambda_{\max}\left(P^{+}Q^{+^{T}}\right)=1$.

On the other hand, consider $\trace\left(Z_{2}^{T}B^{\prime}\right)$. Again by the invariance of the trace under cyclic permutations we have

Hence, it is possible to find a $W$ such that $\Lambda_{\max}\left(W\right)=1$ and such that

Substituting $B^{\prime}$ with $B^{\prime}-B^{+}$ we have that

Since the dual of $\underline{\Omega}$ may not be easy to deal with we bound it by the dual norm of the nuclear norm. By doing so, it will be possible to apply results about the tail of the maximum singular value of a sum of independent random matrices.

For the dual norm of $\underline{\Omega}$ we have that

By Lemma 2.1 in the main text it is known that

Therefore, using Lemma 2.2 in the main text we have for the dual norms that

B.3 Probability bounds for the empirical process

Now we check the assumptions of Theorem 2.1 in the main text.

For the case of matrix completion with uniform sampling we obtain that

where $\iota$ is a $q-$vector consisting of only $1$’s. It follows that

with probability at least $1-\exp(-p\log(p+q))$.

Let $\rho$ be a Lipschitz continuous function with Lipschitz constant $L$. Assume that for a constant $K_{\underline{\Omega}}$

with probability at least $1-\exp(-p\log(p+q))$.

The proof of this lemma is based on the Symmetrization Theorem, on the Contraction Theorem, on the dual norm inequality, and on Bousquet’s Concentration Theorem. We use the notation $\rho(X_{i}B,Y_{i}):=\rho(Y_{i}-\trace(X_{i}B))$ and $\dot{\rho}(X_{i}B,Y_{i}):=\dot{\rho}(Y_{i}-\trace(X_{i}B))$.

The first inequality follows from Theorem B.1, the second from Theorem B.2 and the dual norm inequality. With $S=1/\sqrt{q}$ and $K=\sqrt{1/\log 2}$ in Theorem 2.1 in the main text together with the concavity of the logarithm

we obtain integrating the tail that the expectation of the largest singular value of the sum of masks is bounded by

we have using the Lipschitz continuity of the loss function and the fact that the $X_{i}\in\chi$ are i.i.d.

we obtain from Bousquet’s Concentration Theorem B.3 that for all $t>0$

For the second claim we proceed similarly. We have

The first inequality follows from Theorem B.1, the second from Theorem B.2. Moreover, we have for all $i=1,\dots,n$ that

In view of applying Bousquet’s Concentration Theorem B.3 we have with a similar calculation as for the first claim

To obtain a uniform bound for all $B\in\mathcal{B}$ we use the peeling device given in van de Geer 2000.

We subdivide the set $\mathcal{B}$ as follows for a fixed $B\in\mathcal{B}$

We then consider the set $\left\{B^{\prime}\in\mathcal{B}:1<\underline{\Omega}(B-B^{\prime})\leq 14q\sqrt{pq}\eta\right\}$. We first refine it by choosing $j_{0}$ as the smallest integer such that $j_{0}+1>\log_{2}(14q\sqrt{pq}\eta)$. This leads us to

By the first claim in Lemma B.4 we can conclude that

To keep the notation clean we define the event

Then $\mathcal{C}\subset\bigcup_{j=0}^{j_{0}}\mathcal{C}_{j}$. Therefore, by the union bound

The second claim follows by an analogous reasoning. ∎

The proof of Lemma B.4 requires the Symmetrization Theorem and the Contraction Theorem (or Contraction Principle). The first theorem allows us to reduce the case of non-centred random variables to the case of mean zero random variables. A main tool for this type of reduction is a sequence of Rademacher random variables. The second theorem, the Contraction Principle, is used to compare the limit behaviour of two series of Rademacher random variables with different coefficients. We state these results for the sake of completeness.

Suppose that for all $i=1,\dots,n$ and for all $f\in\mathcal{F}$

Assume further for a constant $D>0$ and for all $f\in\mathcal{F}$ that

B.4 On the distribution of the errors

In this section we discuss the consequences arising from the assumption that the distribution of the errors is symmetric around $0$. Let $F$ be the distribution function of the errors. For purposes of illustration, we discuss the location model. The case of low-rank matrix estimation then follows easily. The location model is as follows

where $\mu^{*}$ is some fixed real number such that $|\mu^{*}|\leq\eta$ and $\varepsilon$ is additive noise with symmetric (around $0$) distribution.

Assume that $f$ is the density with respect to Lebesgue measure of the errors. Suppose further that $f(u)>0$ for all $|u|\leq 2\eta$. Define

We first notice that $\mu^{0}_{1}$ is the median of the distribution of $X$. Since the distribution is continuous and the density positive everywhere the median is unique. Since the distribution of the errors is symmetric around $0$ the distribution of $X$ is symmetric around $\mu^{*}$. This implies that $\mu^{*}$ must be the median. But since the median is unique we must have $\mu^{*}=\mu^{0}_{1}$. ∎

Assume that the distribution function of the errors satisfies

where $\rho_{H}$ is the Huber loss as defined in the main text. Then,

We notice that the first derivative of the Huber loss evaluated in $X-\mu^{0}_{2}$ is given by

We now use the symmetry of the distribution of the errors. This translates to

By changing the variables in the previous integrals we arrive at

The previous integral is always larger than $0$ by assumption unless $\mu^{*}=\mu^{0}_{2}$. ∎

B.5 Proof of Lemma 4.1 in the main text

The proof of this lemma is analogous to the first part of the proof of Corollary 2 in Negahban and Wainwright 2012. Take $\sigma>0$ such that

where $\Lambda_{k}^{0}$ denotes a singular value of the matrix $B^{0}$. Therefore we obtain

References