Estimation of (near) low-rank matrices with noise and high-dimensional scaling
Sahand Negahban, Martin J. Wainwright
Introduction
This initial set of results, though appealing in terms of their simple statements and generality, are somewhat abstractly formulated. Our next contribution is to show that by specializing our main result (Theorem 1) to three classes of models, we can obtain some concrete results based on readily interpretable conditions. In particular, Corollary 3 deals with the case of low-rank multivariate regression, relevant for applications in multitask learning. We show that the random operator satisfies the RSC property for a broad class of observation models, and we use random matrix theory to provide an appropriate choice of the regularization parameter. Our next result, Corollary 4, deals with the case of estimating the matrix of parameters specifying a vector autoregressive (VAR) process . Here we also establish that a suitable RSC property holds with high probability for the random operator , and also specify a suitable choice of the regularization parameter. We note that the technical details here are considerably more subtle than the case of low-rank multivariate regression, due to dependencies introduced by the autoregressive sampling scheme. Accordingly, in addition to terms that involve the size, the matrix dimensions and rank, our bounds also depend on the mixing rate of the VAR process. Finally, we turn to the compressed sensing observation model for low-rank matrix recovery, as introduced by Recht et al. . In this setting, we again establish that the RSC property holds with high probability, specify a suitable choice of the regularization parameter, and thereby obtain a Frobenius error bound for noisy observations (Corollary 5). A technical result that we prove en route—namely, Proposition 1—is of possible independent interest, since it provides a bound on the constrained norm of a random Gaussian operator. In particular, this proposition allows us to obtain a sharp result (Corollary 6) for the problem of recovering a low-rank matrix from perfectly observed random projections, one that removes a logarithmic factor from past work .
The remainder of this paper is organized as follows. Section 2 is devoted to background material, and the set-up of the problem. We present a generic observation model for low-rank matrices, and then illustrate how it captures various cases of interest. We then define the convex program based on nuclear norm regularization that we analyze in this paper. In Section 3, we state our main theoretical results and discuss their consequences for different model classes. Section 4 is devoted to the proofs of our results; in each case, we break down the key steps in a series of lemmas, with more technical details deferred to the appendices. In Section 5, we present the results of various simulations that illustrate excellent agreement between the theoretical bounds and empirical behavior.
Background and problem set-up
We begin with some background on problems and applications in which rank constraints arise, before describing a generic observation model. We then introduce the semidefinite program (SDP) based on nuclear norm regularization that we study in this paper.
where is an i.i.d. sequence of -dimensional zero-mean noise vectors. Given a collection of observations of covariate-output pairs, our goal is to estimate the unknown matrix . This type of model has been used in many applications, including analysis of fMRI image data , analysis of EEG data decoding , neural response modeling and analysis of financial data. This model and closely related ones also arise in the problem of collaborative filtering , in which the goal is to predict users’ preferences for items (such as movies or music) based on their and other users’ ratings of related items. The papers discuss additional instances of low-rank decompositions. In all of these settings, the low-rank condition translates into the existence of a smaller set of “features” that are actually controlling the prediction.
2 A generic observation model
which is specified by the sequence of observation matrices and observation noise . This observation model can be written in a more compact manner using operator-theoretic notation. In particular, let us define the observation vector
Let us illustrate the form of the observation model (5) for some of the applications that we considered earlier.
By re-indexing this collection of observations via the mapping , we recognize multivariate regression as an instance of the observation model (4) with observation matrix and scalar observation .
Recall that a vector autoregressive (VAR) process is defined by the recursion (2), and suppose that observe an -sequence produced by this recursion. Since each is -variate, the scalarized sample size is . Letting index the dimension, we have
In this case, we re-index the collection of observations via the mapping . After doing so, we see that the autoregressive problem can be written in the form (4) with and observation matrix .
3 Regression with nuclear norm regularization
where is a regularization parameter. Note that the optimization problem (9) can be viewed as the analog of the Lasso estimator , tailored to low-rank matrices as opposed to sparse vectors. An important property of the optimization problem (9) is that it can be solved in time polynomial in the sample size and the matrix dimensions and . Indeed, the optimization problem (9) is an instance of a semidefinite program , a class of convex optimization problems that can be solved efficiently by various polynomial-time algorithms . For instance, interior point methods are a classical method for solving semidefinite programs; moreover, as we discuss in Section 5, there are a variety of other methods for solving the semidefinite program (SDP) defining our -estimator
Like in any typical -estimator for statistical inference, the regularization parameter is specified by the statistician. As part of the theoretical results in the next section, we provide suitable choices of this parameter in order for the estimate to behave well, in the sense of being close in Frobenius norm to the unknown matrix .
Main results and some consequences
In this section, we state our main results and discuss some of their consequences. Section 3.1 is devoted to results that apply to generic instances of low-rank problems, whereas Section 3.2 is devoted to the consequences of these results for more specific problem classes, including low-rank multivariate regression, estimation of vector autoregressive processes, and recovery of low-rank matrices from random projections.
We note that analogous conditions have been used to establish error bounds in the context of sparse linear regression , in which case the set corresponded to certain subsets of sparse vectors.
With this notation, we come to the first result of our paper. It is a deterministic result, which specifies two conditions—namely, an RSC condition and a choice of the regularizer—that suffice to guarantee for any solution of the SDP (9) fall within a certain radius.
Suppose that the operator satisfies restricted strong convexity with parameter over the set , and that the regularization parameter is chosen such that . Then any solution to the semidefinite program (9) satisfies
Apart from the tolerance parameter , the two main terms in the bound (14) have a natural interpretation. The first term (involving ) corresponds to estimation error, capturing the difficulty of estimating a rank matrix. The second is an approximation error, in which the projection onto the set describes the gap between the true matrix and the rank approximation.
Let us begin by illustrating the consequences of Theorem 1 when the true matrix has exactly rank , in which case there is a very natural choice of the subspaces represented by and . In particular, we form from the non-zero left singular vectors of , and from its non-zero right singular vectors. Note that this choice of ensures that . For technical reasons to be clarified, it suffices to set in the case of exact rank constraints, and we thus obtain the following result:
Suppose that has rank , and satisfies RSC with respect to . Then as long as , any optimal solution to the SDP (9) satisfies the bound
Like Theorem 1, Corollary 1 is a deterministic statement on the SDP error. It takes a much simpler form since when is exactly low rank, then neither tolerance parameter nor the approximation term are required.
As a more delicate example, suppose instead that is nearly low-rank, an assumption that we can formalize by requiring that its singular value sequence decays quickly enough. In particular, for a parameter and a positive radius , we define the set
Note that the error bound (17) reduces to the exact low rank case (15) when , and . The quantity acts as the “effective rank” in this setting; as clarified by our proof in Section 4.2. This particular choice is designed to provide an optimal trade-off between the approximation and estimation error terms in Theorem 1. Since is chosen to decay to zero as the sample size increases, this effective rank will increase, reflecting the fact that as we obtain more samples, we can afford to estimate more of the smaller singular values of the matrix .
2 Results for specific model classes
As stated, Corollaries 1 and 2 are fairly abstract in nature. More importantly, it is not immediately clear how the key underlying assumption—namely, the RSC condition—can be verified, since it is specified via subspaces that depend on , which is itself the unknown quantity that we are trying to estimate. Nonetheless, we now show how, when specialized to more concrete models, these results yield concrete and readily interpretable results. As will be clear in the proofs of these results, each corollary requires overcoming two main technical obstacles: establishing that the appropriate form of the RSC property holds in a uniform sense (so that a priori knowledge of is not required), and specifying an appropriate choice of the regularization parameter . Each of these two steps is non-trivial, requiring some random matrix theory, but the end results are simply stated upper bounds that hold with high probability.
with probability greater than .
Remarks: Corollary 3 takes a particularly simple form when : then there exists a constant such that |\!|\!|\widehat{\Theta}-\Theta^{*}|\!|\!|_{{F}}^{2}\leq c^{\prime}_{1}\nu^{2}\>R_{q}\;\big{(}\frac{k+p}{n}\big{)}^{1-q/2}. When is exactly low rank—that is, , and has rank —this simplifies even further to
The scaling in this error bound is easily interpretable: naturally, the squared error is proportional to the noise variance , and the quantity counts the number of degrees of freedom of a matrix with rank . Note that if we did not impose any constraints on , then since a matrix has a total of free parameters, we would expect at best to obtain rates of the order . Note that when is low rank—in particular, when —then the nuclear norm estimator achieves substantially faster rates. Finally, we note that as stated, the result requires that tend to infinity in order for the claim to hold with high probability. Although such high-dimensional scaling is the primary focus of this paper, we note that for application to the classical setting of fixed , the same statement (with different constants) holds with replaced by .
Next we turn to the case of estimating the system matrix of an autoregressive (AR) model, as discussed in Example 2.
with probability greater than .
Remarks: Like Corollary 3, the result as stated requires that tend to infinity, but the same bounds hold with replaced by , yielding results suitable for classical (fixed dimension) scaling. Second, the factor , like the analogous termThe term in Corollary 3 has a factor , since the matrix in that case could be non-square in general. in Corollary 3, shows that faster rates are obtained if can be well-approximated by a low rank matrix, namely for choices of the parameter that are closer to zero. Indeed, in the limit , we again reduce to the case of an exact rank constraint , and the corresponding squared error scales as . In contrast to the case of multivariate regression, the error bound (19) also depends on the upper bound on the operator norm of the system matrix . Such dependence is to be expected since the quantity controls the (in)stability and mixing rate of the autoregressive process. As clarified in the proof, the dependence of the sampling in the AR model also presents some technical challenges not present in the setting of multivariate regression.
with probability greater than .
The central challenge in proving this result is in proving an appropriate form of the RSC property. The following result on the random operator may be of independent interest here:
Under the stated conditions, the random operator satisfies
with probability at least .
The proof of this result, provided in Appendix D, makes use of the Gordon-Slepian inequalities for Gaussian processes, and concentration of measure. As we show in Section 4.5, it implies the form of the RSC property needed to establish Corollary 5.
Proposition 1 also implies an interesting property of the null space of the operator ; one that can be used to establish a corollary about recovery of low-rank matrices when the observations are noiseless. In particular, suppose that we are given the noiseless observations for , and that we try to recover the unknown matrix by solving the SDP
a recovery procedure that was studied by Recht et al. . Proposition 1 allows us to obtain a sharp result on recovery using this method:
Suppose that has rank , and that we are given noiseless samples. Then with probability at least , the SDP (22) recovers the matrix exactly.
This result removes some extra logarithmic factors that were included in the earlier work , and provides the appropriate analog to compressed sensing results for sparse vectors . Note that (like in most of our results) we have made little effort to obtain good constants in this result: the important property is that the sample size scales linearly in both and .
Proofs
We now turn to the proofs of Theorem 1, and Corollaries 1 through 6. In each case, we provide the primary steps in the main text, with more technical details stated as lemmas and proved in the Appendix.
By the optimality of for the SDP (9), we have
Defining the error matrix and performing some algebra yields the inequality
By definition of the adjoint and Hölder’s inequality, we have
By the triangle inequality, we have . Substituting this inequality and the bound (24) into the inequality (23) yields
where the second inequality makes use of our choice .
It remains to lower bound the term on the left-hand side, while upper bounding the quantity on the right-hand side. The following technical result allows us to do so. Recall our earlier definition (11) of the sets and associated with a given subspace pair.
Let represent a pair of -dimensional subspaces of left and right singular vectors of . Then there exists a matrix decomposition of the error such that
The matrix satisfies the constraint , and
If , then the nuclear norm of is bounded as
See Appendix A for the proof of this claim. Using Lemma 1, we can complete the proof of the theorem. In particular, from the bound (25) and the RSC assumption, we find that
Using the triangle inequality together with inequality (25), we obtain
From the rank constraint in Lemma 1(a), we have . Putting together the pieces, we find
2 Proof of Corollary 2
On the other hand, we also have , which implies that . From the general error bound with , we obtain
Setting yields that
3 Proof of Corollary 3
For the proof of this corollary, we adopt the following notation. We first define the three matrices
With this notation and using the relation , the SDP objective function (9) can be written as \frac{1}{k}\big{\{}\frac{1}{2n}|\!|\!|Y-X\Theta|\!|\!|_{{F}}^{2}+\lambda_{n}|\!|\!|\Theta|\!|\!|_{{1}}\big{\}}, where we have defined .
In order to establish the RSC property for this model, some algebra shows that we need to establish a lower bound on the quantity
where denotes the minimum eigenvalue. The following lemma follows by adapting known concentration results for random matrices (see the paper for details):
As a consequence, we have with probability at least for all , which establishes that the RSC property holds with .
Next we need upper bound the quantity for this model, so as to verify that the stated choice for is valid. Following some algebra, we find that
The following lemma is proved in Appendix B:
Using these two lemmas, we can complete the proof of Corollary 3. First, recalling the scaling , we see that Lemma 3 implies that the choice satisfies the conditions of Corollary 2 with high probability. Lemma 2 shows that the RSC property holds with , again with high probability. Consequently, Corollary 2 implies that
with probability greater than , as claimed.
4 Proof of Corollary 4
For the proof of this corollary, we adopt the notation
The following lemma provides the lower bound needed to establish RSC for the autoregressive model:
The eigenspectrum of the matrix is well-controlled in terms of the stationary covariance matrix: in particular, as long as , we have
both with probability greater than .
Thus, from the bound (28)(b), we see with the high probability, the RSC property holds with as long as .
As before, in order to verify the choice of , we need to control the quantity . The following inequality, proved in Appendix C.2, yields a suitable upper bound:
There exist constants , independent of etc. such that
From Lemma 5, we see that it suffices to choose . With this choice, Corollary 2 of Theorem 1 yields that
with probability greater than , as claimed.
5 Proof of Corollary 5
Let us now show how Proposition 1 implies the RSC property with an appropriate tolerance parameter. In particular, let us define \delta^{2}\;:=\;R_{q}\,\big{[}\sqrt{\frac{k}{N}}+\sqrt{\frac{p}{N}}\big{]}^{2-q}, so that if we have the inequality , the result of Corollary 5 follows immediately. Therefore, we may take . Now recall from Lemma 1 that the error satisfies the bound (25). Combining these facts, we are guaranteed that , where the set was previously defined (12), and it is sufficient to establish the RSC property over this set.
Observe that the bound (21) implies that for any ,
Following the arguments used in the proofs of Theorem 1 and Corollary 2, we find that
where is a parameter to be chosen. We now set \tau\,=\,\big{(}\sqrt{k}+\sqrt{p}\,\big{)}/\sqrt{N}, and substitute the resulting bound (31) into equation (30), thereby obtaining
If we choose , then we are guaranteed that , which shows that the RSC property holds with .
The next step is to control the quantity , required for specifying a suitable choice of .
If , then
6 Proof of Corollary 6
This corollary follows from a combination of Proposition 1 and Lemma 1. Let be an optimal solution to the SDP (22), and let be the error. Since is optimal and is feasible for the SDP, we have . Using the decomposition from Lemma 1 and applying triangle inequality, we have . From the properties of the decomposition in Lemma 1 (see Appendix A), we find that
Combining the pieces yields that , and hence . By Lemma 1(a), the rank of is at most , so that we obtain .
Note that , since both and agree with the observations. Consequently, from Proposition 1, we have that
where the final inequality follows from the assumption that . We have thus shown that , which implies that as claimed.
Experimental results
In this section, we report the results of various simulations that demonstrate the close agreement between the scaling predicted by our theory, and the actual behavior of the SDP-based -estimator (9) in practice. In all cases, we solved the convex program (9) by using our own implementation in MATLAB of an accelerated gradient descent method which adapts a non-smooth convex optimization procedure to the nuclear-norm . We chose the regularization parameter in the manner suggested by our theoretical results; in doing so, we assumed knowledge of quantities such as the noise variance . (In practice, one would have to estimate such quantities from the data using standard methods.)
Figure 1 shows results for a multivariate regression model with the covariates chosen randomly from a distribution. Panel (a) plots the Frobenius error on a logarithmic scale versus the sample size for three different matrix sizes, . Naturally, in each case, the error decays to zero as increases, but larger matrices require larger sample sizes, as reflected by the rightward shift of the curves as is increased. Panel (b) of Figure 1 shows the exact same set of simulation results, but now with the Frobenius error plotted versus the rescaled sample size . As predicted by Corollary 3, the error plots now are all aligned with one another; the degree of alignment in this particular case is so close that the three plots are now indistinguishable. (The blue curve is the only one visible since it was plotted last by our routine.) Consequently, Figure 1 shows that acts as the effective sample size in this high-dimensional setting.
Figure 2 shows similar results for the autoregressive model discussed in Example 2. As shown in panel (a), the Frobenius error again decays as the sample size is increased, although problems involving larger matrices are shifted to the right. Panel (b) shows the same Frobenius error plotted versus the rescaled sample size ; as predicted by Corollary 4, the errors for different matrix sizes are again quite well-aligned. In this case, we find (both in our theoretical analysis and experimental results) that the dependence in the autoregressive process slows down the rate at which the concentration occurs, so that the results are not as crisp as the low-rank multivariate setting in Figure 1.
Finally, Figure 3 presents the same set of results for the compressed sensing observation model discussed in Example 3. Even though the observation matrices here are qualitatively different (in comparison to the multivariate regression and autoregressive examples), we again see the “stacking” phenomenon of the curves when plotted versus the rescaled sample size , as predicted by Corollary 5.
Discussion
In this paper, we have analyzed the nuclear norm relaxation for a general class of noisy observation models, and obtained non-asymptotic error bounds on the Frobenius norm that hold under high-dimensional scaling. In contrast to most past work, our results are applicable to both exactly and approximately low-rank matrices. We stated a main theorem that provides high-dimensional rates in a fairly general setting, and then showed how by specializing this result to some specific model classes—namely, low-rank multivariate regression, estimation of autoregressive processes, and matrix recovery from random projections—it yields concrete and readily interpretable rates. Lastly, we provided some simulation results that showed excellent agreement with the predictions from our theory.
Acknowledgements
This work was partially supported by a Sloan Foundation Fellowship, AFOSR-09NL184 grant, and an NSF-CCF-0545862 CAREER grant to MJW.
Appendix A Proof of Lemma 1
which establishes Lemma 1(a). Moreover, we note for future reference that by construction of , the nuclear norm satisfies the decomposition
We now turn to the proof of Lemma 1(b). Recall that the error associated with any optimal solution must satisfy the inequality (23), which implies that
Using the triangle inequality and the relation (33), we have
Substituting this inequality into the bound (34), we obtain
Finally, since by assumption, we conclude that
Appendix B Proof of Lemma 3
For positive scalars and , define the (random) quantity
and note that our goal is to upper bound . Note moreover that , a relation which will be useful in the analysis.
Let and denote coverings of and , respectively. We now claim that we have the upper bound
To establish this claim, we note that since the sets and are -covers, for any pair , there exists a pair such that and , with . Consequently, we can write
By construction, we have the bound , and similarly as well as . Substituting these bounds into the decomposition (36) and taking suprema over the left and right-hand sides, we conclude that
We now apply the union bound to control the discrete maximum. It is known (e.g., ) that there exists a covering of and with at most and elements respectively. Consequently, we have
Combining this tail bound with the upper bound (37), we have
Setting , this probability vanishes as long as .
Appendix C Technical details for Corollary 4
In this appendix, we collect the proofs of Lemmas 4 and 5.
Recalling that denotes the unit-norm Euclidean sphere in -dimensions, we first observe that . Our next step is to reduce the supremum to a maximization over a finite set, using a standard covering argument. Let denote a -cover of it. By definition, for any , there is some such that , where . Consequently, for any , the triangle inequality implies that
and hence that . Re-arranging yields the useful inequality
where the last inequality follows from the union bound, and the fact that there exists a -covering of with at most elements.
Moreover, we have . Applying Lemma 8 with , we conclude that
with probability at least , which establishes the upper bound (28)(a).
since . Since |\Psi(\Delta v,v)|\leq\epsilon\,|\!|\!|\Big{(}\frac{1}{n}X^{T}X\Big{)}|\!|\!|_{{\operatorname{op}}}, we obtain the lower bound
By the previously established upper bound(28)(a), have with high probability. Hence, choosing ensures that .
Consequently, it suffices to lower bound the minimum over the covering set. We first establish a concentration result for the function that holds for any fixed . Note that we can write
Note that this bound holds for any fixed . Setting and applying the union bound yields that
which vanishes as long as .
C.2 Proof of Lemma 5
We now apply the union bound to control the discrete maximum. It is known (e.g., ) that there exists a covering of with at most elements. Consequently, we have
For each , let and denote the row of and . Following some simple algebra, we have the decomposition , where
We may now bound each in turn; in doing so, we make repeated use of Lemma 8, which provides concentration bounds for a random variable of the form , where for some matrix .
We begin with , which the easiest to control since (up to scaling by ), it corresponds to the deviation away from the mean of -variable with degrees of freedom. Consequently, applying Lemma 8 with , we obtain
where is the Kronecker delta for the event . As before, by symmetry of , we have , and hence
Morever, we have , so that by applying Lemma 8, we conclude that
which completes the analysis of this term.
Combining the bounds (45), (44) and (46), we conclude that for all ,
Setting and combining with the bound (43), we conclude that
Appendix D Proof of Proposition 1
In particular, our goal is to prove that for any , the lower bound
holds with probability at least . By a standard peeling argument (see Raskutti et al. for details), this lower bound implies the claim (21).
We establish the lower bound (48) using Gaussian comparison inequalities and concentration of measure (see Lemma 7). For each pair , consider the random variable , and note that it is Gaussian with zero mean. For any two pairs and , some calculation yields
We now define a second Gaussian process via
It can be shown that for all pairs , we have
Moreover, equality holds whenever . The conditions of the Gordon-Slepian inequality are satisfied, so that we are guaranteed that
Finally, we need to establish sharp concentration around the mean. Note that the function is Lipschitz with constant , so that Lemma 7 implies that
Appendix E Some useful concentration results
The following lemma is classical , and yields sharp concentration of a Lipschitz function of Gaussian random variables around its mean.
Given a Gaussian random vector , for all , we have
Let be the symmetrix matrix square root, and consider the function . Since it is Lipschitz with constant , Lemma 7 implies that
By integrating this tail bound, we find that the variable satisfies the bound , and hence conclude that
Combining this bound with the tail bound (54), we conclude that
Setting in the bound (56) yields that
Similarly, setting in the tail bound (56) yields that with probability greater than , we have