Estimation of high-dimensional low-rank matrices

Consider the observations $(X_{i},Y_{i})$ satisfying the model

where $X_{i}\in\mathbb{R}^{m\times T}$ are given matrices ($m$ rows, $T$ columns), $A^{*}\in\mathbb{R}^{m\times T}$ is an unknown matrix, $\xi_{i}$ are i.i.d. random errors, $\operatorname{tr}(B)$ denotes the trace of square matrix $B$ and $X^{\prime}$ stands for the transposed of $X$. Our aim is to estimate the matrix $A^{*}$ and to predict the future $Y$-values based on the sample $(X_{i},Y_{i})$, $i=1,\ldots,N$.

We will call model (1) the trace regression model. Clearly, for $T=1$ it reduces to the standard regression model. The “design” matrices $X_{i}$ will be called masks. This name is motivated by the fact that we focus on the applications of trace regression where $X_{i}$ are very sparse, that is, contain only a small percentage of nonzero entries. Therefore, multiplication of $A^{*}$ by $X_{i}$ masks most of the entries of $A^{*}$. The following two examples are of particular interest.

(i) Point masks. For some, typically small, integer $d$ the point masks $X_{i}$ are defined as elements of the set

where $e_{k}(m)$ are the canonical basis vectors of $\mathbb{R}^{m}$. In particular, for $d=1$ the point masks $X_{i}$ are matrices that have only one nonzero entry, which equals to 1. The problem of estimation of $A^{*}$ in this case becomes the problem of matrix completion; the observations $Y_{i}$ are just some selected entries of $A^{*}$ corrupted by noise, and the aim is to reconstruct all the entries of $A$. The problem of matrix completion dates back at least to Srebro, Rennie and Jaakkola (2005), Srebro and Shraibman (2005) and is mainly motivated by applications in recommendation systems. We will analyze the following two special cases of matrix completion:

• –

  USR (Uniform Sampling at Random) matrix completion. The masks $X_{i}$ are independent, uniformly distributed on

  $$\mathcal{X}_{1}=\{e_{k}(m)e_{l}^{\prime}(T)\dvtx 1\leq k\leq m,1\leq l\leq T\},$$

  and independent from $\xi_{1},\dots,\xi_{N}$.

• –

  Collaborative sampling (CS) matrix completion. The masks $X_{i}$ (random or deterministic) belong to $\mathcal{X}_{1}$, are all distinct and independent from $\xi_{1},\ldots,\xi_{N}$.

The CS matrix completion model is natural to describe recommendation systems where every user rates every product only once. The USR matrix completion can be used for transmission of a large-dimensional matrix trough a noisy communication channel; only a chosen small number of entries is transmitted, and nevertheless the original matrix $A^{*}$ can be reconstructed by the receiver. An important feature of the real-world matrix completion problems is that the number of observed entries is much smaller than the size of the matrix: $N\ll mT$, whereas $mT$ can be very large. For example, $mT$ is of the order of hundreds of millions for the Netflix problem.

(ii) Column or row masks. If $X_{i}$ has only a small number $d$ of nonzero columns or rows, it is called column or row mask, respectively. We suppose here that $d$ is much smaller than $m$ and $T$. A remarkable case $d=1$ is covering the problem known in Statistics and Econometrics as longitudinal (or panel, or cross-section) data analysis and in Machine Learning as multi-task learning. In what follows, we will designate this problem as multi-task learning, to avoid ambiguity. In the simplest version of multi-task learning, we have $N=nT$ where $T$ is the number of tasks (for instance, in image detection each task $t$ is associated with a particular type of visual object, e.g., face, car, chair, etc.), and $n$ is the number of observations per task. The tasks are characterized by vectors of parameters $a^{*}_{t}\in\mathbb{R}^{m}$, $t=1,\dots,T$, which constitute the columns of matrix $A^{*}$:

The $X_{i}$ are column masks, each containing only one nonzero column $\mathbf{x}^{(t,s)}\in\mathbb{R}^{m}$ (with the convention that $\mathbf{x}^{(t,s)}$ is the $t$th column):

The column $\mathbf{x}^{(t,s)}$ is interpreted as the vector of predictor variables corresponding to $s$th observation for the $t$th task. Thus, for each $i=1,\dots,N$ there exists a pair $(t,s)$ with $t=1,\dots,T,s=1,\dots,n$, such that

If we denote by $Y^{(t,s)}$ and $\xi^{(t,s)}$ the corresponding values $Y_{i}$ and $\xi_{i}$, then the trace regression model (1) can be written as a collection of $T$ standard regression models:

This is the usual formulation of the multi-task learning model in the literature.

For both examples given above, the matrices $X_{i}$ are sparse in the sense that they have only a small portion of nonzero entries. On the other hand, such a sparsity property is not necessarily granted for the target matrix $A^{*}$. Nevertheless, we can always characterize $A^{*}$ by its rank $r=\operatorname{rank}(A^{*})$, and say that a matrix is sparse if it has small rank; cf. Recht, Fazel and Parrilo (2010). For example, the problem of estimation of a square matrix $A^{*}\in\mathbb{R}^{m\times m}$ is a parametric problem which is formally of dimension $m^{2}$ but it has only $(2m-r)r$ free parameters. If $r$ is small as compared to $m$, then the intrinsic dimension of the problem is of the order $rm$. In other words, the rank sparsity assumption $r\ll m$ is a dimension reduction assumption. This assumption will be crucial for the interpretation of our results. Another sparsity assumption that we will consider is that Schatten-$p$ norm of $A^{*}$ (see the definition in Section 2 below) is small for some $0<p\leq 1$. This is an analog of sparsity expressed in terms of the $\ell_{p}$ norm, $0<p\leq 1$, in vector estimation problems.

Estimation of high-dimensional matrices has been recently studied by several authors in settings different from the ours [cf., e.g., Meinshausen and Bühlmann (2006), Bickel and Levina (2008), Ravikumar et al. (2008), Amini and Wainwright (2009), Cai, Zhang and Zhou (2010) and the references cited therein]. Most of attention was devoted to estimation of a large covariance matrix or its inverse. In these papers, sparsity is characterized by the number of nonzero entries of a matrix.

Candès and Recht (2009), Candès and Tao (2009), Gross (2009), Recht (2009) considered the nonnoisy setting ($\xi_{i}\equiv 0$) of the matrix completion problem under conditions that the singular vectors of $A^{*}$ are sufficiently spread out on the unit sphere or “incoherent.” They focused on exact recovery of $A^{*}$. Until now, the sharpest results are those of Gross (2009) and Recht (2009) who showed that under “incoherence condition” the exact recovery is possible with high probability if $N>Cr(m+T)\log^{2}m$ with some constant $C>0$ when we observe $N$ entries of a matrix $A^{*}\in\mathbb{R}^{m\times T}$ with locations uniformly sampled at random. Candès and Plan (2010a), Keshavan, Montanari and Oh (2009) explored the same setting in the presence of noise, proposed estimators $\hat{A}$ of $A^{*}$ and evaluated their Frobenius norm $\|\hat{A}-A^{*}\|_{F}$. The better bounds are in Keshavan, Montanari and Oh (2009) who suggest $\hat{A}$ such that for $A^{*}\in\mathbb{R}^{m\times T}$ and $T=\alpha m$ with $\alpha>1$ the squared error $\|\hat{A}-A^{*}\|_{F}^{2}$ is of the order $\alpha^{5/2}rm^{3}(\log N)/N$ with probability close to 1 when the noise is i.i.d. Gaussian.

In this paper, we consider the general noisy setting of the trace regression problem. We study a class of Schatten-$p$ estimators $\hat{A}$, that is, the penalized least squares estimators with a penalty proportional to Schatten-$p$ norm; cf. (7). The special case $p=1$ corresponds to the “matrix Lasso.” We study the convergence properties of their prediction error

and of their Schatten-$q$ error. The main contributions of this paper are the following.

1. [(a)]

2. (a)

   For all $0<p\leq 1$, under various assumptions on the masks $X_{i}$ (no assumption, USR matrix completion, CS matrix completionmatrix completionmatrix compl) we obtain different bounds on the prediction error of Schatten-$p$ estimators involving the Schatten-$p$ norm of $A^{*}$.

3. (b)

   For $p$ sufficiently close to 0, under a mild assumption on $X_{i}$, we show that Schatten-$p$ estimators achieve the prediction error rate of convergence $\frac{r\max(m,T)}{N}$, up to a logarithmic factor. This result is valid for matrices $A^{*}$ whose eigenvalues are not exponentially large in $N$. It covers the matrix completion and high-dimensional multi-task learning problems.

4. (c)

   For all $0<p\leq 1$, we obtain upper bounds for the prediction error under the matrix Restricted Isometry (RI) condition on the masks $X_{i}$, which is a rather strong condition, and under the assumption that ${\rm rank}(A^{*})\leq r$. We also derive the bounds for the Schatten-$q$ error of ${\hat{A}}$. The rate in the bounds for the prediction error is $r\max(m,T)/N$ when the RI condition is satisfied with scaling factor 1 (i.e., for the case not related to matrix completion and high-dimensional multi-task learning).

5. (d)

   We prove the lower bounds showing that the rate $r\max(m,T)/N$ is minimax optimal for the prediction error and Schatten-$2$ (i.e., Frobenius) norm estimation error under the RI condition on the class of matrices $A^{*}$ of rank smaller than $r$. Our result is even more general because we prove our lower bound on the intersection of the Schatten-0 ball with the Schatten-$p$ ball for any $0<p\leq 1$, which allows us to show minimax optimality of the upper bounds of (a) as well. Furthermore, we prove minimax lower bounds for collaborative sampling and USR matrix completion problems.

The main point of this paper is to show that the suitably tuned Schatten estimators attain the optimal rate of prediction error up to logarithmic factors. The striking fact is that we can achieve this not only under the very restrictive assumption, such as the RI condition, but also under very mild assumptions on the masks $X_{i}$.

Finally, it is useful to compare the results for matrix estimation when the sparsity is expressed by the rank with those for the high-dimensional vector estimation when the sparsity is expressed by the number of nonzero components of the vector. For the vector estimation, we have the linear model

where $X_{i}\in\mathbb{R}^{p}$, $\beta\in\mathbb{R}^{p}$ and, for example, $\xi_{i}$ are i.i.d. $\mathcal{N}(0,1)$ random variables. Consider the high-dimensional case $p\gg N$. (This is analogous to the assumption $m^{2}\gg N$ in the matrix problem and means that the nominal dimension is much larger than the sample size.) The sparsity assumption for the vector case has the form $s\ll N$, where $s$ is the number of nonzero components, or the intrinsic dimension of $\beta$. Let $\hat{\beta}$ be an estimator of $\beta$. Then the optimal rate of convergence of the prediction risk $N^{-1}\sum_{i=1}^{N}(X_{i}^{\prime}(\hat{\beta}-\beta))^{2}$ on the class of vectors $\beta$ with given $s$ is of the order $s/N$, up to logarithmic factors. This rate is shown to be attained, up to logarithmic factors, for many estimators, such as the BIC, the Lasso, the Dantzig selector, Sparse Exponential Weighting, etc.; cf., for example, Bunea, Tsybakov and Wegkamp (2007), Koltchinskii (2008), Bickel, Ritov and Tsybakov (2009), Dalalyan and Tsybakov (2008). Note that this rate is of the form $\frac{\mathrm{intrinsic\ dimension}}{\mathrm{sample\ size}}=\frac{s}{N}$, up to a logarithmic factor. The general interpretation is therefore completely analogous to that of the matrix case: Assume for simplicity that $A^{*}$ is a square $m\times m$ matrix with $\operatorname{rank}(A^{*})=r$. As mentioned above, the intrinsic dimension (the number of parameters to be estimated to recover $A^{*}$) is then $(2m-r)r$, which is of the order $\sim rm$ if $r\ll m$. An interesting difference is that the logarithmic risk inflation factor is inevitable in the vector case [cf. Donoho et al. (1992), Foster and George (1994)], but not in the matrix problem, as our results reveal.

This paper is organized as follows. In Section 2, we introduce notation, some definitions, basic facts about the Schatten quasi-norms and define the Schatten-$p$ estimators. Section 3 describes elementary steps in their convergence analysis and presents two general approaches to upper bounds on the estimation and prediction error (cf. Theorems 1 and 2) depending on the efficient noise level $\tau$. Our main results are stated in Sections 4, 5 (matrix completion), 6 (multi-task learning). They are obtained from Theorems 1 and 2 by specifying the effective noise level $\tau$ under particular assumptions on the masks $X_{i}$. Concentration bounds for certain random matrices leading to the expressions for the effective noise level are collected in Section 8. Section 7 is devoted to minimax lower bounds. Sections 9 and 10 contain the main proofs. Finally, in Section 11 we establish bounds for the $k$th entropy numbers of the quasi-convex Schatten class embeddings $S_{p}^{M}\hookrightarrow S_{2}^{M}$, $p<1$, which are needed for our proofs and are of independent interest.

In this section, we discuss two schemes of proving upper bounds on the prediction error of $\hat{A}$. The first bound involves only the Schatten-$p$ norm of matrix $A^{*}$. The second involves only the rank of $A^{*}$ but needs the RI condition on the sampling operator.

We start by sketching elementary steps in the convergence analysis of Schatten-$p$ estimators. By the definition of $\hat{A}$,

Recalling that $Y_{i}=\operatorname{tr}(X_{i}^{\prime}A^{*})+\xi_{i}$, we can transform this by a simple algebra to

In the sequel, inequality (8) will be referred to as basic inequality and the random variable $N^{-1}\sum_{i=1}^{N}\xi_{i}\operatorname{tr}((\hat{A}-A^{*})^{\prime}X_{i})$ will be called the stochastic term. The core in the analysis of Schatten-$p$ estimators consists in proving tight bounds for the right-hand side of the basic inequality (8). For this purpose, we first need a control of the stochastic term. Section 8 below demonstrates that such a control strongly depends on the properties of $\mathcal{L}$, that is, of the problem at hand. In summary, Section 8 establishes that, under suitable conditions, for any $0<p\leq 1$ the stochastic term can be bounded for all $\delta>0$ with probability close to 1 as follows:

where $0<\tau<\infty$ depends on $m,T$ and $N$. The quantity $\tau$ plays a crucial role in this bound. We will call $\tau$ the effective noise level. Exact expressions for $\tau$ under various assumptions on the sampling operator $\mathcal{L}$ and on the noise $\xi_{i}$ are derived in Section 8. In Table 1, we present the values of $\tau$ for three important examples under the assumption that $\xi_{i}$ are i.i.d. Gaussian $\mathcal{N}(0,\sigma^{2})$ random variables. In the cases listed in Table 1, inequality (3) holds with probability $1-\varepsilon$, where $\varepsilon=(1/C)\exp(-C(m+T))$ (first and third example) and $\varepsilon=(1/C^{\prime})(\max(m,T)+1)^{-C^{\prime}}$ (second example) with constants $C,C^{\prime}>0$ independent of $N,m,T$.

The following two points will be important to understand the subsequent results:

• •

  In this paper, we will always choose the regularization parameter $\lambda$ in the form $\lambda=4\tau$.

• •

  With this choice of $\lambda$, the effective noise level $\tau$ characterizes the rate of convergence of the Schatten estimator. The smaller is $\tau$, the faster is the rate.

In particular, the first line in Table 1 reveals that when $M=\max(m,T)<N$ the largest $\tau$ corresponds to $p=1$ and it becomes smaller when $p$ decreases to 0. This suggests that choosing Schatten-$p$ estimators with $p<1$ and especially $p$ close to 0 might be advantageous. Note that the assumption of uniform boundedness of $\mathcal{L}$ is very mild. For example, it is trivially satisfied with $c_{0}=1$ for USR matrix completion and with $c_{0}=1/N$ for CS matrix completion. However, in these two cases a specific analysis leads to sharper bounds on the effective noise level (i.e., on the rate of convergence of the estimators); cf. the second and third lines of Table 1.

In this section, we provide two bounds on the prediction error of $\hat{A}$ with a general effective noise level $\tau$. We then detail them in Sections 4–6 for particular values of $\tau$ depending on the assumptions on the $X_{i}$. The first bound involves the Schatten-$p$ norm of matrix $A^{*}$.

From (8) and (3) with $\delta=1/2$ and $\lambda=4\tau$, we get

This and the $p$-norm inequality (3) yield (10).

The bound (10) depends on the magnitude of the elements of $A^{*}$ via $\|A^{*}\|_{S_{p}}$. The next theorem shows that under the RI condition this dependence can be avoided, and only the rank of $A^{*}$ affects the rate of convergence.

Proof of Theorem 2 is given in Section 9. The values $a=a(p)$ and $\delta_{0}(p)$ can be deduced from the proof. In particular, for $p=1$ it is sufficient to take $a=19$.

Discussion about the scaling factor $\nu$. Remark 1 deals with the case $\nu=1$, which seems to be not always appropriate for trace regression models. To our knowledge, the only available examples of matrices $X$ such that the sampling operator $\mathcal{L}$ satisfies the RI condition with $\nu=1$ are complete matrices, that is, matrices with all nonzero entries, which are random and have specific distributions [typically, i.i.d. Rademacher or Gaussian entries; cf. Recht, Fazel and Parrilo (2010)]. Except for degenerate cases [such as $N=mT$, the $X_{i}$ distinct and of the form $\sqrt{N}e_{k}(m)e_{l}(T)^{\prime}$ for $1\leq k\leq m,1\leq l\leq T$] the sampling operator $\mathcal{L}$ defines typically a restricted isometry with $\nu=1$ only if the matrices $X_{i}$ contain a considerable number of (uniformly bounded) nonzero entries.

Let us now specify the form of the RI condition in the context of multi-task learning discussed in the Introduction. Using (2) for a matrix $A=(a_{1}\cdots a_{T})$, we obtain

where $\Psi_{t}=n^{-1}\sum_{s=1}^{n}\mathbf{x}^{(t,s)}(\mathbf{x}^{(t,s)})^{\prime}$ is the Gram matrix of predictors for the $t$th task. These matrices correspond to $T$ separate regression models. The standard assumption is that they are normalized so that all the diagonal elements of each $\Psi_{t}$ are equal to 1. This suggests that the natural RI scaling factor $\nu$ for such model is of the order $\nu\sim\sqrt{T}$. For example, in the simplest case when all the matrices $\Psi_{t}$ are just equal to the $m\times m$ identity matrix, we find $|\mathcal{L}(A)|_{2}^{2}=T^{-1}\sum_{t=1}^{T}a_{t}^{\prime}\Psi_{t}a_{t}=T^{-1}\|A\|_{S_{2}}^{2}.$ Similarly, we get the RI condition with scaling factor $\nu\sim\sqrt{T}$ when the spectra of all the Gram matrices $\Psi_{t},t=1,\dots,T,$ are included in a fixed interval $[a,b]$ with $0<a<b<\infty$. However, this excludes the high-dimensional task regressions, such that the number of parameters $m$ is larger than the sample size, $m>n$. In conclusion, application of the matrix RI techniques in multi-task learning is restricted to low-dimensional regression and the scaling factor is $\nu\sim\sqrt{T}$.

The reason for the failure of the RI approach is that the masks $X_{i}$ are sparse. The sparser are $X_{i}$, the larger is $\nu$. The extreme situation corresponds to matrix completion problems. Indeed, if $N<mT$, then there exists a matrix of rank $1$ in the null-space of the sampling operator $\mathcal{L}$ and hence the RI condition cannot be satisfied. For $N\geq mT$ we can have the RI condition with scaling factor $\nu\sim\sqrt{mT}$, but $N\geq mT$ means that essentially all the entries are observed, so that the very problem of completion does not arise.

The above discussion suggests that Theorem 2 and, in general, the argument based on the restricted isometry or related conditions are not well adapted for several interesting settings. Motivated by this, we propose another approach described in the next theorem, which requires only the comparably mild uniform boundedness condition (5). For simplicity, we focus on Gaussian errors $\xi_{i}$. Set $M=\max(m,T)$.

Inequality (3) holds with probability at least $1-C\exp(-\vartheta M/\break C^{2})$ by Lemma 5. We then use (10) and note that, under our choice of $p$, $\tau\leq c\vartheta M/(Np)$ for some constant $c<\infty$, which does not depend on $M$ and $N$, and

Finally, we give the following theorem quantifying the rates of convergence of the prediction risk in terms of the Schatten norms of $A^{*}$. Its proof is straightforward in view of Theorem 1 and Lemmas 2 and 5.

In Theorem 5 below we show that these rates are optimal in a minimax sense on the corresponding Schatten-$p$ balls for the sampling operators satisfying the RI condition.

As discussed in Section 3, for matrix completion problems the restricted isometry argument as in Theorem 2 is not applicable. We will therefore use Theorems 1 and 3. First, combining Theorem 1 with Lemma 3 of Section 8 we get the following corollary.

Next, combining Theorem 3 with Lemma 5 of Section 8 we get the following corollary.

Note that the bounds of Corollaries 1(i) and 2 achieve the rate $r\max(m,\break T)/N$, up to logarithmic factors under different conditions on the maximal singular value of $A^{*}$. If $\max(m,T)<N<mT$ then the condition in Corollary 2 does not imply more than a polynomial in $\max(m,T)$ growth on $\sigma_{1}(A^{*})$, which is a mild assumption. On the other hand, (17) requires uniform boundedness of $\sigma_{1}(A^{*})$ by some constant to achieve the same rate. However, the estimators of Corollary 2 correspond to nonconvex penalty and are computationally hard.