Optimal selection of reduced rank estimators of high-dimensional matrices

Introduction

where $E$ is a random $m\times n$ matrix, with independent entries with mean zero and variance $\sigma^{2}$.

Standard least squares estimation in (1), under no constraints, is equivalent to regressing each response on the predictors separately. It completely ignores the multivariate nature of the possibly correlated responses, see, for instance, Izenman (2008) for a discussion of this phenomenon. Estimators restricted to have rank equal to a fixed number $k\leq n\wedge p$ were introduced to remedy this drawback. The history of such estimators dates back to the 1950’s, and was initiated by Anderson (1951). Izenman (1975) introduced the term reduced-rank regression for this class of models and provided further study of the estimates. A number of important works followed, including Robinson (1973, 1974) and Rao (1978). The monograph on reduced rank regression by Reinsel and Velu (1998) has an excellent, comprehensive account of more recent developments and extensions of the model. All theoretical results to date for estimators of $A$ constrained to have rank equal to a given value $k$ are of asymptotic nature and are obtained for fixed $p$, independent of the number of observations $m$. Most of them are obtained in a likelihood framework, for Gaussian errors $E_{ij}$. Anderson (1999) relaxed this assumption and derived the asymptotic distribution of the estimate, when $p$ is fixed, the errors have two finite moments, and the rank of $A$ is known. Anderson (2002) continued this work by constructing asymptotic tests for rank selection, valid only for small and fixed values of $p$.

The aim of our work is to develop a non-asymptotic class of methods that yield reduced rank estimators of $A$ that are easy to compute, have rank determined adaptively from the data, and are valid for any values of $m,n$ and $p$, especially when the number of predictors $p$ is large. The resulting estimators can then be used to construct a possibly much smaller number of new transformed predictors or can be used to construct the most important canonical variables based on the original $X$ and $Y$. We refer to Chapter 6 in Izenman (2008) for a historical account of the latter.

We propose to estimate $A$ by minimizing the sum of squares $\|Y-XB\|_{F}^{2}=\sum_{i}\sum_{j}\{Y_{ij}-(XB)_{ij}\}^{2}$ plus a penalty $\mu r(B)$, proportional to the rank $r(B)$, over all matrices $B$. It is immediate to see, using Pythagoras’ theorem, that this is equivalent with computing $\min_{B}\left\{\|PY-XB\|_{F}^{2}+\mu r(B)\right\}$ or $\min_{k}\left\{\min_{B:\ r(B)=k}\|PY-XB\|_{F}^{2}+\mu k\right\}$, with $P$ being the projection matrix onto the column space of $X$. In Section 2.1 we show that the minimizer $\widehat{k}$ of the above expression is the number of singular values $d_{k}(PY)$ of $PY$ that exceed ${\mu}^{1/2}$. This observation reveals the prominent role of the tuning parameter $\mu$ in constructing $\widehat{k}$. The final estimator $\widehat{A}$ of the target matrix $A$ is the minimizer of $\|PY-XB\|_{F}^{2}$ over matrices $B$ of rank $\widehat{k}$, and can be computed efficiently even for large $p$, using the procedure that we describe in detail in Section 2.1 below.

The theoretical analysis of our proposed estimator $\widehat{A}$ is presented in Sections 2.2 – 2.4. The rank of $A$ may not be the most appropriate measure of sparsity in multivariate regression models. For instance, suppose that the rank of $A$ is 100, but only three of its singular values are large and the remaining 97 are nearly zero. This is an extreme example, and in general one needs an objective method for declaring singular values as “large” or “small”. We introduce in Section 2.1 a slightly different notion of sparsity, that of effective rank. The effective rank counts the number of singular values of the signal $XA$ that are above a certain noise level. The relevant notion of noise level turns out to be the largest singular value of $PE$. This is central to our results, and influences the choice of the tuning sequence $\mu$. In Appendix C we prove that the expected value of the largest singular value of $PE$ is bounded by $(q+n)^{1/2}$, where $q\leq m\wedge p$ is the rank of $X$. The effective noise level is at most $(m+n)^{1/2}$, for instance in the model $Y=A+E$, but it can be substantially lower, of order $(q+n)^{1/2}$, in model (1).

In Section 2.2 we give tight conditions under which $\widehat{k}$, the rank of our proposed estimator $\widehat{A}$, coincides with the effective rank. As an immediate corollary we show when $\widehat{k}$ equals the rank of $A$. We give finite sample performance bounds for $\|X\widehat{A}-XA\|_{F}^{2}$ in Section 2.3. These results show that $\widehat{A}$ mimics the behavior of reduced rank estimates based on the ideal effective rank, had this been known prior to estimation. If $X$ has a restricted isometrity property, our estimate is minimax adaptive. In the asymptotic setting, for $n+(m\wedge p)\geq n+q\rightarrow\infty$, all our results hold with probability close to one, for tuning parameter chosen proportionally to the square of the noise level.

We often particularize our main findings to the setting of Gaussian $N(0,\sigma^{2})$ errors $E_{ij}$ in order to obtain sharp, explicit numerical constants for the penalty term. To avoid technicalities, we assume that $\sigma^{2}$ is known in most cases, and we treat the case of unknown $\sigma^{2}$ in Section 2.4.

We contrast our estimator with the penalized least squares estimator $\widetilde{A}$ corresponding to a penalty term $\tau\|B\|_{1}$ proportional to the nuclear norm $\|B\|_{1}=\sum_{j}d_{j}(B)$, the sum of the singular values of $B$. This estimator has been studied by, among others, Yuan et al. (2007) and Lu et al (2010), for model (1). Nuclear norm penalized estimators in general models y=\mbox{\mathcal{X}}(A)+\varepsilon involving linear maps $\mathcal{X}$ have been studied by Candès and Plan (2010) and Negahban and Wainwright (2009). A special case of this model is the challenging matrix completion problem, first investigated theoretically, in the noiseless case, by Candès and Tao (2010). Rohde and Tsybakov (2010) studied a larger class of penalized estimators, that includes the nuclear norm estimator, in the general model y=\mbox{\mathcal{X}}(A)+\varepsilon.

In Section 3 we give bounds on $\|X\widetilde{A}-XA\|_{F}^{2}$ that are similar in spirit to those from Section 2. While the error bounds of the two estimators are comparable, albeit with cleaner results and milder conditions for our proposed estimator, there is one aspect in which the estimates differ in important ways. The nuclear norm penalized estimator is far less parsimonious than the estimate obtained via our rank selection criterion. In Section 3, we offer a correction of the former estimate that yields a correct rank estimate.

Section 4 complements our theoretical results by an extensive simulation study that supports our theoretical findings and suggests strongly that the proposed estimator behaves very well in practice, in most situations is preferable to the nuclear norm penalized estimator and it is always much faster to compute.

Technical results and some intermediate proofs are presented in Appendices A – D.

The Rank Selection Criterion

We propose to estimate $A$ by the penalized least squares estimator

We denote its rank by $\widehat{k}$. The minimization is taken over all $p\times n$ matrices $B$. Here and in what follows $r(B)$ is the rank of $B$ and $\|C\|_{F}=\left(\sum_{i}\sum_{j}C_{ij}^{2}\right)^{1/2}$ denotes the Frobenius norm for any generic matrix $C$. The choice of the tuning parameter $\mu>0$ is discussed in Section 2.2. Since

one needs to compute the restricted rank estimators $\widehat{B}_{k}$ that minimize $\|Y-XB\|_{F}^{2}$ over all matrices $B$ of rank $k$. The following computationally efficient procedure for calculating each $\widehat{B}_{k}$ has been suggested by Reinsel and Velu (1998). Let $M=X^{\prime}X$ be the Gram matrix, $M^{-}$ be its Moore-Penrose inverse and let $P=XM^{-}X^{\prime}$ be the projection matrix onto the column space of $X$.

Compute the eigenvectors $V=[v_{1},v_{2},\cdots,v_{n}]$, corresponding to the ordered eigenvalues arranged from largest to smallest, of the symmetric matrix $Y^{\prime}PY$.

Compute the least squares estimator $\widehat{B}=M^{-}X^{\prime}Y$. Construct $W=\widehat{B}V$ and $G=V^{\prime}$. Form $W_{k}=W[\,,1:k]$ and $G_{k}=G[1:k,\,]$.

Compute the final estimator $\widehat{B}_{k}=W_{k}G_{k}$.

In step 2 above, $W_{k}$ denotes the matrix obtained from $W$ by retaining all its rows and only its first $k$ columns, and $G_{k}$ is obtained from $G$ by retaining its first $k$ rows and all its columns.

Our first result, Proposition 1 below, characterizes the minimizer $\widehat{k}=r(\widehat{A})$ of (3) as the number of eigenvalues of the square matrix $Y^{\prime}PY$ that exceed $\mu$ or, equivalently, as the number of singular values of the matrix $PY$ that exceed ${\mu}^{1/2}$. The final estimator of $A$ is then $\widehat{A}=\widehat{B}_{\widehat{k}}$.

Lemma 14 in Appendix B shows that the fitted matrix $X\widehat{A}$ is equal to $\sum_{j\leq\widehat{k}}d_{j}u_{j}v_{j}^{\prime}$ based on the singular value decomposition $UDV=\sum_{j}d_{j}u_{j}v_{j}^{\prime}$ of the projection $PY$.

Let $\lambda_{1}(Y^{\prime}PY)\geq\lambda_{2}(Y^{\prime}PY)\geq\cdots$ be the ordered eigenvalues of $Y^{\prime}PY$. We have $\widehat{A}=\widehat{B}_{\widehat{k}}$ with

For $\widehat{B}_{k}$ given above, and by the Pythagorean theorem, we have

and we observe that $X\widehat{B}=PY$. By Lemma 14 in Appendix B, we have

where $d_{j}(C)$ denotes the $j$-th largest singular value of a matrix $C$. Then, the penalized least squares criterion reduces to

and we find that $\min_{B}\left\{\|Y-XB\|_{F}^{2}+\mu r(B)\right\}$ equals

It is easy to see that $\sum_{j>k}\left\{\lambda_{j}(Y^{\prime}PY)-\mu\right\}$ is minimized by taking $k$ as the largest index $j$ for which $\lambda_{j}(Y^{\prime}PY)-\mu\geq 0$, since then the sum only consists of negative terms. This concludes our proof. ∎

Remark. The two matrices $W_{\widehat{k}}$ and $G_{\widehat{k}}$, that yield the final solution $\widehat{A}=W_{\widehat{k}}G_{\widehat{k}}$, have the following properties: (i) $G_{\widehat{k}}G_{\widehat{k}}^{\prime}$ is the identity matrix; and (ii) $W_{\widehat{k}}^{\prime}MW_{\widehat{k}}$ is a diagonal matrix. Moreover, the decomposition of $\widehat{A}$ as a product of two matrices with properties (i) and (ii) is unique, see, for instance, Theorem 2.2 in Reinsel and Vélu (1998). As an immediate consequence, one can construct new orthogonal predictors as the columns of $Z=XW_{\widehat{k}}$. If $\widehat{k}$ is much smaller than $p$, this can result in a significant dimension reduction of the predictors’ space.

2 Consistent effective rank estimation

In this section we study the properties of $\widehat{k}=r(\widehat{A})$. We will state simple conditions that guarantee that $\widehat{k}$ equals $r=r(A)$ with high probability. First, we describe in Theorem 2 what $\widehat{k}$ estimates and what quantities need to be controlled for consistent estimation. It turns out that $\widehat{k}$ estimates the number of the singular values of the signal $XA$ above the threshold $\mu^{1/2}$, for any value of the tuning parameter $\mu$. The quality of estimation is controlled by the probability that this threshold level exceeds the largest singular value $d_{1}(PE)$ of the projected noise matrix $PE$. We denote the $j$th singular value of a generic matrix $C$ by $d_{j}(C)$ and we use the convention that the singular values are indexed in decreasing order.

Suppose that there exists an index $s\leq r$ such that

Using the characterization of $\widehat{k}$ given in Proposition 1 we have

Invoke the conditions on $d_{s+1}({XA})$ and $d_{s}({XA})$ to complete the proof. ∎

Theorem 2 indicates that we can consistently estimate the index $s$ provided we use a large enough value for our tuning parameter $\mu$ to guarantee that the probability of the event $\left\{d_{1}(PE)\leq\delta\mu^{1/2}\right\}$ approaches one. We call $s$ the effective rank of $A$ relative to $\mu$, and denote it by $r_{e}=r_{e}(\mu)$.

Let $q=r(X)$ and assume that $E_{ij}$ are independent $N(0,\sigma^{2})$ random variables. Then

In view of this result, we take $\mu^{1/2}>\sigma(n^{1/2}+q^{1/2})$ as our measure of the noise level. The following corollary summarizes the discussion above and lists the main results of this section: the proposed estimator based on the rank selection criterion (RSC) recovers consistently the effective rank $r_{e}$ and, in particular, the rank of $A$.

Assume that $E$ has independent $N(0,\sigma^{2})$ entries. For any $\theta>0$, set

with $\delta$ as in Theorem 2. Then we have, for any $\theta>0$,

In particular, if $d_{r}(XA)>2\mu^{1/2}$ and $\mu^{1/2}=(1+\theta)\sigma(\sqrt{n}+\sqrt{q})$, then

Remark. Corollary 4 holds when $q+n\rightarrow\infty$. If $q+n$ stays bounded, but $m\rightarrow\infty$, the consistency results continue to hold when $q$ is replaced by $q\ln(m)$ in the expression of the tuning parameter $\mu$ given above. Lemma 3 justifies this choice. The same remark applies to all theoretical results in this paper.

for all $x>0$. The conclusion of Corollary 4 then holds for $\mu=C_{0}\Gamma_{E}(n+q)$ with $C_{0}$ large enough. Moreover, all oracle inequalities presented in the next sections remain valid for this choice of the tuning parameter, if $E$ has independent subGaussian entries.

3 Errors bounds for the RSC estimator

In this section we study the performance of $\widehat{A}$ by obtaining bounds for $\|X\widehat{A}-XA\|_{F}^{2}$. First we derive a bound for the fit $\|X\widehat{B}_{k}-XA\|_{F}^{2}$, based on the restricted rank estimator $\widehat{B}_{k}$, for each value of $k$.

Set $c(\theta)=1+2/\theta$. For any $\theta>0$, we have

for all $p\times n$ matrices $B$ of rank $k$. Working out the squares we obtain

for generic $m\times n$ matrices $C$ and $D$. The inner product $<C,D>_{F}$, operator norm $\|C\|_{2}=d_{1}(C)$ and nuclear norm $\|D\|_{1}=\sum_{j}d_{j}(D)$ are related via the inequality $<C,D>_{F}\leq\|C\|_{2}\|D\|_{1}$. As a consequence we find

Using the inequality $2xy\leq x^{2}/a+ay^{2}$ with $a>0$ twice, we obtain that $\|X\widehat{B}_{k}-XA\|_{F}^{2}$ is bounded above by

Hence we obtain, for any $a,b>0$, the inequality

Lemma 14 in the Appendix B states that the minimum of $\|XA-XB\|_{F}^{2}$ over all matrices $B$ of rank $k$ is achieved for the GSVD of $A$ and the minimum equals $\sum_{j>k}d_{j}^{2}(XA)$. The claim follows after choosing $a=(2+\theta)/2$ and $b=\theta/2$. ∎

Assume that $E$ has independent $N(0,\sigma^{2})$ entries. Set $c(\theta)=1+2/\theta$. Then, for any $\theta,\xi>0$, the inequality

holds with probability $1-\exp(-\xi^{2}(n+q)/2)$. In addition,

The symbol $\lesssim$ means that the inequality holds up to multiplicative numerical constants.

Set $t=(1+\xi)^{2}\sigma^{2}(\sqrt{n}+\sqrt{q})^{2}$ for some $\xi>0$. From Lemma 3, it follows that

Theorem 5 bounds the error $\|X\widehat{B}_{k}-XA\|_{F}^{2}$ by an approximation error, $\sum_{j>k}d_{j}^{2}(XA)$, and a stochastic term, $kd_{1}^{2}(PE)$, with probability one. The approximation error is decreasing in $k$ and vanishes for $k>r(XA)$.

We observe that $k(n+q)$ is essentially the number of free parameters of the restricted rank problem. Indeed, our parameter space consists of all $p\times n$ matrices $B$ of rank $k$ and each $XB$ matrix has $k(n+q-k)$ free parameters. Hence we can interpret the bound in Corollary 6 above as the squared bias plus the dimension of the parameter space.

Remark(ii), following Corollary 8 below, shows that $k(n+q)$ is also the minimax lower bound for $\|X\widehat{B}_{k}-XA\|_{F}^{2}$, if the smallest eigenvalue of $X^{\prime}X$ is larger than a strictly positive constant. This means that $X\widehat{B}_{k}$ is a minimax estimator under this assumption.

We now turn to the penalized estimator $\widehat{A}$ and show that it achieves the best (squared) bias-variance trade-off among all rank restricted estimators $\widehat{B}_{k}$ for the appropriate choice of the tuning parameter $\mu$ in the penalty $\text{pen}(B)=\mu r(B)$.

We have, for any $\theta>0$, on the event $(1+\theta)d_{1}^{2}(PE)\leq\mu$,

for any $p\times n$ matrix $B$. In particular, we have, for $\mu\geq(1+\theta)d_{1}^{2}(PE)$

for all $p\times n$ matrices $B$. Working out the squares we obtain

Consequently, using the inequality $2xy\leq x^{2}/a+ay^{2}$ twice, we obtain, for any $a>0$ and $b>0$,

Hence, if $(a+b)d_{1}^{2}(PE)-\mu\leq 0$, we obtain

for any $a>1$ and $b>0$. Lemma 14 in Appendix B evaluates the minimum of $\|XA-XB\|_{F}^{2}$ over all matrices $B$ of rank $k$ and shows that it equals $\sum_{j>k}d_{j}^{2}(XA)$. We conclude our proof by choosing $a=1+\theta/2$ and $b=\theta/2$. ∎

Remark. The first two parts of the theorem show that $\widehat{A}$ achieves the best (squared) bias-variance trade-off among all reduced rank estimators $\widehat{B}_{k}$ if $\mu>d_{1}^{2}(PE)$. Moreover, the index $k$ which minimizes $\sum_{j>k}\{d_{j}^{2}(XA)+\mu k\}$ essentially coincides with the effective rank $r_{e}=r_{e}(\mu)$ defined in the previous section. Therefore, the fit of the selected estimator $X\widehat{A}$ is comparable with that of the estimator $X\widehat{B}_{k}$ with rank $k=r_{e}$. Since the ideal $r_{e}$ depends on the unknown matrix $A$, this ideal estimator cannot be computed. Although our estimator $\widehat{A}$ is constructed independently of $r_{e}$, it mimics the behavior of the ideal estimator $\widehat{B}_{r_{e}}$ and we say that the bound on $\|X\widehat{A}-XA\|_{F}^{2}$ adapts to $r_{e}\leq r$.

The last part of our result is a particular case of the second part, but it is perhaps easier to interpret. Taking the index $k$ equal to the rank $r$, the bias term disappears and the bound reduces to $rd_{1}^{2}(PE)$ up to constants. This shows clearly the important role played by $r$ in the estimation accuracy: the smaller the rank of $A$, the smaller the estimation error.

For Gaussian errors, we have the following precise bounds.

Assume that $E$ has independent $N(0,\sigma^{2})$ entries. Set

with $\theta,\xi>0$ arbitrary. Let $c(\theta)=1+2/\theta$. Then, we have

Apply Lemma 16 in Appendix D to deduce that

Remarks. (i) We note that for $n+q$ large,

(iii) The same type of upper bound as the one of Corollary 8 can be proved if the entries of $E$ are subGaussian: take $\text{pen}(B)=C(n+q)r(B)$ for some $C$ large enough, and invoke Proposition 15 in Appendix C.

(iv) Although the error bounds of $\|X\widehat{A}-XA\|_{F}$ are guaranteed for all $X$ and $A$, the analysis of the estimation performance of $\widehat{A}$ depends on $X$. If $\lambda_{p}(M)\geq\rho>0$, for some constant $\rho$, then, provided $\mu>(1+\theta)d_{1}^{2}(PE)$ with $\theta>0$ arbitrary,

(v) Our results are slightly more general than stated. In fact, our analysis does not require that the postulated multivariate linear model $Y=XA+E$ holds exactly. We denote the expected value of $Y$ by $\Theta$ and write $Y=\Theta+E$. We denote the projection of $\Theta$ onto the column space of $X$ by $XA$, that is, $P\Theta=XA$. Because minimizing $\|Y-XB\|_{F}^{2}+\mu r(B)$ is equivalent with minimizing $\|PY-XB\|_{F}^{2}+\mu r(B)$ by Pythagoras’ theorem, our least squares procedure estimates $XA$, the mean of $PY$. The statements of Theorems 2 and 7 remain unchanged, except that $XA$ is the mean of the projection $PY$ of $Y$, not the mean of $Y$ itself.

4 A data adaptive penalty term

In this section we construct a data adaptive penalty term that employs the unbiased estimator

of $\sigma^{2}$. Set, for any $\theta>0$, $\xi>0$ and $0<\delta<1$,

Notice that the estimator $S^{2}$ requires that $n(m-q)$ be large, which holds whenever $m>>q$ or $m-q\geq 1$ and $n$ is large. The challenging case $m=q<<p$ is left for future research.

Assume that $E$ is an $m\times n$ matrix with independent $N(0,\sigma^{2})$ entries. Using the penalty given above we have, for $c(\theta)=1+2/\theta$,

It remains to bound the expected value of

We split the expectation into two parts: $S^{2}\geq(1-\delta)\sigma^{2}$ and its complement. We observe first that

using Lemma 16 for the last inequality. Next, we observe that

using Lemmas 16 and 17 in Appendix D for the last inequality. This proves the result. ∎

Remark. We see that for large values of $n+q$ and $n(m-q)$,

as the additional terms in the theorem above decrease exponentially fast in $n+q$ and $n(m-q)$. This bound is similar to the one in Corollary 8, obtained for the RSC estimator corresponding to the penalty term that employs the theoretical value of $\sigma^{2}$.

Comparison with nuclear norm penalized estimators

In this section we compare our RSC estimator $\widehat{A}$ with the alternative estimator $\widetilde{A}$ that minimizes

On the event $d_{1}(X^{\prime}E)\leq\tau$, we have, for any $B$,

for all $m\times n$ matrices $B$. Working out the squares we obtain

on the event $d_{1}(X^{\prime}E)\leq{\tau}$, we obtain the claim using the triangle inequality. ∎

We see that $\widetilde{A}$ balances the bias term $\|XA-XB\|_{F}^{2}$ with the penalty term $\tau\|B\|_{1}$, provided $\tau>d_{1}(X^{\prime}E)$. Since $X^{\prime}E=X^{\prime}PE+X^{\prime}(I-P)E=X^{\prime}PE$, we have $d_{1}(X^{\prime}E)\leq d_{1}(X)d_{1}(PE)$. We immediately obtain the following corollary using the results for $d_{1}(PE)$ of Lemma 3.

Assume that $E$ has independent $N(0,\sigma^{2})$ entries. For

The same result, up to constants, can be obtained if the errors $E_{ij}$ are subGaussian, if we replace $\sigma$ in the choice of $\tau$ above by a suitably large constant $C$. The proof of this generalization uses Proposition 15 in Appendix C in lieu of Lemma 3. The same remark applies for all the results in this section.

The next result obtains an oracle inequality for $\widetilde{A}$ that resembles the oracle inequality for the RSC estimator $\widehat{A}$ in Theorem 7. We stress the fact that Theorem 12 below requires that $\lambda_{p}(X^{\prime}X)>0$; this was not required for the derivation of the oracle bound on $\|X\widehat{A}-XA\|_{F}^{2}$ in Theorem 7, which holds for all $X$. We denote the condition number of $M=X^{\prime}X$ by $c_{0}(M)=\lambda_{1}(M)/\lambda_{p}(M)$.

Assume that $E$ has independent $N(0,\sigma^{2})$ entries. For

Both inequalities hold with probability at least $1-\exp\left(-\theta^{2}(n+q)/2\right)$. The symbol $\lesssim$ means that the inequality holds up to multiplicative numerical constants (depending on $\theta$).

To keep the paper self contained, we give a simple proof of this result in Appendix A. Similar results for the NNP estimator of $A$ in the general model y=\mbox{\mathcal{X}}(A)+\varepsilon, where $\mathcal{X}$ is a random linear map, have been obtained by Negahban and Wainwright (2009) and Candès and Plan (2010), each under different sets of assumptions on $\mathcal{X}$. We refer to Rohde and Tsybakov (2010) for more general results on Schatten norm penalized estimators of $A$ in the model y=\mbox{\mathcal{X}}(A)+\varepsilon, and a very thorough discussion on the assumptions on $\mathcal{X}$ under which these results hold.

Theorem 10 shows that the error bounds of the nuclear norm penalized (NNP) estimator $\widetilde{A}$ and the RSC estimator $\widehat{A}$ are comparable, although it is worth pointing out that our bounds for $\widehat{A}$ are much cleaner and obtained under fewer restrictions on the design matrix. However, there is one aspect in which the two estimators differ radically: correct rank recovery. We showed in Section 2.2 that the RSC estimator corresponding to the effective value of the tuning sequence $\mu_{e}$ has the correct rank and achieves the optimal bias-variance trade-off. This is also visible in the left panel of Figure 1 which shows the plots of the MSE and rank of the RSC estimate as we varied the tuning parameter of the procedure over a large grid. The numbers on the vertical axis correspond to the range of values of the rank of the estimator considered in this experiment, 1 to 25. The rank of $A$ is 10. We notice that for the same range of values of the tuning parameter, RSC has both the smallest MSE value and the correct rank. We repeated this experiment for the NNP estimator. The right panel shows that the smallest MSE and the correct rank are not obtained for the same value of the tuning parameter. Therefore, a different strategy for correct rank estimation via NNP is in order.

Rather than taking the rank of $\widetilde{A}$ as the estimator of the rank of $A$, we consider instead, for $M=X^{\prime}X$,

Let $r=r(A)$ and assume that $d_{r}(MA)>4\tau$. Then

If $E$ has independent $N(0,\sigma^{2})$ entries and $\tau=(1+\theta)\sigma d_{1}(X)(\sqrt{n}+\sqrt{q})$, the above probability is bounded by $\exp\left(-\theta^{2}(n+q)/2\right)$.

Empirical Studies

We performed an extensive simulation study to evaluate the performance of the proposed method, RSC, and compare it with the NNP method. The RSC estimator $\widehat{A}$ was computed via the procedure outlined in Section 2.1. This method is computationally efficient in large dimensions. Its computational complexity is the same as that of PCA. Our choice for the tuning parameter $\mu$ was based on our theoretical findings in Section 2. In particular, Corollary 4 and Corollary 8 guarantee good rank selection and prediction performance of RSC provided that $\mu$ is just a little bit larger than $\sigma^{2}(\sqrt{n}+\sqrt{q})^{2}$. Under the assumption that $q<m$, we can estimate $\sigma^{2}$ by $S^{2}$; see Section 2.4 for details. In our simulations we used the adaptive tuning parameter $\mu_{adap}=2S^{2}(n+q)$. We experimented with other constants and found that the constant equal to 2 was optimal; constants slightly larger than 2 gave very similar results.

We compared the RSC estimator with the NNP estimator $\widetilde{A}$ and with the proposed trimmed or calibrated NNP estimator, denoted in what follows by NNP(c). The NNP estimator is the minimizer of the convex criterion $\|Y-XB\|_{F}^{2}+2\tau\|B\|_{1}.$ By the equivalent SDP characterization of the NNP-norm given in Fazel (2002), the original minimization problem is equivalent to the convex optimization problem

In our simulation study we compared the rank selection and the estimation performances of the RSC estimator RSC$|_{adap}$, corresponding to $\mu_{adap}$, with the optimally tuned RSC estimator, and the optimally tuned NNP and NNP(c) estimators. The last three estimators are called RSC$|_{val}$, NNP$|_{val}$ and NNP${}^{(c)}|_{val}$. They correspond to those tuning parameters $\mu_{val}$, $\tau_{val}$ and $\tau_{val}$, respectively, that gave the best prediction accuracy, when prediction was evaluated on a very large independent validation set. This comparison helps us understand the true potential of each method in an ideal situation, and allows us to draw a stable performance comparison between the proposed adaptive RSC estimator and the best possible versions of RSC and NNP.

We considered the following large sample-size set up and large dimensionality set up.

We constructed the matrix of dependent variables $X=[x_{1},x_{2},\cdots,x_{m}]^{\prime}$ by generating its rows $x_{i}$ as i.i.d. realizations from a multivariate normal distribution $\mbox{MVN}(\boldsymbol{0},\Sigma)$, with $\Sigma_{jk}=\rho^{|j-k|}$, $\rho>0$, $1\leq j,k\leq p$. The coefficient matrix $A=bB_{0}B_{1}$, with $b>0$, $B_{0}$ is a $p\times r$ matrix and $B_{1}$ is a $r\times n$ matrix. All entries in $B_{0}$ and $B_{1}$ are i.i.d. $N(0,1)$. Each row in $Y=[y_{1},\cdots,y_{m}]^{\prime}$ is then generated as $y_{i}=x_{i}^{\prime}A+E_{i}$, $1\leq i\leq m$, with $E_{i}$ denoting the $i$-th row of the noise matrix $E$ which has $m\times n$ independent $N(0,1)$ entries $E_{ij}$.

Experiment 2 (p>m(>q)𝑝annotated𝑚absent𝑞p>m(>q))

Each simulated model is characterized by the following control parameters: $m$ (sample size), $p$ (number of independent variables), $n$ (number of response variables), $r$ (rank of $A$), $\rho$ (design correlation), $q$ (rank of the design), and $b$ (signal strength). In Experiment 1, we set $m=100,\,p=25,\,n=25,\,r=10$, and varied the correlation coefficient $\rho=0.1,0.5,0.9$ and signal strength $b=0.1,0.2,0.3,0.4$. All combinations of correlation and signal strength are covered in the simulations. The results are summarized in Table 1. In Experiment 2, we set $m=20$, $p=100$, $n=25$, $q=10$, $r=5$, and varied the correlation $\rho=0.1,\,0.5,\,0.9$ and signal strength $b=0.1,\,0.2,\,0.3$. The corresponding results are reported in Table 2. In both tables, MSE($A$) and MSE($XA$) denote the $40\%$ trimmed-means of $100\cdot\|A-\hat{B}\|_{F}^{2}/(pn)$ and $100\cdot\|XA-X\hat{B}\|_{F}^{2}/(mn)$, respectively. We also report the median rank estimates (RE) and the successful rank recovery percentages (RRP).

(i) We found that the RSC estimator corresponding to the adaptive choice of the tuning parameter $\mu_{adap}=2S^{2}(n+q)$ has excellent performance. It behaves as well as the RSC estimator that uses the parameter $\mu$ tuned on the large validation set or the RSC estimator corresponding to the theoretical $\mu=2\sigma^{2}(n+q)$.

(ii) When the signal-to-noise ratio SNR := ${d_{r}(XA)}/{(\sqrt{q}+\sqrt{n})}$ is moderate or high, with values approximately 1, 1.5 and 2, corresponding to $b=0.2,0.3,0.4$, and for low to moderate correlation between the predictors ($\rho=0.1,0.5$), RSC has excellent behavior in terms of rank selection and means squared errors. Interestingly, NNP does not have optimal behavior in this set-up: its mean squared errors are slightly higher than those of the RSC estimator. When the noise is very large relative to the signal strength, corresponding to $b=0.1$ in Table 1, or when the correlation between some covariates is very high, $\rho=0.9$ in Table 1, NNP may be slightly more accurate than the RSC.

(iii) The NNP does not recover the correct rank, when its regularization parameter is tuned by validation. Both Tables 1 and 2 show that the correct rank $r$ ($r=10$ in Experiment 1 and $r=5$ in Experiment 2) is overestimated by NNP. Our trimmed estimator, NNP(c), provides a successful improvement over NNP in this respect. This supports Theorem 13.

In additional simulations, we found that especially for low or moderate SNRs, the NNP parameter tuning problem is much more challenging than the RSC parameter tuning. NNP cannot accurately estimate $A$ and consistently select the rank at the same time, for the same value of the tuning parameter. This echoes the findings presented in Figure 1, and is to be expected: in NNP regularization, the threshold value $\tau$ also controls the amount of shrinkage, which should be mild for large samples with relatively low contamination. This is the case for moderate SNR and moderate correlation between predictors: the tuned $\tau$ tends to be too small, so it cannot introduce enough sparsity. The same continues to be true for slightly larger values of $\tau$ that compensate for high noise level and very high correlation between predictors. In summary, one may not be able to build an accurate and parsimonious model via the NNP method, without further adjustments.

Overall, RSC is recommended over the NNP estimators, especially when we suspect that the SNR is not very low. With large validation tuning, NNP(c) has the same properties as RSC – they coincide when both methods select the same rank. But in general, the rank estimation via NNP(c) is much more difficult to tune and much more computationally involved than RSC.

For data with low SNR, an immediate extension of the RSC estimator that involves a second penalty term, of ridge-type, may induce the right amount of shrinkage needed to offset the noise in the data. This conjecture will be investigated carefully in future research.

2 Tightness of the rank consistency results

It can be shown, using arguments similar to those used in the proof of Theorem 2, that

On the other hand, the proof of Theorem 2 reveals that

Appendix A Proof of Theorem 12

that holds on the event $d_{1}(X^{\prime}E)\leq\tau$. The inequality can be deduced from the proof of Theorem 10. Then, by Lemmas 3.4 and 2.3 in Recht et al (2007) there exist two matrices $\widetilde{A}_{1}$ and $\widetilde{A}_{2}$ such that

$\widetilde{A}=\widetilde{A}_{1}+\widetilde{A}_{2}$

$\|\widetilde{A}-B\|_{1}=\|\widetilde{A}_{1}-B\|_{1}+\|\widetilde{A}_{2}\|_{1}$

$\|\widetilde{A}-B\|_{F}^{2}=\|\widetilde{A}_{1}-B\|_{F}^{2}+\|\widetilde{A}_{2}\|_{F}^{2}\geq\|\widetilde{A}_{1}-B\|_{F}^{2}$

$\|\widetilde{A}\|_{1}=\|\widetilde{A}_{1}\|_{1}+\|\widetilde{A}_{2}\|_{1}$.

Using $\lambda_{p}(M)\|\widetilde{A}-B\|_{F}^{2}\leq\|X\widetilde{A}-XB\|_{F}^{2}$ and $2xy\leq x^{2}/2+2y^{2}$, we obtain

The proof is complete by choosing the truncated GSVD $B^{\prime}$ under metric $M$, see Lemma 14 below. ∎

Appendix B Generalized singular value decomposition

with $M=X^{\prime}X=NN$ and $B_{0}$ is a fixed $p\times n$ matrix of rank $r$. By the Eckhart-Young theorem, we have the lower bound

for all $p\times n$ matrices $B$ of rank $k$. We now show that this infimum is achieved by the generalized singular value decomposition (GSVD) under metric $M$, limited to its $k$ largest generalized singular values. Following Takane and Hunter (2001, pages 399-400), the GSVD of $B_{0}$ under metric $M$ is $UDV^{\prime}$ where $U$ is an $p\times r$ matrix, $U^{\prime}MU=I_{r}$, $V$ is an $n\times r$ matrix, $V^{\prime}V=I_{r}$ and $D$ is a diagonal $r\times r$ matrix, and $NB_{0}=NUDV^{\prime}.$ It can be computed via the (regular) SVD $\bar{U}\bar{D}\bar{V}^{\prime}$ of $NB_{0}$. From $B_{0}^{\prime}X^{\prime}XB_{0}=VD^{2}V^{\prime}$, the generalized singular values $d_{j}$ are the regular singular values of $NB_{0}$. Let $B_{k}=U_{k}D_{k}V_{k}^{\prime}$ by retaining as usual the first $k$ columns of $U$ and $V$.

Let $B_{k}$ be the GSVD of $B_{0}$ under metric $M$, restricted to the $k$ largest generalized singular values. We have

Since $NB_{0}=NUDV^{\prime}$ and $NB_{k}=NU_{k}D_{k}V_{k}^{\prime}$, we obtain

using the notation $U_{(k)}$ for the $p\times(r-k)$ matrix consisting of the last $r-k$ column vectors of $U$, $D_{(k)}$ is the diagonal $(r-k)\times(r-k)$ matrix based on the last $r-k$ singular values, and $V_{(k)}$ for the $n\times(r-k)$ matrix consisting of the last $r-k$ column vectors of $V$. Finally,

Recall that in the construction of the GSVD, the generalized singular values $d_{j}$ are the singular values of $NB_{0}$. Since

Remark. The rank restricted estimator $\widehat{B}_{k}$ given in Section 2.1 is the GSVD of the least squares estimator $\widehat{B}$ under the metric $M=X^{\prime}X$, see Takane and Hwang (2007).

Appendix C Largest singular values of transformations of subGaussian matrices

We call a random variable $W$ subGaussian with subGaussian moment $\Gamma_{W}$, if

for all $t>0$. Markov’s inequality implies that $W$ has Gaussian type tails:

holds for any $t>0$. Normal $N(0,\sigma^{2})$ random variables are subGaussian with $\Gamma_{W}=\sigma^{2}$. General results on the largest singular values of matrices $E$ with subGaussian entries can be found in the survey paper by Rudelson and Vershynin (2010). The analysis of our estimators require bounds for the largest singular values of $PE$ and $X^{\prime}E$, for which the standard results on $E$ do not apply directly.

Let $E$ be a $m\times n$ matrix with independent subGaussian entries $E_{ij}$ with subGaussian moment $\Gamma_{E}$. Let $X$ be an $m\times p$ matrix of rank $q$ and let $P=X(X^{\prime}X)^{-}X^{\prime}$ be the projection matrix on $R[X]$. Then, for each $x>0$,

with $U=PS^{p-1}=\{u=Ps:\ s\in S^{p-1}\}$. Let $\mathcal{M}$ be a $\delta$-net of $U$ and $\mathcal{N}$ be a $\delta$-net for $S^{n-1}$ with $\delta=1/2$. Since the dimension of $U$ is $q$ and$\|u\|\leq 1$ for each $u\in U$, we need at most $5^{q}$ elements in $\mathcal{M}$ to cover $U$ and $5^{n}$ elements to cover $S^{n-1}$, see Kolmogorov and Tikhomirov (1961). A standard discretization trick, see, for instance, Rudelson and Vershynin (2010, proof of Proposition 2.4), gives

Next, we write $<u,Ev>=\sum_{i=1}^{m}u_{i}<E_{i},v>$ and note that each $<E_{i},v>$ is subGaussian with moment $\Gamma_{E}$, as

It follows that each term in $\sum_{i=1}^{m}u_{i}<E_{i},v>$ is subGaussian, and $<u,Ev>$ is subGaussian with subGaussian moment $\Gamma_{E}\sum_{i=1}^{m}u_{i}^{2}=\Gamma_{E}$. This implies the tail bound

for each fixed $u$ and $v$ and all $t>0$. Combining the previous two steps, we obtain

for all $t>0$. Taking $t^{2}=2\{\ln(5)(n+q)+x\}\Gamma_{E}$ we obtain the first claim. The second claim follows from this tail bound. ∎

Appendix D Auxiliary results

The following string of inequalities are self-evident:

This proves our first claim. The second claim is easily deduced as follows:

Let $Z_{d}$ be a $\chi^{2}_{d}$ random variable with $d$ degrees of freedom. Then

See Cavalier et al (2002, page 857) for the first claim. The second claim follows by taking $x=t(d/2)^{1/2}$. ∎

Acknowledgement. We would like to thank Emmanuel Candès, Angelika Rohde and Sasha Tsybakov for stimulating conversations in Oberwolfach, Tallahassee and Paris, respectively. We also thank the associate editor and the referees for their constructive remarks.

References