High-Dimensional Asymptotics of Prediction: Ridge Regression and Classification

Introduction

There are various enabling hypotheses that allow for successful prediction in high dimensions. These encode domain-specific knowledge and guide model fitting. Popular options include the “sparsity hypothesis”, i.e., that there is a good predictive rule depending only on $w\cdot x$ for some sparse weight vector $w$ (Candès and Tao, 2007; Hastie et al., 2015), the “manifold hypothesis” positing that the $x_{i}$ have useful low-dimensional geometric structure (Rifai et al., 2011; Simard et al., 2000), and several variants of an “independence hypothesis” that rely on independence assumptions for the feature distribution (Bickel and Levina, 2004; Ng and Jordan, 2001). The choice of enabling hypothesis is important from a practical perspective, as it helps choose which predictive method to use, e.g., the lasso with sparsity, neighborhood-based methods under the manifold hypothesis, or naive Bayes given independent features.

There are several applications, however, where the above enabling hypotheses are not known to apply, and where practitioners have achieved accurate high-dimensional prediction using dense—i.e., non-sparse—ridge-regularized linear methods trained on highly correlated features. One striking example is the case of document classification with dictionary-based features of the form “how many times does the $j$-th word in the dictionary appear in the current document.” Even though $p\gg n$, dense ridge-regularized methods reliably work well across a wide range of problem settings (Sutton and McCallum, 2006; Toutanova et al., 2003), and sometimes even achieve state-of-the-art performance on important engineering tasks (Wang and Manning, 2012). As another example, in a recent bioinformatics test of prediction algorithms (Bernau et al., 2014), ridge regression—and a method that was previously proposed by those same authors—performed best, better than lasso regression and boosting.

The goal of this paper is to gain better understanding of when dense, ridge-regularized linear prediction methods can be expected to work well. We focus on a random-effects hypothesis: we assume that the effect size of each feature is drawn independently at random. This can be viewed as an average-case analysis over dense parameters. Our hypothesis is of course very strong; however, it yields a qualitatively different theory for high-dimensional prediction than popular approaches, and thus may motivate future conceptual developments.

Each predictor has a small, independent random effect on the outcome.

This hypothesis is fruitful both conceptually and methodologically. Using random matrix theoretic techniques (see, e.g., Bai and Silverstein, 2010), we derive closed-form expressions for the limiting predictive risk of idge-regularized regression and discriminant analysis, allowing for the features $x$ to have a general covariance structure $\Sigma$. The resulting formulas are pleasingly simple and depend on $\Sigma$ through the Stieltjes transform of the limiting empirical spectral distribution. More prosaically, $\Sigma$ only enters into our formulas through the almost-sure limits of $p^{-1}\,\operatorname{tr}((\widehat{\Sigma}+\lambda I_{p\times p})^{-1})$ and $p^{-1}\operatorname{tr}((\widehat{\Sigma}+\lambda I_{p\times p})^{-2})$, where $\widehat{\Sigma}$ is the sample covariance and $\lambda>0$ the ridge-regularization parameter. Notably, the same mathematical tools can describe the two problems.

From a practical perspective, we identify several high-dimensional regimes where mildly regularized discriminant analysis performs strikingly well. Thus, it appears that the random-effects hypothesis can at least qualitatively reproduce the empirical successes of Bernau et al. (2014), Sutton and McCallum (2006), Toutanova et al. (2003), Wang and Manning (2012), and others. We hope that further work motivated by generalizations of the random-effects hypothesis could yield a new theoretical underpinning for dense high-dimensional prediction.

We begin with an informal overview of our results; in Section 1.4, we switch to a formal and fully rigorous presentation. In this paper we analyze the predictive risk of ridge-regularized regression and classification when $n,\,p\rightarrow\infty$ jointly. We work in a high-dimensional asymptotic regime where $p/n$ converges to a limiting aspect ratio $p/n\rightarrow\gamma>0$. The spectral distribution—i.e., the cumulative distribution function of the eigenvalues—of the feature covariance matrix $\Sigma$ converges weakly to a limiting spectral measure supported on $[0,\,\infty)$. This allows $\Sigma$ to be general, and we will see several examples later. In random matrix theory, this framework goes back to Marchenko and Pastur (1967); see, e.g., Bai and Silverstein (2010). It has been used in statistics and wireless communications by, among others, Couillet and Debbah (2011), Serdobolskii (2007), Tulino and Verdú (2004), and Yao et al. (2015).

The predictive risk of ridge regression, i.e.,

has an almost-sure limit under high-dimensional asymptotics. This limit only depends on the signal strength $\alpha^{2}$, the aspect ratio $\gamma$, the regularization parameter $\lambda$, and the Stieltjes transform of the limiting eigenvalue distribution of $\widehat{\Sigma}$. For the optimal tuning parameter, $\lambda^{*}=\gamma/\alpha^{2}$

where $v$ is the companion Stieltjes transform of the limiting eigenvalue distribution of $\widehat{\Sigma}$, defined in Section 1.4.

The required functionals of the limiting empirical eigenvalue distribution can be written in terms of almost-sure limits of simple quantities. For example, the result (1) can be written as

For general $\lambda$, the limiting error rate depends on the almost sure limits of both $p^{-1}\,\operatorname{tr}((\widehat{\Sigma}+\lambda I_{p\times p})^{-1})$ and $p^{-1}\operatorname{tr}((\widehat{\Sigma}+\lambda I_{p\times p})^{-2})$.

Thanks to the simple form of (1), we can use our results to gain qualitative insights about the behavior of ridge regression. We show that, when the signal-to-noise ratio is high, i.e., $\alpha\gg 1$, the accuracy of ridge regression has a sharp phase transition at $\gamma=1$ regardless of $\Sigma$, essentially validating a conjecture of Liang and Srebro (2010) on the “regimes of learning” problem. We also find that ridge regression obeys an inaccuracy principle, whereby there are no correlation structures $\Sigma$ for which prediction and estimation of $w^{*}$ are both easy. For $\gamma=1$ this simplifies to

this bound is tight for optimally-tuned ridge regression. We refer to Section 2.2 for the general relation. We find the simplicity of the inverse relation remarkable.

In the second part of the paper, we study regularized discriminant analysis in the two-class Gaussian problem

In high dimensions, and in the metric induced by $\Sigma$, the angle between $w^{*}$ and $\hat{w}_{\lambda}$, i.e.,

has an almost-sure limit. The classification error of regularized discriminant analysis converges to an almost-sure limit that depends only on this limiting angle, as well as the limiting Bayes error. The limiting risk can be expressed in terms of $\alpha$, $\gamma$, $\lambda$, as well the Stieltjes transform of the limit eigenvalue distribution of the empirical within-class covariance matrix.

We can again use our result to derive qualitative insights about the behavior of RDA. We find that the limiting angle between $\hat{w}_{\lambda}$ and $w^{*}$ converges to a non-trivial quantity as $\alpha^{2}\to\infty$, implying that our analysis is helpful in understanding the asymptotics of RDA even in a very high signal-to-noise regime. Finally, by studying the limits $\lambda\rightarrow 0$ and $\lambda\rightarrow\infty$, we can recover known high-dimensional asymptotic results about Fisher’s linear discriminant analysis and naive Bayes methods going back to Bickel and Levina (2004), Raudys (1967), Saranadasa (1993), and even early work by Kolmogorov.

Mathematically, our results build on recent advances in random matrix theory. The main difficulty here is finding explicit limits of certain trace functionals involving both the sample and the population covariance matrix. For instance, the Stieltjes transform $m$ of the empirical spectral distribution satisfies $m(-\lambda)=\lim_{p\to\infty}\smash{p^{-1}\operatorname{tr}((\widehat{\Sigma}+\lambda I_{p\times p})^{-1})}$. However, standard random matrix theory does not provide simple expressions for the limits of functionals like $\smash{p^{-1}\operatorname{tr}(\Sigma(\widehat{\Sigma}+\lambda I_{p\times p})^{-1})}$ or $\smash{p^{-1}\operatorname{tr}([\Sigma(\widehat{\Sigma}+\lambda I_{p\times p})]^{-2})}$ that involve both $\Sigma$ and $\widehat{\Sigma}$. For this we leverage and build on recent results, including the work of Chen et al. (2011), Hachem et al. (2007), and Ledoit and Péché (2011). Our contributions include some new explicit formulas, for which we refer to the proofs. These formulas may prove useful for the analysis of other statistical methods under high-dimensional asymptotics, such as principal component regression and kernel regression.

2 A First Example

A key contribution of our theory is a precise understanding of the effect of correlations between the features on regularized discriminant analysis. Correlated features have a non-trivial effect, and cannot be summarized using standard notions such as the condition number of $\Sigma$ or the classification margin. The full eigenvalue spectrum of $\Sigma$ matters. This is in contrast with popular analyses of high-dimensional classification methods, in which the bounds often depend on the operator norm $\|\Sigma\|_{op}$ (see for instance the review Fan et al. (2011)), thus suggesting that existing analyses of many classification methods are not sharp.

Consider the following examples: First, $\Sigma$ has eigenvalues corresponding to evenly-spaced quantiles of the standard Exponential distribution; Second, $\Sigma$ has a depth-$d$ BinaryTree covariance structure that has been used in population genetics to model the correlations between populations whose evolutionary history is described by a balanced binary tree (Pickrell and Pritchard, 2012). In both cases, we set the class means $\mu_{y}$ such as to keep the Bayes error constant across experiments. Figure 1 plots our formulas for the asymptotic error rate along with empirical realizations of the classification error.

Both covariance structures are far from the identity, and have similar condition numbers. However, the Exponential problem is vastly more difficult for RDA than the BinaryTree problem. This example shows that classical notions like the classification margin or the condition number of $\Sigma$ cannot satisfactorily explain the high-dimensional predictive performance of RDA; meanwhile, our asymptotic formulas are accurate even in moderate sample sizes. Our computational results are reproducible and open-source software to do so is available from https://github.com/dobriban/high-dim-risk-experiments/.

3 Related Work

Random matrix theoretic approaches have been used to study regression and classification in high-dimensional statistics (Serdobolskii, 2007; Yao et al., 2015), as well as wireless communications (Couillet and Debbah, 2011; Tulino and Verdú, 2004). Various regression and M-estimation problems have been studied in high dimensions using approximate message passing (Bayati and Montanari, 2012; Donoho and Montanari, 2015) as well as methods inspired by random matrix theory (Bean et al., 2013). We also note a remarkably early random matrix theoretic analysis of regularized discriminant analysis by Serdobolskii (1983). In the wireless communications literature, the estimation properties of ridge regression are well understood; however its prediction error has not been studied.

El Karoui and Kösters (2011) study the geometric sensitivity of random matrix results, and discuss the consequences to ridge regression and regularized discriminant analysis, under weak theoretical assumptions. In contrast, we make stronger assumptions that enable explicit formulas for the limiting risk of both methods, and allow us to uncover several qualitative phenomena. Our use of Ledoit and Péché (2011)’s results simplifies the proof.

We review the literature focusing on ridge regression or RDA specifically in Sections 2.3 and 3.5 respectively. Important references include, among others, Bickel and Levina (2004), Dicker (2014), El Karoui (2013), Fujikoshi et al. (2011), Hsu et al. (2014), Saranadasa (1993), and Zollanvari and Dougherty (2015).

4 Basics and Notation

(high-dimensional asymptotics) The following conditions hold.

The sample size $n\to\infty$ while the dimensionality $p\to\infty$ as well, such that the aspect ratio $p/n\to\gamma>0$.

The spectral distribution $\smash{F_{\Sigma}}$ of $\Sigma$ converges to a limit probability distribution $H$ supported $\smash{[0,\,\infty)}$, called the population spectral distribution (PSD).

Families of covariance matrices $\Sigma$ that fit the setting of this theorem include the identity covariance, BinaryTree, Exponential and the autoregressive AR-1 model with $\Sigma_{ij}=\rho^{|i-j|}$ (see Grenander and Szegő, 1984, for the last one).

Under assumption Assumption A, the spectral distribution $\smash{F_{\widehat{\Sigma}}}$ of the sample covariance matrix $\smash{\widehat{\Sigma}}$ also converges weakly, with probability 1, to a limiting distribution supported on $[0,\,\infty)$.

The limiting distrbution $F$ is called the empirical spectral distribution (ESD), and is determined uniquely by a fixed point equation for its Stieltjes transform, which is defined for any distribution $G$ supported on $[0,\infty)$ as

Given this notation, the Stieltjes transform of the spectral measure of $\smash{\widehat{\Sigma}}$ satisfies

In addition, we write the derivatives asWe will denote by $v^{\prime}(-\lambda)$ the derivative of the Stieltjes transform, $v^{\prime}(z)$, evaluated at $z=-\lambda$; and not the derivative of the function $\lambda\to v(-\lambda)$. $\smash{m^{\prime}(z)=\int_{l=0}^{\infty}{dG(l)}/\left(l-z\right)^{2}}$ and $\smash{v^{\prime}(z)=\gamma(m^{\prime}(z)-z^{-2})+z^{-2}}$. These derivatives can also be understood in terms of empirical observables, through the relation

Finally, our analysis also relies on several more recent formulas for limits of trace functionals involving both $\Sigma$ and $\widehat{\Sigma}$. In particular, we use a formula due to Ledoit and Péché (2011), who in the analysis of eigenvectors of sample covariance matrices showed that, under certain moment conditions:See the supplement for more details about this result.

Predictive Risk of Ridge Regression

Writing $\gamma_{p}=p/n$ and $\lambda_{p}^{*}=\gamma_{p}\alpha^{-2}$, the finite sample predictive risk $r_{\lambda^{*}_{p}}(X)$ converges almost surely

Moreover, for any $\lambda>0$, the predictive risk converges almost surely to the limiting predictive risk $R_{\lambda}(H,\alpha^{2},\gamma)$, where

The choice $\lambda^{*}$ minimizes $R_{\lambda}(H,\alpha^{2},\gamma)$.

The proof of part 1 is sufficiently simple to outline here; see the supplement for part 2. We begin by verifying formula (6):

On the last line we used the choice of $\lambda_{p}^{*}$. From the results of Ledoit and Péché (2011), and from $\gamma_{p}\to\gamma$, it can be verified that $\smash{p^{-1}\operatorname{tr}(\Sigma(\widehat{\Sigma}+\gamma_{p}\alpha^{-2}I_{p\times p})^{-1})}$ converges almost surely to limit in (5), finishing part 1. ∎

This result fully characterizes the first order behavior of the predictive risk of ridge regression under high-dimensional asymptotics. To verify its finite-sample accuracy, we perform a simulation with the BinaryTree and Exponential models. We compute the limit risks using the algorithms in the supplement. The results in Figure 2 show that the formulas given in Theorem 2.1 appear to be accurate, even in small sized problems. In Figure 2, for BinaryTree we train on $n=\gamma^{-1}p$ samples, where $p=2^{4}$; for Exponential on $n=20$. We set the signal strength to $\alpha^{2}=1$ and generated $w$, $X$, and $\varepsilon$ as Gaussian random variables with i.i.d. entries and the desired variance. The results are averaged over 500 simulation runs; we evaluated the empirical prediction error using a test set of size 100.

Intriguingly, Figure 2 shows that the prediction performance of ridge regression is very similar on the two problems. This presents a marked contrast to the RDA example given in the introduction, where the two covariance structures led to very different classification performance. Thus, it appears that the loss function and the spectrum of $\Sigma$ can interact non-trivially.

Meanwhile, in the special case of identity covariance, the quantity $R_{\lambda}-1$ coincides with the normalized estimation error, so we recover known results described in, e.g., Tulino and Verdú (2004). When $\smash{\Sigma=I_{p\times p}}$, we have an explicit expression for the Stieltjes transform (e.g., Bai and Silverstein, 2010, p. 52), valid for $\lambda>0$:

Theorem 2.1 implies that the limit predictive risk of ridge regression for general $\lambda$ equals

which has an explicit form. Furthermore, the optimal risk has a particularly simple form:

See the supplement for details on these derivations.

As an application of Theorem 2.1, we study how the difficulty of ridge regression depends on the signal strength $\alpha^{2}$. Liang and Srebro (2010) call this the regimes of learning problem and argue that, for small $\alpha^{2}$ the complexity of ridge regression should be tightly characterized by dimension-independent Rademacher bounds, while for large $\alpha^{2}$ the error rate should only depend on $\gamma$. Liang and Srebro (2010) justify their claims using generalization bounds for the identity-covariance case $\Sigma=I_{p\times p}$, and conjecture that similar relationships should hold in general. Using our results, we can give a precise characterization of the regimes of learning of optimally-regularized ridge regression with general covariance $\Sigma$.

From Theorem 2.1, we know that given a signal strength $\alpha^{2}$, the asymptotically optimal choice for $\lambda$ is $\smash{\lambda^{*}(\alpha)=\gamma/\alpha^{2}}$, in which case the predictive risk of ridge regression converges to

We now use this formula to examine the two limiting behaviors of the risk, for weak and strong signals. The results below are proved in the supplement, assuming that the population spectral distribution $H$ is supported on a set bounded away from 0 and infinity.

Conversely, the strong-signal limiting behavior of the risk depends the aspect ratio $\gamma$, and experiences a phase transition at $\gamma=1$. When $\gamma<1$, the predictive risk converges to

regardless of $\Sigma$. This quantity is known to be the $n,\,p\rightarrow\infty$, $p/n\to\gamma$ limit of the risk of ordinary least squares (OLS). Thus when $p<n$ and we have a very strong signal, ridge regression cannot outperform OLS, although of course it can do much better with a small $\alpha$.

When $\gamma>1$, the risk $R^{*}(H,\alpha^{2},\gamma)$ can grow unboundedly large with $\alpha$. Moreover, we can verify that

Thus, the limiting error rate depends on the covariance matrix through $v(0)$. In general there is no closed-form expression for $v(0)$, which is instead characterized as the unique $c>0$ for which

In the special case $\Sigma=I_{p\times p}$, however, the limiting expression simplifies to $\smash{1/(\gamma v(0))}$ = $\smash{(\gamma-1)/\gamma}$. In other words, when $p>n$, optimally tuned ridge regression can capture a constant fraction of the signal, and its test-set fraction of explained variance tends to $\gamma^{-1}$.

Finally, in the threshold case $\gamma=1$, the risk $R^{*}(H,\alpha^{2},\gamma)$ scales with $\alpha$:

In summary, we find that for general covariance $\Sigma$, the strong-signal risk $\smash{R^{*}(\alpha^{2},\,\gamma)}$ scales as $\smash{\Theta(1)}$ if $\gamma<1$, as $\smash{\Theta(\alpha)}$ if $\gamma=1$, and as $\smash{\Theta(\alpha^{2})}$ if $\gamma>1$. We illustrate this phenomenon in Figure 3, in the case of the identity covariance $\Sigma=I_{p\times p}$. We see that when $\gamma<1$ the error rate stabilizes, whereas when $\gamma>1$, the error rate eventually gets a slope of 1 on the log-log scale. Finally, when $\gamma=1$, the error rate has a log-log slope of $1/2$.

Thus, thanks to Theorem 2.1, we can derive a complete and exact answer the regimes of learning question posed by Liang and Srebro (2010) in the case of ridge regression. The results (9) and (10) not only show that the scalings found by Liang and Srebro (2010) with $\smash{\Sigma=I_{p\times p}}$ hold for arbitrary $\Sigma$, but make explicit how the slopes depend on the limiting population spectral distribution. The ease with which we were able to read off this scaling from Theorem 2.1 attests to the power of the random matrix approach.

2 An Inaccuracy Principle for Ridge Regression

where $m$ is the Stieltjes transform of the limiting empirical spectral distribution (see, e.g., Tulino and Verdú, 2004, Chapter 3). Combining this result with Theorem 2.1 and (4), we find the following relationship between the limiting predictive risk $R_{P}$ and the limiting estimation risk $R_{E}$.

Under the conditions of Theorem 2.1, the asymptotic predictive and estimation risks of optimally-tuned ridge regression are inversely related,

The equation holds for all limit eigenvalue distributions $H$ of the covariance matrices $\Sigma$. Both sides of the above equation are non-negative: $R_{P}$ cannot fall below the intrinsic noise level \operatorname{Var}\left[Y\,\big{|}\,X\right]=1, while $R_{E}\leq\limsup_{p\to\infty}R_{E,\,n}(\lambda^{*})\leq\limsup_{p\to\infty}R_{E,\,n}(0)=\alpha^{2}$. When $\gamma=1$, we get the even simpler equation

The product of the estimation and prediction risks equals the signal strength. Since this holds for the optimal $\lambda^{*}$, it also implies that for any $\lambda$ we have the lower bound $R_{E}(\lambda)\cdot R_{P}(\lambda)\geq\alpha^{2}$; we find the explicit formula relating the two risks remarkable. The inverse relationship may be somewhat surprising, but it has an intuitive explanation.A similar heuristic was given by Liang and Srebro (2010), without theoretical justification. When the features are highly correlated and $v$ is correspondingly large, prediction is easy because $y$ lies close to the “small” column space of the feature matrix $X$, but estimation of $w$ is hard due to multi-collinearity. As correlation decreases, prediction gets harder but estimation gets easier.

3 Related Work for High-Dimensional Ridge Regression

Regularized Discriminant Analysis

In the second part of the paper, we return to regularized discriminant analysis and the two-class Gaussian discrimination problem (2). For simplicity we will first discuss the case when the population label proportions are balanced. In this case, the Bayes oracle predicts using (Anderson, 2003)

and has an error rate $\text{Err}_{\text{Bayes}}=\Phi\left(-\Delta_{n,p}\right)$, where $\Delta_{n,p}=\sqrt{\delta^{\top}\Sigma^{-1}\delta}$ is half the between-class Mahalanobis distance. The Gaussian classification problem has a rich history, going back to Fisher’s pioneering work on linear discriminant analysis (LDA). When we have the same number of examples from both the positive and negative classes, i.e., $\smash{n_{-1}=n_{+1}=n/2}$, LDA classifies using the linear rule

Here $\widehat{\Sigma}_{c}$ is the centered covariance matrix. In the low-dimensional case where $n$ gets large while $p$ remains fixed, LDA is the natural plug-in rule for Gaussian classification and efficiently converges to the Bayes discrimination function (Anderson, 2003; Efron, 1975). When $p$ is on the same order as $n$, however, the matrix inverse $\smash{\widehat{\Sigma}_{c}^{-1}}$ is unstable and the performance of LDA declines, as discussed among others by Bickel and Levina (2004). Instead, we will study regularized discriminant analysis, the linear rule $\smash{\hat{y}=h_{\hat{w}_{\lambda}}\left(x-(\hat{\mu}_{-1}+\hat{\mu}_{+1})/{2}\right)}$ where $\smash{h_{w}(x)=\operatorname{sign}(w\cdot x)}$ and the weight vector is The notation was chosen to emphasize the similarities between ridge regression and RDA. There will be no possibility for confusion with the ridge regression weight vector, also denoted $\hat{w}_{\lambda}$. $\smash{\hat{w}_{\lambda}=(\widehat{\Sigma}_{c}+\lambda I_{p\times p})^{-1}\hat{\delta}}$ (Friedman, 1989; Serdobolskii, 1983).

(random weights in classification) The following conditions hold:

$\mu_{-1}$ and $\mu_{+1}$ are randomly generated as $\mu_{-1}=\bar{\mu}-\delta$ and $\mu_{+1}=\bar{\mu}+\delta$, where $\delta$ has i.i.d. coordinates with

for some fixed constants $\eta>0$ and $C$.

$\smash{\bar{\mu}=(\mu_{-1}+\mu_{+1})/2}$ is either fixed, or random and independent of $\delta$, $X$ and $y$, and satisfies $\limsup_{p\to\infty}\|\bar{\mu}\|^{2}_{2}/p^{1/2-\zeta}\leq C$ almost surely for some fixed constants $\zeta>0$ and $C$.

and $\tau$, $\eta$, and $\xi$ are determined by the problem parameters $H$ and $\gamma$:

Here, $m=m(-\lambda)$ is the Stieltjes transform of the limit empirical spectral distribution $F$ of the covariance matrix $\smash{\widehat{\Sigma}_{c}}$, and $v=v(-\lambda)$ is the companion Stieltjes transform defined in (4).

The proof of Theorem 3.1 (in the supplement) is similar to Theorem 2.1, but more involved. The main difficulty is to evaluate explicitly the limits of certain functionals of the population and sample covariance matrices. We extend the result of Ledoit and Péché (2011), and rely on additional results and ideas from Hachem et al. (2007) and Chen et al. (2011). In particular, we use a derivative trick for Stieltjes transforms, similar to that employed in a related context by El Karoui and Kösters (2011), Rubio et al. (2012), and Zhang et al. (2013). The limits are then combined with results on concentration of quadratic forms.

Under the conditions of Theorem 3.1, and with unequal sampling, the classification error of RDA converges almost surely:

where the effective classification margins have the form

In this section we assumed Gaussianity, but by a Lindeberg–type argument it should be possible to extend the result to non-Gaussian observations with matching moments. It should also be interesting to study how sensitive the results are to the particular assumptions of our model, such as independence across samples.

It is worth mentioning that the regression and classification problems are very different statistically. In the random effects linear model, ridge regression is a linear Bayes estimator, thus the ridge regularization $\widehat{\Sigma}+\lambda I_{p\times p}$ of the covariance matrix is justified statistically. However, for classification, the ridge regularization is merely a heuristic to help with the ill-conditioned sample covariance. It is thus interesting to know how much this heuristic helps improve upon un-regularized LDA, and how close we get to the Bayes error. We now turn to this problem, which can be studied equivalently from a geometric perspective.

2 The Geometry of RDA

The asymptotics of RDA can be understood in terms of a simple picture. The angle between the Bayes decision boundary hyperplane and the RDA discriminating hyperplane tends to an asymptotically deterministic value in the metric of the covariance matrix, and the limiting risk of RDA can be described in terms of this angle.

Recall that, in the balanced $\pi_{+}=\pi_{-}$ and $n_{+}=n_{-}$ case, the estimated RDA weight vector is $\hat{w}_{\lambda}=(\widehat{\Sigma}_{c}+\lambda I_{p\times p})^{-1}\hat{\delta}$, while the Bayes weight vector is $w^{*}=\Sigma^{-1}\delta$. Writing $\langle a,b\rangle_{\Sigma}=a^{\top}\Sigma b$ for the inner product in the $\Sigma$-metric, the cosine of the angle—in the same metric—between the two is

where $\Theta\left(H,\gamma,\alpha^{2},\lambda\right)$ is the classification margin of RDA with the dependence on each parameter made explicit. The limit of the cosine quantifies the inefficiency of the RDA estimator relative to the Bayes one.

We gain some insight into this angle for two special cases: when $H=\delta_{1}$, and by taking the limit $\alpha^{2}\to\infty$. First, with $H=\delta_{1}$, or equivalently $\Sigma=I_{p\times p}$, we curiously find that the effects of the signal strength $\alpha^{2}$ and regularization rate $\lambda$ decouple completely, as shown in Corollary 3.3 below. See the supplement for a proof of the following result.

Under the conditions of Theorem 3.1, let $\smash{\Sigma=I_{p\times p}}$ for all $p$. Then, the limiting cosine $\Gamma$ of the angle between the Bayes and RDA hyperplanes is

where the Stieltjes transform $m_{I}(-\lambda;\gamma)$ for $\Sigma=I_{p\times p}$ is given in (7). For $\gamma=1$, this expression simplifies further to

Examining $\Gamma(\delta_{1},\gamma,\alpha^{2},\lambda)$, we can attribute the suboptimality to two sources of noise: We need to pay a price $\smash{{\alpha}/{\sqrt{\alpha^{2}+\gamma}}}$ for estimating $\smash{\mu_{\pm 1}}$, while the cost of estimating $\Sigma$ is $\smash{([1+\gamma\lambda\,m_{I}^{2}(-\lambda;\gamma)]/[1+\gamma\,m_{I}(-\lambda;\gamma)])^{1/2}}$. If we knew that $\Sigma=I$, we could send $\lambda\rightarrow\infty$. It is easy to verify that this would eliminate the second term, leading to a loss of efficiency $\alpha/\sqrt{\alpha^{2}+\gamma}$.

In the case of a general covariance matrix $\Sigma$ we get a similar asymptotic decoupling in the strong-signal limit $\alpha^{2}\rightarrow\infty$. The following claim follows immediately from Theorem 3.1.

Under the conditions of Theorem 3.1, the cosine of the angle between the optimal and learned hyperplanes has the limit as $\alpha^{2}\to\infty$:

Thus, RDA is in general inconsistent for the Bayes hyperplane in the case of strong signals. Corollary 3.4 also implies that, in the limit $\alpha\rightarrow\infty$, the optimal $\lambda$ for RDA converges to a non-trivial limit that only depends on the spectral distribution $H$. No such result is true for ridge regression, where $\lambda^{*}=\alpha^{-2}\gamma\rightarrow 0$ as $\alpha\rightarrow\infty$, regardless of $\Sigma$. This contrast arises because ridge regression is a linear Bayes estimator, while RDA is a heuristic which is strictly suboptimal to the Bayes classifier.

We illustrate the behavior of the cosine $\Gamma$ for the AR-1(0.9) model in Figure 4, which displays $\Gamma$ for values of $\alpha$ ranging from $\alpha=0.1$ to $\alpha=2$. We see that the $\Gamma$-curve converges to its large-$\alpha$ limit fairly rapidly. Moreover, somewhat strikingly, we see that the optimal regularization parameter $\lambda^{*}$, i.e., the maximizer of $\Gamma$, increases with the signal strength $\alpha^{2}$.

Finally, we note that Efron (1975) studies the angle $\Gamma$ in detail for low-dimensional asymptotics where $p$ is fixed while $n\rightarrow\infty$; in this case, $\Gamma$ converges in probability to $1$, and the sampling distribution of $n(1-\Gamma)$ converges to a (scaled) $\chi_{p-1}^{2}$ distribution. Establishing the sampling distribution in high dimensions is interesting future work.

3 Do existing theories explain the behavior of RDA?

Theorems 3.1 and 3.2 give precise information about the error rate of RDA. It is of interest to compare this to classical theories, such as Vapnik-Chervonenkis theory or Rademacher bounds, to if they explain the qualitative behavior of RDA. In this section, we study a simple simulation example, and conclude that existing theory does not precisely explain the behavior of RDA.

Existing results give us some intuition about what to expect. Since $n=p$, classical heuristics based on the theory of Vapnik and Chervonenkis (1971) as well as more specialized analyses (Saranadasa, 1993; Bickel and Levina, 2004) predict that unregularized LDA will not work. As we will see, this matches our simulation results. Meanwhile, Bickel and Levina (2004) study worst-case performance of the independence rule relative to the Bayes rule. In our setting, it can be verified that their results imply $\Theta_{\text{IR}}\geq(1-\rho^{2})/(1+\rho^{2})\,\Delta$, where the error rate of the independence rule is $\Phi(-\Theta_{\text{IR}})$. This predicts that independence rules will work better for small correlation $\rho$, which again will match the simulations.

The existing theory, however, is much less helpful for understanding the behavior of RDA for intermediate values of $\lambda$. A learning theoretic analysis based on Rademacher complexity suggests that the generalization performance of RDA should depend on terms that scale like $\smash{\sqrt{\lVert{\hat{w}_{\lambda}}\rVert_{2}^{2}\operatorname{tr}{\Sigma}/n}\asymp\sqrt{\lambda^{-2}p/n}}$ for large values of $\lambda$ (e.g., Bartlett and Mendelson, 2003). In other words, based on a classical approach, we might expect that mildly regularized RDA should not work, but using a large $\lambda$ may help. Rademacher theory is not tight enough to predict what will happen for $\lambda\approx 1$.

Given this background, Figure 5 displays the performance of RDA for different values of $\rho$, along with our theoretically derived error from Theorem 3.1. In the $\smash{\alpha^{2}=1}$ case, we find that—as predicted—unregularized LDA does poorly. However, when $\rho$ is large, mildly regularized RDA does quite well.

4 Linear Discriminant Analysis vs. Independence Rules

Two points along the RDA risk curve that allow for particularly simple analytic expressions occur as $\lambda\to 0$ and $\lambda\to\infty$: the former is just classical linear discriminant analysis while the latter is equivalent to an independence rule (or “naïve Bayes”). In this section, we show how taking these limits we can recover known results about the high-dimensional asymptotics of linear discriminant analysis and naïve Bayes. Further, we compare these two methods over certain parameter classes.

Note that $\lambda\to\infty$ leads to a linear discriminant rule with weight vector $\smash{\hat{\delta}=\hat{\mu}_{+1}-\hat{\mu}_{-1}}$. Usual independence rules take the form $\smash{\operatorname{diag}(\widehat{\Sigma}_{c})^{-1}\hat{\delta}}$. We will assume that all features are normalized to have equal variance, $\Sigma_{ii}=\sigma>0$. In this case the $\lambda\to\infty$ rule corresponds to an independence rule with oracle information about the equality of variances; which we still call “indpendence rule” for simplicity.

Extending our previous notation, we define the asymptotic margin of LDA and independence rules, by taking the limits of $\Theta(\lambda)$ at 0 and $\infty$:

Both limits are well-defined and admit simple expressions, as given below.

Let $H$ be the limit population spectral distribution of the covariance matrices $\Sigma$; and let $T$ be a random variable with distribution $H$. Under the conditions of Theorem 3.1, the margins of LDA and independence rules are equal to

The formula for LDA is valid for $\gamma<1$ while that for IR is valid for any $\gamma$.

This result is proved in the supplement. The formulas are simpler than Theorem 3.1, as they involve the population spectral distribution $H$ directly through its moments. For RDA, the error rate depends on $H$ implicitly through the Stieltjes transform of the ESD $F$.

These formulas are equivalent to known results, some of which were obtained under slightly different parametrization. In particular, Serdobolskii (2007) attributes the IR formula with $H=\delta_{1}$ to unpublished work by Kolmogorov in 1967, while the Raudys and Young (2004) attributes it to Raudys (1967). The LDA formula was derived by Deev (1970) and Raudys (1972); see Section 3.5. Here, our goal was to show how these simple formulas can be recovered from the more powerful Theorem 3.1.

Theorem 3.5 enables us to compare the worst-case performance of LDA and IR over suitable parameter classes of limit spectra. For $0<k_{1}\leq 1\leq k_{2}$ we define the class

The bounds $k_{1},k_{2}$ control the ill-conditioning of the population covariance matrix. We normalize such that the average population eigenvalue is 1, to ensure that the scaling of the problem does not affect the answer. This parameter space is somewhat similar to the one considered by Bickel and Levina (2004). Perhaps surprisingly, however, a direct comparison over these natural problem classes appears to be missing from the literature, and so we provide it below (see the supplement for a proof).

Under the conditions of Theorem 3.5, consider the behavior of LDA and independence rules for $H\in\mathcal{H}(k_{1},k_{2})$.

The least favorable distribution for LDA from the class $\mathcal{H}(k_{1},k_{2})$ is the point mass at 1: $H=\delta_{1}$, i.e., $\Sigma=I$.

The worst-case margin for independence rules is:

If $k_{1}<k_{2}$ the least favorable distribution is the mixture $H=w_{1}\delta_{k_{1}}+w_{2}\delta_{k_{2}}$, where the weights are $w_{1}=(k_{2}-1)/(k_{2}-k_{1})$ and $w_{2}=(1-k_{1})/(k_{2}-k_{1})$; while if $k_{1}=k_{2}=1$, it is the point mass at 1: $H=\delta_{1}$.

This result shows a stark contrast between the worst-case behavior of LDA and independence rules: for fixed signal strength, the worst-case risk of LDA over $\mathcal{H}$ only depends on $\gamma$, and is attained with the limit of identity covariances $\Sigma=I_{p\times p}$ regardless of the values of $k_{1},\,k_{2}$. In contrast, the worst-case behavior of IR occurs for a least favorable distribution $H$ that is as highly spread as possible. This highlights the sensitivity of IR to ill-conditioned covariance matrices. For $0<\gamma<1$, we see that IR are better than LDA in the worst case over $\mathcal{H}$, i.e. $\bar{\Theta}_{\text{LDA}}(\gamma;\alpha^{2})<\bar{\Theta}_{\text{IR}}(\mathcal{H},\gamma;\alpha^{2})$, if and only if

In particular IR performs better than LDA for strong signals $\alpha$; with weaker signals, LDA can sometimes have an edge, particularly if the covariance is poorly conditioned, quantified by a large measure of spread $k_{1}+k_{2}-k_{1}k_{2}=(k_{2}-1)(1-k_{1})+1$.

5 Literature Review for High-Dimensional RDA

There has been substantial work in the former Soviet Union on high-dimensional classification; references on this work include Raudys and Young (2004), Raudys (2012), and Serdobolskii (2007). Raudys (1967)Raudys (1967) is in Russian; see Raudys and Young (2004) for a review. derived the $n,\,p\to\infty$ asymptotic error rate of independence rules in identity-covariance case $\Sigma=I_{p\times p}$, while Deev (1970) and Raudys (1972) obtained the error rate of un-regularized linear discriminant analysis (LDA) for general covariance $\Sigma$, again in the $n,\,p\to\infty$ regime.Serdobolskii (2007) attributes some early results on independence rules with identity covariance to Kolmogorov in 1967, and calls the framework $n,\,p\to\infty$, $p/n\to\gamma$ the Kolmorogov asymptotic regime. However, Kolmogorov apparently never published on the topic. As shown in the previous section, these results can be obtained as special cases of our more general formulas.

For RDA, Serdobolskii (1983) (see also Chapter 5 of Serdobolskii (2007)) considered a more general setting than this paper: classification with a weight vector of the form $\Gamma(\widehat{\Sigma}_{c})^{-1}\hat{\delta}$ instead of just $(\widehat{\Sigma}_{c}+\lambda I_{p\times p})^{-1}\hat{\delta}$, where the scalar function $\Gamma$ admits the integral representation $\Gamma(x)=\int(x+t)^{-1}\,d\eta(t)$ for a suitable measure $\eta$, and is extended to matrices in the usual way. He derived a limiting formula for the error rate of this classifier under high-dimensional asymptotics. However, due to their generality, his results are substantially more involved and much less explicit than ours. In some cases it is unclear to us how one could numerically compute his formulas in practice. Furthermore, his results are proved when $\gamma<1$, and show convergence in probability, not almost surely. We also note the work of Raudys and Skurichina (1995), who derived results about the risk of usual RDA with vanishingly small regularization $\lambda=o(1)$, and for the special case $\gamma<1$.

In another line of work, a Japanese school (e.g., Fujikoshi et al., 2011, and references therein) has studied the error rates of LDA and RDA under high-dimensional asymptotics, with a focus on obtaining accurate higher-order expansions to their risk. For instance Fujikoshi and Seo (1998) obtained asymptotic expansions for the error rate of un-regularized LDA, which can be verified to be equivalent to our results in the $\lambda\rightarrow 0$ limit. More recently, Kubokawa et al. (2013) obtained a second-order expansion of the error rate of RDA with vanishingly small regularization parameter $\lambda=O(1/n)$ in the case $\gamma<1$.

Finally, in the signal processing and pattern recognition literature, Zollanvari et al. (2011) provided asymptotic moments of estimators of the error rate of LDA, under an asymptotic framework where $n,p\to\infty$; however, this paper assumes that the covariance matrix $\Sigma$ is known. More recently, Zollanvari and Dougherty (2015) provide consistent estimators for the error rate of RDA in a doubly asymptotic framework, using deterministic equivalents for random matrices. The goal of our work is rather different from theirs, in that we do not seek empirical estimators of the error rate of RDA, but instead seek simple formulas that help us understand the behavior of RDA.

Acknowledgment

We are grateful to everyone who has provided comments on this manuscript, in particular David Donoho, Jerry Friedman, Iain Johnstone, Percy Liang, Asaf Weinstein, and Charles Zheng. E. D. is supported in part by NSF grant DMS-1418362.

References

Supplement

The supplement is organized as follows: Section 5 describes the efficient computation of the risk formulas, Section 6 has the proofs for ridge regression, and Section 7 has the proofs for regularized discriminant analysis. At several locations we refer to equation numbers from the main text.

Efficient computation of the risk formulas

Consider the spectral distribution of the companion matrix $\smash{\underline{\widehat{\Sigma}}=n^{-1}XX^{\top}}$. Since its spectral distribution $\smash{F_{\underline{\widehat{\Sigma}}}}$ differs from $\smash{F_{\widehat{\Sigma}}}$ by $|n-p|$ zeros, it follows that $\smash{\underline{\widehat{\Sigma}}}$ has a limit ESD $\smash{\underline{F}}$, given by $\smash{\underline{F}=\gamma F+(1-\gamma)I_{[0,\infty)}}$. The companion Stieltjes transform satisfies the Silverstein equation (Silverstein and Combettes, 1992; Silverstein and Choi, 1995):

We now explain how to compute the key quantities $m,v,m^{\prime},v^{\prime}$ that will come up in our risk formulas. On the interval $v\in[0,\infty)$, Silverstein and Choi (1995) prove that the functional inverse of $z\to v:=v(z)$ has the explicit form:

This result enables the efficient computation of the function $z\to v(z)$ for $z<0$. Indeed, assuming one can compute the corresponding integral against $H$, one can tabulate (13) on a dense grid of $v_{i}>0$, to find pairs $(v_{i},z_{i})$, where $z_{i}=z(v_{i})$. Then for the values $z_{i}<0$, the Silverstein and Choi (1995) result shows that $v(z_{i})=v_{i}$. Further, the Silverstein equation can be differentiated with respect to $z$ to obtain an explicit formula for $v^{\prime}$ in terms of $v,H$:

Therefore, once $v(z)$ is computed for a value $z$, the computation of $v^{\prime}(z)$ can be done conveniently in terms of $v(z)$ and $H$, assuming again that the integral involving $H$ can be computed. This is one of the main steps in the Spectrode method for computing the limit ESD (Dobriban, 2015). Finally, $m(z)$ can be computed from $v(z)$ via the equation (4), and $m^{\prime}(z)$ can be computed from $v^{\prime}(z)$ by differentiating (4): $\gamma\left(m^{\prime}(z)-1/z^{2}\right)=v^{\prime}(z)-1/z^{2}$.

Proofs for Ridge Regression

where $\varepsilon_{0}$ is the noise in the new observation. Now $\hat{w}_{\lambda}=(X^{\top}X+\lambda\,n\,I_{p\times p})^{-1}$ $X^{\top}Y$, and $Y=Xw+\varepsilon$, where $\varepsilon$ is the vector of noise terms in the original data. Hence

When we plug this back into the risk formula $r_{\lambda}(X)$, and use that $w,\varepsilon$ are conditionally independent given $X$, we see that the cross-term involving $w,\varepsilon$ cancels. The risk simplifies to

Now using that the components of $w$ and $\varepsilon$ are each uncorrelated conditional on $X$, we obtain the further simplification

Introducing $\widehat{\Sigma}=n^{-1}X^{\top}X$ and $\gamma_{p}=p/n$, and splitting the last term in two by using $X^{\top}X=X^{\top}X+\lambda\,n\,I_{p\times p}-\lambda\,n\,I_{p\times p}$ this yields

For the particular choice $\lambda^{*}=\gamma_{p}\alpha^{-2}$, we obtain the claimed formula $r_{\lambda*}(X)$. Next, we show the convergence of $r_{\lambda}(X)$ for arbitrary fixed $\lambda$. First, by assumption we have $\gamma_{p}\to\gamma$. Therefore, it is enough to show the almost sure convergence of the two functionals:

The convergence of the first one follows directly from the theorem of Ledoit and Péché (2011), given in (5). The second is shown later in the proof of Lemma 7.4 in Section 7.1.4. In that section it is assumed that the eigenvalues of $\Sigma$ are bounded away from 0 an infinity; but one can check that in the proof of Lemma 7.4 only the upper bound is used, and that holds in our case. Therefore the risk $r_{\lambda}(X)$ converges almost surely for each $\lambda$.

Next we find the fomulas for the limits of the two functionals. The limit of $p^{-1}\operatorname{tr}\left(\Sigma\left(\widehat{\Sigma}+\lambda I_{p\times p}\right)^{-1}\right)$ equals $\kappa(\lambda)=\gamma^{-1}(1/[\lambda v(-\lambda)]-1)$ by (5). In the proof of Lemma 7.4 in Section 7.1.4, it is shown that the limit of $p^{-1}\operatorname{tr}\left(\Sigma\left(\widehat{\Sigma}+\lambda I_{p\times p}\right)^{-2}\right)$ is $-\kappa^{\prime}(\lambda)=[v(-\lambda)-\lambda v^{\prime}(-\lambda)]/[\gamma(\lambda v(-\lambda))^{2}]$.

Simplified expression for $R_{\lambda}$: Putting together the results above, we obtain (with $v=v(-\lambda),v^{\prime}=v^{\prime}(-\lambda)$) the desired claim:

Second part: Convergence: First, we note that for $\lambda_{p}^{*}=\gamma_{p}\alpha^{-2}$, the finite sample risk equals by (14)

Introduce the function $k_{p}(\lambda,X)=\frac{1}{p}\operatorname{tr}\left(\Sigma\left(\widehat{\Sigma}+\lambda I_{p\times p}\right)^{-1}\right)$. We need to show $k_{p}(\lambda_{p}^{*},X)\to\kappa(\lambda^{*})$. First, we notice that $\lambda_{p}^{*}\to\lambda^{*}$, and $k_{p}(\lambda,X)\to\kappa(\lambda)$ almost surely. Second, we verify the equicontinuity of $k_{p}$ as a function of $\lambda$, by proving the stronger claim that the derivatives of $k_{p}$ are uniformly bounded:

Therefore, by the equicontinuity of the family $k_{p}(\lambda,X)$ as a function of $\lambda$, we obtain $k_{p}(\lambda_{p}^{*},X)\to\kappa(\lambda^{*})$. Further, the explicit form of $r_{\lambda_{p}^{*}}$ shows that the limit equals $1+\gamma\kappa(\lambda^{*})=1/(\lambda^{*}v(-\lambda^{*}))$ by (5), as desired.

Optimality of $\lambda^{*}$: The limiting risk $R_{\lambda}$ is the same if we assume Gaussian observations. In this case, the finite sample Bayes-optimal choice for $\lambda_{p}$ is $\lambda_{p}^{*}=\gamma_{p}\alpha^{-2}$. We will use the following classical lemma to conclude that the limit of the minimizers $\lambda_{p}^{*}$ is the minimizer of the limit.

Let $f_{n}(x):\mathcal{I}\to\mathcal{I}$ be an equicontinuous family of functions on an interval $\mathcal{I}$, converging pointwise to a continuous function, $f_{n}(x)\to f(x)$. Suppose $y_{n}$ is a minimizer of $f_{n}$ on $\mathcal{I}$, and $y_{n}\to y$. Then $y$ is a minimizer of $f$.

Since $y_{n}$ is a minimizer of $f_{n}$, we have

for all $x\in\mathcal{I}$. But $f_{n}(y_{n})=f_{n}(y_{n})-f_{n}(y)+f_{n}(y)-f(y)+f(y)$. As $n\to\infty$, $f_{n}(y_{n})-f_{n}(y)\to 0$ by the convergence of $y_{n}\to y$ and by the equicontinuity of the family $\{f_{n}\}$; and $f_{n}(x)-f(x)\to 0$ for all $x$ by the convergence of $f_{n}\to f$. Therefore, taking the limit as $n\to\infty$ in (15), we obtain $f(y)\leq f(x)$ for any $x\in\mathcal{I}$, showing that $y$ is a minimizer of $f$. ∎

We use the notation $r_{\lambda,p}(X)=r_{\lambda}(X)$ for the risk, showing that it depends on $p$. Fix an arbitrary sequence of $X$ matrices on the event having probability one where $r_{\lambda,p}(X)\to R_{\lambda}$. By an argument similar to the one given above, the sequence of functions $f_{p}(\lambda)=r_{\lambda,p}(X)$ equicontinuous in $\lambda$ on the set $\lambda\geq 0$. Since $f_{p}(\lambda)$ converges to $R_{\lambda}$ for each $\lambda>0$, and $\lambda_{p}^{*}$ is a sequence of minimizers of $f_{p}(\lambda)$ that converges to $\lambda^{*}$, by Lemma 6.1 $\lambda^{*}$ is a minimizer of $R_{\lambda}$.

2 The Risk of Ridge Regression with Identity Covariance

From Theorem 2.1 it follows that the risk equals the limit $r_{\lambda}(X)$ (14). Since $\Sigma=I$ this simplifies to

By the Marchenko-Pastur theorem, Eq. (3), it follows that $p^{-1}\operatorname{tr}\left(\left(\widehat{\Sigma}+\lambda I_{p\times p}\right)^{-1}\right)$ $\to m_{I}(-\lambda;\gamma)$ is defined in (7).

In the proof of Lemma 7.4 in Section 7.1.4, it is shown that $p^{-1}\operatorname{tr}\left(\left(\widehat{\Sigma}+\lambda I_{p\times p}\right)^{-2}\right)$ $\to-\kappa^{\prime}(\lambda)$. For an identity covariance matrix $\kappa(\lambda)=m_{I}(-\lambda;\gamma)$ by definition of $\kappa(\lambda)$. Therefore, the limit of the second term equals $m^{\prime}_{I}(-\lambda;\gamma)$. We obtain the desired formula: $R_{\lambda}=1+\gamma m_{I}(-\lambda;\gamma)+\lambda\left(\lambda\alpha^{2}-\gamma\right)m^{\prime}_{I}(-\lambda;\gamma)$. For $\lambda^{*}=\gamma\alpha^{-2}$, we obtain $R^{*}=1+\gamma m_{I}(-\lambda^{*};\gamma)$. It is a matter of simple algebra to verify the formula (8) for the risk.

3 Proof of strong-signal limit of ridge

We will first show the results for the strong-signal limit. We start by verifying the following lemma.

Suppose the limit population eigenvalue distribution $H$ has support contained in a compact set bounded away from 0. Let $v(z)$ be the companion Stieltjes transform of the ESD. Then

If $\gamma<1$, $\lim_{\lambda\downarrow 0}\lambda v(-\lambda)=1-\gamma$.

If $\gamma>1$, $\lim_{\lambda\downarrow 0}v(-\lambda)=v(0)$.

Let $\underline{F}$ be the ESD of the companion matrix $n^{-1}XX^{\top}$. It is related to $F$ via $\underline{F}=(1-\gamma)\delta_{0}+\gamma F$. It is well known (e.g., Bai and Silverstein, 2010, Chapter 6), that for $H$ whose support is contained in a compact set bounded away from 0, the following hold for $F$ and $\underline{F}$: if $\gamma<1$, then $F$ has support contained in a compact set bounded away from 0, and $\underline{F}$ has a point mass of $1-\gamma$ at 0; while if $\gamma>1$, then $\underline{F}$ has support contained in a compact set bounded away from 0, and $F$ has a point mass of $1-\gamma^{-1}$ at 0.

If $\gamma<1$, we let $Y$ be distributed according to $F$. Since $Y>c>0$ for some $c$, we have by the dominated convergence theorem

Since $\lambda v(-\lambda)=1-\gamma+\gamma\lambda m(-\lambda)$, this shows $\lim_{\lambda\downarrow 0}\lambda v(-\lambda)=1-\gamma$.

Finally, for $\gamma=1$, the Silverstein equation (12) is equivalent to

Consequently, using the formula for the optimal risk from Theorem 2.1, we have for $\gamma<1$, $\lim_{\alpha^{2}\to\infty}R^{*}(H,\alpha^{2},\gamma)=(1-\gamma)^{-1}$. For $\gamma>1$,

The explicit formula for $R^{*}(\alpha^{2},1)$ is obtained by plugging in the expression (7) into the formula for the optimal risk.

Proofs for Regularized Discriminant Analysis

We will first outline the high-level steps to prove our main result for classification, Theorem 3.1. We break down the proof into several lemmas, whose proof is deferred to later sections. These lemmas are then put together to prove the theorem in the final part of the proof outline.

We start with the well-known finite-sample formula for the expected test error of an arbitrary linear classifier $h_{w,b}(x)=\operatorname{sign}(w\cdot x+b)$ in the Gaussian model (2), conditional on the weight parameters $w,b$ and the means $\mu_{\pm 1}$:

In RDA the weight vector is $\hat{w}_{\lambda}=\left(\widehat{\Sigma}_{c}+\lambda I_{p\times p}\right)^{-1}\hat{\delta}\,$ and the offset is $\hat{b}=\hat{\delta}^{\top}\left(\widehat{\Sigma}_{c}+\lambda I_{p\times p}\right)^{-1}\hat{\mu}$. The first simplification we notice is that $\hat{b}\to_{a.s.}0$.

Under the conditions of Theorem 3.1, we have $\hat{b}\to_{a.s.}0$.

Lemma 7.1 is proved in Section 7.1.1. Since the denominator in the error rate (16) converges almost surely to a fixed, strictly positive constant (see Lemmas 7.4 and 7.5), this will allow us to use the following simpler formula - that does not involve $\hat{b}$ - in evaluating the limit of the error rate.

Recall that $\mu_{-1}=\bar{\mu}-\delta$, $\mu_{+1}=\bar{\mu}+\delta$. The second simplification we notice is that $\hat{w}_{\lambda}^{\top}\bar{\mu}\to_{a.s.}0$.

Under the conditions of Theorem 3.1, we have $\hat{w}_{\lambda}^{\top}\bar{\mu}\to_{a.s.}0$.

Lemma 7.2 is proved in Section 7.1.2. By the same argument as above, this Lemma allows us to use the following even simpler formula - that does not involve $\bar{\mu}$ - in evaluating the limit of the error rate:

To show the convergence of $\smash{\Phi\left(-\hat{w}_{\lambda}^{\top}\delta/\sqrt{\hat{w}_{\lambda}^{\top}\Sigma\hat{w}_{\lambda}}\right)}$, we argue that the linear and quadratic forms involving $\hat{w}_{\lambda}$ concentrate around their means, and then apply random matrix results to find the limits of those means. We start with the numerator.

We have the limit $\hat{w}_{\lambda}^{\top}\delta\rightarrow_{a.s.}\alpha^{2}m(-\lambda)$, where $m(z)$ is the Stieltjes transform of the limit empirical eigenvalue distribution $F$ of the covariance matrix $\widehat{\Sigma}_{c}$.

Lemma 7.3 is proved in Section 7.1.3. To prove the convergence of the denominator, we decompose it as:

where $M:=\left(\widehat{\Sigma}_{c}+\lambda I_{p\times p}\right)^{-1}\Sigma\left(\widehat{\Sigma}_{c}+\lambda I_{p\times p}\right)^{-1}$ and

Lemma 7.4 is proved in Section 7.1.4, using Ledoit and Péché (2011)’s result and a derivative trick similar to that employed in a similar context by El Karoui and Kösters (2011); Rubio et al. (2012); Zhang et al. (2013). Finally, the last statement that we need is:

where $v$ the companion Stieltjes transform of the ESD of the covariance matrix, defined in (4).

Lemma 7.5 is proved in Section 7.1.5, as an application of the results of Hachem et al. (2008) and Chen et al. (2011). With all these results, we can now prove Theorem 3.1.

By the decomposition (19) and Lemmas 7.4 and 7.5, we have the convergence

By Lemma 7.5, the second term is strictly positive. Therefore, combining with Lemma 7.3 and the continuous mapping theorem, we have

Denote by $\Theta$ the parameter on the right hand side. After algebraic simplification, we obtain that $\Theta$ has exactly the form stated in the theorem for the margin of RDA. To finish the proof, we show that the error rate is indeed determined by $\Theta$. From (20) and the continuous mapping theorem, recalling the error rate $\operatorname{Err}_{2}\left(w\right)$ from (18), we have $\operatorname{Err}_{2}\left(\hat{w}_{\lambda}\right)\to_{a.s.}\Phi(-\Theta).$

From Lemma 7.2 and the definition of the error rate $\operatorname{Err}_{1}\left(w\right)$ from (17), we can move from $\operatorname{Err}_{2}$ to $\operatorname{Err}_{1}$: $\operatorname{Err}_{2}\left(\hat{w}_{\lambda}\right)-\operatorname{Err}_{1}\left(\hat{w}_{\lambda}\right)\to_{a.s.}0.$

Finally, from Lemma 7.1 and the definition of the error rate $\operatorname{Err}_{0}\left(w\right)$ in Equation (16), we can discard the offset $\smash{\hat{b}}$, and move from $\operatorname{Err}_{1}$ to $\operatorname{Err}_{0}$: $\operatorname{Err}_{1}\left(\hat{w}_{\lambda}\right)-\operatorname{Err}_{0}\left(\hat{w}_{\lambda},\hat{b}\right)\to_{a.s.}0.$

The last three statements imply that $\smash{\operatorname{Err}_{0}\left(\hat{w}_{\lambda},\hat{b}\right)\to_{a.s.}\Phi(-\Theta)}$, which finishes the proof of Theorem 3.1. ∎

In the proofs of the lemmas we will use the following well-known statement repeatedly:

Lemma 7.6 requires a small proof, which is provided in Section 7.1.6. The rest of this section contains the proofs of the lemmas.

We start by conditioning on the random variables $\bar{\mu},\delta$, or equivalently on $\mu_{\pm 1}$. Conditional on $\bar{\mu},\delta$, we have that $\hat{\mu}_{+1}\sim\mathcal{N}(\mu_{+1},2\Sigma/n)$, independently of $\hat{\mu}_{-1}\sim\mathcal{N}(\mu_{-1},2\Sigma/n)$. Therefore, $\hat{\delta}\sim\mathcal{N}\left(\delta,\,\Sigma/n\right)$, and independently $\hat{\mu}\sim\mathcal{N}\left(\bar{\mu},\,\Sigma/n\right)$. Further - still conditionally on $\bar{\mu},\delta$ - it holds that $\widehat{\Sigma}_{c}$ and $\hat{\mu}_{\pm 1}$ are independent by Gaussianity. This shows that conditionally on $\bar{\mu},\delta$, the random variables $\widehat{\Sigma}_{c},\hat{\mu},\hat{\delta}$ are independent.

Crucially, this representation has the same form regardless of the value of $\delta$, $\bar{\mu}$, therefore the random variables $Z,W$ are unconditionally independent of $\delta,\,\bar{\mu},\,\widehat{\Sigma}_{c}$. The unconditional indepdendence of $\delta,\bar{\mu},\widehat{\Sigma}_{c},Z,W$ will lead to convenient simplifications.

We decompose $\hat{b}$ according to (21) into the four terms that arise from expanding $\hat{\delta},\hat{\mu}$:

where $L=\left(\widehat{\Sigma}_{c}+\lambda I_{p\times p}\right)^{-1}$ and

The proof proceeds by showing that each of the $T_{i}$ converge to zero.

The first term: $T_{1}$. Let us denote $l=L\bar{\mu}$. Then by the independence and the zero-mean property of the coordinates of $\delta$

The second term: $T_{2}$. This term differs from $T_{1}$ because $n^{-1/2}Z$ replaces $\delta$. To show the convergence $T_{1}\to_{a.s.}0$ we only used the properties of the first four moments of $\delta$. The moments of $n^{-1/2}Z$ scale in the same way with $p$ as the moments of $\delta$. Therefore, the same proof shows $T_{2}\to_{a.s.}0$.

The last two terms: $T_{3}$ and $T_{4}$. The convergence of these terms follows directly from a well-known lemma, which we cite from Couillet and Debbah (2011):

While this lemma was originally stated for complex vectors, it holds verbatim for real vectors as well. The lemma applied with $x_{n}=\delta$, $y_{n}=n^{-1/2}W$ and $A_{n}=L\Sigma^{1/2}$ shows convergence of $T_{3}\to_{a.s.}0$; and similarly it shows $T_{4}\to_{a.s.}0$. This finishes the proof of Lemma 7.1.

1.2 Proof of Lemma 7.2

This follows from the proof of Lemma 7.2 by noting $\hat{w}_{\lambda}^{\top}\bar{\mu}=\hat{\delta}^{\top}L\bar{\mu}=T_{1}+T_{2}$, where $L=\left(\widehat{\Sigma}_{c}+\lambda I_{p\times p}\right)^{-1}$, and with $T_{i}$ from Equation (25).

1.3 Proof of Lemma 7.3

We will analyze the two terms separately, and show that the second term converges to 0.

Let $\widehat{\Sigma}_{c}$ be the centered and rescaled covariance matrix used in RDA. Then

Under the conditions of Theorem 3.1, the limit ESD of $\widehat{\Sigma}_{c}$ equals, with probability 1, the limit ESD $F$ of sample covariance matrices $\widehat{\Sigma}=n^{-1}\Sigma^{1/2}V^{\top}V\Sigma^{1/2}$, where $V$ is $n\times p$ with i.i.d. entries of mean 0 and variance 1. Therefore centering the covariance does not change the limit.

Let $u_{i}$ be the centered data points, $u_{i}=x_{i}-\mu_{+1}$ for the positive training examples, and $u_{i}=x_{i}-\mu_{-1}$ for the negative training examples. The $u_{i}$ have mean 0 and covariance matrix $\Sigma$. Let further $\hat{\nu}_{+1}=\hat{\mu}_{+1}-\mu_{+1}$ be the centered mean of the positive training examples; and define $\hat{\nu}_{-1}=\hat{\mu}_{-1}-\mu_{-1}$ analogously. We observe that