arXiv:1607.00743v1 [] 4 Jul 2016

A Residual Bootstrap for High-Dimensional
Regression with Near Low-Rank Designs

Miles E. Lopes
Department of Statistics
University of California, Berkeley
Berkeley, CA 94720
mlopes@stat.berkeley.edu

Abstract
We study the residual bootstrap (RB) method in the context of high-dimensional
linear regression. Specifically, we analyze the distributional approximation of linear contrasts c⊤ (βbρ − β), where βbρ is a ridge-regression estimator. When regression coefficients are estimated via least squares, classical results show that RB
consistently approximates the laws of contrasts, provided that p ≪ n, where the
design matrix is of size n × p. Up to now, relatively little work has considered
how additional structure in the linear model may extend the validity of RB to the
setting where p/n ≍ 1. In this setting, we propose a version of RB that resamples
residuals obtained from ridge regression. Our main structural assumption on the
design matrix is that it is nearly low rank — in the sense that its singular values
decay according to a power-law profile. Under a few extra technical assumptions, we derive a simple criterion for ensuring that RB consistently approximates
the law of a given contrast. We then specialize this result to study confidence
intervals for mean response values Xi⊤ β, where Xi⊤ is the ith row of the design. More precisely, we show that conditionally on a Gaussian design with near
low-rank structure, RB simultaneously approximates all of the laws Xi⊤ (βbρ − β),
i = 1, . . . , n. This result is also notable as it imposes no sparsity assumptions on
β. Furthermore, since our consistency results are formulated in terms of the Mallows (Kantorovich) metric, the existence of a limiting distribution is not required.

1 Introduction
Until recently, much of the emphasis in the theory of high-dimensional statistics has been on “first
order” problems, such as estimation and prediction. As the understanding of these problems has
become more complete, attention has begun to shift increasingly towards “second order” problems,
dealing with hypothesis tests, confidence intervals, and uncertainty quantification [1, 2, 3, 4, 5,
6]. In this direction, much less is understood about the effects of structure, regularization, and
dimensionality — leaving many questions open. One collection of such questions that has attracted
growing interest deals with the operating characteristics of the bootstrap in high dimensions [7, 8, 9].
Due to the fact that bootstrap is among the most widely used tools for approximating the sampling
distributions of test statistics and estimators, there is much practical value in understanding what
factors allow for the bootstrap to succeed in the high-dimensional regime.
The regression model and linear contrasts. In this paper, we focus our attention on highdimensional linear regression, and our aim is to know when the residual bootstrap (RB) method
consistently approximates the laws of linear contrasts. (A review of RB is given in Section 2.)
1

To specify the model, suppose that we observe a response vector Y ∈ Rn , generated according to
Y = Xβ + ε,
n×p

(1)

p

where X ∈ R
is the observed design matrix, β ∈ R is an unknown vector of coefficients, and
the error variables ε = (ε1 , . . . , εn ) are drawn i.i.d. from an unknown distribution F0 , with mean
0 and unknown variance σ 2 < ∞. As is conventional in high-dimensional statistics, we assume
the model (1) is embedded in a sequence of models indexed by n. Hence, we allow X, β, and p to
vary implicitly with n. We will leave p/n unconstrained until Section 3.3, where we will assume
p/n ≍ 1 in Theorem 3, and then in Section 3.4, we will assume further that p/n is bounded strictly
between 0 and 1. The distribution F0 is fixed with respect to n, and none of our results require F0
to have more than four moments.
Although we are primarily interested in cases where the design matrix X is deterministic, we will
also study the performance of the bootstrap conditionally on a Gaussian design. For this reason,
we will use the symbol E[. . . |X] even when the design is non-random so that confusion does not
arise in relating different sections of the paper. Likewise, the symbol E[. . . ] refers to unconditional
expectation over all sources of randomness. Whenever the design is random, we will assume X ⊥
⊥ ε,
denoting the distribution of X by PX , and the distribution of ε by Pε .
Within the context of the regression, we will be focused on linear contrasts c⊤ (βb − β), where c ∈ Rp
is a fixed vector and βb ∈ Rp is an estimate of β. The importance of contrasts arises from the fact
that they unify many questions about a linear model. For instance, testing the significance of the ith
coefficient βi may be addressed by choosing c to be the standard basis vector c⊤ = e⊤
i . Another
important problem is quantifying the uncertainty of point predictions, which may be addressed by
choosing c⊤ = Xi⊤ , i.e. the ith row of the design matrix. In this case, an approximation to the law
of the contrast leads to a confidence interval for the mean response value E[Yi ] = Xi⊤ β. Further
applications of contrasts occur in the broad topic of ANOVA [10].
Intuition for structure and regularization in RB. The following two paragraphs explain the core
conceptual aspects of the paper. To understand the role of regularization in applying RB to highdimensional regression, it is helpful to think of RB in terms of two ideas. First, if βbLS denotes the
ordinary least squares estimator, then it is a simple but important fact that contrasts can be written
as c⊤ (βbLS − β) = a⊤ ε where a⊤:= c⊤ (X ⊤ X)−1 X ⊤ . Hence, if it were possible to sample directly
from F0 , then the law of any such contrast could be easily determined. Since F0 is unknown, the
second key idea is to use the residuals of some estimator βb as a proxy for samples from F0 . When
p ≪ n, the least-squares residuals are a good proxy [11, 12]. However, it is well-known that leastsquares tends to overfit when p/n ≍ 1. When βbLS fits “too well”, this means that its residuals are
b
“too small”, and hence they give a poor proxy for F0 . Therefore, by using a regularized estimator β,
b
overfitting can be avoided, and the residuals of β may offer a better way of obtaining “approximate
samples” from F0 .
The form of regularized regression we will focus on is ridge regression:
βbρ := (X ⊤ X + ρIp×p )−1 X ⊤ Y,

(2)

where ρ > 0 is a user-specificed regularization parameter. As will be seen in Sections 3.2 and 3.3,
the residuals obtained from ridge regression lead to a particularly good approximation of F0 when
the design matrix X is nearly low-rank, in the sense that most of its singular values are close to
0. In essence, this condition is a form of sparsity, since it implies that the rows of X nearly lie
in a low-dimensional subspace of Rp . However, this type of structural condition has a significant
advantage over the the more well-studied assumption that β is sparse. Namely, the assumption that
X is nearly low-rank can be inspected directly in practice — whereas sparsity in β is typically
unverifiable. In fact, our results will impose no conditions on β, other than that kβk2 remains
bounded as (n, p) → ∞. Finally, it is worth noting that the occurrence of near low-rank design
matrices is actually very common in applications, and is often referred to as collinearity [13, ch.
17].
Contributions and outline. The primary contribution of this paper is a complement to the work of
Bickel and Freedman [12] (hereafter B&F 1983) — who showed that in general, the RB method fails
2

to approximate the laws of least-squares contrasts c⊤ (βbLS − β) when p/n ≍ 1. Instead, we develop
an alternative set of results, proving that even when p/n ≍ 1, RB can successfully approximate the
laws of “ridged contrasts” c⊤ (βbρ − β) for many choices of c ∈ Rp , provided that the design matrix
X is nearly low rank. A particularly interesting consequence of our work is that RB successfully
approximates the law c⊤ (βbρ − β) for a certain choice of c that was shown in B&F 1983 to “break”
RB when applied to least-squares. Specifically, such a c can be chosen as one of the rows of X with
a high leverage score (see Section 4). This example corresponds to the practical problem of setting
confidence intervals for mean response values E[Yi ] = Xi⊤ β. (See [12, p. 41], as well as Lemma 2
and Theorem 4 in Section 3.4). Lastly, from a technical point of view, a third notable aspect of our
results is that they are formulated in terms of the Mallows-ℓ2 metric, which frees us from having to
impose conditions that force a limiting distribution to exist.
Apart from B&F 1983, the most closely related works we are aware of are the recent papers [7]
and [8], which also consider RB in the high-dimensional setting. However, these works focus on
role of sparsity in β and do not make use of low-rank structure in the design, whereas our work deals
only with structure in the design and imposes no sparsity assumptions on β.
The remainder of the paper is organized as follows. In Section 2, we formulate the problem of
approximating the laws of contrasts, and describe our proposed methodology for RB based on ridge
regression. Then, in Section 3 we state several results that lay the groundwork for Theorem 4, which
shows that that RB can successfully approximate all of the laws L(Xi⊤ (βbρ − β)|X), i = 1, . . . , n,
conditionally on a Gaussian design. Due to space constraints, all proofs are deferred to material that
will appear in a separate work.
Notation and conventions. If U and V are random variables, then L(U |V ) denotes the law of U ,
conditionally on V . If an and bn are two sequences of real numbers, then the notation an . bn
means that there is an absolute constant κ0 > 0 and an integer n0 ≥ 1 such that an ≤ κ0 bn for all
n ≥ n0 . The notation an ≍ bn means that an . bn and bn . an . For a square matrix A ∈ Rk×k
whose eigenvalues are real, we denote them by λmin (A) = λk (A) ≤ · · · ≤ λ1 (A) = λmax (A).

2 Problem setup and methodology
Problem setup. For any c ∈ Rp , it is clear that conditionally on X, the law of c⊤ (βbρ − β) is
completely determined by F0 , and hence it makes sense to use the notation

 
Ψρ (F0 ; c) := L c⊤ (βbρ − β)
(3)
X .

The problem we aim to solve is to approximate the distribution Ψρ (F0 ; c) for suitable choices of c.

Review of the residual bootstrap (RB) procedure. We briefly explain the steps involved in the
residual bootstrap procedure, applied to the ridge estimator βbρ of β. To proceed somewhat indirectly,
consider the following “bias-variance” decomposition of Ψρ (F0 ; c), conditionally on X,


 
Ψρ (F0 ; c) = L c⊤ βbρ − E[βbρ |X] 
(4)
X + c⊤ E[βbρ |X] − β .
|
{z
}
{z
}
|
=: Φρ (F0 ;c)

=: bias(Φρ (F0 ;c))

Note that the distribution Φ(F0 ; c) has mean zero, and so that the second term on the right side is
the bias of Φρ (F0 ; c) as an estimator of Ψρ (F0 ; c). Furthermore, the distribution Φρ (F0 ; c) may be
viewed as the “variance component” of Ψρ (F0 ; c). We will be interested in situations where the
regularization parameter ρ may be chosen small enough so that the bias component is small. In this
case, one has Ψρ (F0 ; c) ≈ Φρ (F0 ; c), and then it is enough to find an approximation to the law
Φρ (F0 ; c), which is unknown. To this end, a simple manipulation of c⊤ (βbρ − E[βbρ ]) leads to


Φρ (F0 ; c) = L(c⊤ (X ⊤ X + ρIp×p )−1 X ⊤ ε
(5)
X).

Now, to approximate Φρ (F0 ; c), let Fb be any centered estimate of F0 . (Typically, Fb is obtained by
using the centered residuals of some estimator of β, but this is not necessary in general.) Also, let
ε∗ = (ε∗1 , . . . , ε∗n ) ∈ Rn be an i.i.d. sample from Fb . Then, replacing ε with ε∗ in line (5) yields


Φρ (Fb ; c) = L(c⊤ (X ⊤ X + ρIp×p )−1 X ⊤ ε∗ 
(6)
X).
3

At this point, we define the (random) measure Φρ (Fb ; c) to be the RB approximation to Φρ (F0 ; c).
Hence, it is clear that the RB approximation is simply a “plug-in rule”.
A two-stage approach. An important feature of the procedure just described is that we are free
to use any centered estimator Fb of F0 . This fact offers substantial flexibility in approximating
Ψρ (F0 ; c). One way of exploiting this flexibility is to consider a two-stage approach, where a “pilot”
ridge estimator βb̺ is used to first compute residuals whose centered empirical distribution function
is Fb̺ , say. Then, in the second stage, the distribution Fb̺ is used to approximate Φρ (F0 ; c) via the
relation (6). To be more detailed, if (b
e1 (̺), . . . , ebn (̺)) = eb(̺) := Y − X βb̺ are the residuals of βb̺ ,
b
then we define F̺ to be the distribution that places mass 1/n at each of the values ebi (̺) − ē(̺) with
P
ē(̺) := n1 ni=1 ebi (̺). Here, it is important to note that the value ̺ is chosen to optimize Fb̺ as an
approximation to F0 . By contrast, the choice of ρ depends on the relative importance of width and
coverage probability for confidence intervals based on Φρ (Fb̺ ; c). Theorems 1, 3, and 4 will offer
some guidance in selecting ̺ and ρ.
Resampling algorithm. To summarize the discussion above, if B is user-specified number of
bootstrap replicates, our proposed method for approximating Ψρ (F0 ; c) is given below.
1. Select ρ and ̺, and compute the residuals eb(̺) = Y − X βb̺ .

2. Compute the centered distribution function Fb̺ , putting mass 1/n at each ebi (̺) − ē(̺).
3. For j = 1, . . . , B:

• Draw a vector ε∗ ∈ Rn of n i.i.d. samples from Fb̺ .
• Compute zj := c⊤ (X ⊤ X + ρIp×p )−1 X ⊤ ε∗ .

4. Return the empirical distribution of z1 , . . . , zB .

Clearly, as B → ∞, the empirical distribution of z1 , . . . , zB converges weakly to Φρ (Fb̺ ; c), with
probability 1. As is conventional, our theoretical analysis in the next section will ignore Monte Carlo
issues, and address only the performance of Φρ (Fb̺ ; c) as an approximation to Ψρ (F0 ; c).

3 Main results

The following metric will be central to our theoretical results, and has been a standard tool in the
analysis of the bootstrap, beginning with the work of Bickel and Freedman [14].
The Mallows (Kantorovich) metric. For two random vectors U and V in a Euclidean space, the
Mallows-ℓ2 metric is defined by
n h
i
o
(7)
d22 (L(U ), L(V )) := inf E kU − V k22 : (U, V ) ∼ π
π∈Π

where the infimum is over the class Π of joint distributions π whose marginals are L(U ) and L(V ).
It is worth noting that convergence in d2 is strictly stronger than weak convergence, since it also
requires convergence of second moments. Additional details may be found in the paper [14].
3.1 A bias-variance decomposition for bootstrap approximation

To give some notation for analyzing the bias-variance decomposition of Ψρ (F0 ; c) in line (4), we
define the following quantities based upon the ridge estimator βbρ . Namely, the variance is
vρ = vρ (X; c) := var(Ψρ (F0 ; c)|X) = σ 2 kc⊤ (X ⊤ X + ρIp×p )−1 X ⊤ k22 .

To express the bias of Φρ (F0 ; c), we define the vector δ(X) ∈ Rp according to


δ(X) := β − E[βbρ ] = Ip×p − (X ⊤ X + ρIp×p )−1 X ⊤ X β,
4

(8)

and then put
b2ρ = b2ρ (X; c) := bias2 (Φρ (F0 ; c)) = (c⊤ δ(X))2 .

(9)

We will sometimes omit the arguments of vρ and b2ρ to lighten notation. Note that vρ (X; c) does not
depend on β, and b2ρ (X; c) only depends on β through δ(X).
The following result gives a regularized and high-dimensional extension of some lemmas in Freedman’s early work [11] on RB for least squares. The result does not require any structural conditions
on the design matrix, or on the true parameter β.
Theorem 1 (consistency criterion). Suppose X ∈ Rn×p is fixed. Let Fb be any estimator of F0 , and
let c ∈ Rp be any vector such that vρ = vρ (X; c) 6= 0. Then with Pε -probability 1, the following
inequality holds for every n ≥ 1, and every ρ > 0,


b2
(10)
d22 √1vρ Ψρ (F0 ; c), √1vρ Φρ (Fb ; c) ≤ σ12 d22 (F0 , Fb ) + vρρ .

√
Remarks. Observe that the normalization 1/ vρ ensures that the bound is non-trivial, since the
√
distribution Ψρ (F0 ; c)/ vρ has variance equal to 1 for all n (and hence does not become degenerate
for large n). To consider the choice of ρ, it is simple to verify that the ratio b2ρ /vρ decreases monotonically as ρ decreases. Note also that as ρ becomes small, the variance vρ becomes large, and
likewise, confidence intervals based on Φρ (Fb; c) become wider. In other words, there is a trade-off
between the width of the confidence interval and the size of the bound (10).
Sufficient conditions for consistency of RB. An important practical aspect of Theorem 1 is that
for any given contrast c, the variance vρ (X; c) can be easily estimated, since it only requires an
estimate of σ 2 , which can be obtained from Fb. Consequently, whenever theoretical bounds on
d22 (F0 , Fb) and b2ρ (X; c) are available, the right side of line (10) can be controlled. In this way,
Theorem 1 offers a simple route for guaranteeing that RB is consistent. In Sections 3.2 and 3.3 to
follow, we derive a bound on E[d22 (F0 , Fb)|X] in the case where Fb is chosen to be Fb̺ . Later on in
Section 3.4, we study RB consistency in the context of prediction with a Gaussian design, and there
we derive high probability bounds on both vρ (X; c) and b2ρ (X; c) where c is a particular row of X.
3.2 A link between bootstrap consistency and MSPE
If βb is an estimator of β, its mean-squared prediction error (MSPE), conditionally on X, is defined
as


 
mspe(βb |X) := n1 E kX(βb − β)k22 
(11)
X .
The previous subsection showed that in-law approximation of contrasts is closely tied to the approximation of F0 . We now take a second step of showing that if the centered residuals of an estimator
βb are used to approximate F0 , then the quality of this approximation can be bounded naturally in
terms of mspe(βb |X). This result applies to any estimator βb computed from the observations (1).
Theorem 2. Suppose X ∈ Rn×p is fixed. Let βb be any estimator of β, and let Fb be the empirical
b Also, let Fn denote the empirical distribution of n i.i.d.
distribution of the centered residuals of β.
samples from F0 . Then for every n ≥ 1,



2


E d22 (Fb, F0 )X
≤ 2 mspe(βb |X) + 2 E[d22 (Fn , F0 )] + 2σn .
(12)
Remarks. As we will see in the next section, the MSPE of ridge regression can be bounded in a
sharp way when the design matrix is approximately low rank, and there we will analyze mspe(βb̺ |X)
for the pilot estimator. Consequently, when near low-rank structure is available, the only remaining
issue in controlling the right side of line (12) is to bound the quantity E[d22 (Fn , F0 )|X]. The very
recent work of Bobkov and Ledoux [15] provides an in-depth study of this question, and they derive
a variety bounds under different tail conditions on F0 . We summarize one of their results below.
Lemma 1 (Bobkov and Ledoux, 2014). If F0 has a finite fourth moment, then
E[d22 (Fn , F0 )] . log(n)n−1/2 .
5

(13)

Remarks. The fact that the squared distance is bounded at the rate of log(n)n−1/2 is an indication that d2 is a rather strong metric on distributions. For a detailed discussion of this result, see
Corollaries 7.17 and 7.18 in the paper [15]. Although it is possible to obtain faster rates when more
b
stringent tail conditions are placed on F0 , we will only need a fourth moment, since the mspe(β|X)
term in Theorem 2 will often have a slower rate than log(n)n−1/2 , as discussed in the next section.
3.3 Consistency of ridge regression in MSPE for near low rank designs
In this subsection, we show that when the tuning parameter ̺ is set at a suitable rate, the pilot ridge
estimator βb̺ is consistent in MSPE when the design matrix is near low-rank — even when p/n is
large, and without any sparsity constraints on β. We now state some assumptions.
A1. There is a number ν > 0, and absolute constants κ1 , κ2 > 0, such that
b ≤ κ2 i−ν
κ1 i−ν ≤ λi (Σ)

for all i = 1, . . . , n ∧ p.

A2. There are absolute constants θ, γ > 0, such that for every n ≥ 1,

̺
n

= n−θ and

ρ
n

= n−γ .

A3. The vector β ∈ Rp satisfies kβk2 . 1.

Due to Theorem 2, the following bound shows that the residuals of βb̺ may be used to extract a
consistent approximation to F0 . Two other notable features of the bound are that it is non-asymptotic
and dimension-free.
Theorem 3. Suppose that X ∈ Rn×p is fixed and that Assumptions 1–3 hold, with p/n ≍ 1. Assume
1
ν
1
further that θ is chosen as θ = 2ν
3 when ν ∈ (0, 2 ), and θ = ν+1 when ν > 2 . Then,
(
2ν
if ν ∈ (0, 21 ),
n− 3
ν
(14)
mspe(βb̺ |X) .
− ν+1
n
if ν > 12 .
Also, both bounds in (14) are tight in the sense that β can be chosen so that βb̺ attains either rate.

b are observable, they may be used to estimate ν and guide
Remarks. Since the eigenvalues λi (Σ)
the selection of ̺/n = n−θ . However, from a practical point of view, we found it easier to select ̺
via cross-validation in numerical experiments, rather than via an estimate of ν.

A link with Pinsker’s Theorem. In the particular case when F0 is a centered Gaussian distribution, the “prediction problem” of estimating Xβ is very similar to estimating the mean parameters of
a Gaussian sequence model, with error measured in the ℓ2 norm. In the alternative sequence-model
format, the decay condition on the eigenvalues of n1 X ⊤ X translates into an ellipsoid constraint on
the mean parameter sequence [16, 17]. For this reason, Theorem 3 may be viewed as “regression
version” of ℓ2 error bounds for the sequence model under an ellipsoid constraint (cf. Pinsker’s Theorem, [16, 17]). Due to the fact that the latter problem has a very well developed literature, there
may be various “neighboring results” elsewhere. Nevertheless, we could not find a direct reference
for our stated MSPE bound in the current setup. For the purposes of our work in this paper, the more
important point to take away from Theorem 3 is that it can be coupled with Theorem 2 for proving
consistency of RB.
3.4 Confidence intervals for mean responses, conditionally on a Gaussian design
In this section, we consider the situation where the design matrix X has rows Xi⊤ ∈ Rp drawn
i.i.d. from a multivariate normal distribution N (0, Σ), with X ⊥
⊥ ε. (The covariance matrix Σ may
vary with n.) Conditionally on a realization of X, we analyze the RB approximation of the laws
Ψρ (F0 ; Xi ) = L(Xi⊤ (βbρ − β)|X). As discussed in Section 1, this corresponds to the problem of
setting confidence intervals for the mean responses E[Yi ] = Xi⊤ β. Assuming that the population
eigenvalues λi (Σ) obey a decay condition, we show below in Theorem 4 that RB succeeds with high
PX -probability. Moreover, this consistency statement holds for all of the laws Ψρ (F0 ; Xi ) simultaneously. That is, among the n distinct laws Ψρ (F0 ; Xi ), i = 1, . . . , n, even the worst bootstrap
approximation is still consistent. We now state some population-level assumptions.
6

A4. The operator norm of Σ ∈ Rp×p satisfies kΣkop . 1.
Next, we impose a decay condition on the eigenvalues of Σ. This condition also ensures that Σ is
invertible for each fixed p — even though the bottom eigenvalue may become arbitrarily small as p
becomes large. It is important to notice that we now use η for the decay exponent of the population
eigenvalues, whereas we used ν when describing the sample eigenvalues in the previous section.
A5. There is a number η > 0, and absolute constants k1 , k2 > 0, such that for all i = 1, . . . , p,
k1 i−η ≤ λi (Σ) ≤ k2 i−η .
A6. There are absolute constants k3 , k4 ∈ (0, 1) such that for all n ≥ 3, we have the bounds
k3 ≤ np ≤ k4 and p ≤ n − 2.
The following lemma collects most of the effort needed in proving our final result in Theorem 4.
Here it is also helpful to recall the notation ρ/n = n−γ and ̺/n = n−θ from Assumption 2.
Lemma 2. Suppose that the matrix X ∈ Rn×p has rows Xi⊤ drawn i.i.d. from N (0, Σ), and that
Assumptions 2–6 hold. Furthermore, assume that γ chosen so that 0 < γ < min{η, 1}. Then, the
statements below are true.
(i) (bias inequality)
Fix any τ > 0. Then, there is an absolute constant κ0 > 0, such that for all large n, the following
event holds with PX -probability at least 1 − n−τ − ne−n/16 ,
max b2ρ (X; Xi ) ≤ κ0 · n−γ · (τ + 1) log(n + 2).

1≤i≤n

(15)

(ii) (variance inequality)
There are absolute constants κ1 , κ2 > 0 such that for all large n, the following event holds with
γ
PX -probability at least 1 − 4n exp(−κ1 n η ),
1
max
1≤i≤n vρ (X;Xi )

γ

≤ κ2 n1− η .

(16)

(iii) (mspe inequalities)
η
when
Suppose that θ is chosen as θ = 2η/3 when η ∈ (0, 21 ), and that θ is chosen as θ = 1+η
1
η > 2 . Then, there are absolute constants κ3 , κ4 , κ5 , κ6 > 0 such that for all large n,
(
2η
with PX -probability at least 1 − exp(−κ3 n2−4η/3 ), if η ∈ (0, 21 )
κ4 n− 3
η
mspe(βb̺ |X) ≤
2
− η+1
κ6 n
with PX -probability at least 1 − exp(−κ5 n 1+η ),
if η > 21 .

Remarks. Note that the two rates in part (iii) coincide as η approaches 1/2. At a conceptual level,
the entire lemma may be explained in relatively simple terms. Viewing the quantities mspe(βb̺ |X),
b2ρ (X; Xi ) and vρ (X; Xi ) as functionals of a Gaussian matrix, the proof involves deriving concentration bounds for each of them. Indeed, this is plausible given that these quantities are smooth
functionals of X. However, the difficulty of the proof arises from the fact that they are also highly
non-linear functionals of X.
We now combine Lemmas 1 and 2 with Theorems 1 and 2 to show that all of the laws Ψρ (F0 ; Xi )
can be simultaneously approximated via our two-stage RB method.

Theorem 4. Suppose that F0 has a finite fourth moment, Assumptions 2–6 hold, and γ is chosen
η
so that 1+η
< γ < min{η, 1}. Also suppose that θ is chosen as θ = 2η/3 when η ∈ (0, 12 ), and
η
when η > 12 . Then, there is a sequence of positive numbers δn with limn→∞ δn = 0, such
θ = η+1
that the event

h

i


E max d22 √1vρ Ψρ (F0 ; Xi ), √1vρ Φρ (Fb̺ ; Xi ) 
(17)
 X ≤ δn

1≤i≤n

has PX -probability tending to 1 as n → ∞.

Remark. Lemma 2 gives explicit bounds on the numbers δn , as well as the probabilities of the
corresponding events, but we have stated the result in this way for the sake of readability.
7

4 Simulations
In four different settings of n, p, and the decay parameter η, we compared the nominal 90% confidence intervals (CIs) of four methods: “oracle”, “ridge”, “normal”, and “OLS”, to be described
below. In each setting, we generated N1 := 100 random designs X with i.i.d. rows drawn from
N (0, Σ), where λj (Σ) = j −η , j = 1, . . . , p, and the eigenvectors of Σ were drawn randomly by
setting them to be the Q factor in a QR decomposition of a standard p × p Gaussian matrix. Then,
for each realization of X, we generated N2 := 1000 realizations of Y according to the model (1),
where β = 1/k1k2 ∈ Rp , and F0 is the centered t distribution on 5 degrees of freedom, rescaled to
have standard deviation σ = 0.1. For each X, and each corresponding Y , we considered the problem of setting a 90% CI for the mean response value Xi⊤⋆ β, where Xi⊤⋆ is the row with the highest
leverage score, i.e. i⋆ = argmax1≤i≤n Hii and H := X(X ⊤ X)−1 X ⊤ . This problem was shown in
B&F 1983 to be a case where the standard RB method based on least-squares fails when p/n ≍ 1.
Below, we refer to this method as “OLS”.
To describe the other three methods, “ridge” refers to the interval [Xi⊤⋆ βbρ − qb0.95 , Xi⊤⋆ βbρ − qb0.05 ],
where qbα is the α% quantile of the numbers z1 , . . . , zB computed in the proposed algorithm in
Section 2, with B = 1000 and c⊤ = Xi⊤⋆ . To choose the parameters ρ and ̺ for a given X and Y ,
we first computed rb as the value that optimized the MSPE error of a ridge estimator βbr with respect
to 5-fold cross validation; i.e. cross validation was performed for every distinct pair (X, Y ). We then
put ̺ = 5b
r and ρ = 0.1b
r, as we found the prefactors 5 and 0.1 to work adequately across various
settings. (Optimizing ̺ with respect to MSPE is motivated by Theorems 1, 2, and 3. Also, choosing ρ
to be somewhat smaller than ̺ conforms with the constraints on θ and γ in Theorem 4.) The method
“normal” refers to the CI based on the normal approximation L(Xi⊤⋆ (βbρ −β)|X) ≈ N (0, τb2 ), where
r, and σ
b2 is the usual unbiased estimate of σ 2 based
τb2 = σ
b2 kXi⊤⋆ (X ⊤ X +ρIp×p )−1 X ⊤ k22 , ρ = 0.1b
on OLS residuals. The “oracle” method refers to the interval [Xi⊤⋆ βbρ − q̃0.95 , Xi⊤⋆ βbρ − q̃0.05 ], with
ρ = 0.1b
r, and q̃α being the empirical α% quantile of Xi⊤ (βbρ − β) over all 1000 realizations of Y
based on a given X. (This accounts for the randomness in ρ = 0.1b
r.)

Within a given setting of the triplet (n, p, η), we refer to the “coverage” of a method as the fraction of
the N1 ×N2 = 105 instances where the method’s CI contained the parameter Xi⊤⋆ β. Also, we refer to
“width” as the average width of a method’s intervals over all of the 105 instances. The four settings of
(n, p, η) correspond to moderate/high dimension and moderate/fast decay of the eigenvalues λi (Σ).
Even in the moderate case of p/n = 0.45, the results show that the OLS intervals are too narrow
and have coverage noticeably less than 90%. As expected, this effect becomes more pronounced
when p/n = 0.95. The ridge and normal intervals perform reasonably well across settings, with
both performing much better than OLS. However, it should be emphasized that our study of RB
is motivated by the desire to gain insight into the behavior of the bootstrap in high dimensions
— rather than trying to outperform particular methods. In future work, we plan to investigate the
relative merits of the ridge and normal intervals in greater detail.
Table 1: Comparison of nominal 90% confidence intervals

setting 1
n = 100, p = 45, η = 0.5
setting 2
n = 100, p = 95, η = 0.5
setting 3
n = 100, p = 45, η = 1
setting 4
n = 100, p = 95, η = 1

width
coverage
width
coverage
width
coverage
width
coverage

oracle
0.21
0.90
0.22
0.90
0.20
0.90
0.21
0.90

ridge
0.20
0.87
0.26
0.88
0.21
0.90
0.26
0.92

normal
0.23
0.91
0.26
0.88
0.22
0.91
0.23
0.87

OLS
0.16
0.81
0.06
0.42
0.16
0.81
0.06
0.42

Acknowledgements. MEL thanks Prof. Peter J. Bickel for many helpful discussions, and gratefully acknowledges the DOE CSGF under grant DE-FG02-97ER25308, as well as the NSF-GRFP.
8

A

Proof of Theorem 1

Proof. Due to line (4) and Lemma 8.8 in B&F 1981,

d22 (Ψρ (F0 ; c), Φρ (Fb ; c)) = d22 Φρ (F0 ; c), Φρ (Fb ; c) + (c⊤ δ(X))2 .

(18)

If ε∗ ∈ Rn is a random vector whose entries are drawn i.i.d. from Fb , then the definition of Φρ gives
the matching relations


Φρ (F0 ; c) = L(c⊤ (X ⊤ X + ρ)−1 X ⊤ ε
X)
(19)


⊤
⊤
−1
⊤
∗
Φρ (Fb; c) = L(c (X X + ρ) X ε 
X).
To make use of these relations, we apply Lemma 8.9 in B&F 1981, which implies that if w ∈ Rn is
a generic deterministic vector, and if U = (U1 , . . . , Un ) and V = (V1 , . . . , Vn ) are random vectors
with i.i.d. entries, then
d22 (w⊤ U, w⊤ V ) ≤ kwk22 · d2 (U1 , V1 )2 .
Therefore,

d22 Φρ (F0 ), Φρ (Fb) ≤ kc⊤ (X ⊤ X + ρ)−1 X ⊤ k22 · d22 (ε1 , ε∗1 )
=

1
σ2

(20)

· vρ (X; c) · d22 (F0 , Fb ).

Combining this with line (18) and dividing through by vρ proves the claim.

B Proof of Theorem 2
Proof. By the triangle inequality,
d22 (Fb , F0 ) ≤ 2 d22 (Fb, Fn ) + 2 d22 (Fn , F0 ).

(21)

b which places mass 1/n at
Let F̃n be the (uncentered) empirical distribution of the residuals eb of β,
each value ebi , for i = 1, . . . , n. The proofs of Lemmas 2.1 and 2.2 in Freedman 1981, show that



 Pn
2 

 
 
E d2 (Fb, Fn )2 
+ E d2 (F̃n , Fn )2 
X ≤ E n1 i=1 εi
X




(22)
≤ n1 σ 2 + n1 E kb
e − εk22 
X


 
b 2
= n1 σ 2 + n1 E kX(β − β)k
2 X ,
b which holds for any estimator β.
b This
where we have used the algebraic identity eb − ε = X(β − β),
completes the proof.

C

Proof of Theorem 3

Proof. We begin with a simple bias-variance decomposition,
i
h

 2 



X
+
mspe(βb̺ |X) = n1 E X βb̺ − E βb̺ |X


2

1
n



X(E βb̺ |X − β)

2
.
2

(23)

We will handle the bias and variance terms separately. To consider the bias term, note that E[βbρ |X]−
β = Bβ, where
B = (X ⊤ X + ̺Ip×p )−1 X ⊤ X − Ip×p .
Hence,
1
b
n kX(E[β̺ |X]

− β)k22 =

1
2
n kXBβk2

= β ⊤ B( n1 X ⊤ X)Bβ.
9

(24)

If we let li = λi ( n1 X ⊤ X), then the eigenvalues of B( n1 X ⊤ X)B are of the form µi :=
In particular, it is simple to check1 that maxi µi ≍ ̺/n whenever θ ≤ ν, and so

li (̺/n)2
(li +̺/n)2 .

β ⊤ B( n1 X ⊤ X)Bβ . n̺ kβk22 = n−θ kβk22 .
(25)
Note that this bound is tight, since it is achieved whenever β is parallel to the top eigenvector of
B( n1 X ⊤ X)B.

To consider the variance term, note that βb̺ − E[βb̺ |X] = (X ⊤ X + ̺Ip×p )−1 X ⊤ ε, and so

h

 i
 2 
1
b̺ − E[βb̺ |X]) 2 
X = 1 tr X ⊤ X X ⊤ X + ̺Ip×p −1
E
X(
β

n
n
2
2
Pn∧p
li
.
= n1 i=1 li +̺/n
It is natural to decompose the sum in terms of the index set
I(n) := {i ∈ {1, . . . , n ∧ p} : li ≥ ̺/n},

(26)

(27)

θ
ν

which satisfies |I(n)| ≍ n . We will bound the variance term in two complementary cases; either
ν > 1/2 or ν ≤ 1/2. First assume ν > 1/2. Then,
1
n

p
X
i=1

2
li
li +̺/n

=

1
n

X

i∈I(n)

.

1
n

2
li
li +̺/n

|I(n)| +

1
n

Z

+

n∧p

|I(n)|

1
n

X

i6∈I(n)

2
li
li +̺/n

(28)

x−2ν
(̺/n)2 dx

(29)


 θ
using ν >
. n−1 n ν + n2θ · (|I(n)|)1−2ν
 θ

θ
≍ n−1 n ν + n2θ · (n ν )(1−2ν)

1
2

(30)
(31)

θ−ν
ν

= 2n
.
(32)
To see that this upper bound is tight, note that in line (28), we can use the term-wise lower bounds
( 1
if i ∈ I(n)
2
4
li
(33)
≥
l2i
1
li +̺/n
if i 6∈ I(n),
4 (̺/n)2
and then apply an integral approximation from below (which leads to the same rate). Combining the
bias and variance pieces, we have shown that
θ−ν
1
if ν > 1 .
EkX(βb̺ − β)k2 . n ν + n−θ
2

n

2

The bound is optimized when the two terms on the right side have the same rate, which leads to the
ν
.
choice θ = ν+1

In the case where ν ∈ (0, 21 ), the calculation proceeds in the same way up to line (30), where we
obtain the bound
n∧p
 θ

X
2
li
1
. n−1 n ν + n2θ · n1−2ν
(34)
n
li +̺/n
i=1

θ−ν

(35)
= n ν + n2(θ−ν) .
This bound is also tight due to the same reasoning as above. Note that in order for the bound to tend
to 0 as n → ∞, we must choose θ < ν. Furthermore, since we are working under the assumption
ν ∈ (0, 12 ), it follows that the right side of line (35) has rate equal to n2(θ−ν) . Combining the rates
for the bias and variance shows that
1
2
2(θ−ν)
b
+ n−θ
if ν ∈ (0, 21 ).
n EkX(β̺ − β)k2 . n
The bound is optimized when the two terms on the right side have the same rate, which leads to the
choice θ = 2ν
3 .
2

t(̺/n)
Note that if t ∈ R and f (t) := (t+̺/n)
2 , then f is maximized at t = ̺/n. Also, if θ ≤ ν, then there at
least one li that scales at the rate of ̺/n.
1

10

D

Proof of Lemma 2

The proof is split up into three pieces, corresponding to parts (i), (ii), and (iii) in the statement of the
result.
D.1

The bias inequality (15)

We prove inequality (15) by combining Lemmas 3 and 5 below.
Lemma 3. Assume the conditions of Lemma 2. For each i ∈ {1, . . . , n}, there are independent
random vectors ui (X), w(X) ∈ Rp such that the random variable Xi⊤ δ(X) can be represented
algebraically as
bρ (X; Xi ) = Xi⊤ δ(X) = ui (X)⊤ w(X).
Here, the vectors ui (X) can be represented in law as
L

ui (X) =

1
kzk2 Πp (z),

(36)

where z ∈ Rn is a standard Gaussian vector, and Πp (z) := (z1 , . . . , zp ). Also, the vector w(X)
satisfies the bound kw(X)k22 ≤ ρ4 kβk22 almost surely.
Proof. To fix notation, we write X ⊤ = Σ1/2 Z ⊤ where Z ⊤ ∈ Rp×n is a standard Gaussian matrix.
Recall that δ(X) = Bβ, where
B = Ip×p − (X ⊤ X + ρIp×p )−1 X ⊤ X.
Let Z = HLG⊤ be a “signed s.v.d.” for Z, as defined in Appendix E.2, where H ∈ Rn×p , L ∈
Rp×p , and G ∈ Rp×p . Now define ui (X) and w(X) according to
⊤
1/2
Xi⊤ δ(X) = e⊤
Bβ =
i XBβ = ei ZΣ

e⊤
H LG⊤ Σ1/2 Bβ .
{z
}
|i{z } |
=: ui (X)⊤ =: w(X)

From Lemma 17 in Appendix E.2, the rows e⊤
i H can be represented in distribution as
The same lemma also shows that the three matrices H, L, and G are independent.

(37)

1
kzk2 Πp (z).

Hence, to show that ui (X) and w(X) are independent, it suffices to show that w(X) is a function
only of G and L. In turn, it is enough to show that B is a function only of G and L. But this is
simple, because B is a function only of the matrix X ⊤ X, which may be written as
X ⊤ X = Σ1/2 Z ⊤ ZΣ1/2 = Σ1/2 GL2 G⊤ Σ1/2 .

(38)

It remains to show that kw(X)k22 ≤ ρ4 kβk22 almost surely. Combining the definition of w(X) with
line (38) gives

(39)
kw(X)k22 = β ⊤ BX ⊤ XB β.
2

li
1
⊤
The eigenvalues of BX ⊤ XB are of the form µi := n (l(ρ/n)
2 where li = λi ( n X X), and it is
i +ρ/n)
ρ
simple to check that the inequality maxi µi ≤ 4 holds for every realization of X.

Before proceeding to the second portion of the proof of inequality (15), we record some well-known
tail bounds for Gaussian quadratic forms due to Laurent and Massart [18], which will be useful at
various points later on.
Lemma 4 (Laurent & Massart, 2001). Let A ∈ Rn×n be a fixed symmetric matrix, and let z ∈ Rn
be a standard Gaussian vector. Then, for every t > 0,
h
i
√
P z ⊤ Az ≥ tr(A) + 2 |||A|||F t + 2 |||A|||op t ≤ exp(−t)
(40)
and

h
√i
P z ⊤ Az ≤ tr(A) − 2 |||A|||F t ≤ exp(−t).
11

(41)

The next lemma completes the proof of inequality (15).
Lemma 5. Assume the conditions of Lemma 2, and let τ > 0 be a constant. Then for every n ≥ 1,
the following event holds with probability at least 1 − n−τ − ne−n/16 ,
max b2ρ (X; Xi ) ≤ 5kβk22 · n−γ · (τ + 1) log(n + 2).

1≤i≤n

(42)

Proof. Applying the representation for bρ (X; Xi ) given in Lemma 3, there is a standard Gaussian
L

vector z ∈ Rn , such that ui (X) = Πp (z)/kzk2. Consequently,


L
1
⊤
b2ρ (X; Xi ) = kzk
Πp (z),
2 · Πp (z) w(X)w(X)

(43)

The same lemma also implies that for all t′ ∈ (0, 14 ),
i
h
1
1
√
≤ exp(−nt′ ).
P kzk
2 ≥
(1−2 t′ )n

(45)

2

where we may take z and w(X) to be independent by the same lemma. Using Lemma 4 on Gaussian
quadratic forms, as well as the fact that kw(X)k22 ≤ ρ4 kβk22 almost surely, we have for all t > 0,

"
#





√

ρ
⊤
⊤
2
P Πp (z) w(X)w(X) Πp (z) ≥ 4 kβk2 1 + 2 t + 2t)
w(X) ≤ exp(−t).
(44)




2

Now, we combine the bounds by integrating out w(X) in line (44) and choosing t′ = 1/16 in
line (45). Taking a union bound, we conclude that for any t > 0, and any fixed i = 1, . . . , n,
"
#
√

ρ
2
2 1
P bρ (X; Xi ) ≤ n · kβk2 · 2 1 + 2 t + 2t ≥ 1 − e−t − e−n/16 .
(46)
Finally, another union bound shows that the maximum of the bρ (X, Xi⊤ ) satisfies
P

"

max

1≤i≤n

b2ρ (X; Xi )

≤

ρ
n

·

kβk22

·

1
2 (1

#
√
+ 2 t + 2t) ≥ 1 − e−t+log(n) − ne−n/16 ,

(47)

which implies the
√stated result after choosing t = (τ + 1) log(n + 2), and noting that since t ≥ 1,
we have 21 (1 + 2 t + 2t) ≤ 5t = 5(τ + 1) log(n + 2), as well as e−t+log(n) ≤ e−τ log(n+2) ≤ n−τ
for every n ≥ 1.
D.2

The variance inequality (16)

The following “representation lemma” will serve as the basis for controlling the variance
vρ (X; Xi ) = σ 2 kXi⊤ (X ⊤ X + ρIp×p )−1 X ⊤ k22 .
Lemma 6. Assume the conditions of Lemma 2. For each i ∈ {1, . . . , n}, there is a random vector
vi (X) ∈ Rp and a random matrix M (X) ∈ Rp×p that are independent and satisfy the algebraic
relation
kXi⊤ (X ⊤ X + ρIp×p )−1 X ⊤ k22 = vi (X)⊤ M (X)vi (X).

Here, the vector vi (X) can be represented in law as
L

vi (X) =

1
kzk2 Πp (z),

(48)

where z ∈ Rn is a standard Gaussian vector and Πp (z) = (z1 , . . . , zp ). Also, the matrix M (X)
satisfies the algebraic relation
tr(M (X)) = kX(X ⊤X + ρIp×p )−1 X ⊤ k2F .
An explicit formula for M (X) is given below.
12

(49)

Proof. Define the matrix A := (X ⊤ X + ρIp×p )−1 X ⊤ X(X ⊤X + ρIp×p )−1 . Then,
⊤
kXi⊤ (X ⊤ X + ρIp×p )−1 X ⊤ k22 = e⊤
i XAX ei .

(50)

Using the notation in the proof of the previous lemma, let X = ZΣ1/2 where Z ∈ Rn×p is a
standard Gaussian random matrix. Furthermore, let Z = HLG⊤ be a signed s.v.d. for Z, as defined
in Appendix E.2. Then, we define vi (X) and M (X) according to
⊤
⊤
e⊤
LG⊤ Σ1/2{z
AΣ1/2 GL⊤} H ⊤ ei .
i XAX ei = ei H |
| {z }

(51)

=:M(X)

=:vi (X)⊤

Some algebra shows that M (X) satisfies the relation (49). As in the proof of Lemma 3, the argument
is completed using two properties of the signed s.v.d. of a standard Gaussian matrix: The rows of H
can be represented as Πp (z)/kzk2 where z ∈ Rn is a standard Gaussian vector, and the matrices H,
L, and G⊤ are independent. (See Lemma 17 in Appendix E.2.) To show that vi (X) and M (X) are
independent, first note that vi (X) only depends on H. Also, it is simple to check that M (X) only
depends on G and L, because A is a function only of X ⊤ X = Σ1/2 GL2 G⊤ Σ1/2 .
D.2.1 Concentration of the variance and bounds on its expected value
Due to Lemma 6, for each i = 1, . . . , n, we have the representation
L

vρ (X; Xi ) =

1
Π (z)⊤ M (X)Πp (z),
kzk22 p

(52)

where z ∈ Rn is a standard Gaussian vector, independent of M (X). Conditionally on M (X),
the quadratic form Πp (z)⊤ M (X)Πp (z) concentrates around tr(M (X)) by Lemma 4. The same
2
lemma
p also implies that kzk2 concentrates around n. In the next three subsections, we will show
that tr(M (X)) concentrates around its expected value, and obtain upper and lower bounds on the
expected value. We will need two-sided bounds in preparation for Theorem 4.
D.2.2 Concentration of

p
tr(M (X))

Lemma 7. Assume the conditions of Lemma 2. Then for every t > 0, and every n ≥ 1,
#
"
p
p
64 n1−γ t2
tr(M (X)) − E tr(M (X)) ≥ t ≤ 2 exp(− 54
P
kΣkop ).

(53)

p
Proof. We will show that tr(M (X)) is a Lipschitz function of a standard Gaussian matrix. Define
s2
, which satisfies the Lipschitz condition
the function gρ : R+ → [0, 1] by gρ (s) = s2 +(ρ/n)
|gρ (s) − gρ (s′ )| ≤ Ln |s − s′ |,

for all s, s′ ≥ 0, where Ln :=

√
3 3
8

√1 .
ρ/n

If σ(A) = (σ1 (A), . . . , σk (A)) denotes the vector of singular values of a rank k matrix A, then we
define gρ to act on σ(A) component-wise, i.e. gρ (σ(A)) = (gρ (σ1 (A)), . . . , gρ (σk (A))). Recall
from Lemma 6 that
p
(54)
tr(M (X)) = kX(X ⊤X + ρIp×p )−1 X ⊤ kF

and note that the ith singular value of the matrix X(X ⊤ X +ρIp×p )−1 X ⊤ is given by gρ (σi ( √1n X)).
Viewing the Frobenius norm of a matrix as the ℓ2 norm of its singular values, we have
p
tr(M (X)) = kgρ (σ( √1n X))k2 .
(55)
Write X ⊤ = Σ1/2 Z ⊤ for a standard Gaussian matrix Z ∈ Rn×p , and let f : Rn×p → R be defined
according to
p
f (Z) := tr(M (X)).
13

We claim that f is Lipschitz with respect to the Frobenius norm. Let W ⊤ ∈ Rp×n be a generic
matrix, and put A = √1n Σ1/2 Z ⊤ and B = √1n Σ1/2 W ⊤ . Then,
|f (Z) − f (W )| = kgρ (σ(A))k2 − kgρ (σ(B))k2
≤ gρ (σ(A)) − gρ (σ(B))

(56)
(57)

2

≤ Ln kσ(A) − σ(B)k2

(58)

≤ Ln kA − BkF

(59)

√1 Σ1/2
n

= Ln
≤

Ln
√
n

(Weilandt-Hoffman)

Z⊤ − W ⊤ F

p
kΣkop · Z ⊤ − W ⊤

F

,

(60)
(61)

where we have used a version of the Weilandt-Hoffman inequality for singular values [19, p.186],
as well as the inequality kM1 M2 kF ≤ kM1 kop kM2 kF , which holds for any square matrix M1 that
is compatible with M2 . (See Lemma 12 in Appendix E.) The statement of the lemma now follows
from the Gaussian concentration inequality. (See Lemma 15 in Appendix E).
p
D.2.3 Upper bound on E tr(M (X))

Lemma 8. Assume the conditions of Lemma 2. Then, the matrix M (X) satisfies
(
1
p
n(γ−η)+ 2 if η ∈ (0, 21 )
γ
E tr(M (X)) .
if η > 21 .
n 2η

(62)

p
Proof. By Jensen’s inequality, it is enough to bound E[tr(M (X))] from above. Define the unis
variate function ψ : R+ → R+ by ψ(s) := (√s+ρ/n)
2 , and observe that

2 
tr(M (X)) = tr (X ⊤ X + ρIp×p )−1 X ⊤ X
=

p
X

b
λ2i (Σ)
2
b
(λi (Σ)+ρ/n)

i=1

b 2 ))
ψ(λi (Σ

i=1

=

p
X

(63)


b 2) .
= tr ψ(Σ

Here where we use the “operator calculus” notation ψ(A) = U ψ(D)U ⊤ where A is a symmetric
matrix with spectral decomposition A = U DU ⊤ , and ψ(D) is the diagonal matrix whose ith diagb 2 )) is a concave matrix
onal entry is ψ(Dii ). It is simple to check that ψ is a concave, and so tr(ψ(Σ
2
b
functional of Σ by Lemma 13 in Section E.1 of Appendix A. Therefore, Jensen’s inequality implies

b 2 ]))
E[tr(M (X))] ≤ tr(ψ(E[Σ
(64)
P
= pi=1 ψ(λi (S)),
b 2 ]. Since X is Gaussian, Σ
b is a Wishart matrix up to scaling,
where we define the matrix S := E[Σ
and so Lemma 16 in Appendix E shows that this expectation may be evaluated exactly as
S = (1 + n1 )Σ2 +

tr(Σ)
n Σ.

(65)

We will now use this relation to apply an integral approximation to the right side of line (64). Clearly,
the eigenvalues of S are given by
λi (S) = (1 + n1 )λ2i (Σ) +
≍ i−2η +

tr(Σ) −η
.
n i

14

tr(Σ)
n λi (Σ)

(66)

Let r ∈ (0, 1) be a constant to be specified later. On the set of indices 1 ≤ i ≤ ⌈nr ⌉ we use the bound
1
ψ(λi (S)) ≤ 1, and on the set of indices i > ⌈nr ⌉ we use the bound ψ(λi (S)) ≤ (ρ/n)
2 λi (S).
−γ
2
Recalling the assumption ρ/n = n , we may decompose the inequality (64) as
E[tr(M (X))] ≤

⌈nr ⌉

X

ψ(λi (S)) +

i=1

. nr + n2γ

Z

p
X

ψ(λi (S))

(67)

i=⌈nr ⌉+1
p

x−2η +

nr

=: nr + n2γ hn (η, r).

tr(Σ) −η 
dx
n x

(68)
(69)

where the function hn is defined in the last line. The bound is optimized when the two terms on the
right are of the same order; i.e. when r solves the rate equation
nr ≍ n2γ hn (η, r).
Noting that
tr(Σ)
n

≍

 −η
n
n−1

(70)

if η ∈ (0, 1)
if η > 1,

the quantity hn (η, r) may be computed directly as
 1−2η
n
hn (η, r) ≍ nr(1−2η)
 r(1−2η)
n
+ nr(1−η)−1

if η ∈ (0, 12 ),
if η ∈ ( 12 , 1),
if η > 1.

(71)

(72)

If we let r = r∗ (η, γ) denote the solution of the rate equation (70), then some calculation shows that
under the assumption γ ∈ (0, 1),
(
2(γ − η) + 1 if η ∈ (0, 12 ),
(73)
r∗ (η, γ) = γ
if η > 21 .
η
When η ∈ (0, 1) this is straightforward. To show the details for η > 1, note that the rate equation (70) may be written as
nr ≍ n2γ+r(1−2η) + n2γ+r(1−η)−1 ,

(74)

which is the same as

1 ≍ n2(γ−ηr) + n2γ−ηr−1 .
In order for both terms on the right to be O(1), the number r must satisfy the constraints
r ≥ γη ,

r≥

γ
η

+

(75)
(76)

γ−1
η .

(77)

Since Lemma 2 assumes γ ∈ (0, 1), only the first constraint matters. Furthermore, when r ≥ γη ,
the second term in line (75) is o(1), and we are reduced to choosing r so that 1 ≍ n2(γ−ηr) , which
gives r = r∗ (η, γ) = γη . Substituting this value into line (69) completes the proof. (Note from the
discussion preceding line (69) that r must lie in the interval (0, 1), and this requires γ/η < 1, which
explains the assumption γ < min{η, 1} in Lemma 2.)
p
D.2.4 Lower bound on E tr(M (X))

Lemma 9. Assume the conditions of Lemma 2. Then, the matrix M (X) satisfies
p
γ
E tr(M (X)) & n 2η .

(78)

2
Note that if p/n ≍ 1, it is possible that nr > p for small values of n. Since we want nr ≤ p for the
integral in line (68), Lemma 2 is stated for “all large n”.

15

p
Proof. The variable tr(M (X)) may be written as kX ⊤X(X ⊤ X +ρIp×p )−1 kF . Since the Frobenius norm is a convex matrix functional, Jensen’s inequality implies
i
h
p
E tr(M (X)) ≥ E X ⊤ X(X ⊤ X + ρIp×p )−1
F
(79)
i
h

b −1 −1
,
= E Ip×p + nρ Σ
F

b :=
where the last step follows algebraically with Σ
If we define the univariate function
f : R+ → R+ by f (s) = (1 + nρ s)−1 , then last inequality is the same as
p


b −1
.
(80)
E tr(M (X)) ≥ E f Σ
F
1
⊤
n X X.

It is a basic fact that f is operator convex on the domain of positive semidefinite matrices [20,
p.117]. This yields an operator version of Jensen’s inequality with respect to the Loewner ordering
(Lemma 14 in Appendix E):

 −1 

b
b −1  f E Σ
.
(81)
E f Σ

Furthermore, if two matrices satisfy A  B  0, then kAkF ≥ kBkF [21, Corollary 7.7.4]. Using
this fact, as well as the formula for the expected inverse of a Wishart matrix [22, p. 97], we obtain
p
 −1 
b
E tr(M (X)) ≥ f E Σ
F

= f
=

n
−1
n−p−1 Σ

p
X
i=1

=

p
X
i=1

1+

ρ
n

·



F

1


n
−1 ) 2
n−p−1 λi (Σ

λ2i (Σ)
λi (Σ) +

ρ
n

·

2
n
n−p−1

!1/2

!1/2

(82)

.


n
. For any i ∈ J, the ith summand in
Define the index set J = i ∈ {1, . . . , p} : λi (Σ) ≥ nρ n−p−1
the previous line is at least 1/4. Also, assumption A6 that p/n is bounded strictly between 0 and 1,
as well as the decay condition on the λi (Σ), imply that |J| ≍ nγ/η , which completes the proof.
D.2.5 Putting the variance pieces together
Combining Lemmas 7, 8, and 9 with the Gaussian concentration inequality (Lemma 15 in Section E.2 of Appendix A) p
immediately gives the following result. (We choose t to be proportional to
the relevant bound on E[ tr(M (X))] in the Gaussian concentration inequality.)

Lemma 10. Assume the conditions of Lemma 2 and let tr(M (X)) be as in line (49). Then, there are
absolute constants κ1 , κ2 , . . . , κ6 > 0 such that the following upper-tail bounds hold for all large
n,
i
h
(83)
P tr(M (X)) ≥ κ1 n2(γ−η)+1 ≤ exp(−κ2 n2(1−η)+γ ), if η ∈ (0, 21 ),

and

i
h
γ(1−η)
P tr(M (X)) ≥ κ3 nγ/η ≤ exp(−κ4 n1+ η ),

and the following lower-tail bound holds for all large n,
i
h
γ(1−η)
P tr(M (X)) ≤ κ5 nγ/η ≤ exp(−κ6 n1+ η ),

if η > 12 ,

(84)

if η > 0.

(85)

Remarks. Note that in order for the last two probabilities to be small for large values of η > 0,
it is necessary that γ < 1, as assumed in Lemma 2. The next result completes the assembly of the
results in this Subsection D.2. Although the first two bounds in Lemma 11 are not necessary for the
statement of Theorem 4, they show that the variance vρ (X; Xi ) tends 0 as n → ∞ when γ < η,
as assumed in Theorem 4. In other words, we imposed the assumption γ < η so that confidence
intervals based on Φρ (Fb̺ ; Xi ) have width that tends to 0 asymptotically.
16

Lemma 11. Assume the conditions of Theorem 4 and let tr(M (X)) be as in line (49). Assume
γ < min{η, 1}. Then, there are absolute constants k1 , k2 , . . . , k6 > 0 such that the following
upper-tail bounds hold for all large n,
i
h
γ
P max vρ (X; Xi ) ≤ k1 n2(γ−η) ≥ 1 − 4n exp(−k2 n η ), if η ∈ (0, 21 )
(86)
1≤i≤n

and

i
h
γ
γ
P max vρ (X; Xi ) ≤ k3 n η −1 ≥ 1 − 4n exp(−k4 n η ),
1≤i≤n

and

h
P max

1
1≤i≤n vρ (X;Xi )

i
γ
γ
≤ k5 n1− η ≥ 1 − 4n exp(−k6 n η ),

if η > 12 ,

(87)

if η > 0.

(88)

Proof. We only prove the last inequality (88), since the other two inequalities are proven in a similar
way. By Lemma 6, we have
L

vρ (X; Xi ) =

1
Π (z)⊤ M (X)Πp (z)
kzk22 p

(89)

where z ∼ N (0, Ip×p ) and z ⊥
⊥ M (X). To apply the lower-tail
bound for Gaussian quadratic
p
forms, note that Hölder’s inequality implies kM (X)kF ≤
tr(M (X)) since kM (X)kop ≤ 1
almost surely. Therefore, letting t = t′ tr(M (X)) with t′ ∈ (0, 1) in inequality (41) gives

i
h
√



′
(90)
P Πp (z)⊤ M (X)Πp (z) ≥ (1 − 2 t′ ) tr(M (X))

M (X) ≥ 1 − exp − t tr(M (X))
Next, observe that inequality (41) with t = t′ · n for t′ ∈ (0, 1) gives,
h
i
√
P kzk22 ≤ (1 + 4 t′ )n ≥ 1 − exp(−t′ n).
If we define the event

(91)

)

(92)

then the previous two inequalities imply

 

P E1 
M (X) ≥ 1 − exp(−t′ tr(M (X))) − exp(−t′ · n)

(93)

E1 :=

(

1
kzk22

1
≤
Πp (z)⊤ M (X)Πp (z)

√
n
1+4√t′
(1−2 t′ ) tr(M(X))

≥ 1 − 2 exp(−t′ tr(M (X))),

since tr(M (X)) ≤ n almost surely. Next, let κ5 , κ6 > 0 be as in the previous lemma, and define
the event
o
n
γ
n
(94)
E2 := tr(M(X))
≤ κ15 n1− η ,
which has probability P(E2 ) ≥ 1 − exp(−κ6 n1+

γ(1−η)
η

).

We now put these items together. Starting with line (89), if we work on the intersection of E1 and
E2 , then for any fixed i = 1, . . . , n we have
P

h

1
vρ (X;Xi )

≤

√
1− γ
1+4√t′ 1
n η
(1−2 t′ ) κ5

i

h
i
≥ E 1E1 · 1E2

h
i
≥ 1 − E 1E1c + 1E2c

  

c
= 1 − E E 1E1c 
M (X) − P(E2 )
h
i
≥ 1 − 2E exp(−t′ tr(M (X)) − P(E2c)
h

′

i

= 1 − 2E exp(−t tr(M (X))) · (1E2 + 1 ) −
E2c

γ

≥ 1 − exp(−t′ κ5 n η ) − 3P(E2c )
γ

1+

≥ 1 − exp(−t′ κ5 n η ) − 3 exp(−κ6 n
γ

≥ 1 − 4 exp(− min{t′ κ5 , κ6 } · n η )

17

γ(1−η)
η

)

(95)
P(E2c )

where we have used the previous lemma to bound P(E2c ), and also the assumption γ ∈ (0, 1) to
conclude that γη ≤ 1 + γ(1−η)
. Taking a union bound over i = 1, . . . , n, proves the claim.
η
The last component of Lemma 2 is to prove the MSPE inequalities.
D.3

Proof of the MSPE inequalities

The proof of Theorem (3) shows that for any realization of X we have
h
i


mspe(βb̺ |X) := n1 E kX(βb̺ − β)k22 
X . n−θ kβk22 +

1
n

p∧n
X
i=1

2
li
li +̺ ,

(96)

where li = λi ( n1 X ⊤ X). Now observe that the second term on the right side matches the expression
for tr(M (X)) given in line (63) by replacing ρ with ̺ and multiplying by a factor of n1 . Therefore,
using Lemma 10 and recalling ̺/n = n−θ shows that there are absolute constants κ1 , κ2 , κ3 , κ4 > 0
such that for all large n,
and

h
i
P mspe(βb̺ |X) ≥ κ1 n−θ + n2(θ−η) ≤ exp(−κ2 n2(1−η)+θ ), if η ∈ (0, 12 ).
h
θ(1−η)
i
θ
P mspe(βb̺ |X) ≥ κ3 n−θ + n η −1 ≤ exp(−κ4 n1+ η ),

if η > 21 .

(97)
(98)

In line (97), the bound mspe(βb̺ |X) is optimized when n−θ ≍ n2(θ−η) , which explains the choice
θ
in line (98), the bound is optimized when n−θ ≍ n η −1 , which explains the
θ = 2η
3 . Similarly,
η
choice θ = η+1 . Substituting in these values θ yields the stated result.

E Background results
E.1

Results on matrices and convexity

Lemma 12. Let M1 ∈ Rk1 ×k1 and M2 ∈ Rk1 ×k2 . Then,
kM1 M2 kF ≤ kM1 kop kM2 kF .

(99)

Proof. Observe that
kM1 M2 k2F = tr(M2⊤ M1⊤ M1 M2 )

= tr((M1⊤ M1 )(M2 M2⊤ ))

≤

k1
X
i=1

λi (M1⊤ M1 ) · λi (M2 M2⊤ )

≤ kM1 k2op

k1
X

(100)

λi (M2 M2⊤ )

i=1
2
kM1 kop kM2 k2F ,

=
where we have used von Neumann’s trace inequality (also known as Fan’s inequality) [23, p.10] in
the third line.
A result on convex trace functionals. In the following lemma, an interval of the real line refers to
any set of the form (a, b),(a, b],[a, b), or [a, b], where −∞ ≤ a ≤ b ≤ ∞. We also define spec(M )
p×p
to be the set of eigenvalues of a square matrix M . The collection of symmetric
is
Pmatrices in R
p×p
denoted by S
. For a univariate function ϕ, the symbol tr(ϕ(M )) denotes i ϕ(λi (M )).

Lemma 13. Let I ⊂ R be an interval, and let M ⊂ Sp×p be a convex set such that spec(M ) ⊂ I
for all M ∈ M. Let ϕ : I → R be a convex function. Then, the functional
M 7→ tr(ϕ(M ))
(101)
is convex on M.
A proof may be found in the paper [24, Proposition 2].
18

Operator Jensen inequality. A function f : Sp×p → Sp×p is said to be operator convex if for all
λ ∈ [0, 1], and all A, B ∈ Sp×p ,
f (λA + (1 − λ)B)  λf (A) + (1 − λ)f (B),

(102)

where A  B means that B − A is positive semidefinite.

Lemma 14 (Operator Jensen inequality). Suppose f : Sp×p → Sp×p is operator convex, and let A
be a random Sp×p -valued matrix that is integrable. Then,
f (E[A])  E[f (A)].

(103)

Proof. It is enough to show that for all x ∈ Rp ,
x⊤ f (E[A])x ≤ x⊤ E[f (A)]x.

(104)

For any fixed x, consider the function g : Sp×p → R defined by g(A) = x⊤ f (A)x. It is clear that
g is a convex function in the usual sense, and so the ordinary version of Jensen’s inequality implies
g(E[A]) ≤ E[g(A)], which is the same as (104).
E.2

Results on Gaussian vectors and matrices

The following lemma is standard and is often referred to as the Gaussian concentration inequality [25].
Lemma 15. Let Z ∈ Rp be a standard Gaussian vector and let f : Rp → R be an L-Lipschitz
function with respect to the ℓ2 norm. Then for all t > 0,



−t2
(105)
P |f (Z) − E[f (Z)]| ≥ t ≤ 2 exp 2L
2 .
Next, we give a formula for the expected square of a Wishart matrix.

b=
Lemma 16. Let X ∈ Rn×p have rows drawn i.i.d. from N (0, Σ), and let Σ
b =
Proof. Write Σ
then

1
n

Pn

i=1

b 2 ] = (1 + 1 )Σ2 +
E[Σ
n

1
⊤
n X X.

Then,

tr(Σ)
n Σ.

Xi Xi⊤ where Xi⊤ ∈ Rp is the ith row of X. If we put Mi = Xi Xi⊤ ,
b 2 = 12 Pn M 2 + 12 P Mi Mj .
Σ
i
i6=j
i=1
n
n

Clearly, the Mi are independent with E[Mi ] = Σ for all i, and so

b 2 ] = 1 E[M 2 ] + 22 n Σ2 .
E[Σ
1
n
n 2

It remains to compute E[M12 ]. Write Xi = Σ1/2 Zi where Zi is a standard Gaussian vector in Rp ,
and let Σ1/2 = U Λ1/2 U ⊤ be a spectral decomposition of Σ1/2 where U ∈ Rp×p is orthogonal,
and Λ is diagonal with Λii = λi (Σ). By the orthogonal invariance of the normal distribution,
L
Xi = U Λ1/2 Zi , and so
h
i
E[M12 ] = E[X1 X1⊤ X1 X1⊤ ] = U Λ1/2 E Z1 Z1⊤ ΛZ1 Z1⊤ Λ1/2 U ⊤ .
(106)
Define the matrix M̃1 := Z1 Z1⊤ ΛZ1 Z1⊤ . It is straightforward to verify that E[M̃1 ] is diagonal, and
its jth diagonal entry is
E[M̃1 ]jj = tr(Σ) + 2λj .
Therefore,
E[M̃1 ] = tr(Σ)Ip×p + 2Λ,
and combining this with line (106) gives
E[M12 ] = tr(Σ)Σ + 2Σ2 .

19

Signed s.v.d. The following lemma describes the factors of an s.v.d. of a standard Gaussian matrix
Z ∈ Rn×p with n ≥ p. To make the statement of the lemma more concise, we define the term signed
s.v.d. below. This is merely a particular form of the s.v.d. that ensures uniqueness. Specifically, if
Z ∈ Rn×p is a full rank matrix with n ≥ p, then the signed s.v.d. of Z is given by
Z = HLG⊤ ,

(107)

where H ∈ Rn×p has orthonormal columns, the matrix L ∈ Rp×p is diagonal with L11 ≥ L22 ≥
· · · ≥ Lpp > 0, and G ∈ Rp×p is orthogonal with its first row non-negative, i.e. G1i ≥ 0 for all
i = 1, . . . , p. It is a basic fact from linear algebra that the signed s.v.d. of any full rank matrix in
Rn×p exists and is unique [21, Lemma 7.3.1]. This fact applies to Gaussian matrices, since they are
full rank with probability 1.
Lemma 17. Suppose n ≥ p, and let Z ∈ Rn×p be a random matrix with entries drawn i.i.d. from
N (0, 1). Let
Z = HLG⊤ ,

(108)

be the unique signed s.v.d. of Z as defined above. Then, the matrices H, L and G are independent.
Furthermore, if Hi⊤ ∈ Rp denotes the ith row of H, then for each i = 1, . . . , n, the marginal law of
Hi⊤ is given by
L

Hi⊤ =

1
kzk2

· Πp (z),

(109)

where z ∈ Rn is a standard Gaussian vector, and Πp (z) is the projection operator onto the first p
coordinates, i.e. Πp (z) = (z1 , . . . , zp ).
Proof. We first argue that H, L, and G are independent, and then derive the representation for Hi⊤
in the latter portion of the proof.
Due to the fact that the transformation Z 7→ (H, L, G) is invertible, it is possible to obtain the joint
density of (H, L, G) from the density of Z by computing the matrix Jacobian of the factorization
Z = HLG. (See the references [22], [26], and [27] for more background on Jacobians of matrix
factorizations.) To speak in more detail about the joint density, let Vn×p denote the Stiefel manifold
of n×p matrices with orthonormal columns. Also, let Dp×p denote the set of p×p diagonal matrices,
and Op×p the set of orthogonal p × p matrices. The subset of Op×p with non-negative entries in the
n×p
.
first row will be denoted by O+
Let fH,L,G : Vn×p × Dp×p × Op×p
→ [0, ∞) denote the joint density of (H, L, G), where the base
+
measure is the product of Haar measure on Vn×p , Lebesgue measure on Dp×p , and Haar measure
on Op×p restricted to Op×p
+ . (See the book [28] for background on these measures). From lines
8.8-8.10 in the paper [29], it is known that such a density exists, and that fH,L,G factors according
to
fH,L,G(h, l, g) = fH (h) · fL (l) · fG (g),

(110)

where the function fH denotes the density of H with respect to Haar measure on Vn×p , the function
fL denotes the density of L with respect to Lebesgue measure on Dp×p , and the function fG denotes
the density of G with respect to the restriction of Haar measure to Op×p
+ . This proves that H, L, and
G are independent.
We now prove the representation (109). From line 8.10 in the paper [29], it is known that fH is constant on Vn×p . In other words, the matrix H follows the Haar distribution on Vn×p . Consequently,
Theorem 2.2.1(ii) in the book [28] implies that the rows Hi⊤ can be represented as
L

Hi⊤ = Πp (Ji⊤ )

(111)

where Ji⊤ ∈ Rn is the ith row of a Haar-distributed random matrix in On×n . Furthermore, the rows
Ji⊤ are uniformly distributed on the unit sphere on Rn , and hence can be represented as z/kzk2,
where z ∈ Rn is a standard Gaussian vector.
20

References
[1] C.-H. Zhang and S. S. Zhang. Confidence intervals for low dimensional parameters in high
dimensional linear models. Journal of the Royal Statistical Society: Series B, 76(1):217–242,
2014.
[2] A. Javanmard and A. Montanari. Hypothesis testing in high-dimensional regression under the
Gaussian random design model: Asymptotic theory. arXiv preprint arXiv:1301.4240, 2013.
[3] A. Javanmard and A. Montanari. Confidence intervals and hypothesis testing for highdimensional regression. arXiv preprint arXiv:1306.3171, 2013.
[4] P. Bühlmann. Statistical significance in high-dimensional linear models.
19(4):1212–1242, 2013.

Bernoulli,

[5] S. van de Geer, P. Bühlmann, and Y. Ritov. On asymptotically optimal confidence regions and
tests for high-dimensional models. arXiv preprint arXiv:1303.0518, 2013.
[6] J. D. Lee, D. L. Sun, Y. Sun, and J. E. Taylor. Exact inference after model selection via the
lasso. arXiv preprint arXiv:1311.6238, 2013.
[7] A. Chatterjee and S. N. Lahiri. Rates of convergence of the adaptive lasso estimators to the
oracle distribution and higher order refinements by the bootstrap. The Annals of Statistics,
41(3):1232–1259, 2013.
[8] H. Liu and B. Yu. Asymptotic properties of lasso+mls and lasso+ridge in sparse highdimensional linear regression. Electronic Journal of Statistics, 7:3124–3169, 2013.
[9] V. Chernozhukov, D. Chetverikov, and K. Kato. Gaussian approximations and multiplier bootstrap for maxima of sums of high-dimensional random vectors. The Annals of Statistics,
41(6):2786–2819, 2013.
[10] E. L. Lehmann and J. P. Romano. Testing Statistical Hypotheses. Springer, 2005.
[11] D. A. Freedman. Bootstrapping regression models. The Annals of Statistics, 9(6):1218–1228,
1981.
[12] P. J. Bickel and D. A. Freedman. Bootstrapping regression models with many parameters. In
Festschrift for Erich L. Lehmann, pages 28–48. Wadsworth, 1983.
[13] N. R. Draper and H. Smith. Applied Regression Analysis. Wiley-Interscience, 1998.
[14] P. J. Bickel and D. A. Freedman. Some asymptotic theory for the bootstrap. The Annals of
Statistics, pages 1196–1217, 1981.
[15] S. Bobkov and M. Ledoux. One-dimensional empirical measures, order statistics, and Kantorovich transport distances. preprint, 2014.
[16] A. B. Tsybakov. Introduction to Nonparametric Estimation. Springer, 2009.
[17] L. Wasserman. All of Nonparametric Statistics. Springer, 2006.
[18] B. Laurent and P. Massart. Adaptive estimation of a quadratic functional by model selection.
Annals of Statistics, 28(5):1302–1338, 2000.
[19] R. A. Horn and C. R. Johnson. Topics in matrix analysis. Cambridge University Press, 1991.
[20] R. Bhatia. Matrix Analysis. Springer, 1997.
[21] R. A. Horn and C. A. Johnson. Matrix Analysis. Cambridge University Press, 22nd printing
edition, 2009.
[22] R. J. Muirhead. Aspects of Multivariate Statistical Theory. John Wiley & Sons, 1982.
[23] J. M. Borwein and A. S. Lewis. Convex Analysis and Nonlinear Optimization Theory and
Examples. CMS Bookks in Mathematics. Canadian Mathematical Society, 2000.
21

[24] D. Petz. A survey of certain trace inequalities. Functional analysis and operator theory,
30:287–298, 1994.
[25] S. Boucheron, G. Lugosi, and P. Massart. Concentration Inequalities: A Nonasymptotic Theory
of Independence. Oxford University Press, 2013.
[26] A. M. Mathai. Jacobians of Matrix Transformations and Functions of Matrix Argument. World
Scientific, 1997.
[27] A. Edelman and N. R. Rao. Random matrix theory. Acta Numerica, 14:233–297, 2005.
[28] Y. Chikuse. Statistics on Special Manifolds, volume 174. Springer Science & Business Media,
2003.
[29] A. T. James. Normal multivariate analysis and the orthogonal group. The Annals of Mathematical Statistics, pages 40–75, 1954.

22

