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 X\mathfrak{X} 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 X\mathfrak{X}, 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 {Wa}a=1n\{W_{a}\}_{a=1}^{n} is an i.i.d. sequence of kk-dimensional zero-mean noise vectors. Given a collection of observations {Za,Ya}a=1n\{Z_{a},Y_{a}\}_{a=1}^{n} of covariate-output pairs, our goal is to estimate the unknown matrix Θ∗\Theta^{*}. 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 {Xi}i=1N\{X_{i}\}_{i=1}^{N} and observation noise {εi}i=1N\{\varepsilon_{i}\}_{i=1}^{N}. 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 N=nkN=nk observations via the mapping (a,b)↦i=a+(b−1) k(a,b)\mapsto i=a+(b-1)\,k, we recognize multivariate regression as an instance of the observation model (4) with observation matrix Xi=ZaebTX_{i}=Z_{a}e_{b}^{T} and scalar observation yi=Vaby_{i}=V_{ab}.

Recall that a vector autoregressive (VAR) process is defined by the recursion (2), and suppose that observe an nn-sequence {Zt}t=1n\{Z_{t}\}_{t=1}^{n} produced by this recursion. Since each Zt=[Zt1…Ztp]TZ_{t}=\begin{bmatrix}Z_{t1}&\ldots&Z_{tp}\end{bmatrix}^{T} is pp-variate, the scalarized sample size is N=npN=np. Letting b=1,2,…,pb=1,2,\ldots,p index the dimension, we have

In this case, we re-index the collection of N=npN=np observations via the mapping (t,b)↦i=t+(b−1) p(t,b)\mapsto i=t+(b-1)\,p. After doing so, we see that the autoregressive problem can be written in the form (4) with yi=Z(t+1) by_{i}=Z_{(t+1)\,b} and observation matrix Xi=ZtebTX_{i}=Z_{t}e_{b}^{T}.

3 Regression with nuclear norm regularization

where λN>0\lambda_{N}>0 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 NN and the matrix dimensions kk and pp. 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 MM-estimator

Like in any typical MM-estimator for statistical inference, the regularization parameter λN\lambda_{N} 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 Θ^\widehat{\Theta} to behave well, in the sense of being close in Frobenius norm to the unknown matrix Θ∗\Theta^{*}.

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 C\mathcal{C} 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 X\mathfrak{X} satisfies restricted strong convexity with parameter κ(X)>0\kappa(\mathfrak{X})>0 over the set C(r;δ)\mathcal{C}(r;\delta), and that the regularization parameter λN\lambda_{N} is chosen such that λN≥2∣ ⁣∣ ⁣∣X∗(ε⃗)∣ ⁣∣ ⁣∣op⁡/N\lambda_{N}\geq 2|\!|\!|\mathfrak{X}^{*}(\vec{\varepsilon})|\!|\!|_{{\operatorname{op}}}/N. Then any solution Θ^\widehat{\Theta} to the semidefinite program (9) satisfies

Apart from the tolerance parameter δ\delta, the two main terms in the bound (14) have a natural interpretation. The first term (involving r\sqrt{r}) corresponds to estimation error, capturing the difficulty of estimating a rank rr matrix. The second is an approximation error, in which the projection onto the set M⊥(Ur,Vr)\mathcal{M}^{\perp}(U^{r},V^{r}) describes the gap between the true matrix Θ∗\Theta^{*} and the rank rr approximation.

Let us begin by illustrating the consequences of Theorem 1 when the true matrix Θ∗\Theta^{*} has exactly rank rr, in which case there is a very natural choice of the subspaces represented by UU and VV. In particular, we form UU from the rr non-zero left singular vectors of Θ∗\Theta^{*}, and VV from its rr non-zero right singular vectors. Note that this choice of (U,V)(U,V) ensures that ΠM⊥(U,V)(Θ∗)=0\Pi_{\mathcal{M}^{\perp}(U,V)}(\Theta^{*})=0. For technical reasons to be clarified, it suffices to set δ=0\delta=0 in the case of exact rank constraints, and we thus obtain the following result:

Suppose that Θ∗\Theta^{*} has rank rr, and X\mathfrak{X} satisfies RSC with respect to C(r;0)\mathcal{C}(r;0). Then as long as λN≥2∣ ⁣∣ ⁣∣X∗(ε⃗)∣ ⁣∣ ⁣∣op⁡/N\lambda_{N}\geq 2|\!|\!|\mathfrak{X}^{*}(\vec{\varepsilon})|\!|\!|_{{\operatorname{op}}}/N, any optimal solution Θ^\widehat{\Theta} 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 Θ∗\Theta^{*} is exactly low rank, then neither tolerance parameter δ\delta nor the approximation term are required.

As a more delicate example, suppose instead that Θ∗\Theta^{*} is nearly low-rank, an assumption that we can formalize by requiring that its singular value sequence {σi(Θ∗)}i=1min⁡{k,p}\{\sigma_{i}(\Theta^{*})\}_{i=1}^{\min\{k,p\}} decays quickly enough. In particular, for a parameter q∈q\in and a positive radius RqR_{q}, we define the set

Note that the error bound (17) reduces to the exact low rank case (15) when q=0q=0, and δ=0\delta=0. The quantity λN−qRq\lambda_{N}^{-q}R_{q} 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 λN\lambda_{N} is chosen to decay to zero as the sample size NN 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 Θ∗\Theta^{*}.

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 Θ∗\Theta^{*}, 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 Θ∗\Theta^{*} is not required), and specifying an appropriate choice of the regularization parameter λN\lambda_{N}. 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 1−c2exp⁡(−c3(k+p))1-c_{2}\exp(-c_{3}(k+p)).

Remarks: Corollary 3 takes a particularly simple form when Σ=Ip×p\Sigma=I_{p\times p}: then there exists a constant c1′c^{\prime}_{1} 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 Θ∗\Theta^{*} is exactly low rank—that is, q=0q=0, and Θ∗\Theta^{*} has rank r=R0r=R_{0}—this simplifies even further to

The scaling in this error bound is easily interpretable: naturally, the squared error is proportional to the noise variance ν2\nu^{2}, and the quantity r(k+p)r(k+p) counts the number of degrees of freedom of a k×pk\times p matrix with rank rr. Note that if we did not impose any constraints on Θ∗\Theta^{*}, then since a k×pk\times p matrix has a total of kpkp free parameters, we would expect at best to obtain rates of the order ∣ ⁣∣ ⁣∣Θ^−Θ∗∣ ⁣∣ ⁣∣F2=Ω(ν2 k pn)|\!|\!|\widehat{\Theta}-\Theta^{*}|\!|\!|_{{F}}^{2}=\Omega(\frac{\nu^{2}\,k\,p}{n}). Note that when Θ∗\Theta^{*} is low rank—in particular, when r≪min⁡{k,p}r\ll\min\{k,p\}—then the nuclear norm estimator achieves substantially faster rates. Finally, we note that as stated, the result requires that min⁡{k,p}\min\{k,p\} 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 (k,p)(k,p), the same statement (with different constants) holds with k+pk+p replaced by log⁡n\log n.

Next we turn to the case of estimating the system matrix Θ∗\Theta^{*} of an autoregressive (AR) model, as discussed in Example 2.

with probability greater than 1−c2exp⁡(−c3p)1-c_{2}\exp(-c_{3}p).

Remarks: Like Corollary 3, the result as stated requires that pp tend to infinity, but the same bounds hold with pp replaced by log⁡n\log n, yielding results suitable for classical (fixed dimension) scaling. Second, the factor (p/n)1−q/2(p/n)^{1-q/2}, like the analogous termThe term in Corollary 3 has a factor k+pk+p, since the matrix in that case could be non-square in general. in Corollary 3, shows that faster rates are obtained if Θ∗\Theta^{*} can be well-approximated by a low rank matrix, namely for choices of the parameter q∈q\in that are closer to zero. Indeed, in the limit q=0q=0, we again reduce to the case of an exact rank constraint r=R0r=R_{0}, and the corresponding squared error scales as rp/nrp/n. In contrast to the case of multivariate regression, the error bound (19) also depends on the upper bound ∣ ⁣∣ ⁣∣Θ∗∣ ⁣∣ ⁣∣op⁡=γ<1|\!|\!|\Theta^{*}|\!|\!|_{{\operatorname{op}}}=\gamma<1 on the operator norm of the system matrix Θ∗\Theta^{*}. Such dependence is to be expected since the quantity γ\gamma 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 1−c1exp⁡(−c2(k+p))1-c_{1}\exp(-c_{2}(k+p)).

The central challenge in proving this result is in proving an appropriate form of the RSC property. The following result on the random operator X\mathfrak{X} may be of independent interest here:

Under the stated conditions, the random operator X\mathfrak{X} satisfies

with probability at least 1−2exp⁡(−N/32)1-2\exp(-N/32).

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 X\mathfrak{X}; 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 yi=⟨ ⁣⟨Xi,  Θ∗⟩ ⁣⟩y_{i}=\langle\!\langle{X_{i}},\;{\Theta^{*}}\rangle\!\rangle for i=1,…,Ni=1,\ldots,N, and that we try to recover the unknown matrix Θ∗\Theta^{*} 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 Θ∗\Theta^{*} has rank rr, and that we are given N>402r(k+p)N>40^{2}r(k+p) noiseless samples. Then with probability at least 1−2exp⁡(−N/32)1-2\exp(-N/32), the SDP (22) recovers the matrix Θ∗\Theta^{*} 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 NN scales linearly in both rr and k+pk+p.

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 Θ^\widehat{\Theta} for the SDP (9), we have

Defining the error matrix Δ=Θ∗−Θ^\Delta=\Theta^{*}-\widehat{\Theta} 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 ∣ ⁣∣ ⁣∣Θ^+Δ∣ ⁣∣ ⁣∣1−∣ ⁣∣ ⁣∣Θ^∣ ⁣∣ ⁣∣1≤∣ ⁣∣ ⁣∣Δ∣ ⁣∣ ⁣∣1|\!|\!|\widehat{\Theta}+\Delta|\!|\!|_{{1}}-|\!|\!|\widehat{\Theta}|\!|\!|_{{1}}\leq|\!|\!|\Delta|\!|\!|_{{1}}. Substituting this inequality and the bound (24) into the inequality (23) yields

where the second inequality makes use of our choice λN≥2N∣ ⁣∣ ⁣∣X∗(ε⃗)∣ ⁣∣ ⁣∣op⁡\lambda_{N}\geq\frac{2}{N}|\!|\!|\mathfrak{X}^{*}(\vec{\varepsilon})|\!|\!|_{{\operatorname{op}}}.

It remains to lower bound the term on the left-hand side, while upper bounding the quantity ∣ ⁣∣ ⁣∣Δ∣ ⁣∣ ⁣∣1|\!|\!|\Delta|\!|\!|_{{1}} on the right-hand side. The following technical result allows us to do so. Recall our earlier definition (11) of the sets M\mathcal{M} and M⊥\mathcal{M}^{\perp} associated with a given subspace pair.

Let (U~,V~)(\widetilde{U},\widetilde{V}) represent a pair of rr-dimensional subspaces of left and right singular vectors of Θ∗\Theta^{*}. Then there exists a matrix decomposition Δ=Δ′+Δ′′\Delta=\Delta^{\prime}+\Delta^{\prime\prime} of the error Δ\Delta such that

The matrix Δ′\Delta^{\prime} satisfies the constraint rank⁡(Δ′)≤2r\operatorname{rank}(\Delta^{\prime})\leq 2r, and

If λN≥2∣ ⁣∣ ⁣∣X∗(ε⃗)∣ ⁣∣ ⁣∣2/N\lambda_{N}\geq 2|\!|\!|\mathfrak{X}^{*}(\vec{\varepsilon})|\!|\!|_{{2}}/N, then the nuclear norm of Δ′′\Delta^{\prime\prime} 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 ∣ ⁣∣ ⁣∣Δ′∣ ⁣∣ ⁣∣1≤2r∣ ⁣∣ ⁣∣Δ′∣ ⁣∣ ⁣∣F|\!|\!|\Delta^{\prime}|\!|\!|_{{1}}\leq\sqrt{2r}|\!|\!|\Delta^{\prime}|\!|\!|_{{F}}. Putting together the pieces, we find

2 Proof of Corollary 2

On the other hand, we also have Rq ≥  ∑i=1m∣σi(Θ∗)∣q  ≥  ∣K∣ τqR_{q}\,\geq\;\sum_{i=1}^{m}|\sigma_{i}(\Theta^{*})|^{q}\;\geq\;|K|\,\tau^{q}, which implies that ∣K∣≤τ−q Rq|K|\leq\tau^{-q}\,R_{q}. From the general error bound with r=∣K∣r=|K|, we obtain

Setting τ=λN/κ\tau=\lambda_{N}/\kappa 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 N=nkN=nk, 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 λn=λN k\lambda_{n}=\lambda_{N}\,k.

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 σmin⁡\sigma_{\operatorname{min}} 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 σmin⁡(XTX)2n≥σmin⁡(Σ)18\frac{\sigma_{\operatorname{min}}(X^{T}X)}{2n}\geq\frac{\sigma_{\operatorname{min}}(\Sigma)}{18} with probability at least 1−4exp⁡(−n)1-4\exp(-n) for all n≥pn\geq p, which establishes that the RSC property holds with κ(X)=120σmin⁡(Σ)\kappa(\mathfrak{X})=\frac{1}{20}\sigma_{\operatorname{min}}(\Sigma).

Next we need upper bound the quantity ∣ ⁣∣ ⁣∣X∗(ε⃗)∣ ⁣∣ ⁣∣2|\!|\!|\mathfrak{X}^{*}(\vec{\varepsilon})|\!|\!|_{{2}} for this model, so as to verify that the stated choice for λN\lambda_{N} 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 N=knN=kn, we see that Lemma 3 implies that the choice λN=10ν∣ ⁣∣ ⁣∣Σ∣ ⁣∣ ⁣∣op⁡k+pn\lambda_{N}=10\nu\sqrt{|\!|\!|\Sigma|\!|\!|_{{\operatorname{op}}}}\sqrt{\frac{k+p}{n}} satisfies the conditions of Corollary 2 with high probability. Lemma 2 shows that the RSC property holds with κ(X)=σmin⁡(Σ)/20\kappa(\mathfrak{X})=\sigma_{\operatorname{min}}(\Sigma)/20, again with high probability. Consequently, Corollary 2 implies that

with probability greater than 1−c2exp⁡(−c3(k+p))1-c_{2}\exp(-c_{3}(k+p)), 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 XTX/nX^{T}X/n is well-controlled in terms of the stationary covariance matrix: in particular, as long as n>c3pn>c_{3}p, we have

both with probability greater than 1−2c1exp⁡(−c2 p)1-2c_{1}\exp(-c_{2}\,p).

Thus, from the bound (28)(b), we see with the high probability, the RSC property holds with κ(X)=σmin⁡(Σ)/4\kappa(\mathfrak{X})=\sigma_{\operatorname{min}}(\Sigma)/4 as long as n>c3pn>c_{3}p.

As before, in order to verify the choice of λN\lambda_{N}, we need to control the quantity 1n∣ ⁣∣ ⁣∣XTW∣ ⁣∣ ⁣∣op⁡\frac{1}{n}|\!|\!|X^{T}W|\!|\!|_{{\operatorname{op}}}. The following inequality, proved in Appendix C.2, yields a suitable upper bound:

There exist constants ci>0c_{i}>0, independent of n,p,Σn,p,\Sigma etc. such that

From Lemma 5, we see that it suffices to choose λN=80 ∣ ⁣∣ ⁣∣Σ∣ ⁣∣ ⁣∣op⁡1−γpn\lambda_{N}=\frac{80\,|\!|\!|\Sigma|\!|\!|_{{\operatorname{op}}}}{1-\gamma}\sqrt{\frac{p}{n}}. With this choice, Corollary 2 of Theorem 1 yields that

with probability greater than 1−c2exp⁡(−c3p)1-c_{2}\exp(-c_{3}p), 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 ∣ ⁣∣ ⁣∣Δ∣ ⁣∣ ⁣∣F  ≤  δ|\!|\!|\Delta|\!|\!|_{{F}}\;\leq\;\delta, the result of Corollary 5 follows immediately. Therefore, we may take ∣ ⁣∣ ⁣∣Δ∣ ⁣∣ ⁣∣F2 ≥ δ|\!|\!|\Delta|\!|\!|_{{F}}^{2}\,\geq\,\delta. Now recall from Lemma 1 that the error Δ\Delta satisfies the bound (25). Combining these facts, we are guaranteed that Δ∈C(r;δ)\Delta\in\mathcal{C}(r;\delta), where the set C\mathcal{C} 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 Δ∈C\Delta\in\mathcal{C},

Following the arguments used in the proofs of Theorem 1 and Corollary 2, we find that

where τ>0\tau>0 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 N>200 Rq(1−q/2)k pN>200\,R_{q}^{(1-q/2)}k\,p, then we are guaranteed that 14−(4+32)δ≥18\frac{1}{4}-(4+\sqrt{32})\delta\geq\frac{1}{8}, which shows that the RSC property holds with κ(X)=1/8\kappa(\mathfrak{X})=1/8.

The next step is to control the quantity ∥X∗(ε⃗)∥2/N\|\mathfrak{X}^{*}(\vec{\varepsilon})\|_{2}/N, required for specifying a suitable choice of λN\lambda_{N}.

If ∥ε⃗∥2≤2νN\|\vec{\varepsilon}\|_{2}\leq 2\nu\sqrt{N}, then

6 Proof of Corollary 6

This corollary follows from a combination of Proposition 1 and Lemma 1. Let Θ^\widehat{\Theta} be an optimal solution to the SDP (22), and let Δ=Θ^−Θ∗\Delta=\widehat{\Theta}-\Theta^{*} be the error. Since Θ^\widehat{\Theta} is optimal and Θ∗\Theta^{*} is feasible for the SDP, we have ∣ ⁣∣ ⁣∣Θ^∣ ⁣∣ ⁣∣1=∣ ⁣∣ ⁣∣Θ∗+Δ∣ ⁣∣ ⁣∣1≤∣ ⁣∣ ⁣∣Θ∗∣ ⁣∣ ⁣∣1|\!|\!|\widehat{\Theta}|\!|\!|_{{1}}=|\!|\!|\Theta^{*}+\Delta|\!|\!|_{{1}}\leq|\!|\!|\Theta^{*}|\!|\!|_{{1}}. Using the decomposition Δ=Δ′+Δ′′\Delta=\Delta^{\prime}+\Delta^{\prime\prime} from Lemma 1 and applying triangle inequality, we have ∣ ⁣∣ ⁣∣Θ∗+Δ′+Δ′′∣ ⁣∣ ⁣∣1≥∣ ⁣∣ ⁣∣Θ∗+Δ′′∣ ⁣∣ ⁣∣1−∣ ⁣∣ ⁣∣Δ′∣ ⁣∣ ⁣∣1|\!|\!|\Theta^{*}+\Delta^{\prime}+\Delta^{\prime\prime}|\!|\!|_{{1}}\geq|\!|\!|\Theta^{*}+\Delta^{\prime\prime}|\!|\!|_{{1}}-|\!|\!|\Delta^{\prime}|\!|\!|_{{1}}. From the properties of the decomposition in Lemma 1 (see Appendix A), we find that

Combining the pieces yields that ∣ ⁣∣ ⁣∣Δ′′∣ ⁣∣ ⁣∣1≤∣ ⁣∣ ⁣∣Δ′∣ ⁣∣ ⁣∣1|\!|\!|\Delta^{\prime\prime}|\!|\!|_{{1}}\leq|\!|\!|\Delta^{\prime}|\!|\!|_{{1}}, and hence ∣ ⁣∣ ⁣∣Δ∣ ⁣∣ ⁣∣1≤2∣ ⁣∣ ⁣∣Δ′∣ ⁣∣ ⁣∣1|\!|\!|\Delta|\!|\!|_{{1}}\leq 2|\!|\!|\Delta^{\prime}|\!|\!|_{{1}}. By Lemma 1(a), the rank of Δ′\Delta^{\prime} is at most 2r2r, so that we obtain ∣ ⁣∣ ⁣∣Δ∣ ⁣∣ ⁣∣1≤22r∣ ⁣∣ ⁣∣Δ∣ ⁣∣ ⁣∣F≤4r∣ ⁣∣ ⁣∣Δ∣ ⁣∣ ⁣∣F|\!|\!|\Delta|\!|\!|_{{1}}\leq 2\sqrt{2r}|\!|\!|\Delta|\!|\!|_{{F}}\leq 4r|\!|\!|\Delta|\!|\!|_{{F}}.

Note that X(Δ)=0\mathfrak{X}(\Delta)=0, since both Θ^\widehat{\Theta} and Θ∗\Theta^{*} agree with the observations. Consequently, from Proposition 1, we have that

where the final inequality follows from the assumption that N>402r(k+p)N>40^{2}r(k+p). We have thus shown that Δ=0\Delta=0, which implies that Θ^=Θ∗\widehat{\Theta}=\Theta^{*} 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 MM-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 λN\lambda_{N} in the manner suggested by our theoretical results; in doing so, we assumed knowledge of quantities such as the noise variance ν2\nu^{2}. (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 N(0,I)N(0,I) distribution. Panel (a) plots the Frobenius error ∣ ⁣∣ ⁣∣Θ^−Θ∗∣ ⁣∣ ⁣∣F|\!|\!|\widehat{\Theta}-\Theta^{*}|\!|\!|_{{F}} on a logarithmic scale versus the sample size NN for three different matrix sizes, p∈{40,80,160}p\in\{40,80,160\}. Naturally, in each case, the error decays to zero as NN increases, but larger matrices require larger sample sizes, as reflected by the rightward shift of the curves as pp 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 N~: =N/(rp)\widetilde{N}:\,=N/(rp). 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 N/(rp)N/(rp) 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 N/(rp)N/(rp); as predicted by Corollary 4, the errors for different matrix sizes pp 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 XiX_{i} 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 N/rpN/rp, 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 Δ′′\Delta^{\prime\prime}, the nuclear norm satisfies the decomposition

We now turn to the proof of Lemma 1(b). Recall that the error Δ=Θ^−Θ∗\Delta=\widehat{\Theta}-\Theta^{*} 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 ∣ ⁣∣ ⁣∣1NX∗(ε⃗)∣ ⁣∣ ⁣∣op⁡≤λN/2|\!|\!|\frac{1}{N}\mathfrak{X}^{*}(\vec{\varepsilon})|\!|\!|_{{\operatorname{op}}}\leq\lambda_{N}/2 by assumption, we conclude that

Appendix B Proof of Lemma 3

For positive scalars aa and bb, define the (random) quantity

and note that our goal is to upper bound Ψ(1,1)\Psi(1,1). Note moreover that Ψ(a,b)=a b Ψ(1,1)\Psi(a,b)=a\,b\,\Psi(1,1), a relation which will be useful in the analysis.

Let A={u1,…,uA}\mathcal{A}=\{u^{1},\ldots,u^{A}\} and B={v1,…,vB}\mathcal{B}=\{v^{1},\ldots,v^{B}\} denote 1/41/4 coverings of Sk−1S^{k-1} and Sp−1S^{p-1}, respectively. We now claim that we have the upper bound

To establish this claim, we note that since the sets A\mathcal{A} and B\mathcal{B} are 1/41/4-covers, for any pair (u,v)∈Sp−1×Sp−1(u,v)\in S^{p-1}\times S^{p-1}, there exists a pair (ua,vb)∈A×B(u^{a},v^{b})\in\mathcal{A}\times\mathcal{B} such that u=ua+Δuu=u^{a}+\Delta u and v=vb+Δvv=v^{b}+\Delta v, with max⁡{∥Δu∥2,∥Δv∥2}≤1/4\max\{\|\Delta u\|_{2},\|\Delta v\|_{2}\}\leq 1/4. Consequently, we can write

By construction, we have the bound ∣⟨Xvb, WΔu⟩∣≤Ψ(1,1/4)=14Ψ(1,1)|\langle Xv^{b},\,W\Delta u\rangle|\leq\Psi(1,1/4)=\frac{1}{4}\Psi(1,1), and similarly ∣⟨XΔv, Wua⟩∣≤14Ψ(1,1)|\langle X\Delta v,\,Wu^{a}\rangle|\leq\frac{1}{4}\Psi(1,1) as well as ∣⟨XΔv, WΔu⟩∣≤116Ψ(1,1)|\langle X\Delta v,\,W\Delta u\rangle|\leq\frac{1}{16}\Psi(1,1). 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 1/41/4 covering of Sk−1S^{k-1} and Sp−1S^{p-1} with at most A≤8kA\leq 8^{k} and B≤8pB\leq 8^{p} elements respectively. Consequently, we have

Combining this tail bound with the upper bound (37), we have

Setting t2=20ν2∣ ⁣∣ ⁣∣Σ∣ ⁣∣ ⁣∣op⁡k+pnt^{2}=20\nu^{2}|\!|\!|\Sigma|\!|\!|_{{\operatorname{op}}}\frac{k+p}{n}, this probability vanishes as long as n>16(k+p)n>16(k+p).

Appendix C Technical details for Corollary 4

In this appendix, we collect the proofs of Lemmas 4 and 5.

Recalling that Sp−1S^{p-1} denotes the unit-norm Euclidean sphere in pp-dimensions, we first observe that ∣ ⁣∣ ⁣∣X∣ ⁣∣ ⁣∣op⁡=sup⁡u∈Sp−1∥Xu∥2|\!|\!|X|\!|\!|_{{\operatorname{op}}}=\sup_{u\in S^{p-1}}\|Xu\|_{2}. Our next step is to reduce the supremum to a maximization over a finite set, using a standard covering argument. Let A={u1,…,uA}\mathcal{A}=\{u^{1},\ldots,u^{A}\} denote a 1/21/2-cover of it. By definition, for any u∈Sp−1u\in S^{p-1}, there is some ua∈Au^{a}\in\mathcal{A} such that u=ua+Δuu=u^{a}+\Delta u, where ∥Δu∥2≤1/2\|\Delta u\|_{2}\leq 1/2. Consequently, for any u∈Sp−1u\in S^{p-1}, the triangle inequality implies that

and hence that ∣ ⁣∣ ⁣∣X∣ ⁣∣ ⁣∣op⁡≤max⁡ua∈A∥Xua∥2+12∣ ⁣∣ ⁣∣X∣ ⁣∣ ⁣∣op⁡|\!|\!|X|\!|\!|_{{\operatorname{op}}}\leq\max_{u^{a}\in\mathcal{A}}\|Xu^{a}\|_{2}+\frac{1}{2}|\!|\!|X|\!|\!|_{{\operatorname{op}}}. Re-arranging yields the useful inequality

where the last inequality follows from the union bound, and the fact that there exists a 1/21/2-covering of Sp−1S^{p-1} with at most 4p4^{p} elements.

Moreover, we have trace⁡(R)/n=uTΣu≤∣ ⁣∣ ⁣∣Σ∣ ⁣∣ ⁣∣op⁡\operatorname{trace}(R)/n=u^{T}\Sigma u\leq|\!|\!|\Sigma|\!|\!|_{{\operatorname{op}}}. Applying Lemma 8 with t=5pnt=5\sqrt{\frac{p}{n}}, we conclude that

with probability at least 1−c1 exp⁡(−c2 p)1-c_{1}\,\exp(-c_{2}\,p), which establishes the upper bound (28)(a).

since Ψ(Δv,Δv)≥0\Psi(\Delta v,\Delta v)\geq 0. 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 ∣ ⁣∣ ⁣∣1nXTX∣ ⁣∣ ⁣∣op⁡≤24∣ ⁣∣ ⁣∣Σ∣ ⁣∣ ⁣∣op⁡(1−γ)|\!|\!|\frac{1}{n}X^{T}X|\!|\!|_{{\operatorname{op}}}\leq\frac{24|\!|\!|\Sigma|\!|\!|_{{\operatorname{op}}}}{(1-\gamma)} with high probability. Hence, choosing ϵ=(1−γ)σmin⁡(Σ)200∣ ⁣∣ ⁣∣Σ∣ ⁣∣ ⁣∣op⁡\epsilon=\frac{(1-\gamma)\sigma_{\operatorname{min}}(\Sigma)}{200|\!|\!|\Sigma|\!|\!|_{{\operatorname{op}}}} ensures that 2ϵ∣ ⁣∣ ⁣∣1nXTX∣ ⁣∣ ⁣∣op⁡≤σmin⁡(Σ)/42\epsilon|\!|\!|\frac{1}{n}X^{T}X|\!|\!|_{{\operatorname{op}}}\leq\sigma_{\operatorname{min}}(\Sigma)/4.

Consequently, it suffices to lower bound the minimum over the covering set. We first establish a concentration result for the function Ψ(v,v)\Psi(v,v) that holds for any fixed v∈Sp−1v\in S^{p-1}. Note that we can write

Note that this bound holds for any fixed v∈Sp−1v\in S^{p-1}. Setting t∗=(1−γ) σmin⁡(Σ)16∣ ⁣∣ ⁣∣Σ∣ ⁣∣ ⁣∣op⁡t^{*}=\frac{(1-\gamma)\,\sigma_{\operatorname{min}}(\Sigma)}{16|\!|\!|\Sigma|\!|\!|_{{\operatorname{op}}}} and applying the union bound yields that

which vanishes as long as n>4log⁡(4/ϵ)(t∗)2pn>\frac{4\log(4/\epsilon)}{(t^{*})^{2}}p.

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 1/41/4 covering of Sp−1S^{p-1} with at most 8p8^{p} elements. Consequently, we have

For each i=1,…,ni=1,\ldots,n, let XiX_{i} and WiW_{i} denote the ithi^{th} row of XX and WW. Following some simple algebra, we have the decomposition ⟨Xv, Wu⟩n=T1−T2−T3\frac{\langle Xv,\,Wu\rangle}{n}=T_{1}-T_{2}-T_{3}, where

We may now bound each Tj,j=1,2,3T_{j},j=1,2,3 in turn; in doing so, we make repeated use of Lemma 8, which provides concentration bounds for a random variable of the form ∥Y∥22\|Y\|_{2}^{2}, where Y∼N(0,Q)Y\sim N(0,Q) for some matrix Q⪰0Q\succeq 0.

We begin with T2T_{2}, which the easiest to control since (up to scaling by ν\nu), it corresponds to the deviation away from the mean of χ2\chi^{2}-variable with nn degrees of freedom. Consequently, applying Lemma 8 with Q=IQ=I, we obtain

where δij\delta_{ij} is the Kronecker delta for the event {i=j}\{i=j\}. As before, by symmetry of SS, we have ∣ ⁣∣ ⁣∣S∣ ⁣∣ ⁣∣op⁡≤max⁡i=1,…,n∑j=1n∣Sij∣|\!|\!|S|\!|\!|_{{\operatorname{op}}}\leq\max_{i=1,\ldots,n}\sum_{j=1}^{n}|S_{ij}|, and hence

Morever, we have trace⁡(S)/n=ν2+vTΘ∗v\operatorname{trace}(S)/n=\nu^{2}+v^{T}\Theta^{*}v, 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 t>0t>0,

Setting t=10p/nt=10\sqrt{p/n} 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 t≥1t\geq 1, the lower bound

holds with probability at least 1−c1exp⁡(−c2N)1-c_{1}\exp(-c_{2}N). 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 (u,Θ)∈SN−1×R(t)(u,\Theta)\in S^{N-1}\times\mathcal{R}(t), consider the random variable Zu,Θ=⟨u, X(Θ)⟩Z_{u,\Theta}=\langle u,\,\mathfrak{X}(\Theta)\rangle, and note that it is Gaussian with zero mean. For any two pairs (u,Θ)(u,\Theta) and (u′,Θ)(u^{\prime},\Theta), some calculation yields

We now define a second Gaussian process {Yu,Θ ∣ (u,Θ)∈SN−1×R(t)}\{Y_{u,\Theta}\,\mid\,(u,\Theta)\in S^{N-1}\times\mathcal{R}(t)\} via

It can be shown that for all pairs (u,Θ),(u′,Θ′)∈SN−1×R(t)(u,\Theta),(u^{\prime},\Theta^{\prime})\in S^{N-1}\times\mathcal{R}(t), we have

Moreover, equality holds whenever Θ=Θ′\Theta=\Theta^{\prime}. 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 f(X): =inf⁡Θ∈R(t)∥X(Θ)∥2/Nf(\mathfrak{X}):\,=\inf_{\Theta\in\mathcal{R}(t)}\|\mathfrak{X}(\Theta)\|_{2}/\sqrt{N} is Lipschitz with constant 1/N1/\sqrt{N}, 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 Y∼N(0,Q)Y\sim N(0,Q), for all t>2/nt>2/\sqrt{n}, we have

Let Q\sqrt{Q} be the symmetrix matrix square root, and consider the function f(x)=∥Qx∥2/nf(x)=\|\sqrt{Q}x\|_{2}/\sqrt{n}. Since it is Lipschitz with constant ∣ ⁣∣ ⁣∣Q∣ ⁣∣ ⁣∣op⁡/n|\!|\!|\sqrt{Q}|\!|\!|_{{\operatorname{op}}}/\sqrt{n}, Lemma 7 implies that

By integrating this tail bound, we find that the variable Z=∥QX∥2/nZ=\|\sqrt{Q}X\|_{2}/\sqrt{n} satisfies the bound var⁡(Z)≤4∣ ⁣∣ ⁣∣Q∣ ⁣∣ ⁣∣op⁡/n\operatorname{var}(Z)\leq 4|\!|\!|Q|\!|\!|_{{\operatorname{op}}}/n, and hence conclude that

Combining this bound with the tail bound (54), we conclude that

Setting δ=(t−2/n) ∣ ⁣∣ ⁣∣Q∣ ⁣∣ ⁣∣op⁡\delta=(t-2/\sqrt{n})\,\sqrt{|\!|\!|Q|\!|\!|_{{\operatorname{op}}}} in the bound (56) yields that

Similarly, setting δ=∣ ⁣∣ ⁣∣Q∣ ⁣∣ ⁣∣op⁡\delta=\sqrt{|\!|\!|Q|\!|\!|_{{\operatorname{op}}}} in the tail bound (56) yields that with probability greater than 1−2exp⁡(−n/2)1-2\exp(-n/2), we have

References