Year: 2020

Abstract:
We study the following three fundamental problems about ridge regression: (1) what is the structure of the estimator? (2) how to correctly use cross-validation to choose the regularization parameter? and (3) how to accelerate computation without losing too much accuracy? We consider the three problems in a unified large-data linear model. We give a precise representation of ridge regression as a covariance matrix-dependent linear combination of the true parameter and the noise. We study the bias of K-fold cross-validation for choosing the regularization parameter, and propose a simple bias-correction. We analyze the accuracy of primal and dual sketching for ridge regression, showing they are surprisingly accurate. Our results are illustrated by simulations and by analyzing empirical data.

text: We study the following three fundamental problems about ridge regression: (1) what is the structure of the estimator? (2) how to correctly use cross-validation to choose the regularization parameter? and (3) how to accelerate computation without losing too much accuracy? We consider the three problems in a unified large-data linear model. We give a precise representation of ridge regression as a covariance matrix-dependent linear combination of the true parameter and the noise. We study the bias of K-fold cross-validation for choosing the regularization parameter, and propose a simple bias-correction. We analyze the accuracy of primal and dual sketching for ridge regression, showing they are surprisingly accurate. Our results are illustrated by simulations and by analyzing empirical data.

heading: 1 INTRODUCTION
text: Ridge or `2-regularized regression is a widely used method for prediction and estimation when the data dimension p is large compared to the number of datapoints n. This is especially so in problems with many good features, where sparsity assumptions may not be justified. A great deal is known about ridge regression. It is Bayes optimal for any quadratic loss in a Bayesian linear model where the parameters and noise are Gaussian. The asymptotic properties of ridge have been widely studied (e.g., Tulino & Verdú, 2004; Serdobolskii, 2007; Couillet & Debbah, 2011; Dicker, 2016; Dobriban & Wager, 2018, etc). For choosing the regularization parameter in practice, cross-validation (CV) is widely used. In addition, there is an exact shortcut (e.g., Hastie et al., 2009, p. 243), which has good consistency properties (Hastie et al., 2019). There is also a lot of work on fast approximate algorithms for ridge, e.g., using sketching methods (e.g., el Alaoui & Mahoney, 2015; Chen et al., 2015; Wang et al., 2018; Chowdhury et al., 2018, among others).
Here we seek to develop a deeper understanding of ridge regression, going beyond existing work in multiple aspects. We work in linear models under a popular asymptotic regime where n, p→∞ at the same rate (Marchenko & Pastur, 1967; Serdobolskii, 2007; Couillet & Debbah, 2011; Yao et al., 2015). In this framework, we develop a fundamental representation for ridge regression, which shows that it is well approximated by a linear scaling of the true parameters perturbed by noise. The scaling matrices are functions of the population-level covariance of the features. As a consequence, we derive formulas for the training error and bias-variance tradeoff of ridge.
Second, we study commonly used methods for choosing the regularization parameter. Inspired by the observation that CV has a bias for estimating the error rate (e.g., Hastie et al., 2009, p. 243), we study the bias of CV for selecting the regularization parameter. We discover a surprisingly simple form for the bias, and propose a downward scaling bias correction procedure. Third, we study the accuracy loss of a class of randomized sketching algorithms for ridge regression. These algorithms approximate the sample covariance matrix by sketching or random projection. We show they can be surprisingly accurate, e.g., they can sometimes cut computational cost in half, only incurring 5% extra error. Even more, they can sometimes improve the MSE if a suboptimal regularization parameter is originally used.
Our work leverages recent results from asymptotic random matrix theory and free probability theory. One challenge in our analysis is to find the limit of the trace tr (Σ1 + Σ−12 )
−1/p, where Σ1 and Σ2 are p × p independent sample covariance matrices of Gaussian random vectors. The calculation requires nontrivial aspects of freely additive convolutions (e.g., Voiculescu et al., 1992; Nica & Speicher, 2006).
Our work is connected to prior works on ridge regression in high-dimensional statistics (Serdobolskii, 2007) and wireless communications (Tulino & Verdú, 2004; Couillet & Debbah, 2011). Among other related works, El Karoui & Kösters (2011) discuss the implications of the geometric sensitivity of random matrix theory for ridge regression, without considering our problems. El Karoui (2018) and Dicker (2016) study ridge regression estimators, but focus only on the risk for identity covariance. Hastie et al. (2019) study “ridgeless” regression, where the regularization parameter tends to zero.
Sketching is an increasingly popular research topic, see Vempala (2005); Halko et al. (2011); Mahoney (2011); Woodruff (2014); Drineas & Mahoney (2017) and references therein. For sketched ridge regression, Zhang et al. (2013a;b) study the dual problem in a complementary finite-sample setting, and their results are hard to compare. Chen et al. (2015) propose an algorithm combining sparse embedding and the subsampled randomized Hadamard transform (SRHT), proving relative approximation bounds. Wang et al. (2017) study iterative sketching algorithms from an optimization point of view, for both the primal and the dual problems. Dobriban & Liu (2018) study sketching using asymptotic random matrix theory, but only for unregularized linear regression. Chowdhury et al. (2018) propose a data-dependent algorithm in light of the ridge leverage scores. Other related works include Sarlos (2006); Ailon & Chazelle (2006); Drineas et al. (2006; 2011); Dhillon et al. (2013); Ma et al. (2015); Raskutti & Mahoney (2016); Gonen et al. (2016); Thanei et al. (2017); Ahfock et al. (2017); Lopes et al. (2018); Huang (2018).
The structure of the paper is as follows: We state our results on representation, risk, and biasvariance tradeoff in Section 2. We study the bias of cross-validation for choosing the regularization parameter in Section 3. We study the accuracy of randomized primal and dual sketching for both orthogonal and Gaussian sketches in Section 4. We provide proofs and additional simulations in the Appendix. Code reproducing the experiments in the paper are available at https://github. com/liusf15/RidgeRegression.

heading: 2 RIDGE REGRESSION
text: We work in the usual linear regression model Y = Xβ + ε, where each row xi of X ∈ Rn×p is a datapoint in p dimensions, and so there are p features. The corresponding element yi of Y ∈ Rn is its continous response (or outcome). We assume mean zero uncorrelated noise, so Eε = 0, and Cov [ε] = σ2In. We estimate the coefficient β ∈ Rp by ridge regression, solving the optimization problem
β̂ = arg min β∈Rp
1 n ‖Y −Xβ‖22 + λ‖β‖22,
where λ > 0 is a regularization parameter. The solution has the closed form β̂ = ( X>X/n+ λIp )−1 X>Y/n. (1)
We work in a ”big data” asymptotic limit, where both the dimension p and the sample size n tend to infinity, and their aspect ratio converges to a constant, p/n → γ ∈ (0,∞). Our results can be interpreted for any n and p, using γ = p/n as an approximation.
We recall that the empirical spectral distribution (ESD) of a p×p symmetric matrix Σ is the distribution 1p ∑p i=1 δλi where λi, i = 1, . . . , p are the eigenvalues of Σ, and δx is the point mass at x. This is a standard notion in random matrix theory, see e.g., Marchenko & Pastur (1967); Tulino & Verdú (2004); Couillet & Debbah (2011); Yao et al. (2015). The ESD is a convenient tool to summarize all information obtainable from the eigenvalues of a matrix. For instance, the trace of Σ is proportional to the mean of the distribution, while the condition number is related to the range of the support. As is common, we will work in models where there is a sequence of covariance matrices Σ = Σp, and their ESDs converges in distribution to a limiting probability distribution. The results become simpler, because they depend only on the limit.
By extension, we say that the ESD of the n× p matrix X is the ESD of X>X/n. We will consider some very specific models for the data, assuming it is of the form X = UΣ1/2, where U has iid entries of zero mean and unit variance. This means that the datapoints, i.e., the rows of X , have the form xi = Σ1/2ui, i = 1, . . . , p, where ui have iid entries. Then Σ is the ”true” covariance matrix of the features, which is typically not observed. These types of models for the data are very common in random matrix theory, see the references mentioned above.
Under these models, it is possible to characterize precisely the deviations between the empirical covariance matrix Σ̂ = n−1X>X and the population covariance matrix Σ, dating back to the well known classical Marchenko-Pastur law for eigenvectors (Marchenko & Pastur, 1967), extended to more general models and made more precise, including results for eigenvectors (see e.g., Tulino & Verdú, 2004; Couillet & Debbah, 2011; Yao et al., 2015, and references therein). This has been used to study methods for estimating the true covariance matrix, with several applications (e.g., Paul & Aue, 2014; Bun et al., 2017). More recently, such models have been used to study high dimensional statistical learning problems, including classification and regression (e.g., Zollanvari & Genton, 2013; Dobriban & Wager, 2018). Our work falls in this line.
We start by finding a precise representation of the ridge estimator. For random vectors un, vn of growing dimension, we say un and vn are deterministic equivalents, if for any sequence of fixed (or random and independent of un, vn) vectors wn such that lim sup ‖wn‖2 <∞ almost surely, we have |w>n (un − vn)| → 0 almost surely. We denote this by un vn. Thus linear combinations of un are well approximated by those of vn. This is a somewhat non-standard definition, but it turns out that it is precisely the one we need to use prior results from random matrix theory such as from (Rubio & Mestre, 2011).
We extend scalar functions f : R→ R to matrices in the usual way by functional calculus, applying them to the eigenvalues and keeping the eigenvectors. If M = V ΛV > is a spectral decomposition of M , then we define f(M) := V f(Λ)V >, where f(Λ) is the diagonal matrix with entries f(Λii).
For a fixed design matrix X , we can write the estimator as
β̂ = (Σ̂ + λIp) −1Σ̂β + (Σ̂ + λIp)
−1X >ε
n .
However, for a random design, we can find a representation that depends on the true covariance Σ, which may be simpler when Σ is simple, e.g., when Σ = Ip is isotropic. Theorem 2.1 (Representation of ridge estimator). Suppose the data matrix has the form X = UΣ1/2, where U ∈ Rn×p has iid entries of zero mean, unit variance and finite 8 + c-th moment for some c > 0, and Σ = Σp ∈ Rp×p is a deterministic positive definite matrix. Suppose that n, p→∞ with p/n→ γ > 0. Suppose the ESD of the sequence of Σs converges in distribution to a probability measure with compact support bounded away from the origin. Suppose that the noise is Gaussian, and that β = βp is an arbitrary sequence of deterministic vectors, such that lim sup ‖β‖2 <∞. Then the ridge regression estimator is asymptotically equivalent to a random vector with the following representation:
β̂(λ) A(Σ, λ) · β +B(Σ, λ) · σ · Z p1/2 .
Here Z ∼ N (0, Ip) is a random vector that is stochastically dependent only on the noise ε, and A,B are deterministic matrices defined by applying the scalar functions below to Σ:
A(x, λ) = (cpx+ λ) −2(cp + c ′ p)x, B(x, λ) = (cpx+ λ) −1cpx.
Here cp := c(n, p,Σ, λ) is the unique positive solution of the fixed point equation
1− cp = cp n
tr [ Σ(cpΣ + λI) −1] . (2) This result gives a precise representation of the ridge regression estimator. It is a sum of two terms: the true coefficient vector β scaled by the matrixA(Σ, λ), and the noise vectorZ scaled by the matrix B(Σ, λ). The first term captures to what extent ridge regression recovers the ”signal”. Morever, the
noise term Z is directly coupled with the noise in the original regression problem, and thus also the estimator. The result would not hold for an independent noise vector Z.
However, the coefficients are not fully explicit, as they depend on the unknown population covariance matrix Σ, as well as on the fixed-point variable cp.
Some comments are in order:
1. Structure of the proof. The proof is quite non-elementary and relies on random matrix theory. Specifically, it uses the language of the recently developed ”calculus of deterministic equivalents” (Dobriban & Sheng, 2018), and results by (Rubio & Mestre, 2011). A general takeaway is that for n not much larger than p, the empirical covariance matrix Σ̂ is not a good estimator of the true covariance matrix Σ. However, the deviation of linear functionals of Σ̂, can be quantified. In particular, we have
(Σ̂ + λI)−1 (cpΣ + λI)−1,
in the sense that linear combinations of the entries of the two matrices are close (see the proof for more details). 2. Understanding the resolvent bias factor cp. Thus, cp can be viewed as a resolvent bias factor, which tells us by what factor Σ is multiplied when evaluating the resolvent (Σ̂ + λI)−1, and comparing it to its naive counterpart (Σ + λI)−1. It is known that cp is well defined, and this follows by a simple monotonicity argument, see Hachem et al. (2007); Rubio & Mestre (2011). Specifically, the left hand side of (2) is decreasing in cp, while the right hand size is increasing in Also c′p is the derivative of cp, when viewing it as a function of z := −λ. An explicit expression is provided in the proof in Section A.1, but is not crucial right now.
Here we discuss some implications of this representation.
For uncorrelated features, Σ = Ip, A,B reduce to multiplication by scalars. Hence, each coordinate of the ridge regression estimator is simply a scalar multiple of the corresponding coordinate of β. One can use this to find the bias in each individual coordinate.
Training error and optimal regularization parameter. This theorem has implications for understanding the training error, and optimal regularization parameter of ridge regression. As it stands, the theorem itself only characterizes the behavior og linear combinations of the coordinates of the estimator. Thus, it can be directly applied to study the bias Eβ̂(λ) − β of the estimator. However, it cannot directly be used to study the variance; as that would require understanding quadratic functionals of the estimator. This seems to require significant advances in random matrix theory, going beyond the results of Rubio & Mestre (2011). However, we show below that with additional assumptions on the structure of the parameter β, we can derive the MSE of the estimator in other ways.
We work in a random-effects model, where the p-dimensional regression parameter β is random, each coefficient has zero mean Eβi = 0, and is normalized so that Varβi = α2/p. This ensures that the signal strength E‖β‖2 = α2 is fixed for any p. The asymptotically optimal λ in this setting is always λ∗ = γσ2/α2 see e.g., Tulino & Verdú (2004); Dicker (2016); Dobriban & Wager (2018). The ridge regression estimator with λ = pσ2/(nα2) is the posterior mean of β, when β and ε are normal random variables.
For a distribution F , we define the quantities
θi(λ) =
∫ 1
(x+ λ)i dFγ(x),
(i = 1, 2, . . .). These are the moments of the resolvent and its derivatives (up to constants). We use the following loss functions: mean squared estimation error: M(β̂) = E‖β̂ − β‖22, and residual or training error: R(β̂) = E [‖]Y −Xβ̂‖22. Theorem 2.2 (MSE and training error of ridge). Suppose β has iid entries with Eβi = 0, Var [βi] = α2/p, i = 1, . . . , p and β is independent of X and ε. Suppose X is an arbitrary n × p matrix depending on n and p, and the ESD of X converges weakly to a deterministic distribution F as n, p → ∞ and p/n → γ. Then the asymptotic MSE and residual error of the ridge regression estimator β̂(λ) has the form
lim n→∞
M(β̂(λ)) = α2λ2θ2 + γσ 2[θ1 − λθ2], (3)
lim n→∞
R(β̂(λ)) = α2λ2[θ1 − λθ2] + σ2 [ 1− γ(1 + λθ1 − λ2θ2) ] , (4)
Bias-variance tradeoff. Building on this, we can also study the bias-variance tradeoff of ridge regression. Qualitatively, large λ leads to more regularization, and decreases the variance. However, it also increases the bias. Our theory allows us to find the explicit formulas for the bias and variance as a function of λ. See Figure 1 for a plot and Sec. A.3 for the details. As far as we know, this is one of the few examples of high-dimensional asymptotic problems where the precise form of the bias and variance can be evaluated.
Bias-variance tradeoff at optimal λ∗ = γσ2/α2. (see Figure 6) This can be viewed as the ”pure” effect of dimensionality on the problem, keeping all other parameters fixed, and has intriguing properties. The variance first increases, then decreases with γ. In the ”classical” low-dimensional case, most of the risk is due to variance, while in the ”modern” high-dimensional case, most of it is due to bias. This is consistent with other phenomena in proportional-limit asymptotics, e.g., that the map between population and sample eigenvalue distributions is asymptotically deterministic (Marchenko & Pastur, 1967).
Future applications. This fundamental representation may have applications to important statistical inference questions. For instance, inference on the regression coefficient β and the noise variance σ2 are important and challenging problems. Can we use our representation to develop debiasing techniques for this task? This will be interesting to explore in future work.

heading: 3 CROSS-VALIDATION
text: How can we choose the regularization parameter? In practice, cross-validation (CV) is the most popular approach. However, it is well known that CV has a bias for estimating the error rate, because it uses a smaller number of samples than the full data size (e.g., Hastie et al., 2009, p. 243). In this section, we study related questions, proposing a bias-correction method for the optimal regularization parameter. This is closely connected to the previous section, because it relies on the same random-effects theoretical framework. In fact, our conclusions here are a direct consequence of the properties of that framework.
Setup. Suppose we split the n datapoints (samples) into K equal-sized subsets, each containing n0 = n/K samples. We use the k-th subset (Xk, Yk) as the validation set and the other K − 1 subsets (X−k, Y−k), with total sample size n1 = (K − 1)n/K as the training set. We find the ridge
regression estimator β̂−k, i.e. β̂−k(λ) = ( X>−kX−k + n1λIp )−1 X>−kY−k.
The expected cross-validation error is, for isotropic covariance, i.e., Σ = I ,
CV (λ) = EĈV (λ) = E
[ 1
K K∑ k=1 ‖Yk −Xkβ̂−k(λ)‖22/n0
] = σ2 + E [ ‖β̂−k − β‖22 ] .
Bias in CV. When n, p tend to infinity so that p/n → γ > 0, and in the random effects model with Eβi = 0, Varβi = α2/p described above, the minimizer of CV (λ) tends to λ∗k = γ̃σ2/α2, where γ̃ is the limiting aspect ratio of X−k, i.e. γ̃ = γK/(K − 1). Since the aspect ratios of X−k and X differ, the limiting minimizer of the cross-validation estimator of the test error is biased for the limiting minimizer of the actual test error, which is λ∗ = γσ2/α2.
Bias-correction. Suppose we have found λ̂∗k, the minimizer of ĈV (λ). Afterwards, we usually refit ridge regression on the entire dataset, i.e., find
β̂(λ̂∗) = (X>X + λ̂∗nI)−1X>Y.
Based on our bias calculation, we propose to use a bias-corrected parameter
λ̂∗ := λ̂∗k K − 1 K .
So if we use 5 folds, we should multiply the CV-optimal λ by 0.8. We find it surprising that this theoretically justified bias-correction does not depend on any unknown parameters, such as β, α2, σ2.While the bias of CV is widely known, we are not aware that this bias-correction for the regularization parameter has been proposed before.
Numerical examples. Figure 2 shows on two empirical data examples that the debiased estimator gets closer to the optimal λ than the original minimizer of the CV. However, in this case it does not significantly improve the test error. Simulation results in Section A.4 also show that the bias-correction correctly shrinks the regularization parameter and decreases the test error. We also consider examples where p n (i.e., γ 1), because this is a setting where it is known that the bias of CV can be large (Tibshirani & Tibshirani, 2009). However, in this case, we do not see a significant improvement.
Extensions. The same bias-correction idea also applies to train-test validation. In addition, there is a special fast “short-cut” for leave-one-out cross-validation in ridge regression (e.g., Hastie et al.,
2009, p. 243), which has the same cost as one ridge regression. The minimizer converges to λ∗ (Hastie et al., 2019). However, we think that the bias-correction idea is still valuable, as the idea applies beyond ridge regression: CV selects regularization parameters that are too large. See Section A.5 for more details and experiments comparing different ways of choosing the regularization parameter.

heading: 4 SKETCHING
text: A final important question about ridge regression is how to compute it in practice. In this section, we study that problem in the same high-dimensional model used throughout our paper. The computation complexity of ridge regression, O(np min(n, p)), can be intractable in modern large-scale data analysis. Sketching is a popular approach to reducing the time complexity by reducing the sample size and/or dimension, usually by random projection or sampling (e.g. Mahoney, 2011; Woodruff, 2014; Drineas & Mahoney, 2016). Specifically, primal sketching approximates the sample covariance matrix X>X/n by X>L>LX/n, where L is an m× n sketching matrix, and m < n. If L is chosen as a suitable random matrix, then this can still approximate the original sample covariance matrix. Then the primal sketched ridge regression estimator is
β̂p = ( X>L>LX/n+ λIp )−1 X>Y/n. (5)
Dual sketching reduces p instead. An equivalent expression for ridge regression is β̂ = n−1X> ( XX>/n+ λIn )−1 Y . Dual sketched ridge regression reduces the computation cost of the Gram matrix XX>, approximating it by XRR>X> for another sketching matrix R ∈ Rp×d (d < p), so
β̂d = X > (XRR>X>/n+ λIn)−1 Y/n. (6)
The sketching matrices R and L are usually chosen as random matrices with iid entries (e.g., Gaussian ones) or as orthogonal matrices. In this section, we study the asymptotic MSE for both orthogonal (Section 4.1) and Gaussian sketching (Section 4.2). We also mention full sketching, which performs ridge after projecting down both X and Y . In section A.11, we find its MSE. However, the other two methods have better tradeoffs, and we can empirically get better results for the same computational cost.

heading: 4.1 ORTHOGONAL SKETCHING
text: First we consider primal sketching with orthogonal projections. These can be implemented by subsampling, Haar distributed matrices, or subsampled randomized Hadamard transforms (Sarlos, 2006). We recall that the standard Marchenko-Pastur (MP) law is the probability distribution which is the limit of the ESD of X>X/n, when the n× p matrix X has iid standard Gaussian entries, and n, p → ∞ so that p/n → γ > 0, which has an explicit density (Marchenko & Pastur, 1967; Bai & Silverstein, 2010). Theorem 4.1 (Primal orthogonal sketching). Suppose β has iid entries with Eβi = 0, Var [βi] = α2/p, i = 1, . . . , p and β is independent of X and ε. Suppose X has iid standard normal entries.
We compute primal sketched ridge regression (5) with an m × n orthogonal matrix L (m < n, LL> = Im). Let n, p and m tend to infinity with p/n→ γ ∈ (0,∞) and m/n→ ξ ∈ (0, 1). Then the MSE of β̂p(λ) has the limit
M(λ) = α2
[ (λ+ ξ − 1)2 + γ(1− ξ) ] θ2 ( γ ξ , λ ξ ) ξ2 + γσ2 ξθ1 ( γ ξ , λ ξ ) − (λ+ ξ − 1)θ2 ( γ ξ , λ ξ ) ξ2 ,
(7) where θi(γ, λ) = ∫
(x + λ)−idFγ(x) and Fγ is the standard Marchenko-Pastur law with aspect ratio γ.
Structure of the proof. The proof is in Section A.6, with explicit formulas in Section A.6.1. The θi are related to the resolvent of the MP law and its derivatives. In the proof, we decompose the
MSE as the sum of variance and squared bias, both of which further reduce to the traces of certain random matrices, whose limits are determined by the MP law Fγ and λ. The two terms on the RHS of Equation (7) are the limits of squared bias and variance, respectively. There is an additional key step in the proof, which introduces the orthogonal complement L1 of the matrix L such that L>L + L>1 L1 = In, which leads to some Gaussian random variables appearing in the proof, and simplifies calculations.
Simulations. A simulation in Figure 3 (left) shows a good match with our theory. It also shows that sketching does not increase the MSE too much. In this case, by reducing the sample size to half the original one, we only increase the MSE by a factor of 1.05. This shows sketching can be very effective. We also see in Figure 3 (right) that variance is compromised much more than bias.
Robustness to tuning parameter. The reader may wonder how strongly this depends on the choice of the regularization parameter λ. Perhaps ridge regression works poorly with this λ, so sketching cannot worsen it too much? What happens if we take the optimal λ instead of a fixed one? In experiments in Section A.12 we show that the behavior is quite robust to the choice of regularization parameter.
The next theorem states a result for dual orthogonal sketching.
Theorem 4.2 (Dual orthogonal sketching). Under the conditions of Theorem 4.1, we compute the dual sketched ridge regression with an orthogonal p × d sketching matrix R (d 6 p, R>R = Id). Let n, p and d go to infinity with p/n→ γ ∈ (0,∞) and d/n→ ζ ∈ (0, γ). Then the MSE of β̂d(λ) has the limit
α2
γ
[ γ − 1 + (λ− γ + ζ)2θ̄2(ζ, λ) + (γ − ζ)θ̄21(ζ, λ) ] + σ2 [ θ̄1(ζ, λ)− (λ+ ζ − γ)θ̄2(ζ, λ) ] ,
where θ̄i(ζ, λ) = (1− ζ)/λi + ζ ∫ (x+λ)−idFζ(x), and Fζ is the standard Marchenko-Pastur law.
Proof structure and simulations. The proof in Section A.7 follows similar path to the previous one. Here θ̄i comes in because of the companion Stieltjes transform of MP law. The simulation results shown in Figure 11 agrees well with our theory. They are similar to the ones before: sketching has favorable properties, and the bias increases less than the variance.
Optimal tuning parameters. For both primal and dual sketching, the optimal regularization parameter minimizing the MSE seems analytically intractable. Instead, we use a numerical approach in our experiments, based on a binary search. Since this is one-dimensional problem, there are no numerical issues. See Figure 13 in Section A.12.3.

heading: 4.1.1 EXTREME PROJECTION — MARGINAL REGRESSION
text: It is of special interest to investigate extreme projections, where the sketching dimension is much reduced compared to the sample size, so m n. This corresponds to ξ = 0. This can also be viewed as a scaled marginal regression estimator, i.e., β̂ ∝ X>Y . For dual sketching, the same case can be recovered with ζ = 0. Another interest of studying this special case is that the formula for MSE simplifies a lot.
Theorem 4.3 (Marginal regression). Under the same assumption as Theorem 4.1, let ξ = 0. Then the form of the MSE is M(λ) = [α2 [ (λ− 1)2 + γ ] + σ2γ]/λ2. Moreover, the optimal λ∗ that minimizes this equals γσ2/α2 + 1 + γ and the optimal MSE is M(λ∗) = α2 ( 1− α2/[α2(1 + γ) + γσ2] ) .
The proof is in Section A.8. When is the optimal MSE of marginal regression small? Compared to the MSE of the zero estimator α2, it is small when γ(σ2/α2 + 1) + 1 is large. In Figure 4 (left), we compare marginal and ridge regression for different aspect ratios and SNR. When the signal to noise ratio (SNR) α2/σ2 is small or the aspect ratio γ is large, marginal regression does not increase the MSE much. As a concrete example, if we take α2 = σ2 = 1 and γ = 0.7, the marginal MSE is 1− 1/2.4 ≈ 0.58. The optimal ridge MSE is about 0.52, so their ratio is only ca. 0.58/0.52 ≈ 1.1. It seems quite surprising that a simple-minded method like marginal regression can work so well. However, the reason is that when the SNR is small, we cannot expect ridge regression to have good performance. Large γ can also be interpreted as small SNR, where ridge regression works poorly and sketching does not harm performance too much.

heading: 4.2 GAUSSIAN SKETCHING
text: In this section, we study Gaussian sketching. The following theorem states the bias of dual Gaussian sketching. The bias is enough to characterize the performance in the high SNR regime where α/σ → ∞, and we discuss the extension to low SNR after the proof. Theorem 4.4 (Bias of dual Gaussian sketch). Suppose X is an n × p standard Gaussian random matrix. Suppose also that R is a p × d matrix with i.i.d. N (0, 1/d) entries. Then the bias of dual sketch has the expression Bias2(β̂d) = α2 + α2/γ · [m′(z)− 2m(z)] |z=0, where m is a function described below, and m′(z) denotes the derivative of m w.r.t. z. Below, we use the branch of the square root with positive imaginary part.
The function m is characterized by its inverse function, which has the explicit formula m−1(z) = 1/[1 + z/ζ]− [γ + 1− √ (γ − 1)2 + 4λz]/(2z) for complex z with positive imaginary part.
About the proof. The proof is in Section A.9.We mention that the same result holds when the matrices involved have iid non-Gaussian entries, but the proof is more technical. The current proof is based on free probability theory (e.g., Voiculescu et al., 1992; Hiai & Petz, 2006; Couillet & Debbah, 2011). The function m is the Stieltjes transform of the free additive convolution of a standard MP law F1/ξ and a scaled inverse MP law λ/γ · F−11/γ (see the proof).
Numerics. To evaluate the formula, we note that m−1(m(0)) = 0, so m(0) is a root of m−1. Also, dm(0)/dz equals 1/(dm−1(y)/dy|y=m(0)), the reciprocal of the derivative of m−1 evaluated at m(0). We use binary search to find the numerical solution. The theoretical result agrees with the simulation quite well, see Figure 4.
Somewhat unexpectedly, the MSE of dual sketching can be below the MSE of ridge regression, see Figure 4. This can happen when the original regularization parameter is suboptimal. As d grows, the MSE of Gaussian dual sketching converges to that of ridge regression.
We have also found the bias of primal Gaussian sketching. However, stating the result requires free probability theory, and so we present it in the Appendix, see Theorem A.1. To further validate our results, we present additional simulations in Sec. A.12, for both fixed and optimal regularization parameters after sketching. A detailed study of the computational cost for sketching in Sec. A.13 concludes, as expected, that primal sketching can reduce cost when p < n, while dual sketching can reduce it when p > n; and also provides a more detailed analysis.

heading: ACKNOWLEDGMENTS
text: The authors thank Ken Clarkson for helpful discussions and for providing the reference Chen et al. (2015). ED was partially supported by NSF BIGDATA grant IIS 1837992. SL was partially supported by a Tsinghua University Summer Research award. A version of our manuscript is available on arxiv at https://arxiv.org/abs/1910.02373.

heading: A APPENDIX

heading: A.1 PROOF OF THEOREM 2.1
text: If p/n → γ and the spectral distribution of Σ converges to H , we have by the general MarchenkoPastur (MP) theorem of Rubio and Mestre (Rubio & Mestre, 2011), that
(Σ̂ + λI)−1 (cpΣ + λI)−1,
where cp := c(n, p,Σ, λ) is the unique positive solution of the fixed point equation
1− cp = cp n
tr [ Σ(cpΣ + λI) −1] . Here, using the terminology of the calculus of deterministic equivalents (Dobriban & Sheng, 2018), two sequences of (not necessarily symmetric) n × n matrices An, Bn of growing dimensions are equivalent, and we write
An Bn if limn→∞ tr [Cn(An −Bn)] = 0 almost surely, for any sequence Cn of (not necessarily symmetric) n × n deterministic matrices with bounded trace norm, i.e., such that lim sup ‖Cn‖tr < ∞ (Dobriban & Sheng, 2018). Informally, linear combinations of the entries of An can be approximated by the entries of Bn.
We start with
β̂ = ( X>X/n+ λIp )−1 X>Y/n = ( X>X/n+ λIp )−1 X>(Xβ + ε) n
= (Σ̂ + λIp) −1Σ̂β + (Σ̂ + λIp)
−1X >ε
n .
Then, by the general MP law written in the language of the calculus of deterministic equivalents
(Σ̂ + λIp) −1Σ̂ = Ip − λ(Σ̂ + λIp)−1 Ip − λ(cpΣ + λI)−1 = cpΣ(cpΣ + λI)−1.
By the definition of equivalence for vectors,
(Σ̂ + λIp) −1Σ̂β cpΣ(cpΣ + λI)−1β.
We note a subtle point here. The rank of the matrix M := (Σ̂ + λIp)−1Σ̂ is at most n, and so it is not a full rank matrix when n < p. In contrast, cpΣ(cpΣ + λI)−1 can be a full rank matrix. Therefore, for the vectors β in the null space of Σ̂, which is also the null space of X , we certainly have that the two sides are not equal. However, here we assumed that the matrix X is random, and so its null space is a random max(p−n, 0) dimensional linear space. Therefore, for any fixed vector β, the random matrix M will not contain it in its null space with high probability, and so there is no contradiction.
We should also derive an asymptotic equivalent for
(Σ̂ + λIp) −1X
>ε n .
Suppose we have Gaussian noise, and let Z ∼ N (0, Ip). Then we can write
(Σ̂ + λIp) −1X
>ε n =d (Σ̂ + λIp) −1Σ̂1/2 σZ n1/2 .
So the question reduces to finding a deterministic equivalent for h(Σ̂), where h(x) = (x + λ)−2x. Note that
h(x) = (x+ λ)−2x = (x+ λ)−2(x+ λ− λ) = (x+ λ)−1 − λ(x+ λ)−2.
By the calculus of determinstic equivalents: (Σ̂ + λ)−1 (cpΣ + λI)−1. Moreover, fortunately the limit of the second part was recently calculated in (Dobriban & Sheng, 2019). This used the so-called ”differentiation rule” of the calculus of deterministic equivalents to find
(Σ̂ + λ)−2 (cpΣ + λI)−2(I − c′pΣ).
The derivative c′p = dcp/dz has been found in Dobriban & Sheng (2019), in the proof of Theorem 3.1, part 2b. The result is (with γp = p/n,Hp the spectral distribution of Σ, and T a random variable distributed according to Hp)
c′p = γpEHp cpT (cpT−z)2
−1 + γpzEHp T(cpT−z)2 . (8)
So, we find the final answer
(Σ̂ + λIp) −1Σ̂1/2 A(Σ, λ) := (cpΣ + λI)−1 − λ(cpΣ + λI)−2(I − c′pΣ).

heading: A.2 RISK ANALYSIS
text: Figure 5 shows a simulation result. We see a good match between theory and simulation.

heading: A.2.1 PROOF OF THEOREM 2.2
text: Proof. The MSE of β̂ has the form
E‖β̂ − β‖2 = bias2 + δ2,
where
bias2 = E [∥∥∥(X>X/n+ λIp)−1X>X/nβ − β∥∥∥2
2
] ,
δ2 = σ2E [∥∥∥(X>X/n+ λIp)−1 n−1X>∥∥∥2
F
] .
We assume that X has iid entries of zero mean and unit variance, and that Eβ = 0, Var [β] = α2/pIp. As p/n→ γ as n goes to infinity, the ESD of 1nX
>X converges to the MP law Fγ . So we have
bias2 = E [∥∥∥λ (X>X/n+ λIp)−1 β∥∥∥2
2 ] = α2λ2E [ 1
p tr[ ( X>X/n+ λIp )−2 ] ] → α2λ2 ∫ 1 (x+ λ)2 dFγ(x),
and
δ2 = σ2 n2 E [ tr[ ( X>X/n+ λIp )−2 X>X] ] = σ2 n E [ tr[ ( X>X/n+ λIp
)−1 − λ (X>X/n+ λIp)−2]] → σ2γ [∫ 1
x+ λ dFγ(x)− λ
∫ 1
(x+ λ)2 dFγ(x)
] .
Denoting θi(γ, λ) = ∫ 1 (x+λ)i dFγ(x), then
AMSE(β̂) = α2λ2θ2 + γσ 2[θ1 − λθ2]. (9)
For the standard Marchenko-Pastur law (i.e., when Σ = Ip), we have the explicit forms of θ1 and θ2. Specifically,
θ1 =
∫ 1
x+ λ dFγ(x) = −
1
2
[ 2(1 + λ)
λγ + 2 √ γλ z2 ] where
z2 = − 1
2
[ ( √ γ +
1 + λ √ γ ) +
√ ( √ γ +
1 + λ √ γ )2 − 4
] .
It is known that the limiting Stieltjes transform mFγ := mγ of Σ̂ has the explicit form (Marchenko & Pastur, 1967):
mγ(z) = (z + γ − 1) + √ (z + γ − 1)2 − 4zγ −2zγ .
As usual in the area, we use the principal branch of the square root of complex numbers. Hence
θ1 = (−λ+γ−1)+ √ (−λ+γ−1)2+4λγ 2λγ . Also
θ2(γ, λ) =
∫ 1
(x+ λ)2 dFγ(x) = −
∫ d
dλ
1
x+ λ dFγ(x)
= − d dλ θ1 = − 1 γλ2 + 1 √ γ d dλ z2 λ
= − 1 γλ2 + γ + 1 2γλ2 − 1 2 √ γ [ λ+ γ + 1 γλ √
( √ γ + 1+λ√γ )
2 − 4 −
√ ( √ γ + 1+λ√γ ) 2 − 4
λ2 ]
For the residual, E [ 1
n ‖Y −Xβ̂‖22|X
] = α2λ2 1 p tr[ ( X>X/n+ λIp )−1 − λ (X>X/n+ λIp)−2] + σ2 1
n [tr(In)− 2 tr
( X>X/n+ λIp )−1 X>X/n+ tr (( X>X/n+ λIp )−1 X>X/n )2 ].
Next, E [ 1 p tr[ (( X>X/n+ λIp )−1 X>X/n )2 ] ] = E [ 1 p tr[ ( Ip − λ ( X>X/n+ λIp )−1)2 ] ] → 1− 2λθ1 + λ2θ2.
Therefore E [ 1
n ‖Y −Xβ̂‖22
] →α2λ2[θ1 − λθ2] + σ2 [ 1− 2γ(1− λθ1) + γ(1− 2λθ1 + λ2θ2) ] = α2λ2[θ1 − λθ2] + σ2 [ 1− γ(1 + λθ1 − λ2θ2) ] .

heading: A.3 BIAS-VARIANCE TRADEOFF
text: The limiting MSE decomposes into a limiting squared bias and variance. The specific forms of these are
bias2 = α2 ∫
λ2
(x+ λ)2 dFγ(x), var = γσ
2
∫ x
(x+ λ)2 dFγ(x).
See Figure 1 for a plot. We can make several observations.
1. The bias increases with λ, starting out at zero for λ = 0 (linear regression), and increasing to α2 as λ→∞ (zero estimator).
2. The variance decreases with λ, from γσ2 ∫ x−1dFγ(x) to zero.
3. In the setting plotted in the figure, when α2 and σ2 are roughly comparable, there are additional qualitative properties we can investigate. When γ is small, the regularization parameter λ influences the bias more strongly than the variance (i.e., the derivative of the normalized quantities in the range plotted is generally larger for the normalized squared bias). In contrast when γ is large, the variance is influenced more.
Next we consider how bias and variance change with γ at the optimal λ∗ = γσ2/α2. This can be viewed as the ”pure” effects of dimensionality on the problem, keeping all other parameters fixed. Ineed, α2/σ2 can be viewed as the signal-to-noise ratio (SNR), and is fixed. This analysis allows us to study for the best possible estimator (ridge regression, a Bayes estimator), behaves with the dimension. We refer to Figure 6, where we make some specific choices of α and σ.
1. Clearly the overall risk increases, as the problem becomes harder with increasing dimension. This is in line with our intuition.
2. The classical bias-variance tradeoff can be summarized by the equation
bias2(λ) + var(λ) >M∗(α, γ),
where we made explicit the dependence of the bias and variance on λ, and whereM∗(α, γ) is the minimum MSE achievable, also known as the Bayes error, for which there are explicit formulas available (Tulino & Verdú, 2004; Dobriban & Wager, 2018).
3. The variance first increases, then decreases with γ. This shows that in the ”classical” low-dimensional case, most of the risk is due to variance, while in the ”modern” highdimensional case, most of it is due to bias. This observation is consistent with other phenomena in proportional-limit asymptotics, for instance that the map between population and sample eigenvalue distributions is asymptotically deterministic (Marchenko & Pastur, 1967; Bai & Silverstein, 2010).

heading: A.4 SIMULATIONS WITH CROSS-VALIDATION
text: See Figure 7. We consider both small and large γ. Our bias-correction procedure shrinks the λ to the correct direction and decreases the test error. It is also shown that the one-standard-error rule (e.g., Hastie et al., 2009) does not perform well here.
K λ
∗
CV . So the bias-correction decreases the test error
by about 0.003. Right: we take n = 200, p = 1000, γ = 5, α = 3, σ = 1. The bias-correction decreases the test error from 8.92 to 8.89, so it decreases by 0.03.

heading: A.5 CHOOSING THE REGULARIZATION PARAMETER- ADDITIONAL DETAILS
text: Another possible prediction method is to use the average of the ridge estimators computed during cross-validation. Here it is also natural to use the CV-optimal regularization parameters, averaging β̂−k(λ̂ ∗ k), i.e.
β̂avg(λ̂ ∗ k) =
1
K K∑ k=1 β̂−k(λ̂ ∗ k).
This has the advantage that it does not require refitting the ridge regression estimator, and also that we use the optimal regularization parameter.

heading: A.5.1 TRAIN-TEST VALIDATION
text: The same bias in the regularization parameter also applies to train-test validation. Since the number of samples is changed when restricting to the training set, the optimal λ chosen by train-test validation is also biased for the true regularization parameter minimizing the test error. We will later see in simulations (Figure 8) that retraining the ridge regression estimator on the whole data will still significantly improve the performance (this is expected based on our results on CV). For prediction, here we can also use ridge regression on the training set. This effectively reduces sample size n → ntrain, where ntrain is the sample size of the training set. However, if the training set grows such that n/ntrain → 1 while ntrain → ∞, the train-test split has asymptotically optimal performance.

heading: A.5.2 LEAVE-ONE-OUT
text: There is a special “short-cut” for leave-one-out in ridge regression, which saves us from burdensome computation. Write loo(λ) for the leave-one-out estimator of prediction error with parameter λ. Instead of doing ridge regression n times, we can calculate the error explicitly as
loo(λ) = 1
n n∑ i=1
[ Yi −X>i β̂(λ)
1− Sii(λ)
]2 .
where S(λ) = X(X>X + nλI)−1X>. The minimizer of loo(λ) is asymptotically optimal, i.e., it converges to λ∗ (Hastie et al., 2019). However, the computational cost of this shortcut is the same as that of a train-test split. Therefore, the method described above has the same asymptotic performance.
Simulations: Figure 8 shows simulation results comparing different cross-validation methods:
1. kf — k-fold cross-validation by taking the average of the ridge estimators at the CV-optimal regularization parameter.
2. kf refit — k-fold cross-validation by refitting ridge regression on the whole dataset using the CV-optimal regularization parameter. 3. kf bic — k-fold cross-validation by refitting ridge regression on the whole dataset using the CV-optimal regularization parameter, with bias correction. 4. tt — train-test validation, by using the ridge estimator computed on the train data, at the validation-optimal regularization parameter. Note: we expect this to be similar, but worse than the ”kf” estimator. 5. tt refit — train-test validation by refitting ridge regression on the whole dataset, using the validation-optimal regularization parameter. Note: we expect this to be similar, but slightly worse than the ”kf refit” estimator. 6. tt bic — train-test validation by refitting ridge regression on the whole dataset using the CV-optimal regularization parameter, with bias correction. 7. loo — leave-one-out
Figure 8 shows that the naive estimators (kf and tt) can be quite inaccurate without refitting or bias correction. However, if we either refit or bias-correct, the accuracy improves. In this case, there seems to be no significant difference between the various methods.

heading: A.6 PROOF OF THEOREM 4.1
text: Proof. Suppose m/n→ ξ as n goes to infinity. For β̂p, we have bias2 = E [∥∥∥(X>L>LX/n+ λIp)−1X>X/nβ − β∥∥∥2
2
] ,
δ2 = σ2E [∥∥∥(X>L>LX/n+ λIp)−1 n−1X>∥∥∥2
F
] .
Denote M = ( X>L>LX/n+ λIp )−1 , the resolvent of the sketched matrix. We further assume that X has iid N (0, 1) entries and LL> = Im. Let L1 be an orthogonal complementary matrix of L, such that L>L+ L>1 L1 = In. We also denote N = X>L>1 L1X n . Then
MX>X/n =M X>L>LX +X>L>1 L1X
n = Ip − λM +MN.
Therefore, using that Cov [β] = α2/p · Ip, we find the bias as
bias2 = α2 p E [ tr(M − Ip)(M> − Ip) ] = α2
p
{ λ2E [ tr[M2] ] + E [ trM2 (X>L>1 L1X) 2
n2
] − 2λE [ trM2N ]} .
By the properties of Wishart matrices (e.g., Anderson, 2003; Muirhead, 2009), we have
E [N ] = n−m n Ip,
E [ (N)2 ] = 1 n2 E [ Wishart(Ip, n−m)2 ] = 1 n2 [n−m+ p(n−m) + (n−m)2]Ip.
Recalling that m,n→∞ such that m/n→ ξ, and that θi(γ, λ) = ∫ (x+ λ)−idFγ(x),
bias2 = α2
p
[ λ2 + n−m+ p(n−m) + (n−m)2
n2 − 2λn−m n
] E [ tr[M2] ] → α2[(λ+ ξ − 1)2 + γ(1− ξ)]θ2(γ, ξ, λ).
Moreover,
δ2 = σ2 n2 E [ tr[M2X>X] ] = σ2 n · { E [tr[M ]]− λE [ tr[M2] ] + E [ tr[M2N ]
]} → γσ2[θ1(γ, ξ, λ)− λθ2(γ, ξ, λ) + (1− ξ)θ2(γ, ξ, λ)].
Here we used the additional definitions
θi(γ, ξ, λ) =
∫ 1
(ξx+ λ)i dFγ/ξ(x)
θi(γ, λ) = θi(γ, ξ = 1, λ).
Note that these can be connected to the previous definitions by
θ1(γ, ξ, λ) = 1
ξ
∫ 1
x+ λ/ξ dFγ/ξ(x) =
1 ξ θ1
( γ
ξ , λ ξ ) θ2(γ, ξ, λ) = 1
ξ2 θ2
( γ
ξ , λ ξ
) .
Therefore the AMSE of β̂p is
AMSE(β̂p) = α 2[(λ+ ξ − 1)2 + γ(1− ξ)]θ2(γ, ξ, λ) + γσ2[θ1(γ, ξ, λ)− (λ+ ξ − 1)θ2(γ, ξ, λ)]
= α2[(λ+ ξ − 1)2 + γ(1− ξ)] 1 ξ2 θ2
( γ
ξ , λ ξ ) + γσ2 [ 1
ξ θ1
( γ
ξ , λ ξ
) − (λ+ ξ − 1) 1
ξ2 θ2
( γ
ξ , λ ξ
)] . (10)
A.6.1 ISOTROPIC CASE
Consider the special case where Γ = I , that is, X has iid N (0, 1) entries. Then Fγ is the standard MP law, and we have the explicit forms for θi = θi(γ, λ) =
∫ 1
(x+λ)i dFγ :
θ1(γ, λ) = − 1 + λ
γλ +
1
2 √ γλ
[ √ γ + 1 + λ √ γ +
√ ( √ γ +
1 + λ √ γ )2 − 4],
θ2(γ, λ) = − 1 γλ2 + γ + 1 2γλ2 − 1 2 √ γ ( λ+ 1 γ + 1)
1 λ √
( √ γ + 1+λ√γ ) 2 − 4 +
1
2 √ γ
√ ( √ γ +
1 + λ √ γ )2 − 4 1 λ2 ,
θ̄1(ζ, λ) = ζθ1(ζ, λ) + 1− ζ λ , θ̄2(ζ, λ) = ζθ2(ζ, λ) + 1− ζ λ2 ,
The results are obtained by the contour integral formula∫ f(x)dFγ(x) = − 1
4πi ∮ |z|=1 f(|1 + γz|2)(1− z2)2 z2(1 + √ γz)(z + √ γ) dz.
See Proposition 2.10 of Yao et al. (2015).

heading: A.7 PROOF OF THEOREM 4.2
text: Proof. Suppose d/p→ ζ as n goes to infinity. For β̂d, we have bias2 = E [∥∥∥n−1X> (XRR>X>/n+ λIn)−1Xβ − β∥∥∥2
2
] ,
δ2 = σ2 tr[ ( XRR>X>/n+ λIn )−2 XX> n2 ].
Denote M = ( XRR>X>/n+ λIn )−1 . Note that, using that Cov [β] = α2/p · Ip
bias2 = α2 p E [ tr[MXX>/n]2 ] − 2α 2 p E [ tr[MXX>/n] ] + α2 p tr(Ip).
Moreover, lettingR1 to be an orthogonal complementary matrix ofR, such thatRR>+R1R>1 = In, and N = XR1R > 1 X >
n ,
E [ 1
p tr[MXX>/n]
] = 1
p tr[In − λE [tr[M ]] + E [MN ]]
→ 1 γ − λ γ
∫ 1
x+ λ dF̄ζ(x) + γ − ζ γ
∫ 1
x+ λ dF̄ζ(x),
where F̄ζ is the companion MP law, that is, F̄ζ = (1 − γ)δ0 + γFζ . The third term calculated by using that XR and XR1 are independent for a Gaussian random matrix X , so that M,N are independent, and that E [N ] = p−dn In. Thus
E [ 1
p tr[MXX>/n]
] → 1
γ − λ+ ζ − γ γ θ̄1(ζ, λ)
= 1 γ − λ+ ζ − γ γ [ 1− ζ λ + ζθ1(ζ, λ) ] .
Then E [ 1
p tr[MXX>/n]2
] = 1 p E [ tr[In + λ 2M2 +MNMN − 2λM + 2MN − λM2N − λMNM ] .
Note that
E [MNMN |M ] = M [(p− d)(M> + tr(M)In) + (p− d)2M ]/n2
= p− d+ (p− d)2
n2 M2 + p− d n2 tr(M)M,
so
E [ 1
p tr[MXX>/n]2
] → 1
γ [1 + (λ2 − 2λ(γ − ζ) + (γ − ζ)2)θ̄2(ζ, λ)
+ 2(γ − ζ − λ)θ̄1(ζ, λ) + (γ − ζ)θ̄21(ζ, λ)]. Thus we find the following exprssion for the limiting squared bias:
bias2 → α 2
γ [γ − 1 + (λ− γ + ζ)2θ̄2 + (γ − ζ)θ̄21].
With similar calculations (that we omit for brevity), we can find
δ2 → σ2(θ̄1(ζ, λ)− (λ+ ζ − γ)θ̄2(ζ, λ)).
Therefore the AMSE of β̂d is
AMSE = α2
γ [γ − 1 + (λ− γ + ζ)2θ̄2 + (γ − ζ)θ̄21] + σ2[θ̄1(ζ, λ)− (λ+ ζ − γ)θ̄2(ζ, λ)].
(11)

heading: A.8 PROOF OF THEOREM 4.3
text: Proof. Recall that we have m,n → ∞, such that m/n → ξ. Then we need to take ξ → 0. However, we find it more convenient to do the calculation directly from the finite sample results as m,n, p → ∞ with m/n → 0, p/n → γ, It is not hard to check that computing the results in the other way (i.e., interchanging the limits), leads to the same results. Starting from our bias formula for primal sketching, we first get
bias2 = α2
p
[ λ2 + n−m+ p(n−m) + (n−m)2
n2 − 2λn−m n
] E [ tr[ ( X>L>LX/n+ λIp )−2 ] ]
→ α2[(λ− 1)2 + γ]/λ2.
The limit of the trace term is not entirely trivial, but it can be calculated by (1) observing that the m × p sketched data matrix P = LX has iid normal entries (2) thus the operator norm of P>P/n vanishes, (3) and so by a simple matrix perturbation argument the trace concentrates around p/λ2. This gives the rough steps of finding the above limit. Moreover,
δ2 = σ2 n2 E [ tr[ ( X>L>LX/n+ λIp )−2 X>X] ] → γσ2/λ2 · EFγX2 = γσ2/λ2
So the MSE is M(λ) = α2[(λ − 1)2 + γ]/λ2 + σ2 · γ/λ2. From this it is elementary to find the optimal λ and its objective value.

heading: A.9 PROOF OF THEOREM 4.4
text: Proof. Note that the bias can be written as
bias2 = α2
p E
[ tr[ ( XRR>X>
nd + λIn
)−1 XX>
nd ]2
]
− 2α 2 p E [ tr[ ( XRR>X>/n+ λIn )−1 XX>/n] ] + α2.
Write G = XX>. Since RR> ∼ Wp(Ip, d), we have XRR>X> ∼ Wn(G, d). So XRR>X> d = G1/2WG1/2, where W ∼ Wn(In, d).
E [ tr[ ( XRR>X>
nd + λIn
)−1 XX>/n] ] = E [ tr[(G1/2WG1/2/d+ nλIn) −1G] ]
= E [ tr[( W
d + λ(
G n )−1)−1]
] .
So we need to find the law of Wd + λ γ ( G p ) −1. Suppose first that G = XX> ∼ Wn(In, p). Then W and G−1 are asymptotically freely independent. The l.s.d. of W/d is the MP law F1/ξ while the l.s.d. of G/p is the MP law F1/γ . We need to find the additive free convolution W Ḡ, where Ḡ = λγG −1.
Recall that the R-transform of a distribution F is defined by
RF (z) = m −1 F (−z)−
1 z ,
where m−1F (z) is the inverse function of the Stieltjes transform of F (e.g., Voiculescu et al., 1992; Hiai & Petz, 2006; Couillet & Debbah, 2011). We can find the R-transform by solving
mF (RF (z) + 1
z ) = −z.
Note that the R-transform of W/d is
RW (z) = 1
1− z/ξ .
The Stieltjes transform of G−1 is
mG−1(z) =
∫ 1
1/x− z dF1/γ(x) = −
1 z − 1 z2 m1/γ( 1 z )
= −1 z −
1− 1γ − 1 z + √ (1 + 1γ + 1 z )
2 − 4γ 2 zγ
= − 1 + 1γ − 1 z +
√ (1 + 1γ − 1 z )
2 − 4γ 2 zγ .
Then the R-transform of G−1 is
RG−1(z) = − 1 z + γ + 1−
√ (γ + 1)2 − 4γ(z + 1)
2z
= γ − 1− √ (γ − 1)2 − 4γz 2z .
Since we have the property that Raµ(z) = aRµ(az),
RḠ = Rλ γG −1(z) =
γ − 1− √
(γ − 1)2 − 4λz 2z .
Hence we have
RW Ḡ = RW +RḠ = 1 1− z/ξ + γ − 1− √ (γ − 1)2 − 4λz 2z .
Moreover, the Stieltjes transform of µ = W Ḡ satisfies
m−1µ (z) = m −1 W Ḡ(z) = RF (−z)−
1 z =
1 1 + z/ξ + γ − 1− √ (γ − 1)2 + 4λz −2z − 1 z .
Note that
2 α2
p E
[ tr[ ( XRR>X>
nd + λIn
)−1 XX>/n] ] → 2α 2 γ Eµ [ 1 x ] = 2 α2 γ lim z→0 m(z),
α2
p E
[ tr[ ( XRR>X>
nd + λIn
)−1 XX>
nd ]2
] → α 2 γ Eµ [ 1 x2 ] = α2 γ lim z→0 d dz m(z).
So it suffices to find m(z) and ddzm(z) evaluated at zero.
This result can characterize the performance of sketching in the high SNR regime, where α σ. To understand the lower SNR regime, we need to study the variance, and thus we need to calculate
var = σ2 1
n E
[ tr[ ( XRR>X>
nd + λIn
)−2 XX>/n] ] = σ2E [ tr[( W
d + λ γ ( G p )−1)−2G−1] ] where G = XX> ∼ Wn(In, p) is a Wishart distribution, and XRR>X> =d G1/2WG1/2, with W ∼ Wn(In, r). This seems to be quite challenging, and we leave it to future work.

heading: A.10 RESULTS FOR PRIMAL GAUSSIAN SKETCHING
text: The statement requires some notions from free probability, see e.g., Voiculescu et al. (1992); Hiai & Petz (2006); Nica & Speicher (2006); Anderson et al. (2010); Couillet & Debbah (2011) for references .
Theorem A.1 (Bias of primal Gaussian sketch). Suppose X is an n× p standard Gaussian random matrix. Suppose also that L is a d× n matrix with i.i.d. N (0, 1/d) entries. Then the bias of primal sketch has the expressionMSE(β̂p) = α2 + α 2 γ [τ((a+b) −1b(a+b)−1b−1)−2τ((a+b)−1)], where a and b two free random variables, that are freely independent in a non-commutative probability space, and τ is their trace. Specifically, the law of a is the MP law F1/ξ and b = λγ b̃, where the law
of b̃ is the MP law F1/γ .
Proof of Theorem A.1. Note that bias2 = E [∥∥∥(X>L>LX/(nd) + λIp)−1 (X>X/n)β − β∥∥∥2
2
] ,
and ( X>L>LX/(nd) + λIp )−1 X> = X>(L>LXX>/(nd) + λIn) −1. Thus
bias2 = E [∥∥∥∥X>(L>LXX>nd + λIn)−1Xn β − β ∥∥∥∥2
2
]
= α2 + α2
p E[tr[(
XX>L>L
nd + λIn)
−1XX > n ( L>LXX> nd + λIn) −1XX > n ]
− 2 tr[(L >LXX>
nd + λIn)
−1XX >
n ]].
First we find the l.s.d. of (L >LXX> nd + λIn) −1XX> n . Write W = L >L, G = XX>. Then
( L>LXX>
nd + λIn)
−1XX >
n = (
WG
nd + λIn)
−1G
n = G−1(
W
d + λ(
G n )−1)−1G,
which is similar to (Wd + λ( G n ) −1)−1. So it suffices to find the l.s.d. of (Wd + λ γ ( G p ) −1)−1.
By the definition, W ∼ Wn(In, d), G ∼ Wn(In, p), therefore the l.s.d. of W/d converges to the MP law F1/ξ and the l.s.d. of G/p converges to the MP law F1/γ .
Also note that
( XX>L>L
nd + λIn)
−1XX > n ( L>LXX> nd + λIn) −1XX > n = ( W d + λnG−1)−1G−1( W d + λnG−1)−1G.
We write A = Wd , B = λ γ ( G p ) −1. Then it suffices to find
α2 p E [ tr[(A+B)−1B(A+B)−1B−1] ] .
We will find an expression for this using free probability. For this we will need to use some series expansions. There are two cases, depending on whether the operator norm of BA−1 is less than or greater than unity, leading to different series expansions. We will work out below the first case, but the second case is similar and leads to the same answer.
tr[(A+B)−1B(A+B)−1B−1] = tr[A−1(I +BA−1)−1BA−1(I +BA−1)−1B−1]
Since the operator norm of BA−1 is less unity, we have the von Neumann series expansion
[I +BA−1]−1 = ∞∑ i=0 (−BA−1)i,
then we have tr[(A+B)−1B(A+B)−1B−1] = ∑ i,j≥0 (−1)i+j tr[(BA−1)i+j+1B−1A−1]
= ∑ i,j≥0 (−1)i+j tr[(A−1B)i+j+1A−1B−1].
SinceA andB are asymptotically freely independent in the free probability space arising in the limit (e.g., Voiculescu et al., 1992; Hiai & Petz, 2006; Couillet & Debbah, 2011), and the polynomial (a−1b)i+j+1a−1b−1 involves an alternating sequence of a, b, we have
1 n tr[(A−1B)i+j+1A−1B−1]→ τ [(a−1b)i+j+1a−1b−1],
where a and the b are free random variables and τ is their law. Specifically, a is a free random variable with the MP law F1/ξ and b is λγ b̃
−1, where b̃ is a free r.v. with MP law F1/γ . Moreover, they are freely independent.
Hence, we have
1 n tr[(A+B)−1B(A+B)−1B−1]→ τ [ ∑ i≥0 (−1)i(a−1b)i+1 ∑ j≥0 (−1)j(a−1b)ja−1b−1]
= τ [(a−1b)(1 + a−1b)−1(1 + a−1b)−1a−1b−1]
= τ [(a+ b)−1b(a+ b)−1b−1].
Therefore,
bias2 → α2 + α 2
γ
[ 1
n tr[(A+B)−1B(A+B)−1B−1 − 2 tr[A+B]−1] ] = α2 + α2
γ [τ((a+ b)−1b(a+ b)−1b−1)− 2τ((a+ b)−1)].

heading: A.11 RESULTS FOR FULL SKETCHING
text: The full sketch estimator projects down the entire data, and then does ridge regression on the sketched data. It has the form
β̂f = ( X>L>LX/n+ λIp )−1 X>L>LY n .
We have
bias2 = α2λ2 ∫
1
ξx+ λ dFγ/ξ(x) = α
2λ2 1
ξ2 θ2
( γ
ξ , λ ξ ) var = σ2γ [∫ 1
ξx+ λ dFγ/ξ(x)− λ
∫ 1
(ξx+ λ)2 dFγ/ξ(x) ] = σ2γ [ 1
ξ θ1
( γ
ξ , λ ξ
) − λ 1
ξ2 θ2
( γ
ξ , λ ξ
)] ,
therefore
AMSE(β̂f ) = α 2λ2
1
ξ2 θ2
( γ
ξ , λ ξ
) + σ2γ [ 1
ξ θ1
( γ
ξ , λ ξ
) − λ 1
ξ2 θ2
( γ
ξ , λ ξ
)]
The optimal λ for full sketch is always λ∗ = γσ 2
α2 , the same as ridge regression. Some simulation results are shown in Figure 10, and they show the expected shape (e.g., they decrease with ξ).

heading: A.12 NUMERICAL RESULTS

heading: A.12.1 DUAL ORTHOGONAL SKETCHING
text: See Figure 11 for additional simulation results for dual orthogonal sketching.

heading: A.12.2 PERFORMANCE AT A FIXED REGULARIZATION PARAMETER
text: First we fix the regularization parameter at the optimal value for original ridge regression. The results are visualized in Figure 12. On the x axis, we plot the reduction in sample size m/n for primal sketch, and the reduction in dimension d/p for dual sketch. In this case, primal and dual sketch will increase both bias and variance, and empirically in the current case, dual sketch increases them more. So in this particular case, primal sketch is preferred.

heading: A.12.3 PERFORMANCE AT THE OPTIMAL REGULARIZATION PARAMETER
text: We find the optimal regularization parameter λ for primal and dual orthogonal sketching. Then we use the optimal regularization parameter for all settings, see Figure 13. Both primal and dual sketch increase the bias, but decrease the variance. It is interesting to note that, for equal parameters ξ and
ζ, and in our particular case, dual sketch has smaller variance, but larger bias. So primal sketch is preferred bias or MSE is important, but dual sketch is more desired when one wants smaller variance. All in all, dual sketch has larger MSE than primal sketch in the current setting. It can also be seen that in this specific example, the optimal λ for primal sketch is smaller than that of dual sketch. However these results are hard to interpret, because there is no natural correspondence between the two parameters ξ and ζ.

heading: A.13 COMPUTATIONAL COMPLEXITY
text: Since sketching is a method to reduce computational complexity, it is important to discuss how much computational efficiency we gain. Recall our three estimators
β̂ = ( X>X/n+ λIp )−1 X>Y/n = n−1X> ( XX>/n+ λIn )−1 Y,
β̂p = ( X>L>LX/n+ λIp )−1 X>Y/n,
β̂d = n −1X> ( XRR>X>/n+ λIn )−1 Y,
Their computational complexity, when computed in the usual way, is:
• No sketch (Standard ridge): if p < n, computing X>Y and X>X requires O(np) and O(np2) flops, then solving the linear equation (X>X/n + λIp)β̂ = X>Y/n requires O(p3) flops by the LU decomposition. It is O(np2) flops in total.
If p > n, we use the second formula for β̂, and the total flops is O(pn2). • Primal sketch: for the Hadamard sketch (and other sketches based on the FFT), computing LX by FFT requires mp log n, computing (LX)>LX requires mp2, so the total flops is O(p3 +mp(log n+ p)). So the primal sketch can reduce the computation cost only when p < n.
• Dual sketch: computing XRR>X> requires nd (log p + n) flops by FFT, solving (XRR>X>/n+ λIn)
−1 Y requires O(n3) flops, the matrix-vector multiplication of X> and (XRR>X>/n+λIn)−1Y requiresO(np) flops, so the total flops isO(n3+nd(log p+ n)). Dual sketching can reduce the computation cost only when p > n.

Review summaries:
This paper deals with 3 theoretical properties of ridge regression. First, it proves that the ridge regression estimator is equivalent to a specific representation which is useful as for instance it can be used to derive the training error of the ridge estimator. Second, it provides a bias correction mechanism for ridge regression and finally it provides proofs regarding the accuracy of several sketching algorithms for ridge regression.

This paper presents a theoretical study of ridge regression, focusing on the practical problems of correcting for the bias of the cross-validation based estimate of the optimal regularisation parameter, and quantification of the asymptotic risk of sketching algorithms for ridge regression, both in the p / n -> gamma in (0, 1) regime (n = # data points, p = # dimensions). The authors derive most of their results exploiting their (AFAICT) new asymptotic characterisation of the ridge regression estimator which may be of independent interest. The whole study is complemented by a series of numerical experiments.