Estimating covariance and precision matrices along subspaces
Zeljko Kereta Timo Klock

Probing the evolution of heavy-ion collisions using direct photon interferometry
    Klaus Reygers^b
    December 30, 2023
================================================================================




We study the accuracy of estimating the covariance and the precision matrix of a D-variate sub-Gaussian distribution along a prescribed subspace or direction using the finite sample covariance. Our results show that the estimation accuracy depends almost exclusively on the components of the distribution that correspond to desired subspaces or directions. This is relevant and important for problems where the behavior of data along a lower-dimensional space is of specific interest, such as dimension reduction or structured regression problems.  We also show that estimation of precision matrices is almost independent of the condition number of the covariance matrix. The presented applications include direction-sensitive eigenspace perturbation bounds, relative bounds for the smallest eigenvalue, and the estimation of the single-index model. For the latter, a new estimator, derived from the analysis, with strong theoretical guarantees and superior numerical performance is proposed.


Keywords: covariance matrix, precision matrix, finite sample bounds, dimension reduction, rate of convergence, ordinary least squares, single-index model, precision matrix 







§ INTRODUCTION


Estimating the covariance =(X- X)(X- X)^⊤ and the precision matrix ^† of a random vector X ∈^D is a
standard and long standing problem in multivariate statistics with applications in a number of mathematical and applied fields.
Notable examples include any form of dimension reduction, such as principal component analysis, nonlinear dimension reduction,
manifold learning, but also problems ranging from classification, regression, and signal processing to econometrics, brain imaging and social networks.

Given independent copies X_1,…,X_N of X, the most widely used estimator is the sample covariance := 1/N∑_i=1^NX_i X_i^⊤, and the inverse thereof.
The crucial question in estimating covariance and precision matrices is to quantify the minimal number of samples N ensuring that for a desired accuracy ε > 0, and a confidence level u > 0, we have

    - _2 ≤ε S_(X),  respectively,  ^† - ^†_2 ≤ε S_^†(X),

with probability at least 1-exp(-u).
Constants S_(X) and S_^†(X) in (<ref>)
describe the dependence of the error with respect to
the distribution of X and properties of , respectively ^†.

In practice however, we are often not directly interested in matrices  or ^† themselves, but rather in
how they, as (bi)-linear operators, act on specific vectors or matrices.
In terms of concentration inequalities, this can be interpreted as developing bounds for
(- )^⊤ and (^† - ^†)^⊤, where  and  are a given pair of (rectangular) matrices.
Bounds of this type in case of sub-Gaussian distributions are the principal subjects of this work.

For any submultiplicative matrix norm ·  perturbations (- )^⊤ and (^† - ^†)^⊤
are bounded by - and ^† - ^†.
This suggests that standard bounds for - and ^† - ^†  (see the overviews in Sections <ref> and <ref>) suffice so long as the distribution is nearly isotropic, i.e. ≈_D.
However, many modern data analysis tasks explicitly rely on anisotropic distributions
because different spectral modalities of the covariance matrix provide crucial, and complementary, information about
the task at hand.
In this case, using norm submultiplicativity and standard bounds for
- and ^† - ^† overestimates incurred errors
because it decouples  and  from their effect on covariance and precision matrices.
Thus, such bounds cannot capture the true behavior of the estimation error.

A typical example that leverages different modalities of (conditional) covariance matrices
are problems that analyze the structure of point clouds, such as manifold learning.
This is because such methods are often prefaced by a linearization step, where the globally non-linear geometry is locally approximated by tangential spaces.
In such a case the conditional covariance of localized data points is anisotropic, i.e. eigenvalues in tangential directions are notably larger than those in non-tangential directions, and the degree of anistropicity increases as the data is more localized <cit.>, facilitating the linearization.

Anisotropic distributions also play an important role in high-dimensional nonparametric regression problems
that use structural assumptions. Denoting Y as the dependent output variable,
a popular example is the multi-index model [Y|X=] = g(^⊤), for a matrix
 with () ≪ D.
The complexity of the underlying nonparametric regression problem can be significantly
reduced by first identifying () and then performing nonparametric regression in
^(). Many methods in the structural dimension reduction literature <cit.>
propose estimators for  that compare the global covariance matrix  (typically assumed to be isotropic),
with conditional covariance matrices X|Y ∈.
Here ⊂(Y) represents a connected level set of the output.
The rationale behind this approach
is that conditioning breaks the isotropicity and induces different spectral modalities
with respect to directions belonging to () and its orthogonal complement. This is
then leveraged to identify () <cit.>.

We showcase the usability of bounds for (^† - ^†)^⊤
on a concrete example in Section <ref> by analyzing the ordinary
least squares estimator ^†X,Y as an estimator for the index
vector  in the single-index model [Y|X=] = g(^⊤).
Our analysis extends studies <cit.> and shows how is the accuracy of the estimator be affected
by anisotropicity of X. Furthermore, an examination of developed estimation bounds motivates a modification of the ordinary least squares estimator
via an approach based on conditioning and averaging in the spirit of structural dimension reduction techniques <cit.>.
We validate the superiority of this modified estimator theoretically and
numerically, and thereby show how a careful tracking of different modalities of the distribution X helps to develop improved
 methods for common data analysis tasks.

Bounds developed in this work also have a few immediate corollaries, which
might be of independent interest. These include eigenspace perturbation bounds similar to
<cit.>, but which are sensitive to the behavior of X
in the direction corresponding to the eigenspace of interest, and a relative bound
for the smallest eigenvalue of  comparable to <cit.>, but without
the isotropicity assumption.



 §.§ State of the art: covariance matrix estimation


The most common bounds for estimating the covariance matrix from finitely many observations
consider sub-Gaussian <cit.> and bounded <cit.> random vectors. They state that with probability at least 1-exp(-u)

    - _2 ≲X_ψ_2^2√(D+u/N)∨D+u/N, 
       - _F   ≲  C_X^2√(u/N),   provided X_2 ≤ C_X a.s,

where A ≲ B means A≤ CB for some universal constant C.
Besides these two cases, researchers have over the years investigated minimal
moment- or tail conditions on X such that bounds similar to (<ref>)
can be achieved. We refer to papers <cit.>
that consider more general classes of distributions.
The most general setting we are aware of is <cit.>
that considers distributions which for universal C, η > 0 satisfy the tail condition

    ( X_2^2 > t)≤ Ct^-1-η,  for  t > C(),

for every orthogonal projection .
Distributions satisfying this condition include log-concave random variables (e.g. uniform distributions on convex sets)
and product distributions, where the marginals have uniformly bounded 4+s moments
for some s > 0.

Our bounds show that the sample covariance estimator  automatically and implicitly
adapts to (), which serves as a complexity parameter of the estimation problem.
However, in the regime N<() the sample covariance is rank deficient and the estimation is in general 	not possible.
Instead, structural assumptions, such as sparsity or bandedness, are needed to reduce the effective complexity of the problem and allow consistent estimation.
These assumptions can be leveraged by regularized estimation techniques, which include banding <cit.>, thresholding <cit.>, or penalized likelihood estimation <cit.>.
We refer to <cit.> for detailed reviews of existing methods.

We also mention <cit.> that considers
Gaussian random variables taking values in general Banach spaces, and <cit.> which
analyze the concentration of spectral projectors of the covariance matrix, and (bi)-linear forms
thereof, for Gaussian random variables in a separable Hilbert space. This is conceptually related
to our work as the bounds take into account information about relevant spectral projectors, and their
interplay with vectors to which they are applied. This has recently been extended to the
estimation of smooth functionals of covariance operators <cit.>
for Gaussians in separable Hilbert spaces.



 §.§ State of the art: precision matrix estimation


Estimation of the precision matrix is relevant for many problems, ranging from simple tasks
such as data transformations (e.g. standardization Z := ^-1/2X), to  applications
that include linear discriminant analysis, graphical modeling, or
complex data visualization.
Furthermore,  precision matrix encodes information about partial correlations of features of X. Namely, if X follows a Gaussian (or paranormal)
distribution, the ij-th entry of ^† is zero if the i-th and the j-th
feature are conditionally independent.

The inverse ^† of the sample covariance, constructed from N independent copies of a mean zero random vector X ∈^D, is a well-behaved estimator of ^† as N→∞
and D is considered fixed <cit.>.
In such a case bounds for the precision matrix can be obtained by using general perturbation bounds for the Moore-Penrose inverse.
One of the first such bounds <cit.> states that for ∈^d_1 × d_2, and an additively perturbed matrix =  +, we have

    ^† - ^†   ≤ωmax{^†_2^2,^†_2^2}, and 
    ^† - ^†   ≤ω^†_2^†_2,  if () = (),

where · is any unitarily invariant norm, and ω is a small universal constant <cit.>.
Recent studies <cit.> examine the influence of the
perturbation in greater detail, implying  the bound

    ^† - ^†_F    ≤min{^†_2 ^†_F, ^†_2^†_2},
       if () = () = min{d_1,d_2}.

Using these general perturbation bounds it is easy to derive first concentration bounds for the precision matrix.
For example, assume X is sub-Gaussian and that the number of independent data samples satisfies N =C4^k(D+u)X_ψ_2^4/λ_()()^2 for some universal C>0 and k∈.
Equation (<ref>) and Weyl's bound <cit.> imply

    λ_()() ≥λ_()() - - _2 ≥λ_()()(1-2^-k),

and consequently ^†_2 ≤ (1-2^-k)^-1^†_2.
The perturbation bound (<ref>)
and covariance bound (<ref>) give with probability 1-exp(-u)

    ^†-^†_2 ≲^†_2^2 - _2 ≲^†_2^2 X_ψ_2^2√(D+u/N),

where the higher order term in (<ref>) vanishes in the applicable regime
N ≥ C(D+u)X_ψ_2^4/λ_()^2().
While this effectively provides bounds as soon as the covariance
perturbation is bounded, in this work we show that (<ref>) overestimates the error by assigning a quadratic scaling of ^†_2
(Corollary <ref>).

Moreover, we are not aware of precision matrix bounds that take into account the specific nature of
the perturbation
nor of bounds for
(^† - ^†)_2 that are sensitive to objects of interest.

If () grows with N, then  is not a consistent estimator
of  and thus the precision matrix cannot be estimated well by inverting the sample covariance matrix.
Various families of structured precision
matrices have been studied to mitigate these issues, motivated by applications in genomics, finance, and other fields.
Dominant assumptions are sparsity and bandedness, which are exploited
through the use of regularized estimators. Algorithms for estimating ^† based on
regularization include computing ^† column by column through entry-wise Lasso
<cit.>, constrained ℓ_1 minimization <cit.>,
 adaptive ℓ_1 minimization <cit.>, ℓ_1 regularized score matching
 <cit.>, or ridge regressors <cit.>. See <cit.>
 for comprehensive overviews.



 §.§ Overview and contributions


Throughout, X ∈^D is a sub-Gaussian random vector with X := X -  X and = X,
and X_1,…,X_N are independent copies of X.
Sub-Gaussians are a broad class of light-tailed distributions that have received increasing attention in recent years and are used in many branches of probability and statistics.
We define finite sample estimators of X and  by
_X := N^-1∑_i=1^NX_i and := N^-1∑_i=1^N(X_i-_X) (X_i-_X)^⊤.

Let ∈^d_1× D, ∈^d_2× D be matrices determining a direction, subspace, or generally an object of interest.
We can summarize our findings as follows.


  * In Section <ref> we show that with probability at least 1-exp(-u)

    (- )^⊤_2 ≲X_ψ_2X_ψ_2(√(d_A + d_B +u/N)),

where (t) = t ∨ t^2, d_A := (), and d_B := ().
This is similar to <cit.> but replaces the sub-Gaussian norm X_ψ_2
by direction/subspace dependent quantities X_ψ_2 and X_ψ_2.

  * In Section <ref> we show  that with probability at least 1-exp(-u),
we have

    (^† - ^†)^⊤_2 ≲^†X_ψ_2^†X_ψ_2√(()+u/N),

provided N ≥ C √(^†)X_ψ_2^4, for an absolute constant C>0,
where √(^†)X is the standardization of a random variable X, with √(^†)X = _D.

  * In Section <ref> we show stronger bounds for (<ref>) and (<ref>) in case of bounded random vectors.




As we point out in Corollary <ref>,
the bound (<ref>) is  interesting even
when =  = _D, particularly in view of  general perturbation
bounds such as (<ref>).
Consider for example a random vector X for which the sub-Gaussian norm is a good proxy for the variance, i.e.
so that ^†X_ψ_2≈√(^†_2) holds.
(This is the case for example if X ∼(,), or more generally for strict sub-Gaussians <cit.>).
The right hand side of (<ref>) then scales linearly
in ^†_2, whereas (<ref>)
shows a quadratic behavior. This has a significant impact for ill-conditioned
covariance matrices and implies that inverting the sample covariance exhibits better conditioning than inverting a general matrix perturbation.
The same effect is observed if  and  are arbitrary.


Two applications of bounds (<ref>) and (<ref>)
are presented. In Section <ref> we use the covariance bound
(<ref>) to establish a bound for perturbations of eigenspaces of the covariance matrix that is sensitive to the distribution of the random vector in the eigenspace of interest.
This is relevant for example when estimating manifolds from unlabeled point cloud data,
see <cit.>.

In Section <ref> we use
(<ref>) to establish sharp concentration bounds for single-index model estimation.
In this model a response Y ∈ is assumed to follow the regression model
[Y|X] = f(^⊤ X), and the task is to estimate the unknown vector  using a finite data set {(X_i,Y_i):i=1,…,N}.
A common estimator is the normalized ordinary least squares vector, for which we provide
direction-sensitive concentration bounds.
Furthermore, our analysis yields an insight into how the estimator can be improved by a simple and straightforward procedure based on conditioning and averaging. This is presented in Section <ref>.

Most proofs are deferred to the Appendix for the sake of brevity and clarity.


 §.§ General notation


We denote [M] = {1,…,M}, a ∧ b = min{a,b}, a ∨ b = max{a,b},
and we may use the auxiliary function (t) = t ∨ t^2.
For a real, symmetric matrix ∈^d× d we denote by λ_1() ≥…≥λ_d()
the ordered set of its eigenvalues, and by _1(),…,_d() the corresponding eigenvectors.
·_2 denotes the spectral norm of a matrix, and the Euclidean norm of a vector, and ·· is the dot product.
·_F is the Frobenius norm. _D∈^D× D is the identity matrix.
^D-1 is the unit sphere in ^D. For any random vector X
we write X := X -  X. For p ≥ 1 and a random variable X we define the Orlicz norm

    X_ψ_p := inf{s > 0 : exp(X/s^p)≤ 2}.

The definition extends to random vectors X ∈^D by

    X_ψ_p := sup_∈^D-1^⊤ X_ψ_p < ∞.

We only use p=1 (sub-exponential) and p=2 (sub-Gaussian).
If Ω is a finite set, Ω denotes its cardinality. If Ω
is an interval, Ω denotes its length.
Throughout the paper C is a placeholder for a positive universal constant that may have a different value on each occurrence, even in the same line of text.
Furthermore, we sometimes use A≲ B instead of A≤ CB.






§ COVARIANCE MATRIX ESTIMATION


In this section we present bounds for covariance  and eigenspace estimation that are sensitive to the distribution of a given random vector in directions of interest.
The following matrix concentration bound is the fundamental tool of our analysis.


Let ∈^d_1 × D, ∈^d_2 × D and
define d_A=(), and d_B=().
Then for all u>0, with probability at least 1 - exp(-u) we have
for (t) = t ∨ t^2

    (- )^⊤_2 ≲X_ψ_2X_ψ_2(√(d_A+d_B+u/N)).




An analogous result holds for almost surely bounded random vectors by using a different concentration inequality.
In that case ·_ψ_2 can be replaced by a bound for the Euclidean norm of X, and the dimensionality
does not appear in the requirement on N. We will return to this point at the end of Section <ref>.


The proof of Lemma <ref> by and large follows along the lines of traditional concentration results.
Indeed, if =, the result would follow by applying <cit.> to the random vector X, along with some minor adjustments to account for the ranks.
When ≠, a somewhat more careful tracking of the behavior of the random vector X, with respect directions induced by matrices  and  is needed.
In the end, as (<ref>) suggests, the payoff is that the error rate scales only with components of X along those directions.


Some works that consider concentration inequalities for sub-Gaussian random variables do so with respect to the effective rank, defined as r() = ()/_2, see for instance <cit.> or <cit.>.
Effective rank cannot exceed the true rank of a matrix, and unlike (), it is less affected by small eigenvalues.
It can thus be a useful surrogate for approximately low dimensional distributions, and can in those cases be used to provide informative estimation bounds
even if N ≪().

On the other hand, r() is not as useful for precision matrix estimation, since inverting a matrix reverses the ordering of the eigenvalues. Thus, r() should be replaced with r(^†).
Furthermore, our current proof technique for precision matrix estimation
requires () = (), which immediately implies N ≥().
To provide a unified framework, we decided to abstain from using the effective rank in our results.



Applying Lemma <ref> we can reconstruct known error rates in case of low-rank  distributions in high-dimensional ambient spaces.

For all u>0  with probability at least 1 - exp(-u) we have

    - _2 ≲X_ψ_2^2 √(()+u/N),  where (t) = t ∨ t^2.



Lemma <ref> also has an immediate effect on the estimation of eigenvectors and eigenspaces.
Denote by _i,l() := ∑_k=i^l_k()_k()^⊤ the orthoprojector onto the space spanned by i-th to l-th eigenvectors of , and let
_i,l() be the corresponding finite sample estimator.
Moreover, denote _i,l() := _D - _i,l(),
and I_1I_2 := inf_t ∈ I_1, t' ∈ I_2t-t' for I_1,  I_2 ⊂.


Let i ≤ l ∈, and define

    δ_il   = [λ_l(),λ_i()][-∞,λ_l+1()]∪[λ_i-1(),+∞], 
            with λ_0() := ∞, λ_D+1() = -∞.

For any u>0, with probability at least 1 - exp(-u) we have, with (t) = t ∨ t^2,

    _i,l()_i,l()_2 ≲_i,l()X_ψ_2X_ψ_2/δ_il√(()+u/N).




Davis-Kahan Theorem in <cit.> gives

    _i,l()_i,l()_2    ≤π/2_i,l()( - )_i,l()_2/δ_il≤π/2( - )_i,l()_2/δ_il.

The claim now follows by applying Lemma <ref> with = _D and = _i,l().


Typical bounds for eigenspace perturbations _i,l()_i,l()_2 take the specific eigenspace into account only through the denominator, whereas the numerator relies on squared terms of the form X_ψ_2^2 in the sub-Gaussian case, or a bound for X_2^2 in the bounded case.
Expression (<ref>) is thus beneficial if _i,l()X_ψ_2 is smaller than X_ψ_2, as it
provides a sharper estimate.

In order to ensure δ_il > 0, the covariance matrix  must have a population eigengap, that is, δ_il^* := (λ_i-1()-λ_i()) ∧ (λ_l() - λ_l+1()) > 0,
and sufficiently many samples are required in order to stabilize δ_il around δ^*_il.
The latter is typically achieved by first using Weyl’s bound <cit.>,
giving λ_j() -λ_j()≤-_2 for all j ∈ [D],
and then applying a concentration bound for -_2. A consequence however
is that Proposition <ref> is only informative if we have
sufficiently many samples with respect to δ^*_il and X_ψ_2.
Thus, the estimation error is no longer sensitive to the eigenspace of interest.

Preserving the dependence on _il()X_ψ_2, instead of on X_ψ_2,
requires the use of relative eigenvalue bounds. Such bounds have been recently provided in <cit.>
under the assumption that X is strictly sub-Gaussian <cit.>,
i.e. for some K > 0, and for arbitrary matrices ∈^k× D, X satisfies

    X_ψ_2^2 ≤ K X_2.

This asserts that the sub-Gaussian norm is a good proxy
for variances of one-dimensional marginals of the random vector X, and is for instance
satisfied for X ∼(0,) with K = 1. For such distributions we obtain the following
concentration of the sample eigengap.



Assume X satisfies (<ref>) for some K > 0.
Let i,l∈ with 1≤ i ≤ l≤(), and {i,l}≠{1,()}.
Let (t) = t ∨ t^2. There exists a constant C_K>0, depending only on K, so that whenever

    N > C_K(()∨ u)(λ_l+1()/λ_l()-λ_l+1()∨λ_i-1()/λ_i-1()-λ_i()),

(with conventions λ_0() = ∞, λ_()+1() = 0, and ∞/∞ = 0)
for any u>0 with probability at least 1 - 2exp(-u) we have δ_il≥δ_il^*/2.


The case {i,l} = {1,()} is equivalent to asking whether () = ().
This holds for N ≥() if the law of X has a density which is absolutely
continuous with respect to the Lebesgue measure on (), because X_1,…,X_N are almost surely linearly
independent.
It also holds with probability at least 1-exp(-u) if X is sub-Gaussian as soon as N > C_K(() + u),
provided (<ref>) holds, and if (<ref>) does not hold then we need
N > C(() + u)√(^†)X_ψ_2^4.

Lemma <ref> now allows to refine Proposition <ref> with a population eigengap.


Assume X satisfies (<ref>) for some K > 0. Let 1≤ i ≤ l ≤()
with {i,l}≠{1,()}. Let (t) = t ∨ t^2.
There exists C_K>0, depending only on K, so that whenever
N satisfies (<ref>) for any u>0,
with probability at least 1 - exp(-u) we have

    _i,l()_i,l()_2    ≤  C_K_i,l()X_ψ_2X_ψ_2/δ_il^*√(()+u/N)
       ≤ C_K√(Kλ_i())X_ψ_2/δ_il^*√(()+u/N)



The first inequality follows immediately by first conditioning on events in Proposition <ref>, Lemma <ref>,
and then using the union bound (probability 1-3exp(-u) can be adjusted by adjusting C_K). The second inequality in (<ref>)
comes from additionally using (<ref>) and _i,l()X_2 = λ_i().

Recently, <cit.> showed a useful alternative

    _i,l()_i,l()_F ≲(D^1/2-_2 ∧- _F)/δ^*_il.


Thus, with (<ref>) we do not have to stabilize the sample eigengap, and the eigenspace perturbation bound (<ref>) can be
used for arbitrary N≥ 1.
A natural question to ask is whether Proposition <ref>
holds if δ_il is replaced with δ^*_il.
The following example strongly suggests this is not the case.


Assume that Proposition <ref> holds with δ_il^∗ in place of δ_il.
Let X ∼(0,), where = ∑_i=1^D-1_i_i^⊤ + η^2 _D_D^⊤,
for η < 1.
Notice that in this case X_ψ_2 = 1 and _D^⊤ X_ψ_2 = η, by the definition of the sub-Gaussian norm.
For any N≥ 1 with probability at least 1-exp(-u) we would have

    _D,D()_D,D()_2 ≲η/1-η^2(D+u/N).

In particular, we have an increasingly small estimation error as η→ 0 by using only one data sample, i.e. by estimating _D,D() based
on the eigendecomposition of a rank one matrix XX^⊤.







§ PRECISION MATRIX ESTIMATION


In this section we investigate directional estimates of the precision matrix
^† through the empirical precision matrix ^†,
analogously to results in Section <ref>.



Let ∈^d_1× D, ∈^d_2× D.
There exists a uniform constant C>0 such that if N > C(() + u)√(^†)X_ψ_2^4
for any u>0 with probability at least  1-exp(-u) we have

    (^† - ^†)^⊤_2 ≲^†X_ψ_2^†X_ψ_2√(()+u/N).



Let us comment on the implications of Theorem <ref>.
First, we note that √(^†)X_ψ_2 is the sub-Gaussian
norm of the standardization of X.
Provided X is strictly sub-Gaussian for some K > 0, see (<ref>),
 we have √(^†)X_ψ_2^2 ≤ K
√(^†)X_2 = K. It follows that for such distributions
√(^†)X_ψ_2 has a negligible effect on estimating precision matrices.

Second, similar to Lemma <ref>, the bound (<ref>) depends only on components of X induced by  and . Since the sub-Gaussian norm can be interpreted as a proxy for the variance, this improves non-directional bounds
whenever the eigenvalues of ^†^⊤ and ^†^⊤ are
small compared to those of ^†.

Third, as soon as N > C(() + u)√(^†)X_ψ_2^4,
the estimation rate in (<ref>) is similar to
the covariance estimation rate in Lemma <ref>.
That is, assume we are trying to estimate ^† with the sample covariance matrix of the random vector Z = ^†X through iid. copies Z_i=^†X_i.
In that case Lemma <ref> gives precisely the bound (<ref>).
This should come as a bit of a surprise, since it implies that estimating the precision matrix through the inverse of the sample covariance has the same theoretical guarantees as if we had access to a random vector Z  whose covariance is exactly ^†.
To further stress this point, we now compare Theorem <ref>
for =  = _D with the bound

    ^† - ^†_2 ≲^†_2^2 _2√(D+u/N),

which was derived in Section <ref> using general perturbation bounds for the matrix inverse.


There exists a uniform constant C>0 such that if N > C(() + u)√(^†)X_ψ_2^4
for any u > 0 with probability at least 1-exp(-u) we have

    ^† - ^†_2 ≲^†X̃_ψ_2^2 √(()+u/N).



Assuming the squared sub-Gaussian norm is a good proxy for the variance, i.e. if
(<ref>) holds, the right hand side in (<ref>)
becomes K ^†_2√((()+u)/N). Compared to (<ref>),
this shows that using a general perturbation bound overestimates the influence the matrix condition
number _2^†_2 has on precision matrix estimation.
The discrepancy between these two results suggests that the finite sample covariance estimator induces
a specific type of a perturbation that performs a form of regularization when estimating the inverse.

An immediate corollary is a relative bound for the smallest eigenvalue of .

Let  d := ()
and assume X satisfies (<ref>) for some K > 0. There exists a constant C>0 such that provided
N > C(d+u)√(^†)X_ψ_2^4,
we have with probability at least 1-exp(-u)

    |λ_d()/λ_d() - 1|≲ K√(d+u/N).

If (<ref>) does not hold the right hand side of (<ref>) is, up to absolute constants, replaced by λ_d()^†X_ψ_2^2√((d+u)/N).

The bound (<ref>) is similar to <cit.>, which holds for isotropic
random variables with finite fourth order moments. Moreover, it is possible to avoid the dependence on d if the number of samples N
is large compared to
∑_i=1^d-1λ_i()/λ_i() - λ_d(),
by using relative eigenvalue bounds from <cit.>.



  
Numerical validation.
To validate the results of Theorem <ref> we consider
X ∼(0,) for two types of covariance matrices ∈^10× 10:


  * Set =^⊤, for ∈^10×10 sampled uniformly at random from the space of orthonormal
matrices, and =1,1,1,1,1,ν,ν,ν,ν,ν. We consider ν = 10^-j+1 for j ∈ [10], which
implies that the matrix condition number _2^†_2
ranges from 1 to 10^9.

  * Set _i,j = ν^i-j, with ν∈{0.5,0.55,…,0.9,0.95}. This is a common model for distributions
where entries of X correspond to values of a certain feature at different time stamps. It leads to correlated entries
when the time stamps are close by, i.e. when i-j is small.

We choose matrices  and  by sampling an orthoprojector
 uniformly at random with () = 3 and setting
= _D -. We repeat each experiment 100 times and report the averaged relative errors in Figure <ref>.
Different lines of the same color correspond to different values of ν.






The estimation error shows the expected N^-1/2 rate.
Moreover, directional errors show  a clear dependence on the directional spectral norm.
On the other hand, since we are sampling Gaussians the squared sub-Gaussian norm is a proxy for the variance of the given random vector, and the results confirm that the error term in Theorem <ref> indeed scales with the corresponding (directional) sub-Gaussian norm.
Lastly,  although the condition number of  depends on ν, the specific value of ν does not affect how accurate ^† is, as predicted by the theory.



  
Bounded random vectors.
As mentioned in Remark <ref>, a stronger form of direction dependent covariance and precision matrix estimation bounds hold for bounded random vectors.
The proofs for the bounded case follow along similar lines as for the sub-Gaussian case, except for the use of a slightly different probabilistic argument.
We now state only the results for the estimation of covariance and precision matrix, since the remaining bounds follow by analogy.


Let ∈^d_1 × D, ∈^d_2 × D.
Assume  X_2≤ C_A, X_2≤ C_B
almost surely.
Then for any u>0 with probability at least 1 - exp(-u) we have

    (- ) ^⊤_F ≲ C_A C_B√(1+u/N).

Assume
^†X_2≤ C_A^†, ^†X_2≤ C_B^†,
√(^†)X_2^2≤Θ almost surely.
There exists C>0 such that provided N > C(1+u)Θ^2,  for any u>0 with probability at least 1-exp(-u) we have

    (^† - ^†)^⊤_F ≲ C_A^† C_B^†√(1+u/N).








§ APPLICATION TO SINGLE-INDEX MODEL ESTIMATION


In this section we use the results of Sections <ref> and <ref> to establish concentration bounds
for estimating the index vector  in the single-index model [Y|X=x] = g(^⊤ x).
Moreover, directional terms that arise from using directional dependent matrix concentration bounds provide an insight into how to further improve the performance of the standard estimator.
The second half of the section is thus devoted to describing this strategy (which is based on splitting up the data, conditioning and averaging), proving the error estimates, and providing numerical evidence to show and examine the claimed performance gains.

Throughout the section, we use a directional sub-Gaussian condition number

    κ(, X) := ^†X_ψ_2^2X_ψ_2^2,    is an orthoprojector.

This quantity can be seen as a restricted condition number ^†_2 _2
for the matrix = X. Indeed, for strict sub-Gaussians, i.e. those satisfying (<ref>) for some K,
the two are equal up to a factor depending on K.




 §.§ Ordinary least squares for the single-index model


Single-index model (SIM) is a popular semi-parametric regression model that poses the relationship between the features X ∈^D  and responses Y ∈
as [Y|X] = f(^⊤ X), where f is an unknown link function and ∈^D-1 is an unknown index vector.
SIM was developed in the 80s and 90s <cit.> as an extension of generalized linear regression that does not specify the link function,
and which could thus avoid errors incurred by model misspecification.
Common applications are in econometrics <cit.> and signal processing under
sparsity assumptions on the index vector <cit.>.
It has been shown, e.g. in  <cit.>, that (in certain scenarios) the minimax estimation rate of SIM equals that of nonparametric univariate regression.

Methods for estimating the SIM from a finite data set {(X_i,Y_i) : i ∈ [N]} often first construct an approximate index vector
, and then use nonparametric regression on {(^⊤ X_i, Y_i) : i ∈ [N]} to estimate the link function.
With such an approach the generalization error of the resulting estimator depends
largely on the error incurred by estimating the index vector. Thus, the construction of  becomes the critical point.

An efficient approach, which first appeared in
<cit.>, and later in modified forms in <cit.>,
is to solve the ordinary least squares (OLS) problem

    (ĉ, ) = _c ∈, ∈^D∑_i=1^N(Y_i - c - ^⊤(X_i- _X))^2,  where _X = ∑_i=1^N X_i/N,

and then set := /_2.
It was shown in <cit.> that √(N)(- ) is  asymptotically
normal with mean zero, provided X has an elliptical distribution, f is a non-decreasing function that is strictly increasing
on some non-empty sub-interval of the support of ^⊤ X.
Assuming only ellipticity of X and statistical independence
of Y and X given ^⊤ X, the population vector :=^†X,Y
is still contained in {}, see e.g. <cit.>.
Provided ≠ 0, in these cases the direction := /_2 equals the index vector  up to sign.
Thus,  is in many cases a consistent estimator of the index vector , and under certain conditions it converges with an N^-1/2 rate.

The minimal ·_2-norm solution of (<ref>) admits a closed form

    ĉ = _Y,  = ^†, with _Y :=∑_i=1^NY_i/N,  :=∑_i=1^N(X_i - _X)(Y_i - _Y)/N.

Using the results of the previous two sections we can show a direction
dependent concentration bound for the vector .


Let Y ∈ be sub-Gaussian. Denote  =^⊤, := _D -, and κ_ = κ(,X) ∨κ(,X).
There exists C>0 such that provided N > C(()+u)√(^†)X_ψ_2^4,
for any u>0, with probability at least 1 - exp(-u) we have

    ( - )_2 ≲Y_ψ_2^†X_ψ_2√(κ_)√(()+u/N),
       ( - )_2 ≲Y_ψ_2^†X_ψ_2√(κ_)√(()+u/N).


For the normalized vector, respectively the directions, the following bound holds.


Assume the setting of Lemma <ref>.
There exists a constant C>0 such that
for any u>0, with probability at least 1-exp(-u)  we have

    - _2 ≲Y_ψ_2^†X_ψ_2√(κ_)/_2√(()+u/N),  provided that 
       N > C(()+u)(√(^†)X_ψ_2^4 ∨Y_ψ_2^2^†X_ψ_2^2κ_/^2_2).



As mentioned before, under certain conditions we have =, where  is the index vector in the given SIM.
In such cases Corollary <ref> confirms that  is a consistent estimator of  and achieves a N^-1/2 convergence rate, which has been observed in previous works.
To provide some context for our result we give a brief
description of the two most popular strategies for estimating .

The first group of methods can be distinguished by their simplicity and efficiency
as they estimate the index vector using empirical estimates of
first and second order moments of random variables X and Y. The studied OLS estimator (sometimes
also referred to as the average derivative estimator <cit.>)
belongs to this group. Inverse regression based techniques also
fall into this category, see e.g. <cit.> or the review <cit.>,
and their convergence rate usually equals N^-1/2 and is thus on par with OLS.
A drawback caused essentially by the simplicity of these methods is that theoretical guarantees typically require X to be an elliptical distribution, and
sometimes a form of monotonicity is needed to avoid issues that arise when dealing with symmetric functions such as ↦ (^⊤)^2.

The second, more sophisticated but also computationally heavier, class of methods is based on
solving more complicated optimal programs for the recovery of , to which a closed form solution does not exist. This includes non-parametric methods for sufficient dimension reduction, e.g. <cit.>, aiming
at estimating the gradients of f(^⊤ X) at all training samples X_i; methods based on combining
index estimation and monotonic regression, e.g. <cit.>;
and methods that simultaneously estimate the index vector and use spline interpolation <cit.>, to name just a few.
For these methods convergence rates up to N^-1/2 (sometimes slower) can typically be proven, and
the theory requires less stringent assumptions than those exhibited by the first class of methods.


A new insight in Corollary <ref>, which has
not been emphasized before, is the directionally-dependent influence of the spectrum of  on the estimation guarantee.
To make this more precise, we now consider a special case when X is strictly sub-Gaussian,
the link function f is Lipschitz-smooth, and the index  is an eigenvector of . We note that similar results hold
if  is an approximate eigenvector of , in the sense that
^⊤^†≈ (^⊤)^-1.



Let Y be sub-Gaussian, X be strictly sub-Gaussian (i.e. satisfying (<ref>) for some K > 0)
and assume Y = f(^⊤ X) + ζ for ζ = 0 with σ_ζ:=ζ_ψ_2<∞.
Assume consistent estimation, i.e. =, and that = σ_^2.
Then there exists C_K>0, depending only on K, so that if

    N > C_K(()+u)(L+σ_ζσ_^-1)^2κ(,X)/_2^2,

for any u>0 we have with probability at least 1-exp(-u)

    - _2 ≲(Lσ_ + σ_ζ)√(^†_2κ(,X))/_2√(()+u/N).



Corollary <ref> shows that the variance in the direction of the index vector,
σ_^2 = ^⊤ X, and the variance in orthogonal directions, quantified through the spectrum of ^†, influence index vector estimation in a significantly different manner.
Namely, as long as σ_≫σ_ζ,
smaller σ_ has a provably beneficial effect on estimation accuracy, whereas small
non-zero eigenvalues of  (i.e. large eigenvalues of ^†), corresponding to eigenvectors in (), can only worsen the accuracy.
A similar observation can be made when inspecting the asymptotic covariance
of √(N)(- ), which after using = σ_^2 and <cit.>, is given by

    σ_^2/f(^⊤ X),^⊤ X_=:T_1^†(Ỹ - X̃^⊤)X̃^†_=:T_2.

The scalar factor T_1 is comparable to _2^2,
and the matrix T_2 shows that the variance of the estimator grows with large eigenvalues
of ^†.
On the other hand, the variance of the residual Ỹ - X̃^⊤ decreases
if the accuracy of the linear fit of X̃^⊤≈ f(^⊤ X)- f(^⊤ X) improves,
which can typically be observed under monotonic links f and decreasing σ_^2 = ^⊤ X.


We will now use this observation as a guiding principle for developing a modified estimator
that splits the data into subsets with small σ_'s, computes corresponding OLS vectors on each subset,
and then combines local estimators into a global estimator by weighted averaging.



 §.§ Averaged conditional least squares for the single-index model




We now study the just discussed alternative procedure, where we first split the data into
subsets, aiming to reduce the variance of the data distribution in the direction of the index vector, and then compute and average out the estimators
from each subset. Since we have no  a priori knowledge about the index vector, constructing such a partition seems challenging.
However, when the link function is monotonic the partitioning is induced by a decomposition
of (Y) into equisized intervals, see Figure <ref>.
For the sake of simplicity, in the following we assume Y ∈ [0,1) holds almost surely.

For ℓ∈[J] let _J,ℓ := [ℓ-1/J,ℓ/J) denote equisized regions partitioning [0,1).
Furthermore, define herein called level-sets _J,ℓ = {(X_i, Y_i) : Y_i ∈_J,ℓ},
which induce a partition of the data set into J subsets based on the responses. Then estimate
 according to the following algorithm.


  *  Solve (<ref>) on each subset _J,ℓ, by computing ĉ_J,ℓ∈ and _J,ℓ∈^D.

  *  Define the empirical density ρ̂_J,ℓ := _J,ℓ/N, set the thresholding parameter α > 0, and compute the averaged outer product matrix

    _J = ∑_ℓ = 1^J1_[α J^-1, 1](ρ̂_J,ℓ)ρ̂_J,ℓ_J,ℓ_J,ℓ^⊤.


  *  Use the eigenvector corresponding to the largest eigenvalue of _J, denoted as _1(_J), as an approximation of the index vector .


The parameter α is used to promote numerical stability by suppressing the contributions of sparsely populated subsets.
In other words, we only keep those level-sets whose empirical mass behaves as if Y were uniformly distributed over (Y) (which can in some problems be achieved by a suitable transformation of the responses).
The parameter J on the other hand defines the number of sets we use in the partition of the given data set,
and dictates the trade-off between ^⊤ X|Y ∈_J,ℓ
and the number of samples |_J,ℓ| in a given level set.

For a random vector Z we denote a conditional random vector Z_J,ℓ:=Z|Y∈_J,ℓ.
Note that Z_J,ℓ inherits sub-Gaussianity of Z provided
(Y ∈_J,ℓ) > 0, see Lemma <ref>.
We now  analyze the approach under the following assumption:


  *  _J,ℓ:=X|Y ∈_J,ℓ^†X,Y|Y ∈_J,ℓ∈{} for all ℓ∈ [J].

Assumption <ref> is not particularly restrictive.
For example, it can be shown that <ref>
holds if X is elliptically symmetric, which is a standard assumption when using the OLS
functional (<ref>), and if the function noise Y-[Y|X] is statistically independent of X given ^⊤ X <cit.>.



Assume <ref> holds and that Y ∈ [0,1) almost surely. Let J > 0, α > 0, and assume we
are given N iid. copies of (X,Y). Denote := ^⊤, := _D -, _J,ℓ := X|Y ∈_J,ℓ,
and

    κ_J,ℓ:=κ(,X_J,ℓ) ∨κ(, X_J,ℓ) , and  K_J,ℓ :=√(_J,ℓ^†)X_J,ℓ_ψ_2.

Let I_J:= {ℓ∈ [J] : ρ̂_J,ℓ >α J^-1} be the index set containing the indices of active level-sets
and set K_J := max_ℓ∈ I_JK_J,ℓ.
There exists C > 0 such that, if

    N > CK_J^4J(()+log(J)+u)/α,

there exists a sign s ∈{-1,1} so that for any u>0, with probability at least 1-exp(-u) we have

    s _1(_J) -_2^2 ≤ε_N,J,u∑_ℓ∈ I_Jρ̂_J,ℓκ_J,ℓ_J,ℓ^†X_J,ℓ_ψ_2^2/∑_ℓ∈ I_Jρ̂_J,ℓ(b_J,ℓ_2^2 -ε_N,J,uκ_J,ℓ_J,ℓ^†X_J,ℓ_ψ_2^2),

where 
    ε_N,J,u := C() + log(J) + u/α J N.

The same guarantee holds when replacing the denominator in (<ref>) with λ_1(_J).

The error rate in Theorem <ref> is dictated by ε_N,J,u
and there are two ways to interpret Theorem <ref>.
First, for a fixed parameter J all terms in (<ref>) and (<ref>) (except N and ε_N,J,u) are constants, and we obtain a  N^-1/2 convergence rate
for the non-squared error as soon as the sample size is sufficiently large. This is
the same rate as the one achieved with the standard OLS estimator (<ref>).

Parameter J however provides additional flexibility and an intriguing option is to select J as a function that grows with N while still complying with (<ref>).
Namely, using u = log(J) and J≍ N/log(N), and plugging into (<ref>), implies that
a √(log(N/log(N))log(N))N^-1≤log(N)N^-1 rate is achievable, so long as the remaining terms involving parameter J are balanced.

To give an illustrative example of such a case, we consider the following idealized
setting. We assume the noise-free regime Y = f(^⊤ X) with strictly monotonic
link in the sense that

    ^⊤ X, f(^⊤ X)| Y ∈_J,ℓ > L Var(^⊤ X|Y∈_J,ℓ)  for some   L > 0.

Furthermore, conditional random variables
X_J,ℓ are assumed strictly sub-Gaussian, where the strict sub-Gaussianity constant
K > 0 in (<ref>) is independent of J and ℓ,
and we model local covariance matrices as

    X̃_J,ℓ = C_ℓ J^-2 +    with  c_1 ≤ C_ℓ≤ c_2,

where c_1,c_2 are universal constants independent of J.
Condition (<ref>) implies that partitioning the data into J level-sets results in a reduction of the variance along the direction of the index
vector, while not affecting the variance in directions orthogonal to the index vector.

With strict sub-Gaussianity
and (<ref>), we have K_J ≤√(K), κ_J,ℓ≤ (c_2/c_1) K^2,
_J,ℓ^†X_J,ℓ_ψ_2^2 ≤ KJ^2/c_1,
and _J,ℓ^†X_J,ℓ_ψ_2^2 ≤ K.
Then, using u = log(J) and
J = τ N/log(N) for some τ > 0,
the result (<ref>) can be written as

    s _1(_J) - _2^2 ≤C_K ()/ατ∑_ℓ∈ I_Jρ̂_J,ℓ(b_J,ℓ_2^2 - τC_K(()+1)/α)log^2(N)/N^2,

where C_K > 0 depends only on K, c_1, and c_2. Furthermore, using
Assumption <ref>, and (<ref>), we have
X,Y|Y ∈_J,ℓ = a^⊤ X, Y | Y ∈_J,ℓa,
and with the asserted monotonicity (<ref>) we get

    b_J,ℓ_2    ≥X,Y|Y ∈_J,ℓ^⊤^†_J,ℓX,Y|Y ∈_J,ℓ/X,Y|Y ∈_J,ℓ_2
       > LVar(a^⊤ X|Y∈_J,ℓ) a^⊤^†_J,ℓ a = L.

Plugging this bound into (<ref>), and choosing τ∈ (0,α L^2/(2C_K(()+ 1))),
we thus obtain a log^2(N)/N^2 convergence rate for the squared error.

The described idealized setting can be a reasonable approximation of reality
for strictly monotonic link functions f in the noise-free regime, since
(<ref>) can be motivated by the observation that in this case slicing the distribution of X based on Y
reduces the variance of the data along the direction of the index vector,
as shown in Figure <ref>. Thus, monotonicity is a key requirement
for establishing a faster rate. In the case of noisy responses and strictly monotone links,
the matter is more delicate. Namely, by the minimax rate of linear regression, which equals Var(Y-f(^⊤ X))/N, result (<ref>)
can only hold for all N ∈, if Var(Y-f(^⊤ X)) ∈(log^2(N)/N)
as N →∞.

The same observation carries over to
more general strictly monotonic links, because, in the presence of noise, condition (<ref>)
and the lower bound (<ref>) do not hold for all J ∈. For instance, if J^2 exceeds 1/Var(Y-f(^⊤ X))
or becomes even larger, the conditional distribution Y|Y ∈_J,ℓ has roughly as much correlation with label noise as it does with the labels, f(^⊤ X).
This suggests b_J,ℓ_2 ≈ 0, and that further slicing the distribution of X based on Y
will not decrease the variance in index vector direction, as required per (<ref>).
For J^-2≫Var(Y-f(^⊤ X)) however, Condition (<ref>) and
(<ref>) can still hold approximately, and our experiments below suggest that
the proposed estimator benefits from an aggresive splitting of the data
in the case of nonlinear, strictly monotonic links.
We now confer the implications of Theorem <ref>,
with synthetic examples, under a data-driven parameter choice rule for J.






  
Numerical setup and parameter choice.
We sample X ∼(0,_D), with D=10, and let = (1,0,…,0)^⊤ (the specific choice of  is irrelevant for the results due to the rotational invariance of X).
Responses are generated by Y = f(^⊤ X) + ζ, where ζ∼(0,σ_ζ^2Var(f(^⊤ X))).
We set α = 0.05 and  additionally exclude subsets with _J,ℓ < 2D, for the sake of numerical stability.
Our goal is twofold.
First, to compare estimation of the index vector using the standard OLS approach (<ref>) with the strategy proposed in this section, and second, to empirically confer the observations described after Theorem <ref>.

A critical step of the modified approach is the selection of J as a function of N.
As mentioned in Theorem <ref>, the denominator effectively
lower bounds λ_1(_J), and the corresponding concentration bound can be written as

    s _1(_J) - _2^2 ≲() + log(J) + u/α NJ∑_ℓ∈ I_Jρ̂_J,ℓκ_J,ℓ_J,ℓ^†X_J,ℓ_ψ_2^2/λ_1(_J).

Neglecting the dependence of κ_J,ℓ and _J,ℓ^†X_J,ℓ_ψ_2^2 on J
and log(J)-terms, originating from union bounds in the proofs, the error is minimized
when maximizing Jλ_1(_J).
We thus propose to adaptively choose J by

    J^* := max{J=⌈ (1.5)^k⌉ : k ∈_0,  Jλ_1(_J) > J' λ_1(_J')  ∀     J' < J},

where we use an exponential grid in  to decrease the computational demand. Note also that
λ_1(_J) ≥∑_ℓ∈ I_Jρ̂_J,ℓb̂_J,ℓ_2^2 ≈∑_ℓ∈ I_Jρ̂_J,ℓb_J,ℓ_2^2,
so, using parameter choice (<ref>), we expect no further increase of J if b_J,ℓ_2's decay rapidly.



  
Numerical experiments.

Figure <ref> shows the errors - _2 and the corresponding optimized choices J^∗ of the number of level-sets J, for different link functions f.
Solid lines correspond to the estimator (<ref>), the dashed lines to _1(_J^*), and different colors represent different noise levels σ_ζ.
In Figures <ref> - <ref> we consider monotonic link functions and we can see that
the approach presented in this section performs substantially better than the standard OLS approach.
On the other hand, in Figures <ref>, <ref>
and <ref> the two approaches achieve similar performance.
In case of <ref> this is because (<ref>) is indeed optimal,
according to the Gauss-Markov Theorem, whereas link functions in <ref>, <ref>,
are not monotonic.
In the latter case, the plots of J^* confirm that λ_1(_J) decays rapidly as a function of J, leaving J^* essentially constant
as a function of N.

Let us examine the results for monotonic functions in more detail.
The plots for J^∗ in Figures <ref> - <ref> show that the
number of level-sets J^* indeed grows as a function of N, and it does so up to a level dictated
by the noise level σ_ζ.
This shows that λ_1(_J) does not decay faster than J^-1 over a reasonably large range of N values, i.e. until  J^-1≈σ_ζ.
As a consequence, our approach achieves a faster, N^-1 estimation rate for the index vector.
Specifically, in the noise-free case this holds asymptotically as N→∞.
On the other hand, in the case of corrupted Y's, we first have a faster N^-1 convergence,
and then a sharp transition into the usual N^-1/2 rate.
The number of points N at which this transition occurs depends inversely on the level of noise.



§ ACKNOWLEDGEMENTS

We like to thank the anonymous referees for their feedback that helped to greatly improve earlier
versions of the manuscript.
This work is supported by the grant Function-driven Data Learning in High Dimensions
by the Research Council of Norway.






§ APPENDIX



 §.§ Additional technical results




Let X,Z be sub-Gaussian random vectors in ^D, let = X, and ∈^D× D.
Then we have


  *  X -  X_ψ_2≲X_ψ_2 (holds also for sub-exponential random variables),

  *  X̃_ψ_2≤_2 _(^⊤)X̃_ψ_2,

  *  ^⊤_2 ≲X̃_ψ_2^2,

  *  X^⊤ Z is sub-exponential with X̃^⊤Z̃ - [X̃^⊤Z̃]_ψ_1≲X̃_ψ_2Z̃_ψ_2.




Property <ref> is shown in <cit.> for sub-Gaussian random variables.
Applying the definition of the sub-Gaussian norm (<ref>) the claim follows,
since ^⊤X is a sub-Gaussian random variable for every ∈^D-1.
The same line of arguments holds for sub-exponential vectors.

For <ref> we compute for arbitrary ∈^D-1 and ^⊤∈(^⊤)

    ^⊤X̃_ψ_2   = ^⊤_2 (^⊤/^⊤_2)^⊤X̃_ψ_2
       ≤^⊤_2 sup_∈(^⊤) ∩^D-1^⊤X̃_ψ_2≤_2 _(^⊤)X̃_ψ_2.


For <ref> we first note that <cit.> implies
u≲ũ_ψ_2^2 for any sub-Gaussian u, where ũ=u-[u].
Thus, for every ∈^D-1 we have

    ^⊤ X = ^⊤ X≲^⊤X̃_ψ_2^2.

Taking the supremum with respect to ∈^D-1, the result follows since X is positive semidefinite.

Property <ref> follows from the centering property <ref> and <cit.>.




Let ∈^D× D be positive semidefinite, _1∈^d_1× D, _2∈^d_2× D.
For ∈^d_1, ∈^d_2 we have ^⊤_1_2^⊤≤√(^⊤_1 _1^⊤)√(^⊤_2 _2^⊤).
 Moreover, _1 _2_2^2≤_1_1^⊤_2_2 _2^⊤_2.


Applying the Cauchy-Schwartz inequality we have

    ^⊤_1 _2^⊤ = ⟨^1/2_1^⊤,^1/2_2^⊤⟩≤^1/2_1^⊤_2^1/2_2^⊤_2.

By the same line of argument we have
^1/2_1^⊤_2^2 = ^⊤_1 _1^⊤
and ^1/2_2^⊤_2^2 = ^⊤_2 _2^⊤, giving the first statement.
Considering now _2=_2=1, we have

    sup__2=1,
    _2=1^⊤_1_2^⊤≤√(sup__2=1^⊤_1 _1^⊤)√(sup__2=1^⊤_2 _2^⊤).

Notice that since  is positive semidefinite, then _1 _1^⊤ and _2 _2^⊤ are positive semidefinite. Therefore,

    sup__2=1^⊤_1 _1^⊤=sup__2=1(_1 _1^⊤)^1/2_2^2=_1 _1^⊤_2.

An analogous expression holds for the other term. Identifying the quadratic form on the left hand side in (<ref>) as the operator norm of _1_2, the conclusion follows.



The following is a standard concentration bound for sub-Gaussian and sub-exponential random vectors around their mean.


Let {X_i : i ∈ [N]} be independent copies of a centered random vector X ∈^D.
Denote := N^-1∑_i=1^NX_i and (t) = t ∨ t^2.
For any u>0 the following hold with probability at least 1-exp(-u).


  * If X_ψ_2 < ∞ we have
_2 ≲X_ψ_2√(D+u/N).

  * If X_ψ_1 < ∞ we have
_2 ≲X_ψ_1(√(D+u/N)).



The argument for the two bounds follows along analogous lines. Let δ < 1/4, and  be a δ-net
of ^D-1. We first use <cit.> to rewrite

    _2 = sup_∈^D-1^⊤= sup_∈^D-1 N^-1∑_i=1^N^⊤ X_i
     ≤ 2 sup_∈ N^-1∑_i=1^N^⊤ X_i.

The term N^-1∑_i=1^N^⊤ X_i is a sum of either sub-Gaussian or sub-exponential random variables, for any ∈.
In the former case we have  ^⊤ X_i_ψ_1≤X_ψ_1, and ^⊤ X_i_ψ_2≤X_ψ_2 in the latter.
Hoeffding's inequality <cit.>, in the sub-Gaussian case, or Bernstein's inequality
<cit.>, in the sub-exponential case, now yield
concentration bounds for the sums.
The claim follows by applying the union bound over ∈, where the number of events is bounded by ≤ 12^D,
see for instance <cit.>.




 §.§ Proofs for Section 2



Let X = X -  X.
Let _∈^d_A× D and _∈^d_B× D be matrices whose rows contain the orthonormal basis for () and (), respectively.
Since  and  are symmetric, and ()⊆(), we have -=(-), where  is the orthogonal projection onto ().
We thus have

    (-)^⊤_2 = (-)^⊤_2,

for =_∈^d_A× D and =_∈^d_B× D.
Denote now = N^-1∑_i=1^NX_iX_i^⊤. Notice that , compared to the empirical covariance , uses the true, instead of the empirical mean of X.
We then have

    (- )^⊤_2 ≤(- )^⊤_2 + ∑_i=1^NX_i/N_2∑_i=1^NX_i/N_2.

By Lemma <ref> the second term
is always of higher order, and the resulting error can be absorbed into a corresponding upper
bound for the first term. Thus, in the following we focus on (- )^⊤_2.

We closely follow the proof of <cit.>.
Let δ < 1/4, and by ,  denote δ-nets of
 ^d_A-1 and ^d_B-1. From <cit.> we have

    (- ) ^⊤_2    ≤ (1-2δ)^-1sup_∈
    ∈⟨(- ) ^⊤, ⟩
       ≤
    2sup_∈
    ∈⟨ (- )  ^⊤, ^⊤⟩.

Consider now any pair (,)∈×, and write

    ⟨^⊤,^⊤⟩ = ∑_i=1^N ⟨X_i (X_i^⊤^⊤),^⊤⟩/N =
    ∑_i=1^N⟨X_i,⟩⟨X_i,⟩/N.

Since ⟨X_i,⟩ and ⟨X_i,⟩ are sub-Gaussian, their product is sub-exponential, and from
<cit.> we have

    ⟨X_i, ⟩⟨X_i,⟩_ψ_1   ≤⟨X_i, ⟩_ψ_2⟨X_i,⟩_ψ_2≤X̃_ψ_2X̃_ψ_2,

which we denote by σ_:= X̃_ψ_2X̃_ψ_2 for short.
Since []=, by Lemma <ref>
we have for any u > 0, with (t) = t∨ t^2,

    (|⟨^⊤, ^⊤⟩ - ⟨^⊤, ^⊤⟩|≥σ_(√(1+u/N))) ≤exp(-u).

The size of the nets can be bounded as ≤ 12^d_A and ≤ 12^d_B, see <cit.>.
Thus, considering all pair (,)∈× and using the union bound we get, for some universal constant C (bounded e.g. by log(12) + 1),

    (|sup_∈, ∈⟨^⊤,^⊤⟩-⟨^⊤,^⊤⟩|≥ C σ_(√(u+(d_A + d_B)/N))) 
                           ≤exp-u - (d_A + d_B)log(12)
                           ≤exp((d_A + d_B)log(12)-u-(d_A + d_B)log(12)).


Using the δ-net approximation bound (<ref>) we thus get

    ((- ) ^⊤_2≥ 2C σ_(√(u+(d_A + d_B)/N))) 
       ≤(|sup_∈
    ∈⟨^⊤,^⊤⟩-⟨^⊤,^⊤⟩|≥ C σ_(√(u+(d_A + d_B)/N))).


The result now follows from X_ψ_2=X_ψ_2 and X_ψ_2=X_ψ_2
and thus σ_ = X_ψ_2X_ψ_2.



We use the shorthand notation λ_i := λ_i() and λ̂_i := λ̂_i().
By the definition of δ_il we have

    δ_il   = λ̂_i-1 - λ_i∧λ_l - λ̂_l+1≥ (λ̂_i-1 - λ_i) ∧ (λ_l - λ̂_l+1)
       ≥((λ_i-1-λ_i) + (λ̂_i-1 - λ_i-1)) ∧((λ_l-λ_l+1) + (λ_l+1 - λ̂_l+1)),

so it suffices to show λ̂_i-1 - λ_i-1≥ -(λ_i-1-λ_i)/2
and λ_l+1 - λ̂_l+1≥ - (λ_l-λ_l+1)/2.
Note that if i=1 or l = (), the corresponding gap is trivial and
no concentration is required. Hence, we can focus on the case i > 1 and l < ().
Using <cit.>, we obtain that for any u>0 a relative eigenvalue bound

    λ̂_i-1 - λ_i-1 < -λ_i-1-λ_i/2,

holds with probability at most exp(-u), provided

    N > C_K(λ_i-1/λ_i-1-λ_i∨ 1)(∑_j=1^i-1λ_j/λ_j - λ_i-1 + λ_i-1-λ_i/2∨ uλ_i-1/λ_i-1-λ_i).

Furthermore, since λ_i-1-(λ_i-1-λ_i)/2 < λ_j for all j ≥ i-1,
the summands in the sum in (<ref>) are increasing, and we thus get

    ∑_j=1^i-1λ_j/λ_j - λ_i-1 + λ_i-1-λ_i/2≤ (i-1)2λ_i-1/λ_i-1-λ_i≲() λ_i-1/λ_i-1-λ_i.

Thus, (<ref>) is implied by the requirement on N in the statement. Similarly, using <cit.> we get that

    λ̂_l+1 - λ_l+1 > λ_l-λ_l+1/2

holds with probability at most exp(-u), provided

    N > C_K(λ_l+1/λ_l-λ_l+1∨ 1)(∑_j=l+1^Dλ_j/λ_l+1 - λ_j + λ_l-λ_l+1/2∨ uλ_l+1/λ_l-λ_l+1).

As before, we have λ_j - (λ_l-λ_l+1)/2 < λ_l+1 for all j ≥ l+1,
so that

    ∑_j=l+1^Dλ_j/λ_l+1 - λ_j + λ_l-λ_l+1/2≤ (()-l)2λ_ℓ+1/λ_l-λ_l+1≲() λ_l+1/λ_l-λ_l+1.

It follows that (<ref>) is implied by the requirement on N in the statement.
Taking the complementary event of the union of events (<ref>)
and (<ref>), which holds with probability 1-2exp(-u),
we get λ̂_i-1 - λ_i-1≥ -(λ_i-1-λ_i)/2
and λ_l+1 - λ̂_l+1≥ - (λ_l-λ_l+1)/2
and thus the result is proven.




 §.§ Proofs for Section <ref>




If
√(^†) (-) √(^†)_2 < 1,
then () ∩() = {0}.


We prove the statement by contraction.
Assume there exists a nonzero ∈() ∩(). Since ∈() we have ∈(√(^†)) implying the existence of a ∈^D-1 with = √(^†) (possibly after scaling ).  Using the properties of  and , we get a contradiction by

    1 = ^⊤ = ^⊤  (-)   = ^⊤√(^†)  (-)  √(^†)≤√(^†)  (-)  √(^†)_2<1.




Assume for now √(^†) (-) √(^†)_2 < 1/2 which will be
ensured later by concentration arguments. Conditioned on this event, Lemma <ref> implies () ∩() = {0}. Since () ⊂() holds almost surely, this implies () = () almost surely.
Denote now := -.
We want to use a von Neumann series argument to express ^† in terms
of ^† and the corresponding finite sample perturbation. Recall that for general symmetric matrices  and  with  () = (), we have  ()^† = ^†^†  <cit.>. Letting
 be the orthogonal projection onto (), and =_D-, we can thus write

    ^†   = ( + Δ)^† = (^1/2( + √(^†)Δ√(^†))^1/2)^†
       =√(^†)( + √(^†)Δ√(^†))^†√(^†)
       =√(^†)(_D -  +  √(^†)Δ√(^†))^†√(^†)
       =√(^†)(_D +  √(^†)Δ√(^†))^-1√(^†),

where we used
the additivity of the Moore-Penrose inverse for matrices with orthogonal non-trivial eigenspaces in the last equality. Using a von Neumann series expansion for the second
term, we obtain

    ^†   = √(^†)(_D +  √(^†)Δ√(^†))^†√(^†) =√(^†)∑_k=0^∞(-√(^†)Δ√(^†))^k√(^†)
       = ^† + √(^†)∑_k=1^∞(-√(^†)Δ√(^†))^k√(^†),

where we used √(^†)√(^†) = ^† in the last equality.
Subtracting ^†, and multiplying by  from the right and ^⊤ from the left, it follows that

    (^† - ^†)^⊤ = √(^†)∑_k=1^∞(-√(^†)Δ√(^†))^k√(^†)^⊤
       =  -^†Δ^†^⊤ + ^†Δ√(^†)∑_k=0^∞(-√(^†)Δ√(^†))^k√(^†)Δ^†^⊤.

The rest of the proof requires to apply Lemma <ref> multiple times to concentrate Δ when multiplied by different matrices from left and right. More specifically, since √(^†)X_ψ_2^2 ≳√(^†)X_2≳ 1,
we can assume the regime N > () + u by the assumption on N in Theorem <ref>.
Then, we have by Lemma <ref> with probability 1-4exp(-u)

    ^†^†^⊤_2    ≲^†X_ψ_2^†X_ψ_2√(()+u/N), 
    ^†√(^†)_2    ≲^†X_ψ_2√(^†)X_ψ_2√(()+u/N),
    √(^†)^†^⊤_2    ≲^†X_ψ_2/√(^†)X_ψ_2,  and √(^†)√(^†)≤1/2,

whenever N > C(() + u)√(^†)X_ψ_2^4 (we used properties of the sub-Gaussian norm in Lemma <ref>
to get √(^†)X_ψ_2^2 ≥ CCov(√(^†)X)_2  = C)
and thus √(^†)X_ψ_2^2≲√(^†)X_ψ_2^4).
Using identity (<ref>), combining the above bounds and adjusting the uniform constant C, the desired result holds with probability 1-exp(-u).



For the covariance estimation bound denote first = N^-1∑_i=1^NX_i X_i^⊤, where X_i = X_i -  X,
and decompose the error, as in the proof of Lemma <ref> into

    (- )^⊤_F   = (- )^⊤_F 
       ≤(- )^⊤_F + 1/N∑_i=1^NX_i_F1/N∑_i=1^NX_i_F.

The second term is again of higher order, with high probability, and can thus be disregarded. For the first term, define random matrices _i=X_iX_i^⊤^⊤-^⊤.
First note [_i]=0, and

    _i_F   ≤X_iX_i^⊤^⊤_F +^⊤^⊤_F
       =X_i_2X_i_2+XX^⊤^⊤_F≤2C_AC_B.

Moreover, _i take values in the Hilbert space (^d_A× d_B, ·_F).
Thus, by <cit.>, for every u>0 we have

    1/N∑_i=1^N _i_F≤2√(2)C_AC_Bu/√(N),

with probability at least 1-exp(-u^2).
The conclusion follows after adjusting the confidence level u to account for the higher order term and the union bound.

The first steps for establishing the bound for precision matrices are the same as in the proof of Theorem <ref>.
Indeed, the only differences are in the following lines of inequalities.
Instead of (<ref>), we have

    ^†^†^⊤_F    ≲ C_A^† C_B^†√(1+u/N),  and √(^†)^†^⊤_F ≲ C_B^†√(Θ)√(1+u/N);
    ^†√(^†)_F    ≲C_A^†/√(Θ),  and √(^†)√(^†)_F≤1/2,

 using the just proven covariance bounds.
The claim now follows by using this together with (<ref>), as in the proof of Theorem <ref>.




 §.§ Proofs for Section <ref>


We first need a bound for = X,Y, and a concentration around the finite
sample estimate = N^-1∑_i=1^N(X_i-_X)(Y_i - _Y).


Let ∈^k× D. If X ∈^D,Y ∈ are sub-Gaussian, we have
_2 ≲X_ψ_2Y_ψ_2.
Furthermore, with probability at least 1 - exp(-u) we have

    ( - )_2 ≲X_ψ_2Y_ψ_2√(k+u/N)∨k+u/N.



By Lemma <ref> we have
^⊤  X≲X_ψ_2^2, and Y≲Y_ψ_2^2, which implies

    _2    =  X,Y_2 = sup_∈^k-1^⊤ X,Y =
    sup_∈^k-1^⊤  X,Y
       ≤sup_∈^k-1√(^⊤  X)√(Y)≲X_ψ_2Y_ψ_2.

Define the random variable Z_i := X_iY_i - X,Y with Z_i = 0.
We write

    ( - ) = 1/N∑_i=1^N  Z_i - (1/N∑_i=1^NX_i)(1/N∑_i=1^NY_i).

As in the proof of Lemma <ref>, by applying Lemma <ref> it follows that the second term is of higher order.
On the other hand, the first term is an empirical mean of a sub-exponential centered variable Z_i,

    Z_i_ψ_1   = X_iY_i -  X,Y_ψ_1≲X_iY_i_ψ_1
       ≲sup_∈^k-1^⊤X_iY_i_ψ_1≲X_i_ψ_2Y_i_ψ_2,

where we use the centering property of the sub-exponential norm, and the bound for the sub-exponential norm by the
product of sub-Gaussian norms (Lemma <ref>).
Applying now Lemma <ref>, the statement follows.


We first note that √(^†)X_ψ_2^2 ≳√(^†)X_2 ≥ 1,
so by choosing C in the statement large enough, we assume the regime N > (() + u) in the following.
We begin with a bound for (^† - ^†)_2.

    (^† - ^†)_2 ≤(^†-^†)_2 +(^†-^†)_2
                           +^†(- )_2
    + ^†(- )_2
         ≤(^†-^†)_2_2
    + (^†-^†)_2_2 
                           + ^†_2(- )_2 +^†_2(- )_2.

By Lemma <ref> we have _2 ≤X_ψ_2Ỹ_ψ_2, _2 ≤X_ψ_2Ỹ_ψ_2.
Furthermore, since ∈(), we can rewrite ( - )_2
= ( - )_2 and ( - )_2
= ( - )_2, where the rows of ∈^d_Q× D and
 ∈^d_P× D contain orthonormal bases for () and (), respectively.
 Using Lemma <ref> and d_Q ∨ d_P ≤(),
we have with probability at least 1 - 2exp(-u)

    ( - )_2 ≲X_ψ_2Ỹ_ψ_2√(() + u/N), 
    ( - )_2 ≲X_ψ_2Ỹ_ψ_2√(() + u/N).

Furthermore, by Theorem <ref> whenever
N > C (()+u)√(^†)X_ψ_2^4 we get with probability 1-2exp(-u)
for each ∈{,}

    (^†-^†)_2 ≲^†X_ψ_2^†X_ψ_2√(() + u/N).

Conditioned on all four events we now have

    (^† - ^†)_2_2
       ≲Ỹ_ψ_2^†X_ψ_2√(κ(,X))√(() + u/N),
    (^† - ^†)_2_2
       ≲Ỹ_ψ_2^†X_ψ_2√(κ(,X))√(() + u/N).

Moreover, since N > () + u as stated in the beginning, we get

    ^†_2(- )_2    ≤(^†_2+ (^†-^†)_2)(- )_2
       ≲Y_ψ_2^†X_ψ_2√(κ(, X))√(() + u/N),

and

    ^†_2(- )_2    ≤(^†_2+ (^†-^†)_2)(- )_2
       ≲Y_ψ_2^†X_ψ_2√(κ(, X))√(() + u/N),

where we use property <ref> in Lemma <ref>, exploiting the identity ^†=^†^†, and using Lemma <ref> in the second line.
Combining the bounds and conditioning on the events above, we have that with probability at least 1-4exp(-u)

    (^† - ^†)_2 ≲Y_ψ_2^†X_ψ_2√(κ(,X) ∨κ(,X))√(() + u/N),

if N> C (() + u)(√(^†)X_ψ_2^4).
Adjusting the constant C to account for the change in the probability constant and the requirement on N, the claim follows.
The proof for the bound on ^†_2 follows analogous lines or argument.



Assume for the moment ^⊤ > 0. In this case  and  are co-linear since
=⟨, ⟩. Therefore, =/_2 = /_2,
which implies

    (- )_2 =
    _2 |1/_2 - 1/_2|≤|_2 - _2/_2|≤_2/_2.

Since =0 we now have

    - _2 = √((- )_2^2 + _2^2)≤√(2)_2/_2≤√(2)_2/_2 - (- )_2.

Therefore, it suffices to ensure (- )_2 < 1/2_2, which would also give ^⊤ > 0.
To that end, Lemma <ref> gives that (- )_2 ≤1/2_2 holds with probability at least 1-exp(-u)
provided

    N > C(()+u)(√(^†)X_ψ_2^4 ∨Y_ψ_2^2^†X_ψ_2^2κ_/^2_2).

The result follows by bounding _2 through Lemma <ref>.



We first note 
    Y_ψ_2≤f(^⊤ X) -  f(^⊤ X)_ψ_2 + σ_ζ.

Furthermore, we can use a bounded difference inequality for unbounded spaces,
see <cit.>, to get
f(^⊤ X) -  f(^⊤ X)_ψ_2≲ L^⊤ X_ψ_2≲√(K) Lσ_,
where we used the strict sub-Gaussianity in (<ref>) and ^⊤ = σ_^2.
Therefore,

    Y_ψ_2   ≲ KLσ_ + σ_ζ.

Using κ(,X) ≲ K^2, √(^†)X_ψ_2 = √(K), and plugging in the bounds on Y_ψ_2, as well as
^†X_ψ_2 = √(K)/σ_ and ^†X_ψ_2 = √(K^†)_2,
the claim follows.



If Z is sub-Gaussian and E an event with (E)>0, then Z| E is also sub-Gaussian


Assume without loss of generality that Z∈. The argument for vectors then follows by the definition.
We use the characterization of sub-Gaussianity by the moment bound <cit.>.
Let p≥ 1. By the law of total expectation it follows

    [| Z|^p] = [| Z|^p| E](E) +[| Z|^p| E^𝖼](E^𝖼)≥[| Z|^p| E](E).

 Dividing by (E) and using the monotonicity of the p-th root yields
 
    ([| Z|^p| E])^1/p≤([| Z|^p])^1/p/(E)^1/p≲Z_ψ_2√(p)/(E),

where we used (E)≤ 1 and the sub-Gaussianity of Z in the last inequality.



Using 0 ≤max_s = ± 1⟨ s _1(_J), ⟩≤ 1, we first compute

    min_s =  ± 1 s_1(_J) - _2^2    = min_s = ± 1 2(1-⟨ s _1(_J), ⟩)
    
       ≤ 2(1-⟨_1(_J), ⟩^2) ≤ 2_1(_J)_2^2.

The Davis-Kahan Theorem <cit.> for _J and = ^⊤ then gives

    _1(_J)_2^2 =_1(_J)_1(_J)^⊤_F^2 ≤( - _J)_F^2/λ_1(_J)^2 =
    _J_F^2/λ_1(_J)^2.

It remains to find a concentration bound for the matrix _J around zero and
a lower bound for the leading eigenvalue λ_1(_J).
Now let π : [I_J] → I_J be bijective and introduce the matrix

    _J = [√(ρ̂_J,π(1))_J,π(1)  |…|√(ρ̂_J,π(I_J))_J,π(I_J)] ∈^D×I_J,

satisfying _J = _J  _J^⊤, and thus λ_1(_J) = λ_1(_J  _J^⊤).
Using _F ≤_2 _F and _F=^⊤_F, which hold for arbitrary  and , yields

    _J_F^2 = _J  _J^⊤_F^2 ≤_J_2^2 _J_F^2 ≤λ_1(_J)_J_F^2.

By Lemma <ref>, conditioning on a sub-Gaussian random vector gives a sub-Gaussian random vector. Thus, by (<ref>) we have
simultaneously for all ℓ∈ I_J

    _J,ℓ_2^2≲Y_J,ℓ_ψ_2^2 _J,ℓ^†X_J,ℓ_ψ_2^2 κ_J,ℓ(_J,ℓ) +  log(I_J) + u/N_J,ℓ,

with probability at least 1 - exp(-u), provided N_J,ℓ > C((_J,ℓ) + log(I_J)+u)K_J,ℓ^4
for all ℓ∈ I_J.
By definition of I_J, we have N_J,ℓ > α NJ^-1, and thus the previous condition
is satisfied whenever for all ℓ∈ I_J

    N/J   > C((_J,ℓ) + log(J) + u)K_J,ℓ^4/α,

which is implied by N > CJK_J^4(()+log(J)+u)/α.
Under the same conditions and with similar probability we obtain

    _J_F^2    ≤λ_1(_J) ∑_ℓ∈ I_Jρ̂_J,ℓ_J,ℓ_2^2
    
       ≲λ_1(_J) ∑_ℓ∈ I_Jρ̂_J,ℓY_J,ℓ_ψ_2^2_J,ℓ^†X_J,ℓ_ψ_2^2 κ_J,ℓ() + log(J) + u/N_J,ℓ
       ≲λ_1(_J)(() +log(J) + u)/α N J∑_ℓ∈ I_Jρ̂_J,ℓ_J,ℓ^†X_J,ℓ_ψ_2^2 κ_J,ℓ,

where we used N_J,ℓ > α N J^-1 and  Y_J,ℓ_ψ_2≲_J,ℓ = J^-1
(definition of _J,ℓ) in the last inequality.
To bound λ_1(_J) from below, we note that on the same event where
bound (<ref>) holds, we also have for all ℓ∈ I_J

    (_J,ℓ -_J,ℓ)_2^2    ≲Y_J,ℓ_ψ_2^2 _J,ℓ^†X_J,ℓ_ψ_2^2 κ_J,ℓ(_J,ℓ) + log(J) + u)/N_J,ℓ
       ≲() + log(J) + u/α N J_J,ℓ^†X_J,ℓ_ψ_2^2 κ_J,ℓ

Thus, multiplying _J from left and right by  and using <ref>, we obtain

    λ_1(_J)≥^⊤_J=∑_ℓ∈ I_Jρ̂_J,ℓ⟨, b̂_J,ℓ⟩^2
    = ∑_ℓ∈ I_Jρ̂_J,ℓ(b_J,ℓ_2^2 - (b_J,ℓ-b̂_J,ℓ)_2^2)
           ≥∑_ℓ∈ I_Jρ̂_J,ℓ(b_J,ℓ_2^2 - C(() + log(J) + u)/α N J_J,ℓ^†X_J,ℓ_ψ_2^2 κ_J,ℓ)

for some universal constant C > 0. The result follows by plugging the bounds on _J_F^2
and λ_1(_J) into (<ref>), and using the definition of ε_N,J,u.





<><#>1#1   

adamczak2010quantitative
Adamczak, R., Litvak, A., Pajor, A., and
  Tomczak-Jaegermann, N. (2010).
Quantitative estimates of the convergence of the empirical covariance
  matrix in log-concave ensembles.
Journal of the American Mathematical Society 23, 2,
  535–561.

adragni2009sufficient
Adragni, K. P. and Cook, R. D. (2009).
Sufficient dimension reduction and prediction in regression.
Philosophical Transactions of the Royal Society A:
  Mathematical, Physical and Engineering Sciences 367, 1906,
  4385–4405.

AndersonStatAn
Anderson, T. W. (2003).
An Introduction to Multivariate Statistical Analysis, 3rd
  ed.
John Wiley & Sons.

arbel2019strictsubgaussian
Arbel, J., Marchal, O., and Nguyen, H. D.
  (2019).
On strict sub-gaussianity, optimal proxy variance and symmetry for
  bounded random variables.
to appear in ESAIM: Probability and statistics
  arxiv:1901.09188.

balabdaoui2019score
Balabdaoui, F., Groeneboom, P., and
  Hendrickx, K. (2019).
Score estimation in the monotone single-index model.
Scandinavian Journal of Statistics 46, 2, 517–544.

rosasco2010learningwithintegraloperator
Belkin, M., Rosasco, L., and Vito, E. D.
  (2010).
On learning with integral operators.
Journal of Machine Learning Research 11, 905–934.

bhatia2013matrix
Bhatia, R. (2013).
Matrix analysis.  Vol. 169.
Springer Science & Business Media.

BickelLevinaThresholding
Bickel, P. and Levina, E. (2008a).
Covariance Regularization by thresholding.
The Annals of Statistics 36, 6, 2577–2604.

BickelLevinaBanding
Bickel, P. and Levina, E. (2008b).
Regularized estimation of large covariance matrices.
The Annals of Statistics 36, 1, 199–227.

brillinger1982generalized
Brillinger, D. (1982).
A generalized linear model with "gaussian" regressor variables.
A Festschrift For Erich L. Lehmann, 97.

TCaiCLIME
Cai, T., Liu, W., and Luo, X. (2011).
A constrained ℓ_1 minimization approach to sparse precision
  matrix estimation.
Journal of the American Statistical Association 106,
  594–607.

ACLIME
Cai, T., Liu, W., and Zhou, H. (2016).
Estimating sparse precision matrix: optimal rates of convergence and
  adaptive estimation.
Annals of Statistics 44, 455–488.

cai2016estimating
Cai, T., Ren, Z., and Zhou, H. H. (2016).
Estimating structured high-dimensional covariance and precision
  matrices: Optimal rates and adaptive estimation.
Electronic Journal of Statistics 10, 1, 1–59.

dalalyan2008new
Dalalyan, A. S., Juditsky, A., and Spokoiny,
  V. (2008).
A new algorithm for estimating the effective dimension-reduction
  subspace.
Journal of Machine Learning Research 9, Aug,
  1647–1678.

dennis2000save
Dennis Cook, R. (2000).
SAVE: A method for dimension reduction and graphics in regression.
Communications in Statistics - Theory and
  Methods 29, 9-10, 2109–2121.

dobson2008introduction
Dobson, A. and Barnett, A. (2008).
An introduction to generalized linear models.
Chapman and Hall/CRC.

fan2016overview
Fan, J., Liao, Y., and Liu, H. (2016).
An overview of the estimation of large covariance and precision
  matrices.
The Econometrics Journal 19, 1, C1–C32.

FHTGLasso
Friedman, J., Hastie, T., and Tibshirani,
  R. (2008).
Sparse inverse covariance estimation with the graphical lasso.
Biostatistics 9, 3, 432–441.

gyorfi2006distribution
Györfi, L., Kohler, M., Krzyzak, A.,
  and Walk, H. (2006).
A distribution-free theory of nonparametric regression.
Springer Science & Business Media.

hristache2001direct
Hristache, M., Juditsky, A., and Spokoiny,
  V. (2001).
Direct estimation of the index coefficient in a single-index model.
The Annals of Statistics 29, 3, 595–623.

PenalisedLikelihoodCovs
Huang, J., Liu, N., Pourahmadi, M., and
  Liu, L. (2006).
Covariance matrix selection and estimation via penalised normal
  likelihood.
Biometrika 93, 1, 85–98.

ichimura1993semiparametric
Ichimura, H. (1993).
Semiparametric least squares (SLS) and weighted SLS estimation of
  single-index models.
Journal of Econometrics 58, 1-2, 71–120.

kakade2011efficient
Kakade, S. M., Kanade, V., Shamir, O., and
  Kalai, A. (2011).
Efficient learning of generalized linear and single index models with
  isotonic regression.
In Advances in Neural Information Processing Systems.
  927–935.

kalai2009isotron
Kalai, A. T. and Sastry, R. (2009).
The isotron algorithm: High-dimensional isotonic regression.
In COLT. Citeseer.

klasing2009comparison
Klasing, K., Althoff, D., Wollherr, D., and
  Buss, M. (2009).
Comparison of surface normal estimation methods for range sensing
  applications.
In 2009 IEEE International Conference on Robotics and
  Automation. IEEE, 3206–3211.

klock2020estimating
Klock, T., Lanteri, A., and Vigogna, S.
  (2020).
Estimating multi-index models with response-conditional least
  squares.
arXiv preprint arXiv:2003.04788.

koltchinskii2017asymptotically
Koltchinskii, V. (2017).
Asymptotically efficient estimation of smooth functionals of
  covariance operators.
arXiv preprint arXiv:1710.09072.

koltchinskii2020efficient
Koltchinskii, V., Löffler, M., Nickl, R.,
  and others. (2020).
Efficient estimation of linear functionals of principal components.
The Annals of Statistics 48, 1, 464–490.

koltchinskii2017concentration
Koltchinskii, V., Lounici, K., and others.
  (2017a).
Concentration inequalities and moment bounds for sample covariance
  operators.
Bernoulli 23, 1, 110–133.

koltchinskii2017normal
Koltchinskii, V., Lounici, K., and others.
  (2017b).
Normal approximation and concentration of spectral projectors of
  sample covariance.
The Annals of Statistics 45, 1, 121–157.

kontorovich2014concentration
Kontorovich, A. (2014).
Concentration in unbounded metric spaces and algorithmic stability.
Proceedings of Machine Learning Research, Vol. 32.
  Proceedings of Machine Learning Research, 28–36.

kuchibhotla2020efficient
Kuchibhotla, A. K., Patra, R. K., and
  others. (2020).
Efficient estimation in single index models through smoothing
  splines.
Bernoulli 26, 2, 1587–1618.

li2007directional
Li, B. and Wang, S. (2007).
On directional regression for dimension reduction.
Journal of the American Statistical
  Association 102, 479, 997–1008.

li1991sliced
Li, K.-C. (1991).
Sliced inverse regression for dimension reduction.
Journal of the American Statistical
  Association 86, 414, 316–327.

li1989regression
Li, K.-C. and Duan, N. (1989).
Regression analysis under link violation.
The Annals of Statistics 17, 3, 1009–1052.

li2018some
Li, W., Chen, Y., Vong, S., and
  Luo, Q. (2018).
Some refined bounds for the perturbation of the orthogonal projection
  and the generalized inverse.
Numerical Algorithms 79, 2, 1–21.

little15multiscaleSVD
Little, A. V., Maggioni, M., and Rosasco,
  L. (2017).
Multiscale geometric methods for data sets I: Multiscale SVD, noise
  and curvature.
Applied and Computational Harmonic Analysis 43, 3,
  504–567.

FASP
Liu, W. and Luo, X. (2015).
Fast and adaptive sparse precision matrix estimation in high
  dimensions.
Journal of Multivariate Analysis 135, 153–162.

ma2013review
Ma, Y. and Zhu, L. (2013).
A review on dimension reduction.
International Statistical Review 81, 1, 134–150.

mcaleer2008single
McAleer, M. and Da Veiga, B. (2008).
Single-index and portfolio models for forecasting value-at-risk
  thresholds.
Journal of Forecasting 27, 3, 217–235.

meinshausenGLASSO
Meinshausen, N. and Bühlmann, P. (2006).
High-dimensional graphs and variable selection with the lasso.
The Annals of Statistics 34, 3, 1436–1462.

meng2010optimal
Meng, L. and Zheng, B. (2010).
The optimal perturbation bounds of the Moore-Penrose inverse under
  the Frobenius norm.
Linear Algebra and its Applications 432, 4,
  956–963.

mitra2003estimating
Mitra, N. J. and Nguyen, A. (2003).
Estimating surface normals in noisy point cloud data.
In Proceedings of the nineteenth annual symposium on
  Computational geometry. ACM, 322–328.

ouyang2005normal
OuYang, D. and Feng, H.-Y. (2005).
On the normal vector estimation for point cloud data from smooth
  surfaces.
Computer-Aided Design 37, 10, 1071–1079.

plan2016generalized
Plan, Y. and Vershynin, R. (2016).
The generalized lasso with non-linear observations.
IEEE Transactions on information theory 62, 3,
  1528–1537.

plan2016high
Plan, Y., Vershynin, R., and Yudovina, E.
  (2016).
High-dimensional estimation with geometric constraints.
Information and Inference: A Journal of the
  IMA 6, 1, 1–40.

radchenko2015high
Radchenko, P. (2015).
High dimensional single index models.
Journal of Multivariate Analysis 139, 266–282.

reiss2020nonasymptotic
Reiß, M. and Wahl, M. (2020).
Nonasymptotic upper bounds for the reconstruction error of PCA.
The Annals of Statistics 48, 2, 1098–1123.

srivastava2013covariance
Srivastava, N. and Vershynin, R. (2013).
Covariance estimation for distributions with 2+ε moments.
The Annals of Probability 41, 5, 3081–3111.

tian2004some
Tian, Y. and Cheng, S. (2004).
Some identities for Moore–Penrose inverses of matrix products.
Linear and Multilinear Algebra 52, 6, 405–420.

WesselRIDGE
van Wieringen, W. N. (2019).
The Generalized Ridge Estimator of the Inverse Covariance Matrix.
Journal of Computational and Graphical
  Statistics 28, 4, 932–942.

vershynin2012close
Vershynin, R. (2012a).
How close is the sample covariance matrix to the actual covariance
  matrix?
Journal of Theoretical Probability 25, 3, 655–686.

vershynin2010introduction
Vershynin, R. (2012b).
Introduction to the non-asymptotic analysis of random matrices.
210–268.

vershynin2018high
Vershynin, R. (2018).
High-dimensional probability: An introduction with
  applications in data science.  Vol. 47.
Cambridge University Press.

wedin1973perturbation
Wedin, P.-Å. (1973).
Perturbation theory for pseudo-inverses.
BIT Numerical Mathematics 13, 2, 217–232.

weyl1912asymptotische
Weyl, H. (1912).
Das asymptotische Verteilungsgesetz der Eigenwerte linearer
  partieller Differentialgleichungen (mit einer Anwendung auf die Theorie
  der Hohlraumstrahlung).
Mathematische Annalen 71, 4, 441–479.

xu2018perturbation
Xu, X. (2018).
On the perturbation of the Moore-Penrose inverse of a matrix.
arXiv preprint arXiv:1809.00203.

yaskov2014lower
Yaskov, P. (2014).
Lower bounds on the smallest eigenvalue of a sample covariance
  matrix.
Electronic Communications in Probability 19.

yu2014useful
Yu, Y., Wang, T., and Samworth, R. (2014).
A useful variant of the Davis-Kahan theorem for statisticians.
Biometrika 102, 2, 315–323.




