# Inference for Projection Parameters in Linear Regression:
beyond $d=o(n^{1/2})$

Woonyoung Chang  Department of Statistics & Data Science, Carnegie Mellon
University  Arun Kumar Kuchibhotla  Department of Statistics & Data Science,
Carnegie Mellon University  Alessandro Rinaldo  Department of Statistics &
Data Science, The University of Texas at Austin

###### Abstract

We consider the problem of inference for projection parameters in linear
regression with increasing dimensions. This problem has been studied under a
variety of assumptions in the literature. The classical asymptotic normality
result for the least squares estimator of the projection parameter only holds
when the dimension $d$ of the covariates is of a smaller order than $n^{1/2}$,
where $n$ is the sample size. Traditional sandwich estimator-based Wald
intervals are asymptotically valid in this regime. In this work, we propose a
bias correction for the least squares estimator and prove the asymptotic
normality of the resulting debiased estimator. Precisely, we provide an
explicit finite sample Berry Esseen bound on the Normal approximation to the
law of the linear contrasts of the proposed estimator normalized by the
sandwich standard error estimate. Our bound, under only finite moment
conditions on covariates and errors, tends to 0 as long as $d=o(n^{2/3})$ up
to the polylogarithmic factors. Furthermore, we leverage recent methods of
statistical inference that do not require an estimator of the variance to
perform asymptotically valid statistical inference and that leads to a sharper
miscoverage control compared to Wald’s. We provide a discussion of how our
techniques can be generalized to increase the allowable range of $d$ even
further.

##  1 Introduction

Linear regression is a fundamental statistical tool that has been widely used
in various fields of research. The classical literature on linear regression
studies the ordinary least square (OLS) estimation has focused primarily on
the well-specified case, where the underlying truth postulates the linear
relation between the response variable and covariates. As elucidated in works
in an assumption-lean framework (Berk et al., 2019; Vansteelandt and Dukes,
2022), although the model assumption sometimes takes account of the pre-
knowledge, the usage of model assumptions is dishonest when used because of
the mathematical convenience. As real data often possess a highly nonlinear
structure, relying on model assumptions and considering them as representing
ground truth in inference may be problematic. Compared to the popularity of
parametric regression in both practical studies in statistics and
econometrics, the theoretically established properties of regression
coefficients in an assumption-lean framework have just started to gain
attraction. To this effort, some approaches start with a traditional estimator
of a parameter indexing a parametric regression model and then characterize
what estimand the estimator converges to, without assuming that the model is
true (Kuchibhotla et al., 2020; Berk et al., 2019). In particular, we focus on
misspecified linear regression models comprised of a $d$-dimensional random
vector of covariates $X$ and a scalar random variable $Y$. If
$\mathbb{P}_{X,Y}$ admits the second moment, then the conditional expectation
$\mathbb{E}[Y|X]$, which is not necessarily linear, is the best $L^{2}$
approximation to $Y$ among functions of $X$. It is well known that the best
linear $L_{2}$ approximation to $Y$ is the linear function $\beta^{\top}X$ is
well-defined where the coefficients $\beta=\beta(\mathbb{P}_{X,Y})$ is given
by

$\displaystyle\beta=\operatorname*{arg\,min}_{\theta\in\mathbb{R}^{d}}\mathbb{E}[(Y-\theta^{\top}X)^{2}].$

Provided that the population Gram matrix $\Sigma=\mathbb{E}_{X}[XX^{\top}]$ is
invertible, the solution is unique and given by the vector of _projection_
parameters (Buja et al., 2019a, b),

$\beta=\Sigma^{-1}\Gamma,$

where $\Gamma=\mathbb{E}[XY]$. The projection parameter is traditionally
estimated using the ordinary least square estimator (OLSE). Suppose that we
observe a sample of $n$ i.i.d observations
$(X_{1},Y_{1}),\ldots,(X_{n},Y_{n})$ from $\mathbb{P}_{X,Y}$. Then, the OLSE
is defined as

$\displaystyle\widehat{\beta}=\operatorname*{arg\,min}_{\theta\in\mathbb{R}^{d}}\frac{1}{n}\sum_{i=1}^{n}(Y_{i}-\theta^{\top}X_{i})^{2}.$

Provided that the sample Gram matrix
$\widehat{\Sigma}:=n^{-1}\sum_{i=1}^{n}X_{i}X_{i}^{\top}$ is invertible with
probability 1, then the OLSE is well-defined and can be expressed as

$\displaystyle\widehat{\beta}=\widehat{\Sigma}^{-1}\widehat{\Gamma},$

where $\widehat{\Gamma}=n^{-1}\sum_{i=1}^{n}X_{i}Y_{i}$.

###  1.1 Related Works

##### Asymptotic Normality and Berry Esseen bound for the Least Square
Estimator

In fixed-dimensional settings, the OLSE has been the subject of large sample
theory and conventional Berry-Esseen bounds, demonstrated by Van der Vaart
(2000) and Pfanzagl (1973), respectively. With increasing dimensions, Bickel
and Freedman (1983) established the asymptotic normality of the least squares
and the consistency of bootstrap, requiring $d=o(\sqrt{n})$ under a well-
specified linear model with fixed covariates and homogeneous errors.
Additionally, a series of papers by Portnoy (1984, 1985, 1986, 1988) showcased
various consistency and asymptotic normality properties for estimators,
including M-estimators and maximum likelihood estimators. Translating these
findings into our context mandates $d=o(\sqrt{n})$ to ensure the
$\sqrt{n}$-consistency of the OLSE within a well-specified model under a set
of arguably strong assumptions on the data-generating process. In a more
general setting, albeit postulating at least a partially linear model, the
asymptotic normality of (non-)studentized OLSE has been investigated with
additional regularity conditions on covariates and errors (Donald and Newey,
1994; Cattaneo et al., 2015, 2018a, 2018b). This reference list is far from
complete. Some of these works pertain to the regime $d=o(n)$ and occasionally
allow the boundary case of $d=O(n)$.

An initial step toward the misspecified linear model with minimal assumptions
on covariates and errors was made in Mammen (1993). They were concerned with
linear contrasts of the regression coefficients and asserted the consistency
of both the resampling bootstrap and the wild bootstrap in estimating the
regression coefficient with increasing dimensions. They claimed that
$d=o(n^{3/4})$ is sufficient for consistent resampling bootstrap and that
$d=o(\sqrt{n})$ for the wild bootstrap. While their work provides useful
intuition regarding the asymptotic behavior of the OLSE under model
misspecification, their findings rely on heuristics and do not explicitly
present the convergence rate. Rinaldo et al. (2019) formulated the Berry
Esseen bound for general nonlinear statistics under the misspecified setting.
Their analysis, specifically applied to projection parameters, established a
uniform Berry Esseen bound for entry-wise asymptotic normality of the Ordinary
Least Square Estimator (OLSE), converging to $0$ under $d=o(n^{1/5})$,
disregarding polylogarithmic factors (refer to Theorem 2 therein). More
recently, Kuchibhotla et al. (2020) derived a novel finite sample bound for
the studentized OLSE under finite moment assumptions on covariates and errors.
They obtained a uniform Berry Esseen bound for entry-wise asymptotic normality
of the OLSE, with their bound scaling as $(d/n)^{1/3}+d/\sqrt{n}$,
disregarding polylogarithmic factors and under sufficient moment conditions on
covariates and errors; See Theorem 10 of Kuchibhotla et al. (2020).
Consequently, they required $d=o(\sqrt{n})$ for the bound to vanish,
facilitating the construction of simultaneous confidence intervals for the
projection parameter coefficients under such dimension range. Notably, these
intervals can exhibit a width of the order $1/\sqrt{n}$. It should be
emphasized that the parameters of interest in this article are linear
contrasts of the projection parameters which are different from those
addressed in Rinaldo et al. (2019) and Kuchibhotla et al. (2020).

##### De-biasing of the Least Square Estimator

Under misspecified linear models, extensive attention has focused on scenarios
where $d=o(\sqrt{n})$, with limited exploration beyond this scaling despite
the OLSE remaining well-defined for $d\leq n$. In cases where $d\gg\sqrt{n}$,
the OLSE demonstrates a bias of order $d/\sqrt{n}$ (Mammen, 1993; Cattaneo et
al., 2019), impeding valid inference on regression parameters without de-
biasing. To tackle this, Cattaneo et al. (2019) proposed the bias-corrected
$M$-estimator via the Jackknife method, ensuring consistency and asymptotic
normality of the two-step point estimate under $d=O(\sqrt{n})$, presenting
bootstrap-based inferential methods as well.

The term “de-biasing” typically refers to a correction applied to the original
regularized estimator, often necessary to resist the curse of dimensionality
or to enhance prior knowledge about the geometric/intrinsic structure of data.
Common examples in linear regression include the LASSO, ridge, and SLOP
estimators. It is worth emphasizing that “de-biasing” here aims to address the
“bias” mainly induced by model misspecification at the population level rather
than a specific regularization technique.

##### Variance Estimation

In the context of linear regression, the inflating variance of regression
coefficient estimates along with increasing dimension has prompted efforts to
propose robust covariance estimates capable of accommodating dimensionality
and potential heteroscedasticity. The degree-of-freedom-corrected covariance
estimator, introduced by MacKinnon and White (1985), has served as a
foundation for subsequent modifications within well-specified linear models,
known as HC-class variance estimators. For comprehensive reviews of these
estimators, we refer our readers to Long and Ervin (2000) and MacKinnon
(2012). Recent additions to the HC class include the HCK estimator proposed by
Cattaneo et al. (2018b) and the HCA estimator proposed by Jochmans (2022).
However, in both cases, consistency of variance estimate requires a
dimensionality constraint of $d=O(n^{1/2-\epsilon})$ in misspecified linear
models, albeit under different settings. Under a very similar setting to ours,
Kuchibhotla et al. (2020) noted that a dimensionality requirement of
$d=o(n^{1/2})$ appears to be unavoidable for the consistency of the sandwich
estimator for the OLSE variance (see Lemma 8 of Kuchibhotla et al. (2020)).

###  1.2 Contributions

In this paper, we employ an assumption-lean framework to offer a finite-sample
distribution approximation for the projection parameter, imposing minimal
assumptions. Specifically, our focus centers on estimating and inferring a
scalar parameter $\theta:=c^{\top}\beta$ in a predetermined direction
$c\in\mathbb{R}^{d}$. The main contributions of this study are summarized as
follows.

  * •

We propose a bias-corrected least square estimator for the projection
parameter. Theorem 3 presents the Berry Esseen bound for the approximation to
the adjusted Normal distribution of the unnormalized linear contrast of the
proposed estimator, and the bound is sharp and roughly scales like

$\displaystyle\frac{1}{\sqrt{n}}+\frac{d^{3/2}}{n},$ (1)

provided that $q_{x}\geq 12$ and $q\geq 4$, where $q_{x}$ and $q$ represent
the finite moments of the covariate $X$ and the error $Y-X^{\top}\beta$,
respectively; See also Corollary 1. Consequently, achieving the vanishing
bound necessitates the scaling $d=o(n^{2/3})$, strictly embracing the
traditional requirement of $d=o(n^{1/2})$. Our result, to the best of our
knowledge, presents the sharpest bound under the (possibly) misspecified
linear model while relying on the weakest assumptions regarding the true data
distribution.

  * •

While the Berry Esseen bound presented in Theorem 3 involves the unknown
parameter for the asymptotic variance of the proposed estimator, Section 4
outlines at least two inference methods for constructing (non-)asymptotically
valid confidence intervals. These methods do not necessitate estimating the
asymptotic variance, thereby enabling us to utilize the sharp Berry Esseen
bound from Theorem 3. Our approach involves inferential techniques based on
resampling and sample splitting, including HulC (Kuchibhotla et al., 2021),
$t-$statistic based inference (Ibragimov and Müller, 2010; Lam, 2022), and
wild bootstrap (Wu, 1986; Mammen, 1993). We also illustrate these proposed
methods through empirical examples in Section 5.

  * •

Our main result, Theorem 5, establishes a finite-sample Berry-Esseen bound for
a studentized linear contrast of the bias-corrected estimator under the
assumptions of finite moments for both the covariates and the errors. Our
bound, disregarding the polylogarithmic factors, roughly scales as follows:

$\displaystyle\frac{1}{n^{1/3}}+\sqrt{\frac{d}{n}}+\frac{d^{3/2}}{n},$ (2)

provided that $q_{x}\geq 12$ and $q\geq 4$; See also Corollary 2. Our bound
immediately leads to the validity of the Wald confidence interval for a linear
contrast of the projection parameters as outlined in Corollary 3.

The slower convergence rate outlined in (2) compared to that in (1) is solely
due to the additional convergence rate introduced by the sandwich variance
estimator. However, it is noteworthy that the consistency of the sandwich
variance estimator does not demand any additional dimension scaling beyond
$d=o(n^{2/3})$. Moreover, as per Theorem 4, the sandwich variance estimator
attains consistency if
$d=o(\min\\{n^{1-1/q_{x}-1/q},n^{1-2/q_{x}-2/(3q)}\\})$. To the best of our
knowledge, our result on the consistency of the sandwich variance estimator
allows the broadest range of dimension scaling in assumption lean linear
regression literature.

  * •

One of the key quantities in our analysis is
$\mathcal{D}_{\Sigma}=\|\Sigma^{-1/2}\widehat{\Sigma}\Sigma^{-1/2}-I_{d}\|_{\rm
op}$, representing the difference between the population and the sample gram
matrix in the operator norm. We introduce a new concentration inequality for
$\mathcal{D}_{\Sigma}$ under finite moment assumptions on covariates
($q_{x}\geq 4$), detailed in Proposition 22. This result is inspired by recent
findings on the universal spectral properties of random matrices, relating
them to Gaussian random matrices with identical mean and covariance in their
entries (Brailovskaya and van Handel, 2023). Our bound for
$\mathcal{D}_{\Sigma}$ matches the sharpness (up to polylogarithmic factors)
of the one provided in Tikhomirov (2018), regarded as the sharpest to date.
Notably, our bound avoids the necessity for a ‘large’ $d$ to ensure a ‘high’
probability concentration of the sample Gram matrix. Section B.3 outlines the
historical context of $\mathcal{D}_{\Sigma}$ and offers a comparison of
concentration inequalities in this regard.

###  1.3 Organization

The remainder of this paper is organized as follows: Section 2.1 introduces
the problem setup and notations, while Section 2.2 outlines the distributional
assumptions on the data generating process. Our primary result on the Berry
Esseen bound for the distribution approximation of the law of the studentized
linear contrast of the proposed estimator is presented in Section 3. Section 4
details several inferential methods for the projection parameter based on the
bootstrap and sample splitting. Numerical results are provided in Section 5.
The Appendix includes all proofs, technical lemmas, and additional numerical
results.

##  2 Problem setup and Assumptions

###  2.1 Projection parameters

Let
$(X_{1},Y_{1}),\ldots,(X_{n},Y_{n})\sim(X,Y)\in\mathbb{R}^{d}\times\mathbb{R}$
be an independent sample of $n$ observations. The _projection parameter_
$\beta$ is defined as

$\beta=\beta_{n}:=\operatorname*{arg\,min}_{\theta\in\mathbb{R}^{d}}\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}[(Y_{i}-X_{i}^{\top}\theta)^{2}].$
(3)

The minimizer $\beta$ in (3) solves the linear equation that

$\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}[X_{i}(Y_{i}-X_{i}^{\top}\beta)]=0.$ (4)

If the population gram matrix
$\Sigma=\Sigma_{n}:=n^{-1}\sum_{i=1}^{n}\mathbb{E}[X_{i}X_{i}^{\top}]$ is
invertible, the projection parameter is well-defined and can be uniquely
expressed as $\beta=\Sigma^{-1}\Gamma$ where
$\Gamma=\Gamma_{n}:=n^{-1}\sum_{i=1}^{n}\mathbb{E}[X_{i}Y_{i}]$. In the case
where the observations follow the linear model $Y=X^{\top}\beta^{*}+\epsilon$
with $\mathbb{E}[\epsilon|X]=0$, the projection parameter corresponds to the
model parameter, i.e., $\beta=\beta^{*}$. On the contrary, if the underlying
truth exhibits a possibly non-linear structure, then the projection parameter
gives the best linear approximation $X^{\top}\beta$ to $Y$ with respect to the
joint distribution of $(X,Y)$. For a detailed discussion on the projection
parameter and its interpretation, see Buja et al. (2019a, b).

The projection parameter is estimated using the ordinary least square
estimator (OLSE) which follows the traditional definition that

$\widehat{\beta}=\widehat{\beta}_{n}:=\operatorname*{arg\,min}_{\theta\in\mathbb{R}^{d}}\frac{1}{n}\sum_{i=1}^{n}(Y_{i}-X_{i}^{\top}\theta)^{2}.$

If the sample population matrix
$\widehat{\Sigma}=\widehat{\Sigma}_{n}:=n^{-1}\sum_{i=1}^{n}X_{i}X_{i}^{\top}$
is invertible with probability 1, then the OLSE is well defined and can be
written as $\widehat{\beta}=\widehat{\Sigma}^{-1}\widehat{\Gamma}$ where
$\widehat{\Gamma}=\widehat{\Gamma}_{n}:=n^{-1}\sum_{i=1}^{n}X_{i}Y_{i}$.

##### Additional Notations

For two real numbers $a$ and $b$, let $a\vee b=\max\\{a,b\\}$, $a\wedge
b=\min\\{a,b\\}$ and $a_{+}=\max\\{a,0\\}$. Let $I_{p}$ be the $p\times p$
identity matrix. We let $1_{p}$ be the $p$-dimensional vector of ones and
$0_{p}$ be of zeros. For any $x\in\mathbb{R}^{d}$, we write
$\|x\|_{2}=\sqrt{x^{\top}x}$. In addition with a positive definite matrix
$A\in\mathbb{R}^{d\times d}$, we denote the scaled Euclidean norm as
$\|x\|_{A}=\sqrt{x^{\top}Ax}$. We let $\|\cdot\|_{\rm op}$ be the operator
norm of the matrices. The unit sphere in $\mathbb{R}^{d}$ is
$\mathbb{S}^{d-1}=\left\\{\theta\in\mathbb{R}^{d}:\|\theta\|_{2}=1\right\\}$.
The $j$:th canonical basis of $\mathbb{R}^{d}$ is written as $\mathbf{e}_{j}$
for all $j=1,\ldots,d$. Convergence in distribution is denoted by
$\xrightarrow{D}$ and convergence in probability by $\xrightarrow{P}$. We
write the smallest eigenvalue and the largest eigenvalue of the square
matrices as $\lambda_{\min}(\cdot)$ and $\lambda_{\max}(\cdot)$, respectively.
The cumulative distribution function of the standard Normal distribution is
denoted as $\Phi(\cdot)$ and its $k$:th derivative as $\Phi^{(k)}(\cdot)$ for
$k\in\mathbb{N}$. Finally, we adopt the conventional asymptotic notations.
Given two non-negative sequences of reals $\\{a_{n}\\}$ and $\\{b_{n}\\}$, we
write $a_{n}=O(b_{n})$ if there exists a constant $C\in(0,\infty)$ and
$N\in\mathbb{N}$ such that for any $n\geq N$, $a_{n}\leq Cb_{n}$, and
$a_{n}=o(b_{n})$ if $a_{n}/b_{n}\to 0$ as $n\to\infty$. We sometimes write
$a_{n}\lesssim b_{n}$ for $a_{n}=O(b_{n})$, and $a_{n}\ll b_{n}$ for
$a_{n}=o(b_{n})$. Furthermore, for a sequence of random variables
$\left\\{\gamma_{n}\right\\}$ and a non-zero sequence $\\{c_{n}\\}$, we write
$\gamma_{n}=o_{P}(c_{n})$ if $|\gamma_{n}/c_{n}|\xrightarrow{P}0$ as
$n\to\infty$, and $\gamma_{n}=O_{P}(c_{n})$ if for any $\epsilon>0$, there
exists $C_{\epsilon}\in(0,\infty)$ and $N_{\epsilon}\in\mathbb{N}$ such that
$\mathbb{P}(|\gamma_{n}/c_{n}|>C_{\epsilon})<\epsilon$ for all $n\geq
N_{\epsilon}$.

###  2.2 Assumptions

Our main results are built upon the following assumptions, primarily
emphasizing minimal assumptions on the data generating process. These include
moment conditions on both the response and covariates, thereby encompassing a
wide range of distributions, such as heavy-tailed distributions and discrete
distributions.

######  Assumption 1.

There exists some $q\geq 2$ and a constant $K_{y}>0$ satisfying that

$\left(\mathbb{E}[|Y_{i}-X_{i}^{\top}\beta|^{q}]\right)^{1/q}\leq K_{y},$

for $i=1,\ldots,n$.

Assumption 1 imposes a moment condition solely on the error, allowing for
heteroscedastic errors that may arbitrarily depend on the covariates and even
accommodate heavy-tailed errors. This assumption is prevalent in recent
literature on heteroskedastic linear regression (Mourtada et al., 2022;
Kuchibhotla and Patra, 2022; Han and Wellner, 2019; Mendelson, 2016).
Additionally, it is often assumed that the conditional variance of the error
is uniformly bounded almost everywhere, described by

$0<\inf_{x}v^{2}(x)\leq\sup_{x}v^{2}(x)<\infty,$ (5)

where $v^{2}(x)=\mathbb{E}[(Y-X^{\top}\beta)^{2}|X=x]$. Here, $\inf_{x}$ and
$\sup_{x}$ are essential infimum and essential supremum, respectively, with
respect to the marginal distribution of $X$. Given (5), our results in the
following section remain valid even under the relaxed moment assumptions on
the covariates.

######  Assumption 2.a.

There exists some $q_{x}\geq 2$ and a constant $K_{x}\geq 1$ satisfying that

$\left(\mathbb{E}[|u^{\top}\Sigma^{-1/2}X_{i}|^{q_{x}}]\right)^{1/q_{x}}\leq
K_{x},$

for $i=1,\ldots,n$ and $u\in\mathbb{S}^{d-1}$.

######  Assumption 2.b.

There exists a constant $K_{x}>0$ such that

$\mathbb{E}\left[\exp\left(\frac{|u^{\top}\Sigma^{-1}X_{i}|^{2}}{K_{x}^{2}}\right)\right]\leq
2,$

for all $i=1,\ldots,n$ and $u\in\mathbb{S}^{d-1}$.

######  Assumption 2.c.

There exists some $q_{x}\geq 2$ and a constant $K_{x}\geq 1$ satisfying that

$\left(\mathbb{E}[|u^{\top}\Sigma^{-1/2}X_{i}|^{q_{x}}]\right)^{1/q_{x}}\leq
K_{x},$

for $i=1,\ldots,n$ and $u\in\mathbb{S}^{d-1}$. Furthermore,
$\Sigma^{-1/2}X_{i}\in\mathbb{R}^{d}$, $i=1,\ldots,n$, have $d$ independent
entries.

Assumption 2.a implies $L^{q_{x}}-L^{2}$ equivalence for one-dimensional
marginals of covariates. This is a significant weakening of sub-Gaussianity in
Assumption 2.b. It is common to impose $q_{x}\geq 4$ when studying the linear
regression with increasing dimensions (Oliveira, 2016; Catoni, 2016; Mourtada
et al., 2022; Mourtada, 2022; Vaškevičius and Zhivotovskiy, 2023) even though
we sometimes require a higher moment condition to analyze a higher order
expansion of the least square estimate.

Furthermore, Assumption 2.a implies the bounded-ness of $q_{x}$:th moment of
the $\Sigma^{-1}$-norm of the random vector $X_{i}$. To illustrate, note that
$\|X_{i}\|_{\Sigma^{-1}}^{2}=\sum_{j=1}^{d}(\mathbf{e}_{j}^{\top}\Sigma^{-1/2}X_{i})^{2}$
for all $i=1,\ldots,n$. By applying Jensen’s inequality, we have:

$\left(\frac{\|X_{i}\|_{\Sigma^{-1}}^{2}}{d}\right)^{q_{x}/2}\leq\frac{1}{d}\sum_{j=1}^{d}\left(\mathbf{e}_{j}^{\top}\Sigma^{-1/2}X_{i}\right)^{q_{x}},$

which holds deterministically as long as $q_{x}\geq 2$. Consequently, it leads
to

$\mathbb{E}\left[\left(\frac{\|X_{i}\|_{\Sigma^{-1}}}{\sqrt{d}}\right)^{q_{x}}\right]\leq\frac{1}{d}\sum_{j=1}^{d}\mathbb{E}\left[(\mathbf{e}_{j}^{\top}\Sigma^{-1/2}X_{i})^{q_{x}}\right]\leq\max_{1\leq
j\leq
d}\mathbb{E}\left[(\mathbf{e}_{j}^{\top}\Sigma^{-1/2}X_{i})^{q_{x}}\right]\leq
K_{x}^{q_{x}},$

where the last inequality follows from Assumption 2.a. This also emphasizes
that the appropriate scaling of $\|X_{i}\|_{\Sigma^{-1}}$ is proportional to
$\sqrt{d}$ as the dimension $d$ grows.

Assumption 2.c assumes that $\Sigma^{-1/2}$-normalized covariates have
independent entries, in addition to the moment conditions. Under Assumption
2.c, we can achieve a more precise control over the $L^{2}$-norm of covariates
(Rio, 2017) and the sample Gram matrix (Srivastava and Vershynin, 2013;
Mendelson and Paouris, 2014).

######  Assumption 3.

Let

$V:={\rm
Var}\left[\frac{1}{\sqrt{n}}\sum_{i=1}^{n}X_{i}(Y_{i}-X_{i}^{\top}\beta)\right]={\rm
Var}\left[X_{1}(Y_{1}-X_{1}^{\top}\beta)\right].$

There exist constants $0<\underline{\lambda}\leq\overline{\lambda}<\infty$
satisfying that

$\underline{\lambda}\leq\lambda_{\rm
min}(\Sigma^{-1/2}V\Sigma^{-1/2})\leq\lambda_{\rm
max}(\Sigma^{-1/2}V\Sigma^{-1/2})\leq\overline{\lambda}.$

Assumption 3 ensures that both $V$ and $\Sigma$ scale in a similar order,
serving as the only technical assumption in our paper. This condition appears
inevitable to prevent the asymptotic variance of the least square estimator
from inflating or vanishing. Specifically, noting that
$V=\mathbb{E}[XX^{\top}(Y-X^{\top}\beta)^{2}]=\mathbb{E}[XX^{\top}\mathbb{E}[(Y-X^{\top}\beta)^{2}|X]]$
the condition (5) is sufficient for Assumption 3 to be held.

##  3 Methods and Results

###  3.1 Approximation of the Sample Gram Matrix and Bias Characterization

In this section, we present an approximation of the distribution of
$\sqrt{n}c^{\top}(\widehat{\beta}-\beta)$ for any
$c\in\mathbb{R}^{d}\setminus\\{0_{d}\\}$. To be presented results are mostly
invariant to the scaling of $c\in\mathbb{R}^{d}$. However, we sometimes
consider those with $\|c\|_{\Sigma^{-1}}=1$, which gives a proper
normalization,
$\sqrt{n}c^{\top}(\widehat{\beta}-\beta)=\sqrt{n}\tilde{c}^{\top}\Sigma^{1/2}(\widehat{\beta}-\beta)$
for some $\tilde{c}\in\mathbb{S}^{d-1}$.

Our analysis starts with the least square estimate $\widehat{\beta}$,
particularly by approximating the sample Gram matrix $\widehat{\Sigma}$ using
Taylor series expansion. We note

$\displaystyle\widehat{\Sigma}^{-1}$ $\displaystyle=$
$\displaystyle\Sigma^{-1}+\widehat{\Sigma}^{-1}(\Sigma-\widehat{\Sigma})\Sigma^{-1}$
$\displaystyle\approx$
$\displaystyle\Sigma^{-1}+\Sigma^{-1}(\Sigma-\widehat{\Sigma})\Sigma^{-1}.$

(See Lemma 33 of Kuchibhotla et al. (2020).) Consequently, we get the
following approximation of $\widehat{\beta}-\beta$:

$\displaystyle\widehat{\beta}-\beta$ $\displaystyle=$
$\displaystyle\widehat{\Sigma}^{-1}\frac{1}{n}\sum_{i=1}^{n}X_{i}(Y_{i}-X_{i}^{\top}\beta)$
(7) $\displaystyle\approx$
$\displaystyle\Sigma^{-1}\frac{1}{n}\sum_{i=1}^{n}X_{i}(Y_{i}-X_{i}^{\top}\beta)$
$\displaystyle+\,\Sigma^{-1}(\Sigma-\widehat{\Sigma})\Sigma^{-1}\frac{1}{n}\sum_{i=1}^{n}X_{i}(Y_{i}-X_{i}^{\top}\beta).$

The term (7) is the first-order approximation and is the average of influence
function—with respect to the square loss—evaluated at independent data points.
Using the conventional large sample approximation for the case when
$d=o(\sqrt{n})$, the variability of $\widehat{\beta}$ can be well-represented
only through the term (7), and the second order term (7) or higher order terms
are negligible. In particular, (7) is responsible for the asymptotic normality
of $\sqrt{n}(\widehat{\beta}-\beta)$ when $d=o(\sqrt{n})$. Here, however, the
term (7) may create a non-vanishing bias when $d\gg\sqrt{n}$. To see this, we
note that

$\displaystyle\Sigma^{-1}(\Sigma-\widehat{\Sigma})\Sigma^{-1}\frac{1}{n}\sum_{i=1}^{n}X_{i}(Y_{i}-X_{i}^{\top}\beta)$
$\displaystyle=$
$\displaystyle\left(\frac{1}{n}\sum_{i=1}^{n}\Sigma^{-1}(\Sigma-
X_{i}X_{i}^{\top})\Sigma^{-1}\right)\left(\frac{1}{n}\sum_{i=1}^{n}X_{i}(Y_{i}-X_{i}^{\top}\beta)\right)$
(8) $\displaystyle=$
$\displaystyle\frac{1}{n^{2}}\sum_{i=1}^{n}\Sigma^{-1}X_{i}(Y_{i}-X_{i}^{\top}\beta)$
(11)
$\displaystyle-\frac{1}{n^{2}}\sum_{i=1}^{n}\Sigma^{-1}X_{i}(Y_{i}-X_{i}^{\top}\beta)\|X_{i}\|_{\Sigma^{-1}}^{2}$
$\displaystyle+\frac{1}{n^{2}}\sum_{1\leq i\neq j\leq n}\Sigma^{-1}(\Sigma-
X_{i}X_{i}^{\top})\Sigma^{-1}X_{j}(Y_{j}-X_{j}^{\top}\beta).$

The quantities in (11) and (11) are resulted from multiplying the individual
terms in (8) with the same indices, while the quantity (11) corresponds to the
cross-index product. It is noteworthy that (11) can be absorbed into the
average of influence functions in (7) as it has an additional scaling factor
$1/n$. The quantity (11) is a mean-zero degenerate $U$-statistic of order $2$
up to a scaling factor that converges to 1 as $n\to\infty$. Therefore, the
only quantity that may have a non-zero mean is (11) and we denote this by

$\mathcal{B}:=-\frac{1}{n^{2}}\sum_{i=1}^{n}\Sigma^{-1}X_{i}(Y_{i}-X_{i}^{\top}\beta)\|X_{i}\|_{\Sigma^{-1}}^{2}.$
(12)

Under the control of $c^{\top}\Sigma^{-1}X_{i}$ and $Y_{i}-X_{i}^{\top}\beta$
using Assumption 2.a and 1, respectively, the quantity $c^{\top}\mathcal{B}$
scales like $O(d/n)$ due the $\|X_{i}\|_{\Sigma^{-1}}^{2}$ factor.
Consequently, $\sqrt{n}c^{\top}\mathcal{B}$ yields a non-degenerate bias when
$d\gg\sqrt{n}$. By manually removing the bias, a linear contrast
$c^{\top}(\widehat{\beta}-\mathcal{B}-\beta)$ can be approximated by the sum
of the average of influence functions and the degenerate second-order
$U$-statistics as

$c^{\top}(\widehat{\beta}-\mathcal{B}-\beta)~{}\approx~{}\frac{1}{n}\sum_{i=1}^{n}c^{\top}\psi(X_{i},Y_{i})+\frac{1}{n(n-1)}\sum_{1\leq
i\neq j\leq n}c^{\top}\phi(X_{i},Y_{i},X_{j},Y_{j}),$

where

$\psi(x,y)=\left(1+\frac{1}{n}\right)\Sigma^{-1}x(y-x^{\top}\beta)\quad\mbox{and}\quad\phi(x,y,x^{\prime},y^{\prime})=\left(1-\frac{1}{n}\right)\Sigma^{-1}(\Sigma-
xx^{\top})\Sigma^{-1}x^{\prime}(y^{\prime}-{x^{\prime}}^{\top}\beta).$ (13)

It is well known that $U$-statistics of order $k\geq 2$ are asymptotically
normally distributed, and there has been a vast literature related to Normal
approximations and the rates of convergence for $k$-order $U$-statistics. We
refer the readers to Bentkus et al. (2009). In order to state our first
theorem, we define

$\displaystyle\sigma_{c}^{2}$ $\displaystyle=$ $\displaystyle{\rm
Var}\left[c^{\top}\Sigma^{-1}X(Y-X^{\top}\beta)\right]=c^{\top}\Sigma^{-1}V\Sigma^{-1}c,$
(14) $\displaystyle\kappa_{c}$ $\displaystyle=$ $\displaystyle
n^{-5/2}\binom{n}{2}\mathbb{E}\left[\\{c^{\top}\psi(X_{1},Y_{1})\\}\\{c^{\top}\psi(X_{2},Y_{2})\\}\\{c^{\top}\phi(X_{1},Y_{1},X_{2},Y_{2})\\}\right]/\sigma_{c}^{3},$
(15)

for any $c\in\mathbb{R}^{d}\setminus\\{0_{d}\\}$. It should be emphasized that
the parameter $\kappa_{c}$ is invariant to the scaling of
$c\in\mathbb{R}^{d}\setminus\\{0_{d}\\}$. That is, $\kappa_{c}=\kappa_{\alpha
c}$ for $\alpha\neq 0$. In addition, recall that $\Phi(x)$ is the cumulative
distribution function of the standard Normal distribution and $\Phi^{(3)}(x)$
is the third order derivative of $\Phi(x)$.

Theorem 1 states a Berry-Essen bound for the Normal approximation of a linear
contrast $c^{\top}(\widehat{\beta}-\mathcal{B}-\beta)$ with known bias
$\mathcal{B}$.

######  Theorem 1.

Suppose that Assumption 1 holds for $q>3$ and that Assumption 3 holds.

  * (i)

Suppose further that Assumption 2.a holds for $q_{x}\geq 4$ and
$s:=(1/q+1/q_{x})^{-1}\geq 3$, and assume $d+2\log(2n)\leq n/(18K_{x})^{2}$.
Then, there exists a constant
$C=C(q,q_{x},K_{x},K_{y},\overline{\lambda},\underline{\lambda})$ such that

$\displaystyle\sup_{c\in\mathbb{R}^{d}\setminus\\{0_{d}\\}}\sup_{t\in\mathbb{R}}\left|\mathbb{P}\left[\sqrt{n}\left\\{c^{\top}(\widehat{\beta}-\mathcal{B}-\beta)\right\\}\leq
t\right]-\left\\{\Phi\left(\frac{t}{\sigma_{c}}\right)+\kappa_{c}\Phi^{(3)}\left(\frac{t}{\sigma_{c}}\right)\right\\}\right|$
$\displaystyle\leq$
$\displaystyle\frac{C}{\sqrt{n}}+C\left\\{\left(\frac{d\log^{3}(2n)}{n^{4/5-8/(5q_{x})}}\right)^{5q_{x}/(2q_{x}+8)}+\frac{d^{3/2}\log(2n)}{n}\right\\}\left(1\vee\sqrt{\frac{\log(2n)}{d}}\right).$

  * (ii)

Suppose further that Assumption 2.b holds. Then, there exists a constant
$C=C(q,K_{x},K_{y},\overline{\lambda},\underline{\lambda})$ such that

$\displaystyle\sup_{c\in\mathbb{R}^{d}\setminus\\{0_{d}\\}}\sup_{t\in\mathbb{R}}\left|\mathbb{P}\left[\sqrt{n}\left\\{c^{\top}(\widehat{\beta}-\mathcal{B}-\beta)\right\\}\leq
t\right]-\left\\{\Phi\left(\frac{t}{\sigma_{c}}\right)+\kappa_{c}\Phi^{(3)}\left(\frac{t}{\sigma_{c}}\right)\right\\}\right|$
$\displaystyle\leq$ $\displaystyle
C\left(\frac{1}{\sqrt{n}}+\frac{d^{3/2}}{n}+\frac{d\log(2n)}{n}\right).$

  * (iii)

Suppose further that Assumption 2.c holds for $q_{x}>4$ and $s\geq 3$, and
assume $d+2\log(2n)\leq n/(18K_{x})^{2}$. Then, there exists a constant
$C=C(q,q_{x},K_{x},K_{y},\overline{\lambda},\underline{\lambda})$ such that

$\displaystyle\sup_{c\in\mathbb{R}^{d}\setminus\\{0_{d}\\}}\sup_{t\in\mathbb{R}}\left|\mathbb{P}\left[\sqrt{n}\left\\{c^{\top}(\widehat{\beta}-\mathcal{B}-\beta)\right\\}\leq
t\right]-\left\\{\Phi\left(\frac{t}{\sigma_{c}}\right)+\kappa_{c}\Phi^{(3)}\left(\frac{t}{\sigma_{c}}\right)\right\\}\right|$
$\displaystyle\leq$
$\displaystyle\frac{C}{\sqrt{n}}+C\left\\{\left(\frac{d^{3/2}\log^{3}(2n)}{n}\right)^{q_{x}/(q_{x}+4)}+\frac{d^{3/2}\log(2n)}{n}\right\\}\left(1\vee\sqrt{\frac{\log(2n)}{d}}\right).$

Furthermore, if $(3/q_{x}+1/q)^{-1}\geq 2$, then there exists a constant
$C=C(q,q_{x},K_{x},K_{y},\overline{\lambda},\underline{\lambda})$ such that

$\displaystyle\sup_{c\in\mathbb{R}^{d}\setminus\\{0_{d}\\}}\sup_{t\in\mathbb{R}}\left|\mathbb{P}\left[\sqrt{n}\left\\{c^{\top}(\widehat{\beta}-\beta)\right\\}\leq
t\right]-\left\\{\Phi\left(\frac{t}{\sigma_{c}}\right)+\kappa_{c}\Phi^{(3)}\left(\frac{t}{\sigma_{c}}\right)\right\\}\right|$
$\displaystyle\leq$ $\displaystyle
C\sqrt{\frac{d}{n}}+C\left\\{\left(\frac{d^{3/2}\log^{3}(2n)}{n}\right)^{q_{x}/(q_{x}+4)}+\frac{d^{3/2}\log(2n)}{n}\right\\}\left(1\vee\sqrt{\frac{\log(2n)}{d}}\right).$

######  Remark 1.

The conventional Berry Esseen type inequality implies the proximity between
the distribution of the law of the estimator and the Normal distribution. In
contrast, we show a distribution approximation to the “adjusted” Normal
distribution, allowing a more precise representation of the approximation. The
degree of adjustment is determined by the parameter $\kappa_{c}$, which turns
out to scale as $O(\sqrt{d/n})$ (see Lemma 16). Hence, we can enhance the
Berry-Esseen bound for the Normal approximation by incorporating an additional
term of $O(\sqrt{d/n})$ to the existing bounds in Theorem 1, but yielding a
slower convergence rate. More importantly, under the well-specified linear
model, i.e., $\mathbb{E}[Y|X]=X^{\top}\beta$, then it follows that
$\kappa_{c}=0$ from its definition in (15). Consequently, in this case, the
presented Berry-Esseen bound in Theorem 1 can be employed for the Normal
approximation as well.

######  Remark 2.

The last part of Theorem 1 says that the OLSE does not necessarily require de-
biasing under Assumption 2.c with a sufficiently large finite moment of
covariates and errors. To provide an intuition for this, let us consider the
bias term $c^{\top}\mathcal{B}$ in (12), which scales with its expected value
$\mathbb{E}[c^{\top}\mathcal{B}]$. By examining the definition of the
projection parameter, we can express $\mathbb{E}[c^{\top}\mathcal{B}]$ as the
covariance between two random variables:

$\displaystyle\mathbb{E}[c^{\top}\mathcal{B}]$ $\displaystyle=$
$\displaystyle-\frac{1}{n}\mathbb{E}[c^{\top}\Sigma^{-1}X(Y-X^{\top}\beta)\|X\|_{\Sigma^{-1}}^{2}]$
$\displaystyle=$ $\displaystyle-\frac{1}{n}{\rm
Cov}\left(c^{\top}\Sigma^{-1}X(Y-X^{\top}\beta),\|X\|_{\Sigma^{-1}}^{2}\right).$

An application of Cauchy Schwarz’s inequality yields that

$\left|\mathbb{E}[c^{\top}\mathcal{B}]\right|\leq\frac{1}{n}\left({\rm
Var}\left[c^{\top}\Sigma^{-1}X(Y-X^{\top}\beta)\right]\right)^{1/2}\left({\rm
Var}\left[\|X\|_{\Sigma^{-1}}^{2}\right]\right)^{1/2}.$ (16)

The leading variance term on the right-hand side of (16), under the finite
moment assumptions of covariates and errors, is $O(1)$. Meanwhile, the
quantity ${\rm Var}[\|X\|_{\Sigma^{-1}}^{2}]$ generally scales as
$\mathbb{E}[\|X\|_{\Sigma^{-1}}^{4}]=O(d^{2})$, resulting in
$|\mathbb{E}[c^{\top}\mathcal{B}]|=O(d/n)$. However, when an additional
assumption of independent entries of covariates, i.e., Assumption 2.c, is
satisfied, $\|X\|_{\Sigma^{-1}}^{2}$ roughly follows $\chi^{2}_{d}$
distribution, thereby leading to ${\rm Var}[\|X\|_{\Sigma^{-1}}^{2}]=O(d)$.
This significantly reduces the magnitude of the right-hand side of (16) to
$O(\sqrt{d}/n)$. Consequently, the expected bias (when scaled by $\sqrt{n}$)
is negligible if $d=o(n)$.

######  Remark 3.

Regarding the results under Assumptions 2.a and 2.c, we conjecture that the
dependence on logarithmic factors of our bound is suboptimal. The presence of
seemingly spurious logarithmic factors mainly stems from the concentration
inequality we employ for the sample Gram matrix, described in Proposition 22
in Section B.3. Although Theorem 1.1 of Tikhomirov (2018) partially addresses
our query when $q_{x}>4$, its application to our context (see Proposition 23)
requires further supposition: $d\to\infty$ as $n\to\infty$.

Theorem 1 gives the Berry Essen bound for the linear contrast of the
_properly-debiased_ projection parameters under three different distributional
assumptions on covariates. The first part of Theorem 1 regards the most
general case which relies only on finite-moment assumptions. Ignoring the
polylogarithmic factors, the Berry-Esseen bound tends to zero as long as

$d=o\left(n^{\min\\{2/3,4/5-8/(5q_{x})\\}}\right),\quad n\to\infty.$

If $q_{x}\geq 12$, this reduces to $d=o(n^{2/3})$ and meets the dimension
requirement for the vanishing Berry Esseen bound under Assumption 2.b of sub-
Gaussianity. Moreover, under Assumption 2.c, the bound tends to 0 as
$d=o(n^{2/3})$ for all $q_{x}>4$, disregarding the polylogarithmic factors.

###  3.2 Bias Estimation

The result given in Theorem 1 is mostly of theoretical importance as the
inequality yet involves unknown parameters through bias $\mathcal{B}$ and
asymptotic variance $\sigma_{c}^{2}$. Therefore, to do inference on $\beta$,
for e.g. building confidence region, it is necessary to incorporate the
estimates of bias and variance in the Berry Esseen bound given in Theorem 1.
To this effect, we consider a method of moment estimator
$\widehat{\mathcal{B}}$ for $\mathcal{B}$ defined as

$\widehat{\mathcal{B}}=-\frac{1}{n^{2}}\sum_{i=1}^{n}\widehat{\Sigma}^{-1}X_{i}(Y_{i}-X_{i}^{\top}\widehat{\beta})\|X_{i}\|_{\widehat{\Sigma}^{-1}}^{2}.$
(17)

The next result provides the consistency rate for the bias estimate in high
dimensions and under mild moment conditions. As elucidated in Remark 2, the
OLSE does not necessitate a debiasing procedure under Assumption 2.c with a
sufficient number of moments of covariates and errors. Therefore, the
following results will consider the cases under Assumption 2.a and 2.b only.

######  Theorem 2.

Suppose that Assumption 1 holds for $q>3$ and that Assumption 3 holds. Fix
$c\in\mathbb{R}^{d}$ with $\|c\|_{\Sigma^{-1}}=1$, and assume $d+2\log(2n)\leq
n/(18K_{x})^{2}$.

  * (i)

Suppose further that Assumption 2.a holds for $q_{x}\geq 8$,
$(1/q_{x}+1/q)^{-1}\geq 3$, and $(3/q_{x}+1/q)^{-1}\geq 2$. Then there exists
a constant $C=C(q,q_{x},\overline{\lambda},K_{x},K_{y})$ such that for all
$\eta\in(0,1]$,

$\displaystyle\mathbb{P}\left(\sqrt{n}|c^{\top}(\widehat{\mathcal{B}}-\mathcal{B})|\geq\delta_{n,d}(\eta)\right)\leq\frac{1}{\sqrt{n}}+\delta_{n,d}(\eta)+\eta,$

where

$\delta_{n,d}(\eta):=C\left(\frac{d^{2}\log^{3/2}(2d/\eta)}{\eta^{2/q_{x}}n^{3/2-2/q_{x}}}+\frac{d^{3}\log^{3}(2d/\eta)}{\eta^{4/q_{x}}n^{5/2-4/q_{x}}}+\frac{d^{3/2}\log(2n/\eta)}{n}\right).$

  * (ii)

Suppose further that Assumption 2.b holds. Then, there exists a constant
$C=C(q,\overline{\lambda},K_{x},K_{y})$ such that

$\displaystyle\mathbb{P}\left(\sqrt{n}|c^{\top}(\widehat{\mathcal{B}}-\mathcal{B})|\geq\delta_{n,d}\right)\leq\frac{1}{\sqrt{n}}+\delta_{n,d},$

where

$\delta_{n,d}:=C\log(2d)\left[\frac{d^{3/2}}{n}+\frac{d\log(2n)}{n}\right].$

Theorem 2 provides the finite sample concentration inequalities for the bias
estimate. Inspecting the result under Assumption 2.a, with the choice of

$\eta=\left(\frac{d}{n^{3/4-1/q_{x}}}\right)^{(2q_{x})/(2+q_{x})}\vee\left(\frac{d}{n^{5/6-4/(3q_{x})}}\right)^{(3q_{x})/(4+q_{x})},$

which was intended to minimize the tail probability, the bias estimate is
$\sqrt{n}-$consistent if $d=o(n^{\min\\{2/3,3/4-1/q_{x}\\}})$, ignoring the
polylogarithmic factors. If $q_{x}\geq 12$, then the dimension requirement
reduces to $d=o(n^{2/3})$. If Assumption 2.a is replaced with Assumption 2.b,
then the dimension requirement for the consistency becomes $d=o(n^{2/3})$
again ignoring the polylogarithmic factors.

###  3.3 Berry Esseen Bound for the Unnormalized Bias-corrected Estimator

Combining the error bound for the bias estimate from Theorem 2 with the Berry-
Esseen bound in Theorem 1, we get to the next result, a uniform Berry-Esseen
bound for the unnormalized bias-corrected OLS estimator. We denote the bias-
corrected estimator for the projection parameters as

$\displaystyle\widehat{\beta}_{\rm
bc}:=\widehat{\beta}-\widehat{\mathcal{B}}.$

######  Theorem 3.

  * (i)

Suppose that assumptions made in Theorem 2(i) hold. Then, there exists a
constant $C=C(q,q_{x},K_{x},K_{y},\overline{\lambda},\underline{\lambda})$
such that

$\displaystyle\sup_{c\in\mathbb{R}^{d}\setminus\\{0_{d}\\}}\sup_{t\in\mathbb{R}}\left|\mathbb{P}\left[\sqrt{n}c^{\top}\left(\widehat{\beta}_{\rm
bc}-\beta\right)\leq
t\right]-\left\\{\Phi\left(\frac{t}{\sigma_{c}}\right)+\kappa_{c}\Phi^{(3)}\left(\frac{t}{\sigma_{c}}\right)\right\\}\right|\leq\epsilon_{n,d},$

where

$\displaystyle\epsilon_{n,d}$ $\displaystyle=$ $\displaystyle
C\left\\{\left(\frac{d\log^{3}(2n)}{n^{4/5-8/(5q_{x})}}\right)^{5q_{x}/(2q_{x}+8)}+\left(\frac{d\log^{3/4}(2n)}{n^{3/4-1/q_{x}}}\right)^{(2q_{x})/(2+q_{x})}+\left(\frac{d\log(2n)}{n^{5/6-4/(3q_{x})}}\right)^{(3q_{x})/(4+q_{x})}\right\\}$
$\displaystyle+C\left(\frac{1}{\sqrt{n}}+\frac{d^{3/2}\log(2n)}{n}\right).$

  * (ii)

Suppose that the assumptions made in Theorem 2(ii) hold. Then, there exists a
constant $C=C(q,K_{x},K_{y},\overline{\lambda},\underline{\lambda})$ such that

$\displaystyle\sup_{c\in\mathbb{R}^{d}\setminus\\{0_{d}\\}}\sup_{t\in\mathbb{R}}\left|\mathbb{P}\left[\sqrt{n}c^{\top}\left(\widehat{\beta}_{\rm
bc}-\beta\right)\leq
t\right]-\left\\{\Phi\left(\frac{t}{\sigma_{c}}\right)+\kappa_{c}\Phi^{(3)}\left(\frac{t}{\sigma_{c}}\right)\right\\}\right|$
$\displaystyle\leq$ $\displaystyle
C\left(\frac{1}{\sqrt{n}}+\frac{d^{3/2}\log(2n)}{n}\right).$

Theorem 3 presents the Berry Esseen bound for approximating the law of the
unnormalized linear contrast of the bias-corrected LSE — with respect to any
query vector $c\in\mathbb{R}^{d}\setminus\\{0_{d}\\}$— to the adjusted Normal
distribution. To the best of our knowledge, our bound stands as the sharpest
among Berry Esseen bounds applicable to existing estimators in the assumption-
lean linear regression framework. Moreover, it accommodates the widest range
of dimensions for vanishing bound.

Under the moment condition of Assumption 2.a, our bound converges to zero if
$d=o(n^{\min\\{2/3,4/5-8/q_{x}\\}})$, disregarding the polylogarithmic
factors. This reduces to $d=o(n^{2/3})$ when $q_{x}\geq 12$. Furthermore, if
the covariates are sub-Gaussian, the Berry Esseen bound tends to zero if
$d=o(n^{2/3})$, again ignoring the polylogarithmic factors. Comparing the
Berry Essen bound of the bias-corrected estimator to the bound given in
Theorem 1 when the bias is known, we observe that using the bias estimate does
not impose a stricter dimension requirement. This is due to the fact that the
dimension requirement for the $\sqrt{n}$-consistency of the bias estimate in
Theorem 2 is less strict than that given in Theorem 1, as shown by the
inequality: ${\rm min}\\{2/3,4/5-8/(5q_{x})\\}\leq{\rm
min}\\{2/3,3/4-1/q_{x}\\}$ for $q_{x}\geq 4$. If $q_{x}\geq 12$, then the
quantities on both sides coincide with 2/3, and thus, the bound converges to
zero if $d=o(n^{2/3})$.

The following corollary presents a simpler Berry Essen bound for the bias-
corrected OLS estimator. We use $A\lesssim B$ as a shorthand for the
inequality $A\leq C_{n,d}B$, where $C_{n,d}$ involves only the constants
defined in the assumptions and log-polynomial factors of $n$ and $d$.

######  Corollary 1.

Suppose that the assumptions made in Theorem 3(i) hold. Furthermore, if
$q_{x}\geq 12$ and $d=O(n^{2/3})$, then

$\sup_{c\in\mathbb{R}^{d}\setminus\\{0_{d}\\}}\sup_{t\in\mathbb{R}}\left|\mathbb{P}\left[\sqrt{n}c^{\top}\left(\widehat{\beta}_{\rm
bc}-\beta\right)\leq
t\right]-\left\\{\Phi\left(\frac{t}{\sigma_{c}}\right)+\kappa_{c}\Phi^{(3)}\left(\frac{t}{\sigma_{c}}\right)\right\\}\right|\lesssim\frac{1}{\sqrt{n}}+\frac{d^{3/2}}{n}.$
(18)

In addition, the approximation to the Normal distribution follows:

$\sup_{c\in\mathbb{R}^{d}\setminus\\{0_{d}\\}}\sup_{t\in\mathbb{R}}\left|\mathbb{P}\left[\sqrt{n}c^{\top}\left(\widehat{\beta}_{\rm
bc}-\beta\right)\leq
t\right]-\Phi\left(\frac{t}{\sigma_{c}}\right)\right|\lesssim\sqrt{\frac{d}{n}}+\frac{d^{3/2}}{n}.$
(19)

A slower rate in (19) for the Normal approximation in comparison to that in
(18) stems from the incorporation of the adjustment term
$\sup_{c}\sup_{t}\kappa_{c}\Phi^{(3)}(t/\sigma_{c})$, which approximately
scales at $O(\sqrt{d/n})$ (See Lemma 16), into the bound. As discussed
previously in Remark 1, however, if the linear model is correctly specified so
that $\kappa_{c}=0$, then the bound provided in (18) can also be utilized for
the Normal approximation.

The Berry Esseen bound in Theorem 3, despite its reliance on the unknown
distribution parameter $\sigma_{c}^{2}$, does not hinder our ability to
conduct inference on the projection parameter $\beta$. We leverage recent
statistical methodologies based on re-sampling or sample splitting to enable
this inference, which remains unaffected by the uncertainty of
$\sigma_{c}^{2}$. In Section 4, a series of distinct methods are introduced,
illustrating and exemplifying this approach.

###  3.4 Consistency of the Sandwich Variance Estimator

Berry Esseen bound for the studentized bias-corrected LSE remains pertinent,
notably in the construction of Wald-type confidence intervals. As previously
discussed, accomplishing this entails devising an estimator for the asymptotic
variance that ideally maintains consistency across a wide range of dimensions
argued in Theorem 3. For a specified query vector
$c\in\mathbb{R}^{d}\setminus\\{0_{d}\\}$, the asymptotic variance of
$\sqrt{n}c^{\top}(\widehat{\beta}_{\rm bc}-\beta)$ can be approximated through
$\sigma_{c}^{2}=c^{\top}\Sigma^{-1}V\Sigma^{-1}c$, as asserted in Theorem 3.
Notably, the $d\times d$ matrix $\Sigma^{-1}V\Sigma^{-1}$ is commonly referred
to as the ‘sandwich variance’ for OLS in a misspecified linear regression
setting (White, 1980; Buja et al., 2019a). A natural estimator for the
sandwich variance is its plug-in counterpart
$\widehat{\Sigma}^{-1}\widehat{V}\widehat{\Sigma}^{-1}$, termed the ‘sandwich
variance estimator’, where $\widehat{\Sigma}$ is the sample Gram matrix and
$\widehat{V}$ is expressed as:

$\widehat{V}:=\frac{1}{n}\sum_{i=1}^{n}X_{i}X_{i}^{\top}(Y_{i}-X_{i}^{\top}\widehat{\beta})^{2}.$

To this end, we define
$\widehat{\sigma}_{c}^{2}:=c^{\top}\widehat{\Sigma}^{-1}\widehat{V}\widehat{\Sigma}^{-1}c$
for $c\in\mathbb{R}^{d}\setminus\\{0_{d}\\}$.

The consistency of the sandwich variance estimator has been studied under
various model- and distributional- assumptions and has often been a stand-
alone line of research. A recent significant finding in a closely related
context is detailed in Lemma 7 of Kuchibhotla et al. (2020). In particular,
they prove that
$\|\widehat{\Sigma}^{-1}\widehat{V}\widehat{\Sigma}^{-1}-\Sigma^{-1}V\Sigma^{-1}\|_{\rm
op}=o_{p}(1)$ under conditions where
$d=o(\min\\{n^{1/2},n^{2(1-2/q_{x})/3}\\})$, up to polylogarithmic factors,
provided that $(1/q_{x}+1/q)^{-1}\geq 4$. Hence, directly applying their
results to our scenario may be challenging, as we are considering situations
beyond $d=o(\sqrt{n})$. Nevertheless, we have found that along a specific
query vector $c\in\mathbb{R}^{d}$, the estimator $\widehat{\sigma}_{c}^{2}$
maintains consistency across a broader range of dimensions, as elucidated in
the following theorem.

######  Theorem 4.

Suppose that Assumption 1,2.a, and 3 holds. Fix
$c\in\mathbb{R}^{d}\setminus\\{0_{d}\\}$, and assume $d+4\log(2n)\leq
n/(18K_{x})^{2}$. If $(1/q_{x}+1/q)^{-1}\geq 4$ and $(3/q_{x}+1/q)^{-1}\geq
2$, then there exists a constant
$C=C(q,q_{x},K_{x},K_{y},\overline{\lambda},\underline{\lambda})$ such that

$\mathbb{P}\left(\left|\frac{\widehat{\sigma}_{c}^{2}}{\sigma_{c}^{2}}-1\right|\geq\eta_{n,d}\right)\leq\eta_{n,d},$

where

$\displaystyle\eta_{n,d}$ $\displaystyle=$ $\displaystyle
C\left(\frac{1}{n^{1/3}}+\sqrt{\frac{d+\log(2n)}{n}}+\frac{d\log(2n)}{n^{1-1/q_{x}}}\right)$
$\displaystyle+C\left[\left(\frac{d\log^{3/2}(2n)}{n^{1-2/q_{x}}}\right)^{\frac{1}{1+2/q_{x}}}+\left(\frac{d\log^{5/4}(2n)}{n^{1-1/q_{x}-1/q}}\right)^{\frac{2}{1+2/q_{x}+2/q}}+\left(\frac{d\log^{3/2}(2n)}{n^{1-2/q_{x}-2/(3q)}}\right)^{\frac{3}{1+3/q_{x}+1/q}}\right]$
$\displaystyle+C\left[\left(\frac{d\log^{5/4}(2n)}{n^{1-2/q_{x}}}\right)^{\frac{2}{1+4/q_{x}}}+\left(\frac{d\log^{3/4}(2n)}{n^{1-2/q_{x}}}\right)^{\frac{2+8/q_{x}+4/q}{1+6/q_{x}+2/q}}+\left(\frac{d\log^{1/2}(2n)}{n^{1-2/q_{x}-2/(3q)}}\right)^{\frac{3+12/q_{x}+6/q}{1+6/q_{x}+2/q}}\right].$

######  Remark 4.

The proof of Theorem 4 is extensive due to the complex nature of the sandwich
variance estimator, encompassed in Section A.7. Theorem 4 primarily focuses on
the concentration inequality for
$|\widehat{\sigma}_{c}^{2}/\sigma_{c}^{2}-1|$, with the pivotal deterministic
inequality established in Lemma 7. This lemma argues that the difference
between the plug-in sandwich variance estimate and its population counterparts
can be deterministically bounded using several factors, including
$\|\Sigma^{-1/2}\widehat{\Sigma}\Sigma^{-1/2}-I_{d}\|_{\rm op}$,
$\|\widehat{\beta}-\beta\|_{\Sigma}$, and $\max_{1\leq i\leq
n}|X_{i}^{\top}(\widehat{\beta}-\beta)|$, among others mentioned.
Consequently, the proof relies on the concentration inequality of these
quantities, which in themselves hold potential as areas of independent
interest. The rate presented in Theorem 4 is a result of leveraging these
concentration inequalities, aiming to achieve the widest possible scaling of
$d$ in terms of $n$. However, alternative choices remain viable.

Theorem 4 shows that for a given $c\in\mathbb{R}^{d}$, the variance estimate
$\widehat{\sigma}_{c}^{2}$ is consistent for the asymptotic variance
$\sigma_{c}^{2}$ as long as

$d=o\left(n^{\min\\{1-2/q_{x},1-1/q_{x}-1/q,1-2/q_{x}-2/(3q)\\}}\right),$ (20)

ignoring the polylogarithmic factors. Under the assumptions that
$(1/q_{x}+1/q)^{-1}\geq 4$ and $(3/q_{x}+1/q)^{-1}\geq 2$, it reduces to
$d=o(n^{2/3})$. This strictly embraces the one obtained in Theorem 3, but it
should be emphasized that Theorem 3 relies on the moment conditions $q_{x}\geq
8$ and $(1/q_{x}+1/q)^{-1}\geq 3$.

###  3.5 Berry Esseen Bound for the Studentized Bias-corrected Estimator

The error bounds established for the sandwich variance estimator in Theorem 4
and the Berry Esseen bound derived in Theorem 3 help to derive the Berry
Esseen bound for the studentized debiased LSE, leading to our main result
Theorem 5. For easier comprehension, we present a simpler bound in Corollary
2.

######  Theorem 5.

Suppose that Assumption 1,2.a, and 3 holds with $q_{x}\geq 8$,
$(1/q_{x}+1/q)^{-1}\geq 4$, and $(3/q_{x}+1/q)^{-1}\geq 2$. Assume
$d+4\log(2n)\leq n/(18K_{x})^{2}$, and let $\epsilon_{n,d}$ and $\eta_{n,d}$
be those defined in Theorem 3 and Theorem 4, respectively, with possibly
different choices of constants. Then,

$\displaystyle\sup_{c\in\mathbb{R}^{d}\setminus\\{0_{d}\\}}\sup_{t\in\mathbb{R}}\left|\mathbb{P}\left[\sqrt{n}\frac{c^{\top}(\widehat{\beta}_{\rm
bc}-\beta)}{\widehat{\sigma}_{c}}\leq
t\right]-\left\\{\Phi(t)+\kappa_{c}\Phi^{(3)}(t)\right\\}\right|\leq\epsilon_{n,d}+\eta_{n,d}.$

######  Corollary 2.

Suppose that the assumptions made in Theorem 5 hold. Furthermore, if
$q_{x}\geq 12$, $q\geq 6$, and $d=O(n^{2/3})$, then

$\sup_{c\in\mathbb{R}^{d}\setminus\\{0_{d}\\}}\sup_{t\in\mathbb{R}}\left|\mathbb{P}\left[\sqrt{n}\frac{c^{\top}(\widehat{\beta}_{\rm
bc}-\beta)}{\widehat{\sigma}_{c}}\leq
t\right]-\Phi(t)\right|\lesssim\frac{1}{n^{1/3}}+\sqrt{\frac{d}{n}}+\frac{d^{3/2}}{n}.$
(21)

Here, $\lesssim$ is a shorthand for the inequality up to constants and log-
polynomial factors of $n$ and $d$.

######  Remark 5.

The bound provided in (21) consists of three terms, none of which entirely
dominates the others. Specifically, when the dimension is relatively ‘small’,
that is, if $d\lesssim n^{1/3}$, the dominant factor in the bound becomes
$n^{-1/3}$. In a moderately dimensional scenario, i.e., when $n^{1/3}\lesssim
d\lesssim n^{1/2}$, $\sqrt{d/n}$ becomes the leading factor. However, in a
‘large’ dimension scenario where $d\gtrsim n^{1/2}$, the prevailing factor in
the bound is $d^{3/2}/n$.

One of the immediate consequences of Theorem 5 is the asymptotic validity of
the Wald confidence interval for the linear contrast of the projection
parameters. Define a level-$\alpha$ Wald confidence interval with respect to a
predetermined $c\in\mathbb{R}^{d}\setminus\\{0_{d}\\}$ as

${\rm CI}^{\rm Wald}_{c,\alpha}=\left[c^{\top}\widehat{\beta}_{\rm
bc}-\Phi^{-1}\left(1-\frac{\alpha}{2}\right)\frac{\widehat{\sigma}_{c}}{\sqrt{n}},c^{\top}\widehat{\beta}_{\rm
bc}-\Phi^{-1}\left(\frac{\alpha}{2}\right)\frac{\widehat{\sigma}_{c}}{\sqrt{n}}\right].$

######  Corollary 3.

Suppose that assumptions made in Theorem 6 hold, and let $\epsilon_{n,d}$ and
$\eta_{n,d}$ be as described therein. Then, the false coverage rate of Wald CI
can be controlled as

$\sup_{c\in\mathbb{R}^{d}\setminus\\{0_{d}\\}}\left|\mathbb{P}\left(c^{\top}\beta\notin{\rm
CI}^{\rm
Wald}_{c,\alpha}\right)-\alpha\right|\leq\epsilon_{n,d}+\eta_{n,d}+C\sqrt{\frac{d}{n}},$

for some $C=C(K_{x},K_{y},\underline{\lambda})$. If $q_{x}\geq 12$, $q\geq 6$,
and $d=O(n^{2/3})$, then

$\sup_{c\in\mathbb{R}^{d}\setminus\\{0_{d}\\}}\left|\mathbb{P}\left(c^{\top}\beta\notin{\rm
CI}^{\rm
Wald}_{c,\alpha}\right)-\alpha\right|\lesssim\frac{1}{n^{1/3}}+\sqrt{\frac{d}{n}}+\frac{d^{3/2}}{n}.$

##  4 Confidence Region for the Projection Parameter

In the previous section, we have shown that the bias-corrected least square
estimator is asymptotically Normal distributed. In this section, we present
three methods that do not require a consistent estimate of the asymptotic
variance but still enable us to get valid confidence intervals for the linear
contrasts of the projection parameters.

###  4.1 HulC

The HulC (Kuchibhotla et al., 2021) is a statistical method that computes
confidence regions for parameters by constructing convex hulls around a set of
estimates. Suppose that we split a sample of $n$ i.i.d observations
$(X_{1},Y_{1}),\ldots,(X_{n},Y_{n})$ into $B$ batches where each batch
contains at least $\lfloor n/B\rfloor$ observations***Kuchibhotla et al.
(2021) noted that the batches need not be equal-sized, but having
approximately equal sizes yields good width properties.. Denote the bias-
corrected least square estimators obtained from $b$:th batch as
$\widehat{\beta}_{\rm bc}^{(b)}$ for $b=1,\ldots,B$. For any
$c\in\mathbb{R}^{d}$, define the maximum median bias (Kuchibhotla et al.,
2021) for $\beta$ as

$\Delta_{c}:=\max_{1\leq b\leq
B}\left(\frac{1}{2}-\min\left\\{\mathbb{P}(c^{\top}(\widehat{\beta}_{\rm
bc}^{(b)}-\beta)>0),\mathbb{P}(c^{\top}(\widehat{\beta}_{\rm
bc}^{(b)}-\beta)<0)\right\\}\right).$

A key idea of the HulC is that the event

$\min_{1\leq b\leq B}c^{\top}\widehat{\beta}_{\rm bc}^{(b)}\leq
c^{\top}\beta\leq\max_{1\leq b\leq B}c^{\top}\widehat{\beta}_{\rm bc}^{(b)}$

occurs unless either of the following happens: (1)
$c^{\top}(\widehat{\beta}_{\rm bc}^{(b)}-\beta)>0$ for all $1\leq b\leq B$ or
(2) $c^{\top}(\widehat{\beta}_{\rm bc}^{(b)}-\beta)<0$ for all $1\leq b\leq
B$. Consequently, we get

$\mathbb{P}\left(\min_{1\leq b\leq B}c^{\top}\widehat{\beta}_{\rm
bc}^{(b)}\leq c^{\top}\beta\leq\max_{1\leq b\leq
B}c^{\top}\widehat{\beta}_{\rm bc}^{(b)}\right)\geq
1-\left(\frac{1}{2}-\Delta_{c}\right)^{B}-\left(\frac{1}{2}+\Delta_{c}\right)^{B}.$

Since the bias-corrected OLS estimator is asymptotically Normal as shown in
Corollory 1, $\Delta_{c}$ converges to the asymptotic median, which is $0$.
The detailed procedure for constructing confidence interval at level $\alpha$
is described in Algorithm 1.

######  Theorem 6.

Suppose that the assumptions made in Theorem 3(i) hold, and denote the Berry
Esseen bound given in Theorem 3(i) as $\epsilon_{n,d}$. For any
$\alpha\in(0,1)$ and $c\in\mathbb{R}^{d}$ with $\|c\|_{\Sigma^{-1}}=1$, let
$B=\lceil\log(2/\alpha)\rceil$ and assume $d\leq n/B$. Let ${\rm CI}^{\rm
HulC}_{c,\alpha}$ be the confidence interval returned by Algorithm 1. Then,
the false coverage rate of HulC CI can be controlled as

$\sup_{c\in\mathbb{R}^{d}\setminus\\{0_{d}\\}}\left\\{\mathbb{P}\left(c^{\top}\beta\notin{\rm
CI}^{\rm
HulC}_{c,\alpha}\right)-\alpha\right\\}\leq\alpha\left(B^{2}\epsilon_{\frac{n}{B},d}^{2}+C\frac{B^{3}d}{n}\right),$

for some constant $C=C(K_{x},K_{y},\underline{\lambda})$. If $q_{x}\geq 12$
and $d=O(n^{2/3})$, then

$\sup_{c\in\mathbb{R}^{d}\setminus\\{0_{d}\\}}\left\\{\mathbb{P}\left(c^{\top}\beta\notin{\rm
CI}^{\rm
HulC}_{c,\alpha}\right)-\alpha\right\\}\lesssim\alpha\left(\frac{B^{3}d}{n}+\frac{B^{4}d^{3}}{n^{2}}\right).$

######  Remark 6 (Comparison with Wald Confidence Interval).

The Wald CI and HulC CI are both (asymptotically) valid as demonstrated in
Corollary 3 and Theorem 6, respectively. A notable distinction between them
lies in the fact that the HulC CI does not require an estimate for the
asymptotic variance. Furthermore, the HulC CI shows a sharper miscoverage rate
due to its reliance on the ‘squared’ Berry Esseen bound in Theorem 3. Despite
this, the HulC CI tends to be wider due to its robust nature, which can be
empirically demonstrated in Section 5. Importantly, the expected width ratio
of the HulC CI to the Wald CI approximately equals
$\sqrt{\log_{2}(\log_{2}(2/\alpha))}$ (Kuchibhotla et al., 2021), slowly
diverging as $\alpha\to 0+$, and reaching 1.55 for the conventional choice of
$\alpha=0.05$. A more comprehensive comparison can be found in Section 2.3 of
Kuchibhotla et al. (2021).

Algorithm 1 HulC for the projection parameter

Require: data $(X_{1},Y_{1}),\ldots,(X_{n},Y_{n})$, target coverage
$1-\alpha$, and a query vector $c\in\mathbb{R}^{d}\setminus\\{0_{d}\\}$.
Ensure: A confidence interval ${\rm CI}^{\rm HulC}_{c,\alpha}$ for
$c^{\top}\beta$

1: Let $B=\lceil\log_{2}(2/\alpha)\rceil$ and $\tau=2^{B-1}\alpha-1$. Note
that $\tau\in[0,1)$.

2: Draw a random variable $U$ from the uniform distribution on the interval
$[0,1]$. Let $B^{*}=B-1$ if $U\leq\tau$ and let $B^{*}=B$ otherwise.

3: Split the data $(X_{1},Y_{1}),\ldots,(X_{n},Y_{n})$ into $B^{*}$ disjoint
batches where each batch contains at least $\lfloor n/B^{*}\rfloor$
observations.

4: Compute $B^{*}$ bias-corrected LS estimators $\widehat{\beta}_{\rm
bc}^{(1)},\ldots,\widehat{\beta}_{\rm bc}^{(B^{*})}$.

5: Return the level-$\alpha$ confidence interval  ${\rm CI}^{\rm
HulC}_{c,\alpha}=\left[\min_{1\leq b\leq B^{*}}c^{\top}\widehat{\beta}_{\rm
bc}^{(b)},\max_{1\leq b\leq B^{*}}c^{\top}\widehat{\beta}_{\rm
bc}^{(b)}\right].$

###  4.2 $t$-statistic based Inference

As in the previous section, consider $B$ partitions of $n$ observations. Here,
the number of batches $B$ does not necessarily depend on the target coverage
$1-\alpha$, and suppose that each batch has at least $\lfloor n/B\rfloor$
observations. Denote by $\widehat{\beta}_{\rm bc}^{(b)}$ the bias-corrected
estimator of $\beta$ using observations in batch $b$, $1\leq b\leq B$. Since
estimators are obtained independently, this amounts to the convergence of the
joint distribution;

$\sqrt{\frac{n}{B}}\left(\frac{c^{\top}(\widehat{\beta}_{\rm
bc}^{(1)}-\beta)}{\sigma_{c}},\ldots,\frac{c^{\top}(\widehat{\beta}_{\rm
bc}^{(B)}-\beta)}{\sigma_{c}}\right)^{\top}\xrightarrow{D}\mathcal{N}(0_{B},I_{B}),$

as $n\to\infty$ if the dimension requirement in Theorem 3 are met. Hence, we
can form an asymptotically valid confidence region based on the $t$-statistic,

$T=\sqrt{B}\frac{c^{\top}(\overline{\widehat{\beta}}_{\rm
bc}-\beta)}{s_{\widehat{\beta}}},$

where $\overline{\widehat{\beta}}_{\rm
bc}=B^{-1}\sum_{b=1}^{B}\widehat{\beta}_{\rm bc}^{(b)}$ and
$s_{\widehat{\beta}}^{2}=(B-1)^{-1}\sum_{b=1}^{B}\\{c^{\top}(\widehat{\beta}_{\rm
bc}^{(b)}-\overline{\widehat{\beta}}_{\rm bc})\\}^{2}$. By the continuous
mapping theorem (Billingsley, 2013, Theorem 2.7), the given statistic $T$ is
asymptotically $t$-distributed with $(B-1)$-degrees of freedom. By
implication, the confidence interval for $c^{\top}\beta$ is given by

${\rm CI}^{\rm T}_{c,\alpha}:=\left[c^{\top}\overline{\widehat{\beta}}_{\rm
bc}-t_{\alpha/2,B-1}\frac{s_{\widehat{\beta}}}{\sqrt{B}},c^{\top}\overline{\widehat{\beta}}_{\rm
bc}+t_{\alpha/2,B-1}\frac{s_{\widehat{\beta}}}{\sqrt{B}}\right],$

where $t_{q,\nu}$ denotes the upper $q$th quantile of $t$-distribution with
$\nu$ degrees of freedom.

###  4.3 Bootstrap Inference

This section presents a method for carrying out post-bias-correction inference
in our scenario, based on a bootstrap approach. Specifically, we utilize the
wild bootstrap technique (Wu, 1986) to derive a computationally efficient
approximation of the distribution of the bias-corrected OLS estimator. The
idea is, as the residual bootstrap, to leave the covariates at their sample
value, but to resample the response variable based on the residual values.
First, we compute from the original dataset that the OLS estimator
$\widehat{\beta}$ and the residuals
$\widehat{\epsilon}_{i}:=Y_{i}-X_{i}^{\top}\widehat{\beta}$ for
$i=1,\ldots,n$. For $i=1,\ldots,n$, let $\xi_{i}^{*}$ be i.i.d. bootstrap
weights such that $\mathbb{E}\xi_{i}^{*}=0$, $\mathbb{E}{\xi_{i}^{*}}^{2}=1$,
and $\mathbb{E}{\xi_{i}^{*}}^{3}=1$. Then, we generate a bootstrap sample of
$n$ observations $(X_{1},Y_{1}^{*}),\ldots,(X_{n},Y_{n}^{*})$ by letting

$Y_{i}^{*}=X_{i}^{\top}\widehat{\beta}+\widehat{\epsilon}_{i}\xi_{i}^{*},$

for $i=1,\ldots,n$. Although Bootstrap weights are often drawn from the
standard Normal distribution, we made use of the asymmetric weights in Section
5, which was proposed in Mammen (1993);
$\mathbb{P}(\xi^{*}=-(\sqrt{5}-1)/2)=(\sqrt{5}+1)/(2\sqrt{5})$ and
$\mathbb{P}(\xi^{*}=(\sqrt{5}+1)/2)=(\sqrt{5}-1)/(2\sqrt{5})$. Let estimator
$\widehat{\beta}^{*}_{\rm bc}$ be the bias-corrected least square estimator
based on bootstrap observations,

$\displaystyle\widehat{\beta}^{*}_{\rm
bc}=\widehat{\beta}^{*}-\widehat{\mathcal{B}}^{*},$

where $\widehat{\beta}^{*}$ and $\widehat{\mathcal{B}}^{*}$ are bootstrap
counterparts of the OLS estimator and estimated bias, respectively;

$\displaystyle\widehat{\beta}^{*}=\widehat{\Sigma}^{-1}\frac{1}{n}\sum_{i=1}^{n}X_{i}Y_{i}^{*}\quad\mbox{and}\quad\widehat{\mathcal{B}}^{*}=-\frac{1}{n^{2}}\sum_{i=1}^{n}\widehat{\Sigma}^{-1}X_{i}(Y_{i}^{*}-X_{i}^{\top}\widehat{\beta}^{*})\|X_{i}\|_{\widehat{\Sigma}^{-1}}^{2}.$

To construct the confidence interval, let $T_{b}^{*}=\widehat{\beta}_{\rm
bc}^{*,(b)}-\widehat{\beta}$ where $\widehat{\beta}_{\rm bc}^{*,(b)}$ is the
bootstrap bias-corrected estimator in $b$:th simulation. The wild bootstrap
level $\alpha$ bias-corrected confidence interval for $c^{\top}\beta$ is given
by

$\displaystyle{\rm CI}_{\alpha}^{\rm WBS}=\left[\widehat{\beta}_{\rm
bc}-\widehat{q}^{*}_{1-\alpha/2},\widehat{\beta}_{\rm
bc}-\widehat{q}^{*}_{\alpha/2}\right],$

where
$\widehat{q}^{*}_{\alpha}=\inf\\{t\in\mathbb{R}:\widehat{F}^{*}(t)\geq\alpha\\}$
is the empirical upper $\alpha$:th quantile of $\\{T_{b}^{*}:1\leq b\leq B\\}$
with $\widehat{F}^{*}(t)=\frac{1}{B}\sum_{b=1}^{B}\mathbbm{1}(T_{b}^{*}\leq
t)$.

######  Remark 7.

The validity of the above wild bootstrap inference may require additional
dimension requirements. Theorem 1 of Mammen (1993) shows that wild bootstrap
provides a consistent estimation of the unconditional distribution of
$\sqrt{n}(\widehat{\beta}-\beta)$ in a well-specified linear regression model,
without any bias correction, under the condition that $d=o(\sqrt{n})$. More
recently, Theorem 11 of Kuchibhotla et al. (2020) establishes a finite sample
approximation of the conditional bootstrap distribution of the law of OLSE in
a similar setting to ours. In particular, the proof utilizes the Gaussian
comparison bound (Lopes, 2020) to compare the conditional distribution of the
bootstrap estimate with the Normal distribution. However, the estimation of
asymptotic variance appears to be inevitable in the process, and this adds an
additional dimension requirement of $d=o(\sqrt{n})$, particularly when
$q_{x}\geq 8$ and $s\geq 4$ (see Lemma 10 of Kuchibhotla et al. (2020)).
Furthermore, although not discussed in a linear regression context, Lemma 4.1
of Zhilova (2020) gives a sequence of random variables that justifies the
necessity of the condition $d=o(\sqrt{n})$ for the consistency of the wild
bootstrap when estimating sample mean.

##  5 Numerical Studies

In this section, we compare the empirical coverages and widths of the 95%
confidence intervals obtained from three inferential methods discussed in
Section 4. Furthermore, we also compare the bootstrap confidence interval
based on the jackknife-debiased OLS proposed in Cattaneo et al. (2019). In our
specific simulation settings, which will be described shortly, we discovered
that when leveraged with the bootstrap inferential methods, the jackknife-
based debiased OLSE slightly outperforms the OLSE and the proposed bias-
corrected OLSE. A large set of implementations, including both resampling
bootstrap- and wild bootstrap-based confidence intervals using the OLS and the
proposed debiased OLS, are presented in Appendix D.

###  5.1 Well-specified Linear Model

This section concerns the well-specified linear model. The simulation setting
is as follows: for $n\in\left\\{1000,2000\right\\}$ and $d\in\left\\{20k:1\leq
k\leq 24\right\\}$, independent observations
$(Y_{i},X_{i})\in\mathbb{R}^{d}\times\mathbb{R}$, $1\leq i\leq n$, are
generated as

$\displaystyle
X_{i}\sim\mathcal{N}(0_{d},I_{d}),\epsilon_{i}\sim\mathcal{N}(0,1),\mbox{ and
}Y_{i}=2X_{i}(1)+\epsilon_{i},$

where $X_{i}(1)$ is the first coordinate of $X_{i}$, that is,
$X_{i}(1)=\mathbf{e}_{1}^{\top}X_{i}$. Thus, the response $Y_{i}$ relies only
on the first dimension of $X_{i}$ and independent error $\epsilon_{i}$. As the
given model is well-specified, the projection parameter is given by
$\beta=2\mathbf{e}_{1}\in\mathbb{R}^{d}$.

Figure 1 compares the empirical coverage and the width of the 95% confidence
intervals of the first coordinate of $\beta$, which is 2, obtained from
various methods. It is noteworthy that our bias-corrected estimator when
incorporated with HULC achieves the target coverage of 0.95. On the contrary,
the $t$-statistic-based confidence interval seems to be fairly conservative
while the wild bootstrap inference yields a seemingly less conservative
confidence interval.

Figure 1:  Comparison of coverages and widths of confidence intervals; wild
bootstrap based on Jackknife estimators, wild bootstrap based on proposed
bias-corrected estimators, HulC using proposed bias-corrected estimators, and
$t$-statistic based inference using proposed bias-corrected estimators, and
Wald confidence interval based on proposed estimators. The empirical coverages
and widths of CI are computed based on 1000 replications. To attain a
0.95-level HULC confidence interval, Algortihm 1 requires data to be split
into at most six subsets. This only allows the dimension to increase by
$\approx$320 in order to ensure that the least square estimator is well-
defined within a subset of data. Furthermore, $t$-statistics-based inference
also requires data splitting, and we used the same data split as employed in
HulC CI for the sake of simplicity.

###  5.2 Misspecified Model

In this section, we concern about a misspecified non-linear model. The data-
generating process for individual observation is as follows. We first
independently generate two $d$-dimensional Gaussian random vector

$\displaystyle Z\sim\mathcal{N}(0_{d},I_{d})\mbox{ and
}W\sim\mathcal{N}(0_{d},\Sigma_{\rho}).$

Here, the covariance matrix of $W$, $\Sigma_{\rho}$, is a compound symmetry
where the diagonal entries are all $1$ and the off-diagonal elements are all
$\rho$ for $\rho\in[0,1)$. That is, $\Sigma_{\rho}=(1-\rho)I_{d}+\rho
1_{d}1_{d}^{\top}$. We define the covariate vector as $X=Z\odot W$ where
$\odot$ denotes the entry-wise product. Then, the entries of $X$ are
uncorrelated, i.e., $\mathbb{E}[XX^{\top}]=I_{d}$, but not necessarily
independent except the case when $\rho=0$. For a $d$-dimensional parameter
$\theta\in\mathbb{R}^{d}$, we let

$Y=(X^{\top}\theta)^{3}+\epsilon,$ (22)

where $\epsilon\sim\mathcal{N}(0,1)$ and $\epsilon$ is independent with $Z$
and $W$.

Under the aforementioned data-generating process, we conducted comprehensive
experiments with the following sample sizes, dimensions, $\theta$, and $\rho$;

$\displaystyle n$ $\displaystyle\in$
$\displaystyle\left\\{1000,2000,5000,10000,20000\right\\},$ $\displaystyle d$
$\displaystyle\in$ $\displaystyle\left\\{50k:1\leq k\leq 20\right\\},$
$\displaystyle\theta$ $\displaystyle\in$
$\displaystyle\left\\{\mathbf{e}_{1},1_{d}/\sqrt{d}\right\\}\subset\mathbb{R}^{d},$
$\displaystyle\rho$ $\displaystyle\in$
$\displaystyle\left\\{0,0.2,0.5\right\\}.$

Under the model (22), the projection parameter $\beta$ has the closed-form
representation and turns out to be determined by $\theta$ and $\rho$ as,

$\displaystyle\beta=3(1+2\rho^{2})\|\theta\|^{2}_{2}\theta+6(1-\rho^{2})\theta^{\odot
3},$

where $\theta^{\odot 3}=\theta\odot\theta\odot\theta$ (see Lemma 28). The
linear contrast of interest was set to $\theta^{\top}\beta$ which becomes the
first coordinate of $\beta$ when $\theta=\mathbf{e}_{1}$ and becomes the
scaled average of coefficients of $\beta$ when $\theta=1_{d}/\sqrt{d}$.
Furthermore, we can show that

$\displaystyle\theta^{\top}\beta=\begin{cases}9,&\mbox{ if
}\theta=\mathbf{e}_{1},\\\ 3(1+2\rho^{2})+6(1-\rho^{2})/d,&\mbox{ if
}\theta=1_{d}^{\top}/\sqrt{d};\end{cases}.$

Figure 2 compares the coverage and length of the 95% confidence intervals for
$\theta^{\top}\beta$ under the setting where $n=20000$,
$\theta=\mathbf{e}_{1}$, and $\rho=0$ or $0.5$. Results under different
combinations of sample size, dimension, $\theta$, and $\rho$ are contained in
Appendix D.

Figure 2: Comparison of coverages and widths of confidence intervals obtained
from 5 different methods under the misspecified model with $n=20000$,
$\theta=\mathbf{e}_{1}$. The model parameter $\rho$ is set to $0.0$ (top) and
$0.5$ (bottom). See Figure 1 for the abbreviations of methods. The empirical
coverages and widths of CI are computed based on 1000. It is noteworthy that
$n^{1/2}\approx 141$ and $n^{2/3}\approx 737$.

##  6 Discussion

We provide an estimator for the bias of the ordinary least squares estimator
that is consistent as long as the dimension $d$ grows slower than $n^{2/3}$,
where $n$ is the sample size. The resulting debiased least squares estimator,
after proper normalization, is asymptotically Normal at an $n^{1/2}$-rate as
long as $d=o(n^{2/3})$. We also provide valid inference along any arbitrary
direction for the projection parameter without having to estimate the
variance. We achieve this inferential goal by leveraging the methods such as
HulC. We believe the results of this paper can be extended to further expand
the allowed growth rate of dimension. This extension would involve further
expansion of the least squares estimator into higher order $U$-statistics and
removing the bias.

##### Acknowledgements.

A. K. Kuchibhotla and A. Rinaldo were partially supported by NSF DMS-2113611.

## References

  * Abdalla (2023) Pedro Abdalla.  Covariance estimation under missing observations and $l\\_4-l\\_2$ moment equivalence.  _arXiv preprint arXiv:2305.12981_ , 2023.
  * Adamczak et al. (2010) Radosław Adamczak, Alexander Litvak, Alain Pajor, and Nicole Tomczak-Jaegermann.  Quantitative estimates of the convergence of the empirical covariance matrix in log-concave ensembles.  _Journal of the American Mathematical Society_ , 23(2):535–561, 2010.
  * Bentkus et al. (2009) Vidmantas Bentkus, Bing-Yi Jing, and Wang Zhou.  On normal approximations to U-statistics.  _The Annals of Probability_ , 37(6):2174 – 2199, 2009.
  * Berk et al. (2019) Richard Berk, Andreas Buja, Lawrence Brown, Edward George, Arun Kumar Kuchibhotla, Weijie Su, and Linda Zhao.  Assumption lean regression.  _The American Statistician_ , 2019.
  * Bickel and Freedman (1983) Peter J Bickel and David A Freedman.  Bootstrapping regression models with many parameters.  _Festschrift for Erich L. Lehmann_ , pages 28–48, 1983.
  * Billingsley (2013) Patrick Billingsley.  _Convergence of probability measures_.  John Wiley & Sons, 2013.
  * Brailovskaya and van Handel (2023) Tatiana Brailovskaya and Ramon van Handel.  Universality and sharp matrix concentration inequalities, 2023.
  * Buja et al. (2019a) Andreas Buja, Lawrence Brown, Richard Berk, Edward George, Emil Pitkin, Mikhail Traskin, Kai Zhang, and Linda Zhao.  Models as approximations I: Consequences illustrated with linear regression.  _Statistical Science_ , 34(4):523–544, 2019a.
  * Buja et al. (2019b) Andreas Buja, Lawrence Brown, Arun Kumar Kuchibhotla, Richard Berk, Edward George, and Linda Zhao.  Models as approximations II: A model-free theory of parametric regression.  _Statistical Science_ , 34(4):545–565, 2019b.
  * Catoni (2016) Olivier Catoni.  Pac-bayesian bounds for the gram matrix and least squares regression with a random design.  _arXiv preprint arXiv:1603.05229_ , 2016.
  * Cattaneo et al. (2015) Matias D Cattaneo, Michael Jansson, and Whitney K Newey.  Treatment effects with many covariates and heteroskedasticity.  Technical report, cemmap working paper, 2015.
  * Cattaneo et al. (2018a) Matias D Cattaneo, Michael Jansson, and Whitney K Newey.  Alternative asymptotics and the partially linear model with many regressors.  _Econometric Theory_ , 34(2):277–301, 2018a.
  * Cattaneo et al. (2018b) Matias D Cattaneo, Michael Jansson, and Whitney K Newey.  Inference in linear regression models with many covariates and heteroscedasticity.  _Journal of the American Statistical Association_ , 113(523):1350–1361, 2018b.
  * Cattaneo et al. (2019) Matias D Cattaneo, Michael Jansson, and Xinwei Ma.  Two-step estimation and inference with possibly many included covariates.  _The Review of Economic Studies_ , 86(3):1095–1122, 2019.
  * Donald and Newey (1994) S.G. Donald and W.K. Newey.  Series estimation of semilinear models.  _Journal of Multivariate Analysis_ , 50(1):30–40, 1994.  ISSN 0047-259X.  doi: https://doi.org/10.1006/jmva.1994.1032.  URL https://www.sciencedirect.com/science/article/pii/S0047259X84710323.
  * Einmahl and Li (2008) Uwe Einmahl and Deli Li.  Characterization of LIL behavior in Banach space.  _Transactions of the American Mathematical Society_ , 360(12):6677–6693, 2008.
  * Guédon and Rudelson (2007) Olivier Guédon and Mark Rudelson.  Lp-moments of random vectors via majorizing measures.  _Advances in Mathematics_ , 208(2):798–823, 2007\.
  * Guédon et al. (2017) Olivier Guédon, Alexander E Litvak, Alain Pajor, and Nicole Tomczak-Jaegermann.  On the interval of fluctuation of the singular values of random matrices.  _Journal of the European Mathematical Society (EMS Publishing)_ , 19(5), 2017.
  * Han and Wellner (2019) Qiyang Han and Jon A Wellner.  Convergence rates of least squares regression estimators with heavy-tailed errors.  _The Annals of Statistics_ , 47(4):2286–2319, 2019.
  * Ibragimov and Müller (2010) Rustam Ibragimov and Ulrich K Müller.  t-statistic based correlation and heterogeneity robust inference.  _Journal of Business & Economic Statistics_, 28(4):453–468, 2010.
  * Jochmans (2022) Koen Jochmans.  Heteroscedasticity-robust inference in linear regression models with many covariates.  _Journal of the American Statistical Association_ , 117(538):887–896, 2022.
  * Koltchinskii and Lounici (2017) Vladimir Koltchinskii and Karim Lounici.  Concentration inequalities and moment bounds for sample covariance operators.  _Bernoulli_ , 23(1):110–133, 2017.
  * Kuchibhotla and Patra (2022) Arun K Kuchibhotla and Rohit K Patra.  On least squares estimation under heteroscedastic and heavy-tailed errors.  _The Annals of Statistics_ , 50(1):277–302, 2022\.
  * Kuchibhotla et al. (2020) Arun Kumar Kuchibhotla, Alessandro Rinaldo, and Larry Wasserman.  Berry-Esseen bounds for projection parameters and partial correlations with increasing dimension.  _arXiv preprint arXiv:2007.09751_ , 2020.
  * Kuchibhotla et al. (2021) Arun Kumar Kuchibhotla, Sivaraman Balakrishnan, and Larry Wasserman.  The HulC: Confidence regions from convex hulls.  _arXiv preprint arXiv:2105.14577_ , 2021.
  * Lam (2022) Henry Lam.  A cheap bootstrap method for fast inference.  _arXiv preprint arXiv:2202.00090_ , 2022.
  * Long and Ervin (2000) J Scott Long and Laurie H Ervin.  Using heteroscedasticity consistent standard errors in the linear regression model.  _The American Statistician_ , 54(3):217–224, 2000\.
  * Lopes (2020) Miles E Lopes.  Central limit theorem and bootstrap approximation in high dimensions with near $1/\sqrt{n}$ rates.  _arXiv preprint arXiv:2009.06004_ , 2020.
  * MacKinnon (2012) James G MacKinnon.  Thirty years of heteroskedasticity-robust inference.  In _Recent advances and future directions in causality, prediction, and specification analysis: Essays in honor of Halbert L. White Jr_ , pages 437–461. Springer, 2012.
  * MacKinnon and White (1985) James G MacKinnon and Halbert White.  Some heteroskedasticity-consistent covariance matrix estimators with improved finite sample properties.  _Journal of econometrics_ , 29(3):305–325, 1985\.
  * Mammen (1993) Enno Mammen.  Bootstrap and wild bootstrap for high dimensional linear models.  _The annals of statistics_ , 21(1):255–285, 1993\.
  * Mendelson (2016) Shahar Mendelson.  Upper bounds on product and multiplier empirical processes.  _Stochastic Processes and their Applications_ , 126(12):3652–3680, 2016.
  * Mendelson and Paouris (2012) Shahar Mendelson and Grigoris Paouris.  On generic chaining and the smallest singular value of random matrices with heavy tails.  _Journal of Functional Analysis_ , 262(9):3775–3811, 2012.
  * Mendelson and Paouris (2014) Shahar Mendelson and Grigoris Paouris.  On the singular values of random matrices.  _Journal of the European Mathematical Society (EMS Publishing)_ , 16(4), 2014.
  * Mourtada (2022) Jaouad Mourtada.  Exact minimax risk for linear least squares, and the lower tail of sample covariance matrices.  _The Annals of Statistics_ , 50(4):2157–2178, 2022.
  * Mourtada et al. (2022) Jaouad Mourtada, Tomas Vaškevičius, and Nikita Zhivotovskiy.  Distribution-free robust linear regression.  _Mathematical Statistics and Learning_ , 4(3):253–292, 2022.
  * Oliveira (2016) Roberto Imbuzeiro Oliveira.  The lower tail of random quadratic forms with applications to ordinary least squares.  _Probability Theory and Related Fields_ , 166:1175–1194, 2016.
  * Pfanzagl (1973) Johann Pfanzagl.  The accuracy of the normal approximation for estimates of vector parameters.  _Zeitschrift für Wahrscheinlichkeitstheorie und verwandte Gebiete_ , 25(3):171–198, 1973.
  * Portnoy (1984) Stephen Portnoy.  Asymptotic behavior of $m$-estimators of $p$ regression parameters when $p^{2}/n$ is large. i. consistency.  _The Annals of Statistics_ , 12(4):1298–1309, 1984.
  * Portnoy (1985) Stephen Portnoy.  Asymptotic behavior of $m$ estimators of $p$ regression parameters when $p^{2}/n$ is large; ii. normal approximation.  _The Annals of Statistics_ , 13(4):1403–1417, 1985.
  * Portnoy (1986) Stephen Portnoy.  Asymptotic behavior of the empiric distribution of m-estimated residuals from a regression model with many parameters.  _The Annals of Statistics_ , 14(3):1152–1170, 1986.
  * Portnoy (1988) Stephen Portnoy.  Asymptotic behavior of likelihood methods for exponential families when the number of parameters tends to infinity.  _The Annals of Statistics_ , pages 356–366, 1988.
  * Puchkin et al. (2023) Nikita Puchkin, Fedor Noskov, and Vladimir Spokoiny.  Sharper dimension-free bounds on the frobenius distance between sample covariance and its expectation.  _arXiv preprint arXiv:2308.14739_ , 2023.
  * Rinaldo et al. (2019) Alessandro Rinaldo, Larry Wasserman, and Max G’Sell.  Bootstrapping and sample splitting for high-dimensional, assumption-lean inference.  _The Annals of Statistics_ , 47(6):3438–3469, 2019.
  * Rio (2017) Emmanuel Rio.  About the constants in the fuk-nagaev inequalities.  _HAL_ , 2017, 2017.
  * Rudelson (1999) Mark Rudelson.  Random vectors in the isotropic position.  _Journal of Functional Analysis_ , 164(1):60–72, 1999.
  * Srivastava and Vershynin (2013) Nikhil Srivastava and Roman Vershynin.  Covariance estimation for distributions with 2+ $\varepsilon$ moments.  _Annals of probability: An official journal of the Institute of Mathematical Statistics_ , 41(5):3081–3111, 2013.
  * Tikhomirov (2018) Konstantin Tikhomirov.  Sample covariance matrices of heavy-tailed distributions.  _International Mathematics Research Notices_ , 2018(20):6254–6289, 2018.
  * Tropp (2016) Joel A Tropp.  The expected norm of a sum of independent random matrices: An elementary approach.  In _High Dimensional Probability VII: The Cargèse Volume_ , pages 173–202. Springer, 2016.
  * Tropp et al. (2015) Joel A Tropp et al.  An introduction to matrix concentration inequalities.  _Foundations and Trends® in Machine Learning_ , 8(1-2):1–230, 2015.
  * Van der Vaart (2000) Aad W Van der Vaart.  _Asymptotic statistics_ , volume 3.  Cambridge university press, 2000.
  * Vansteelandt and Dukes (2022) Stijn Vansteelandt and Oliver Dukes.  Assumption-lean inference for generalised linear model parameters.  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 84(3):657–685, 2022.
  * Vaškevičius and Zhivotovskiy (2023) Tomas Vaškevičius and Nikita Zhivotovskiy.  Suboptimality of constrained least squares and improvements via non-linear predictors.  _Bernoulli_ , 29(1):473–495, 2023.
  * Vershynin (2012) Roman Vershynin.  How close is the sample covariance matrix to the actual covariance matrix?  _Journal of Theoretical Probability_ , 25(3):655–686, 2012.
  * Vershynin (2018) Roman Vershynin.  _High-dimensional probability: An introduction with applications in data science_ , volume 47.  Cambridge university press, 2018.
  * White (1980) Halbert White.  A heteroskedasticity-consistent covariance matrix estimator and a direct test for heteroskedasticity.  _Econometrica: journal of the Econometric Society_ , pages 817–838, 1980.
  * Wu (1986) Chien-Fu Jeff Wu.  Jackknife, bootstrap and other resampling methods in regression analysis.  _the Annals of Statistics_ , 14(4):1261–1295, 1986.
  * Yang and Kuchibhotla (2021) Yachong Yang and Arun Kumar Kuchibhotla.  Finite-sample efficient conformal prediction.  _arXiv preprint arXiv:2104.13871_ , 2021.
  * Zhilova (2020) M Zhilova.  Nonclassical Berry–Esseen inequalities and accuracy of the bootstrap.  _Annals of Statistics_ , 48(4):1922–1939, 2020\.

##  Appendix A Proofs of Theorems and Corollaries

###  A.1 Proof of Theorem 1

Recall that the OLS estimator can be expressed as
$\widehat{\beta}=\beta+\widehat{\Sigma}^{-1}\frac{1}{n}\sum_{i=1}^{n}X_{i}(Y_{i}-X_{i}^{\top}\beta)$.
We write Taylor series expansion (up to the second order) of the sample Gram
matrix $\widehat{\Sigma}$ at $\Sigma$ as

$\displaystyle\widehat{\Sigma}^{-1}-\Sigma^{-1}$ $\displaystyle=$
$\displaystyle\widehat{\Sigma}^{-1}(\Sigma-\widehat{\Sigma})\Sigma^{-1}$
$\displaystyle=$
$\displaystyle\Sigma^{-1}(\Sigma-\widehat{\Sigma})\Sigma^{-1}+(\widehat{\Sigma}^{-1}-\Sigma^{-1})(\Sigma-\widehat{\Sigma})\Sigma^{-1}$
$\displaystyle=$
$\displaystyle\Sigma^{-1}(\Sigma-\widehat{\Sigma})\Sigma^{-1}+\widehat{\Sigma}^{-1}(\Sigma-\widehat{\Sigma})\Sigma^{-1}(\Sigma-\widehat{\Sigma})\Sigma^{-1}.$

Plugging the expansion in the OLS expression gives that

$\displaystyle\widehat{\beta}-\beta$ $\displaystyle=$
$\displaystyle\Sigma^{-1}\frac{1}{n}\sum_{i=1}^{n}X_{i}(Y_{i}-X_{i}\beta)+\Sigma^{-1}(\Sigma-\widehat{\Sigma})\Sigma^{-1}\frac{1}{n}\sum_{i=1}^{n}X_{i}(Y_{i}-X_{i}^{\top}\beta)$
(24)
$\displaystyle+\,\widehat{\Sigma}^{-1}(\Sigma-\widehat{\Sigma})\Sigma^{-1}(\Sigma-\widehat{\Sigma})\Sigma^{-1}\frac{1}{n}\sum_{i=1}^{n}X_{i}(Y_{i}-X_{i}^{\top}\beta).$

As in (14), the low-order approximation (24) of $\widehat{\beta}-\beta$ can be
decomposed into two components: $\mathcal{B}$ and $\mathcal{U}$. In this
decomposition, $\mathcal{B}$ represents the bias, which is precisely defined
in (12), while $\mathcal{U}$ is a combination of first and second-order
$U$-statistics, expressed as:

$\displaystyle\mathcal{U}=\left(1+\frac{1}{n}\right)\sum_{i=1}^{n}\psi(X_{i},Y_{i})+\frac{1}{n(n-1)}\sum_{1\leq
i\neq j\leq n}\phi(X_{i},Y_{i},X_{j},Y_{j}),$

and $\psi$ and $\phi$ are defined in (13). Finally, we denote the
approximation error in (24) as $\mathcal{R}$.

The proof involves two steps of approximations; (1) an approximation of the
distribution of $\sqrt{n}c^{\top}(\widehat{\beta}-\beta-\mathcal{B})$ to that
of $\sqrt{n}\,c^{\top}\mathcal{U}$ and (2) an approximation of the
distribution of $\sqrt{n}\,c^{\top}\mathcal{U}$ to the adjusted Normal
distribution. For any $x\in\mathbb{R}$ and $\epsilon>0$,

$\displaystyle\mathbb{P}(\sqrt{n}c^{\top}(\widehat{\beta}-\beta-\mathcal{B})\leq
x)$ $\displaystyle=$
$\displaystyle\mathbb{P}(\sqrt{n}c^{\top}(\mathcal{U}+\mathcal{R})\leq x)$
(25) $\displaystyle\leq$
$\displaystyle\mathbb{P}(\sqrt{n}\,c^{\top}\mathcal{U}\leq
x+\epsilon)+\mathbb{P}(\left|\sqrt{n}c^{\top}\mathcal{R}\right|>\epsilon).$

The first term on the right-hand side can be approximated to the adjusted
Normal percentile by Lemma 16 as it says that

$\displaystyle\sup_{x\in\mathbb{R}}\left|\mathbb{P}(\sqrt{n}\,c^{\top}\mathcal{U}\leq
x+\epsilon)-N_{c}\left(\frac{x+\epsilon}{\sigma_{c}}\right)\right|\leq\Delta_{n},$
(26)

where we write $N_{c}(x)=\Phi(x)-\kappa_{c}\Phi^{(3)}(x)$ (see (14) for the
definition of $\kappa_{c}$). In addition, we have

$\Delta_{n}=C\left(\frac{\bar{\lambda}^{3/2}K_{x}^{3}K_{y}^{3}}{\sqrt{n}}+\frac{\bar{\lambda}K_{x}^{6}K_{y}^{2}d}{2n}\right)=O\left(\frac{1}{\sqrt{n}}+\frac{d}{n}\right),$

for some absolute constant $C>0$. Combining (25) and (26) leads to

$\displaystyle\mathbb{P}(\sqrt{n}c^{\top}(\widehat{\beta}-\beta-\mathcal{B})\leq
x)-N_{c}\left(\frac{x}{\sigma_{c}}\right)$ $\displaystyle\leq$ $\displaystyle
N_{c}\left(\frac{x+\epsilon}{\sigma_{c}}\right)-N_{c}\left(\frac{x}{\sigma_{c}}\right)+\Delta_{n}+\mathbb{P}(\left|\sqrt{n}c^{\top}\mathcal{R}\right|>\epsilon)$
(27) $\displaystyle\leq$
$\displaystyle\frac{\epsilon}{\sigma_{c}}\sup_{t\in\mathbb{R}}\left|N_{c}^{\prime}(t)\right|+\Delta_{n}+\mathbb{P}(\left|\sqrt{n}c^{\top}\mathcal{R}\right|>\epsilon).$

Instead of (25), we can similarly argue that

$\mathbb{P}(\sqrt{n}c^{\top}(\widehat{\beta}-\beta-\mathcal{B})>x)\leq\mathbb{P}(\sqrt{n}\,c^{\top}\mathcal{U}>x+\epsilon)+\mathbb{P}(\left|\sqrt{n}c^{\top}\mathcal{R}\right|>\epsilon),$

for any $x\in\mathbb{R}$ and $\epsilon>0$. This leads to an analog of (27),
and consequently, we have

$\left|\mathbb{P}(\sqrt{n}c^{\top}(\widehat{\beta}-\beta-\mathcal{B})\leq
x)-N_{c}\left(\frac{x}{\sigma_{c}}\right)\right|\leq\frac{\epsilon}{\sigma_{c}}\sup_{t\in\mathbb{R}}\left|N_{c}^{\prime}(t)\right|+\Delta_{n}+\mathbb{P}(\left|\sqrt{n}c^{\top}\mathcal{R}\right|>\epsilon).$

It follows from the definition of $N_{c}$ that
$\sup_{t\in\mathbb{R}}\left|N_{c}^{\prime}(t)\right|=(1+3\kappa_{c})/\sqrt{2\pi}$,
and $\kappa_{c}$ can be further bounded through Lemma 16 essentially leading
to

$\displaystyle\frac{1}{\sigma_{c}}\sup_{t\in\mathbb{R}}\left|N_{c}^{\prime}(t)\right|\leq\frac{1}{\sqrt{2\pi}\sigma_{c}}\left(1+\frac{3K_{x}^{3}K_{y}}{2\sigma_{c}}\sqrt{\frac{d}{n}}\right)$
$\displaystyle\leq$
$\displaystyle\frac{1}{\sqrt{2\pi}\sigma_{c}}\left(1+\frac{3K_{x}^{3}K_{y}}{2\sigma_{c}}\right)$
(28) $\displaystyle\leq$
$\displaystyle\sqrt{\frac{1}{2\pi\underline{\lambda}}}\left(1+\frac{3K_{x}^{3}K_{y}}{2\underline{\lambda}^{1/2}}\right).$

We used $d\leq n$ for the penultimate inequality, and the last inequality is
due to that $\sigma_{c}^{2}\geq\underline{\lambda}$ which follows from
Assumption 3. Hence, we have

$\left|\mathbb{P}(\sqrt{n}c^{\top}(\widehat{\beta}-\beta-\mathcal{B})\leq
x)-N_{c}\left(\frac{x}{\sigma_{c}}\right)\right|\leq\Delta_{n}+C\epsilon+\mathbb{P}\left(|\sqrt{n}c^{\top}\mathcal{R}|>\epsilon\right),$
(29)

for $C=C(\underline{\lambda},K_{x},K_{y})$ in (28). The next and final step
involves harnessing the value of $\epsilon$ on the right-hand side of (29),
which will serve as our Berry Esseen bound. It is worth highlighting that
Lemma 13—15 establishes the tail bound for the quantity
$\sqrt{n}c^{\top}\mathcal{R}$, contingent upon three different moment
assumptions on covariates (Assumption 2.a—2.c). In particular, their results
have the form of

$\mathbb{P}\left(|\sqrt{n}c^{\top}\mathcal{R}|>\epsilon_{1}(\delta)\right)\leq\delta.$

Hence, each of these scenarios necessitates distinct selections for $\delta$
which were chosen to minimize $\epsilon_{1}(\delta)+\delta$. We note that the
choices we make for $\delta$ are such that $1/\delta$ is no larger than a
polynomial of $n$, so that $\log(1/\delta)$ can be bounded with $\log(2n)$ up
to constants.

##### Under Assumption 2.a

We use Lemma 13. The choice of

$\displaystyle\delta=\left(\frac{d\log^{3}(2n)}{n^{4/5-8/(5q_{x})}}\right)^{5q_{x}/(2q_{x}+8)},$

yields the result.

##### Under Assumption 2.b

The result of Lemma 14 with $\delta=1/n$ yields the result.

##### Under Assumption 2.c

In Lemma 15, taking the right-hand side of (52) as $\epsilon=\epsilon(\delta)$
with the choice of

$\displaystyle\delta=\left(\frac{d^{1/2+4/q_{x}}\log^{3}(2d)}{n^{2-4/q_{x}}}\right)^{q_{x}/(q_{x}+4)}\leq\left(\frac{d^{3/2}\log^{3}(2d)}{n}\right)^{q_{x}/(q_{x}+4)},$

yields the result.

Moreover, in this scenario, we claimed that the bias $\mathcal{B}$ degenerates
with the additional condition that $(3/q_{x}+1/q)^{-1}\geq 2$. To see this, we
revisit to inequality in (25) and instead write,

$\displaystyle\mathbb{P}(\sqrt{n}c^{\top}(\widehat{\beta}-\beta)\leq x)$
$\displaystyle=$
$\displaystyle\mathbb{P}(\sqrt{n}c^{\top}(\mathcal{U}+\mathcal{B}+\mathcal{R})\leq
x)$ $\displaystyle\leq$
$\displaystyle\mathbb{P}(\sqrt{n}\,c^{\top}\mathcal{U}\leq
x+\epsilon+\epsilon^{\prime})+\mathbb{P}(\left|\sqrt{n}c^{\top}\mathcal{R}\right|>\epsilon)++\mathbb{P}(\left|\sqrt{n}c^{\top}\mathcal{B}\right|>\epsilon^{\prime}).$

This leads to

$\displaystyle\left|\mathbb{P}(\sqrt{n}c^{\top}(\widehat{\beta}-\beta-\mathcal{B})\leq
x)-N_{c}\left(\frac{x}{\sigma_{c}}\right)\right|$ $\displaystyle\leq$
$\displaystyle\Delta_{n}+C\epsilon+\mathbb{P}\left(|\sqrt{n}c^{\top}\mathcal{R}|>\epsilon\right)$
(31)
$\displaystyle+C\epsilon^{\prime}+\mathbb{P}\left(|\sqrt{n}c^{\top}\mathcal{B}|>\epsilon^{\prime}\right).$

It is noteworthy that (31) can still be controlled in the aforementioned way.
The quantity in (31) can be controlled using Lemma 17. Let
$\epsilon^{\prime}=\epsilon^{\prime}(\delta)$ be the right-hand side in the
probability notation in (58), then

$\displaystyle
C\epsilon^{\prime}+\mathbb{P}(|\sqrt{n}c^{\top}\mathcal{B}|>\epsilon^{\prime})\leq
CK_{x}^{3}K_{y}\sqrt{\frac{d}{n}}\left(1+\frac{d}{n\sqrt{\delta}}\right)+\delta.$

Finally, choosing $\delta=d/n$ yields the desired result.

###  A.2 Proof of Theorem 2

Throughout the proof, we fix a contrast vector $c\in\mathbb{R}^{d}$ such that
$\|c\|_{\Sigma^{-1}}=1$. To control
$|c^{\top}(\widehat{\mathcal{B}}-\mathcal{B})|$, we begin by writing,

$\displaystyle
c^{\top}\widehat{\Sigma}^{-1}X_{i}=c^{\top}\Sigma^{-1}X_{i}+R_{1,i},$
$\displaystyle R_{1,i}=c^{\top}(\widehat{\Sigma}^{-1}-\Sigma^{-1})X_{i},$
$\displaystyle(Y_{i}-X_{i}^{\top}\widehat{\beta})=(Y_{i}-X_{i}^{\top}\beta)+R_{2,i},$
$\displaystyle R_{2,i}=X_{i}^{\top}(\beta-\widehat{\beta}),$
$\displaystyle\|X_{i}\|_{\widehat{\Sigma}^{-1}}^{2}=\|X_{i}\|_{\Sigma^{-1}}^{2}(1+R_{3,i}),$
$\displaystyle
R_{3,i}=\frac{\|X_{i}\|_{\widehat{\Sigma}^{-1}}^{2}}{\|X_{i}\|_{\Sigma^{-1}}^{2}}-1,$

for $i=1,\ldots,n$. Then, it follows from the definition of
$\widehat{\mathcal{B}}$ in (17) that
$|c^{\top}(\widehat{\mathcal{B}}-\mathcal{B})|$ can be bounded by seven
distinct quantities, denoted as $\mathcal{S}_{k}$ for $k=1,\ldots,7$ with
respect to the presented order.

$\displaystyle\frac{n}{d}\left|c^{\top}(\widehat{\mathcal{B}}-\mathcal{B})\right|\leq$
$\displaystyle\left|\frac{1}{n}\sum_{i=1}^{n}R_{1,i}(Y_{i}-X_{i}^{\top}\beta)\frac{\|X_{i}\|_{\Sigma^{-1}}^{2}}{d}\right|+\left|\frac{1}{n}\sum_{i=1}^{n}c^{\top}\Sigma^{-1}X_{i}R_{2,i}\frac{\|X_{i}\|_{\Sigma^{-1}}^{2}}{d}\right|$
$\displaystyle+$
$\displaystyle\left|\frac{1}{n}\sum_{i=1}^{n}c^{\top}\Sigma^{-1}X_{i}(Y_{i}-X_{i}^{\top}\beta)R_{3,i}\frac{\|X_{i}\|_{\Sigma^{-1}}^{2}}{d}\right|+\left|\frac{1}{n}\sum_{i=1}^{n}R_{1,i}R_{2,i}\frac{\|X_{i}\|_{\Sigma^{-1}}^{2}}{d}\right|$
$\displaystyle+$
$\displaystyle\left|\frac{1}{n}\sum_{i=1}^{n}c^{\top}\Sigma^{-1}X_{i}R_{2,i}R_{3,i}\frac{\|X_{i}\|_{\Sigma^{-1}}^{2}}{d}\right|+\left|\frac{1}{n}\sum_{i=1}^{n}R_{1,i}(Y_{i}-X_{i}^{\top}\beta)R_{3,i}\frac{\|X_{i}\|_{\Sigma^{-1}}^{2}}{d}\right|$
$\displaystyle+$
$\displaystyle\left|\frac{1}{n}\sum_{i=1}^{n}R_{1,i}R_{2,i}R_{3,i}\frac{\|X_{i}\|_{\Sigma^{-1}}^{2}}{d}\right|$
$\displaystyle=:$ $\displaystyle\sum_{k=1}^{7}\mathcal{S}_{k}.$

Before analyzing individual remainder terms, we define

$\mathcal{D}_{\Sigma}=\left\|\Sigma^{-1/2}\widehat{\Sigma}\Sigma^{-1/2}-I_{d}\right\|_{\rm
op}.$ (32)

It is noteworthy that the quantity $R_{3,i}$ can be stochastically bounded by
in terms of $\mathcal{D}_{\Sigma}$. To see this, write the eigenvalue
decomposition of $\Sigma^{-1/2}\widehat{\Sigma}\Sigma^{-1/2}$ as
$\Sigma^{-1/2}\widehat{\Sigma}\Sigma^{-1/2}=U\Lambda U^{\top}$ where $U$ is a
$d\times d$ orthogonal matrix and $\Lambda$ is a $d\times d$ diagonal matrix
with decreasing diagonal entries. Then, the quantity $\mathcal{D}_{\Sigma}$ is
determined by two extreme eigenvalues as

$\mathcal{D}_{\Sigma}=\max\left\\{\left|\lambda_{\rm
max}(\Sigma^{-1/2}\widehat{\Sigma}\Sigma^{-1/2})-1\right|,\left|\lambda_{\rm
min}(\Sigma^{-1/2}\widehat{\Sigma}\Sigma^{-1/2})-1\right|\right\\}.$

Now, we note the deterministic relation that

$\displaystyle\sup_{\theta\in\mathbb{R}^{d}}\left|\frac{\theta^{\top}\widehat{\Sigma}^{-1}\theta}{\theta^{\top}\Sigma^{-1}\theta}-1\right|$
$\displaystyle=$
$\displaystyle\sup_{\omega\in\mathbb{R}^{d}}\left|\frac{\omega^{\top}(\Sigma^{1/2}\widehat{\Sigma}^{-1}\Sigma^{1/2}-I_{d})\omega}{\omega^{\top}\omega}\right|$
$\displaystyle=$
$\displaystyle\sup_{\omega\in\mathbb{R}^{d}}\left|\frac{\omega^{\top}(\Lambda^{-1}-I_{d})\omega}{\omega^{\top}\omega}\right|$
$\displaystyle=$ $\displaystyle\max\left\\{\left|\frac{1-\lambda_{\rm
max}(\Sigma^{-1/2}\widehat{\Sigma}\Sigma^{-1/2})}{\lambda_{\rm
max}(\Sigma^{-1/2}\widehat{\Sigma}\Sigma^{-1/2})}\right|,\left|\frac{1-\lambda_{\rm
min}(\Sigma^{-1/2}\widehat{\Sigma}\Sigma^{-1/2})}{\lambda_{\rm
min}(\Sigma^{-1/2}\widehat{\Sigma}\Sigma^{-1/2})}\right|\right\\}$
$\displaystyle\leq$ $\displaystyle\frac{\mathcal{D}_{\Sigma}}{\lambda_{\rm
min}(\Sigma^{-1/2}\widehat{\Sigma}\Sigma^{-1/2})}.$

Define the event $\mathcal{A}:=\left\\{\lambda_{\rm
min}(\Sigma^{-1/2}\widehat{\Sigma}\Sigma^{-1/2})\geq 1/2\right\\}$. On
$\mathcal{A}$, we have

$\displaystyle R_{3,i}\leq\|\Sigma^{1/2}\Sigma^{-1}\Sigma^{1/2}-I\|_{\rm
op}\leq 2\mathcal{D}_{\Sigma},\quad i=1,\ldots,n.$

Moreover, Proposition 24 states that

$\mathbb{P}\left(\lambda_{\rm
min}(\Sigma^{-1/2}\widehat{\Sigma}\Sigma^{-1/2})\geq
1-9K_{x}\sqrt{\frac{d+2\log(2/\delta)}{n}}\right)\leq 1-\delta,$ (33)

for any $\delta\in(0,1)$. Therefore, $\mathbb{P}(\mathcal{A})\leq 1-1/n$
provided that $d+2\log(2n)\leq n/(18K_{x})^{2}$.

By inspecting the relationship between the remainders, we can reduce the
number of remainders that we need to handle. On $\mathcal{A}$, we have

$\displaystyle\mathcal{S}_{5}\leq
2\mathcal{D}_{\Sigma}\mathcal{S}_{2},\quad\mathcal{S}_{6}\leq
2\mathcal{D}_{\Sigma}\mathcal{S}_{1},\quad\mathcal{S}_{7}\leq
2\mathcal{D}_{\Sigma}\mathcal{S}_{4}.$

Combining all, we have that with probability at least $1-1/n$ that

$\frac{n}{d}\left|c^{\top}(\widehat{\mathcal{B}}-\mathcal{B})\right|\leq\left(1+2\mathcal{D}_{\Sigma}\right)\left[\mathcal{S}_{1}+\mathcal{S}_{2}+\mathcal{S}_{4}\right]+\mathcal{S}_{3}.$
(34)

Hence, it suffices to shift our attention to the remainder $\mathcal{S}_{k}$
for $k=1,2,3,4$, and these four remainders have been carefully analyzed in
Lemmas 18—21, respectively. We provide a brief summary of the key results for
each remainder term: define the events with some proper constants
$\left\\{C_{i}:i=1,2,3,4\right\\}$ which only depend on $K_{x}$ and $K_{y}$
that

$\displaystyle\mathcal{M}_{1}$ $\displaystyle=$
$\displaystyle\left\\{\mathcal{S}_{1}\leq\frac{C_{1}\mathcal{D}_{\Sigma}}{\lambda_{\rm
min}(\Sigma^{-1/2}\widehat{\Sigma}\Sigma^{-1/2})}\right\\}\Rightarrow\mathcal{M}_{1}\cap\mathcal{A}\subseteq\left\\{\mathcal{S}_{1}\leq
2C_{1}\mathcal{D}_{\Sigma}\right\\},$ $\displaystyle\mathcal{M}_{2}$
$\displaystyle=$ $\displaystyle\left\\{\mathcal{S}_{2}\leq
C_{2}\|\widehat{\beta}-\beta\|_{\Sigma}\right\\},$
$\displaystyle\mathcal{M}_{3}$ $\displaystyle=$
$\displaystyle\left\\{\mathcal{S}_{3}\leq\frac{C_{3}\mathcal{D}_{\Sigma}}{\lambda_{\rm
min}(\Sigma^{-1/2}\widehat{\Sigma}\Sigma^{-1/2})}\right\\}\Rightarrow\mathcal{M}_{3}\cap\mathcal{A}\subseteq\left\\{\mathcal{S}_{3}\leq
2C_{3}\mathcal{D}_{\Sigma}\right\\}$ $\displaystyle\mathcal{M}_{4}(\delta)$
$\displaystyle=$ $\displaystyle\left\\{\mathcal{S}_{4}\leq
C_{4}\frac{\mathcal{D}_{\Sigma}\|\widehat{\beta}-\beta\|_{\Sigma}}{\lambda_{\rm
min}(\Sigma^{-1/2}\widehat{\Sigma}\Sigma^{-1/2})}\left(1+\frac{d\log(2d)}{\delta\sqrt{n}}\right)\right\\},\quad\delta\in(0,1]$
$\displaystyle\Rightarrow$
$\displaystyle\mathcal{M}_{4}(\delta)\cap\mathcal{A}\subseteq\left\\{\mathcal{S}_{4}\leq
2C_{4}\mathcal{D}_{\Sigma}\|\widehat{\beta}-\beta\|_{\Sigma}\left(1+\frac{d\log(2d)}{\delta\sqrt{n}}\right)\right\\}$

Lemma 18—21 prove that under Assumption 2.a the events $\mathcal{M}_{i}$,
$i=1,\ldots,4$ occur with high probabilities as,

$\mathbb{P}(\mathcal{M}_{1})\geq
1-\frac{d}{n},\quad\mathbb{P}(\mathcal{M}_{2})\geq
1-\frac{d}{n},\quad\mathbb{P}(\mathcal{M}_{3})\geq
1-\frac{1}{n},\quad\mathbb{P}(\mathcal{M}_{4})\geq 1-\delta.$

It follows from (34) and the definitions of the events $\mathcal{M}_{i}$,
$i=1,2,3,4$ that there exists a constant $C=C()$ such that on the event
$\mathcal{A}\cap\mathcal{M}_{1}\cap\mathcal{M}_{2}\cap\mathcal{M}_{3}\cap\mathcal{M}_{4}(\delta)$,
it holds that

$\frac{n}{d}|c^{\top}(\widehat{\mathcal{B}}-\mathcal{B})|\leq
C\left(1\vee\mathcal{D}_{\Sigma}\right)\left(\mathcal{D}_{\Sigma}+\|\widehat{\beta}-\beta\|_{\Sigma}+\mathcal{D}_{\Sigma}\|\widehat{\beta}-\beta\|_{\Sigma}\left(1+\frac{d\log(2d)}{\delta\sqrt{n}}\right)\right)$
(35)

Consequently, to manage $\lvert
c^{\top}(\widehat{\mathcal{B}}-\mathcal{B})\rvert$, it is sufficient to derive
upper bounds for $\mathcal{D}_{\Sigma}$ and
$\lVert\widehat{\beta}-\beta\rVert_{\Sigma}$ with high probabilities, which
are described in Propositions 22 and 27, respectively. We will separately
control these quantities based on the assumptions pertaining to the
covariates.

##### Under Assumption 2.a.

By Proposition 22 and Proposition 27, we get

$\displaystyle\mathbb{P}(\mathcal{D}_{\Sigma}\leq C\mathbf{A}_{n,d}(\eta))\geq
1-\eta,\quad\mathbb{P}(\|\widehat{\beta}-\beta\|_{\Sigma}\leq
C^{\prime}\mathbf{B}_{n,d})\geq 1-\frac{1}{n^{\min\\{1,s/2-1\\}}},$

where $C$ and $C^{\prime}$ are some constants which depend on
$q_{x},q,K_{x},K_{y},$ and $\overline{\lambda}$, and

$\mathbf{A}_{n,d}(\eta)=\frac{d\log^{3/2}(2d/\eta)}{\eta^{2/q_{x}}n^{1-2/q_{x}}}+\sqrt{\frac{d\log(2d/\eta)}{n}},\quad\mathbf{B}_{n,d}=\sqrt{\frac{d+\log(2n)}{n}}.$

Denote the events $\mathcal{C}_{1}(\eta):=\\{\mathcal{D}_{\Sigma}\leq
C\mathbf{A}_{n,d}(\eta)\\}$ and
$\mathcal{C}_{2}:=\\{\|\widehat{\beta}-\beta\|_{\Sigma}\leq
C^{\prime}\mathbf{B}_{n,d}\\}$, and we let

$\delta^{*}:=\delta^{*}(\eta)=\left[\left(\mathbf{A}_{n,d}^{1/2}(\eta)\vee\mathbf{A}_{n,d}(\eta)\right)\mathbf{B}_{n,d}^{1/2}\frac{\log^{1/2}(2d)d}{\sqrt{n}}\right]\wedge
1.$

Equipped with such choice of $\delta^{*}$, we have from (35) that

$\displaystyle\sqrt{n}|c^{\top}(\widehat{\mathcal{B}}-\mathcal{B})|$
$\displaystyle\leq$ $\displaystyle
C\left[\frac{d}{\sqrt{n}}(1\vee\mathbf{A}_{n,d}(\eta))\left(\mathbf{A}_{n,d}(\eta)+\mathbf{B}_{n,d}+\mathbf{A}_{n,d}(\eta)\mathbf{B}_{n,d}\right)+\frac{\log^{1/2}(2d)d}{\sqrt{n}}(\mathbf{A}_{n,d}^{1/2}(\eta)\vee\mathbf{A}_{n,d}(\eta))\mathbf{B}_{n,d}^{1/2}\right],$
$\displaystyle\leq$ $\displaystyle
C\left[\frac{d}{\sqrt{n}}(\mathbf{A}_{n,d}(\eta)\vee\mathbf{A}_{n,d}^{2}(\eta))+\frac{\log^{1/2}(2d)d}{\sqrt{n}}(\mathbf{A}_{n,d}^{1/2}(\eta)\vee\mathbf{A}_{n,d}(\eta))\mathbf{B}_{n,d}^{1/2}\right]$
$\displaystyle\leq$ $\displaystyle
C\frac{d}{\sqrt{n}}\left[(\mathbf{A}_{n,d}(\eta)\vee\mathbf{A}_{n,d}^{2}(\eta))+\mathbf{B}_{n,d}\log(2d)\right],$

holds with probability at least

$\mathbb{P}\left(\left(\mathcal{A}\cap\mathcal{M}_{1}\cap\mathcal{M}_{2}\cap\mathcal{M}_{3}\cap\mathcal{M}_{4}(\delta^{*})\cap\mathcal{C}_{1}(\eta)\cap\mathcal{C}_{2}\right)^{\complement}\right)\geq\frac{2d}{n}+\frac{1}{n}+\frac{1}{n^{\min\\{1,s/2-1\\}}}+\delta^{*}+\eta.$

##### Under Assumption 2.b

If the covariates are drawn from a sub-Gaussian distribution, a more robust
result regarding the concentration of the sample covariance matrix can be
obtained, as demonstrated in Proposition 22. Specifically, we can establish
that:

$\mathbb{P}(\mathcal{D}_{\Sigma}\leq C\mathbf{B}_{n,d})\geq 1-\frac{1}{n}.$

The remaining step of the proof is identical to the one presented under
Assumption 2.a, with the only difference being the utilization of
$\mathbf{B}_{n,d}$ instead of $\mathbf{A}_{n,d}$ as a stochastic bound for
$\mathcal{D}_{\Sigma}$.

###  A.3 Proof of Theorem 3

For a given $c\in\mathbb{R}^{d}$ with $\|c\|_{\Sigma^{-1}}=1$, write
$\widehat{S}_{n}=\sqrt{n}c^{\top}(\widehat{\beta}-\widehat{\mathcal{B}}-\beta)$
and $\tilde{S}_{n}=\sqrt{n}c^{\top}(\widehat{\beta}-\mathcal{B}-\beta)$. For
any $x\in\mathbb{R}$ and $\epsilon>0$, we have

$\displaystyle\mathbb{P}(\widehat{S}_{n}\leq x)$ $\displaystyle\leq$
$\displaystyle\mathbb{P}(\tilde{S}_{n}\leq
x+\epsilon)+\mathbb{P}(|\tilde{S}_{n}-\widehat{S}_{n}|\geq\epsilon),$
$\displaystyle\mathbb{P}(\widehat{S}_{n}>x)$ $\displaystyle\leq$
$\displaystyle\mathbb{P}(\tilde{S}_{n}>x+\epsilon)+\mathbb{P}(|\tilde{S}_{n}-\widehat{S}_{n}|\geq\epsilon).$
(36)

Define
$\tilde{\Delta}_{n,d}:=\sup_{x\in\mathbb{R}}|\mathbb{P}(\tilde{S}_{n}\leq
x)-N_{c}(x/\sigma_{c})|$ and
$\widehat{\Delta}_{n,d}:=\sup_{x\in\mathbb{R}}|\mathbb{P}(\widehat{S}_{n}\leq
x)-N_{c}(x/\sigma_{c})|$. Then, it follows from (A.3) that

$\displaystyle\widehat{\Delta}_{n}$ $\displaystyle\leq$
$\displaystyle\tilde{\Delta}_{n}+\sigma_{c}^{-1}\|N_{c}^{\prime}\|_{\infty}\epsilon+\mathbb{P}\left(\sqrt{n}|\widehat{\mathcal{B}}-\mathcal{B}|\geq\epsilon\right)$
(37) $\displaystyle\leq$
$\displaystyle\tilde{\Delta}_{n}+C\epsilon+\mathbb{P}\left(\sqrt{n}|\widehat{\mathcal{B}}-\mathcal{B}|\geq\epsilon\right),$

where $C=C(K_{x},K_{y},\underline{\lambda})$ is defined in (28). Note that the
quantity $\tilde{\Delta}_{n}$ is bounded in Theorem 1. Moreover, the tail
probability of $\widehat{\mathcal{B}}$ is controlled in Theorem 2.

##### Under Assumption 2.a

We have from Theorem 2(i) that

$\mathbb{P}\left(\sqrt{n}|\widehat{\mathcal{B}}-\mathcal{B}|\geq\delta_{n,d}(\eta)\right)\leq\frac{1}{\sqrt{n}}+\delta_{n,d}(\eta).$

Taking $\epsilon=\delta_{n,d}(\eta^{*})$ in (37) where

$\eta^{*}=\left(\frac{d\log^{3/4}(2n)}{n^{3/4-1/q_{x}}}\right)^{(2q_{x})/(2+q_{x})}\vee\left(\frac{d\log(2n)}{n^{5/6-4/(3q_{x})}}\right)^{(3q_{x})/(4+q_{x})},$

yields the result.

##### Under Assumption 2.b

Taking $\epsilon=\delta_{n,d}$, which is defined in Theorem 2(ii), gives the
desired result.

###  A.4 Proof of Corollary 1

We will show that Berry Esseen bound in Theorem 3(i) is dominated by
$1/\sqrt{n}+d^{3/2}/n$ up to polylogarithmic factors. Since $q_{x}\geq 12$ and
$d\lesssim n^{2/3}$, we observe that

$\frac{d\log^{3}(2n)}{n^{4/5-8/(5q_{x})}},\frac{d\log^{3/4}(2n)}{n^{3/4-1/q_{x}}},\frac{d\log(2n)}{n^{5/6-4/(3q_{x})}}\lesssim\frac{d}{n^{2/3}}\lesssim
1.$

Furthermore, since the exponents of the above terms appeared in Berry Esseen
bound is larger than $3/2$ as long as $q_{x}\geq 12$, the proof is completed.

###  A.5 Proof of Theorem 5

Throughout, we fix $c\in\mathbb{R}^{d}$ with $\|c\|_{\Sigma^{-1}}=1$. Theorem
3 proves that

$\sup_{t\in\mathbb{R}}\left|\mathbb{P}\left(\sqrt{n}\frac{c^{\top}(\widehat{\beta}_{\rm
bc}-\beta)}{\sigma_{c}}\leq t\right)-N_{c}(t)\right|\leq\epsilon_{n,d}.$

We observe that for a small $\delta>0$,

$\mathbb{P}\left(\sqrt{n}\frac{c^{\top}(\widehat{\beta}_{\rm
bc}-\beta)}{\widehat{\sigma}_{c}}\leq
t\right)\leq\mathbb{P}\left(\sqrt{n}\frac{c^{\top}(\widehat{\beta}_{\rm
bc}-\beta)}{\sigma_{c}}\leq
t(1+\delta)\right)+\mathbb{P}\left(\frac{\widehat{\sigma}_{c}}{\sigma_{c}}>1+\delta\right).$

Consequently,

$\displaystyle\mathbb{P}\left(\sqrt{n}\frac{c^{\top}(\widehat{\beta}_{\rm
bc}-\beta)}{\widehat{\sigma}_{c}}\leq t\right)-N_{c}(t)$ $\displaystyle\leq$
$\displaystyle\mathbb{P}\left(\sqrt{n}\frac{c^{\top}(\widehat{\beta}_{\rm
bc}-\beta)}{\sigma_{c}}\leq
t(1+\delta)\right)-N_{c}(t)+\mathbb{P}\left(\frac{\widehat{\sigma}_{c}}{\sigma_{c}}>1+\delta\right)$
$\displaystyle\leq$
$\displaystyle\epsilon_{n,d}+N_{c}\left((1+\delta)t\right)-N_{c}(t)+\mathbb{P}\left(\frac{\widehat{\sigma}_{c}}{\sigma_{c}}>1+\delta\right).$

Similarly, we can obtain

$\displaystyle\mathbb{P}\left(\sqrt{n}\frac{c^{\top}(\widehat{\beta}_{\rm
bc}-\beta)}{\widehat{\sigma}_{c}}>t\right)-(1-N_{c}(t))$ $\displaystyle\leq$
$\displaystyle\mathbb{P}\left(\sqrt{n}\frac{c^{\top}(\widehat{\beta}_{\rm
bc}-\beta)}{\sigma_{c}}>t(1-\delta)\right)-(1-N_{c}(t))+\mathbb{P}\left(\frac{\widehat{\sigma}_{c}}{\sigma_{c}}\leq
1-\delta\right)$ $\displaystyle\leq$
$\displaystyle\epsilon_{n,d}+N_{c}(t)-N_{c}\left((1-\delta)t\right)+\mathbb{P}\left(\frac{\widehat{\sigma}_{c}}{\sigma_{c}}\leq
1-\delta\right).$

Combining these leads to

$\left|\mathbb{P}\left(\sqrt{n}\frac{c^{\top}(\widehat{\beta}_{\rm
bc}-\beta)}{\widehat{\sigma}_{c}}\leq
t\right)-N_{c}(t)\right|\leq\epsilon_{n,d}+N_{c}((1+\delta)t)-N_{c}((1-\delta)t)+\mathbb{P}\left(\left|\frac{\widehat{\sigma}_{c}}{\sigma_{c}}-1\right|>\delta\right).$
(38)

To bound the middle term on the right-hand side in (38), we note that for some
constant $c>0$,

$\displaystyle\left|N_{c}((1+\delta)t)-N_{c}((1-\delta)t)\right|$
$\displaystyle=$
$\displaystyle\left|\int_{1-\delta}^{1+\delta}N_{c}^{\prime}(ut)t\,du\right|$
$\displaystyle\leq$
$\displaystyle\int_{1-\delta}^{1+\delta}\left|N_{c}^{\prime}(ut)t\,du\right|\leq\frac{1}{1-\delta}\int_{1-\delta}^{1+\delta}\left|N_{c}^{\prime}(ut)ut\,du\right|\leq
c\frac{2\delta}{1-\delta}.$

For the last inequality, we use the fact that
$xN_{c}^{\prime}(x)=x[\Phi^{(1)}(x)-\kappa_{c}\Phi^{(4)}(x)]$ is uniformly
bounded (see Lemma 16). For the rightmost term in (38), define the event

$\mathcal{E}_{c,n,d}=\left\\{\left|\frac{\widehat{\sigma}_{c}^{2}}{\sigma_{c}^{2}}-1\right|\leq\eta_{n,d}\right\\},$

where $\eta_{n,d}$ is defined in Theorem 4 and satisfying that
$\mathbb{P}(\mathcal{E}_{c,n,d})\geq 1-\eta_{n,d}$. On $\mathcal{E}_{c,n,d}$,
we have

$1-(\eta_{n,d}\wedge 1)\leq\sqrt{1-(\eta_{n,d}\wedge
1)}\leq\frac{\widehat{\sigma}_{c}}{\sigma_{c}}\leq\sqrt{1+\eta_{n,d}}\leq
1+\frac{1}{2}\eta_{n,d}.$

Hence, taking $\delta=\eta_{n,d}\wedge 1/2$ in (38) yields

$\left|\mathbb{P}\left(\sqrt{n}\frac{c^{\top}(\widehat{\beta}_{\rm
bc}-\beta)}{\widehat{\sigma}_{c}}\leq
t\right)-N_{c}(t)\right|\leq\epsilon_{n,d}+C\eta_{n,d}.$

###  A.6 Proof of Theorem 6

For $c\in\mathbb{R}^{d}\setminus\\{0_{d}\\}$, Theorem 3 proves that

$\displaystyle\sup_{t\in\mathbb{R}}\left|\mathbb{P}\left[\sqrt{n}c^{\top}(\widehat{\beta}_{\rm
bc}-\beta)\leq
t)\right]-N_{c}\left(\frac{t}{\sigma_{c}}\right)\right|\leq\epsilon_{n,d},$

for some rate $\epsilon_{n,d}>0$. Hence, the median bias of
$c^{\top}(\widehat{\beta}_{\rm bc}-\beta)$ is bounded as

$\displaystyle{\rm Medbias}\left(c^{\top}(\widehat{\beta}_{\rm
bc}-\beta)\right)$ $\displaystyle:=$
$\displaystyle\frac{1}{2}-\min\left\\{\mathbb{P}(c^{\top}(\widehat{\beta}_{\rm
bc}^{(b)}-\beta)>0),\mathbb{P}(c^{\top}(\widehat{\beta}_{\rm
bc}^{(b)}-\beta)<0)\right\\}$ (39) $\displaystyle\leq$
$\displaystyle\epsilon_{n,d}+\left|N_{c}(0)-\frac{1}{2}\right|=\epsilon_{n,d}+\kappa_{c}$
$\displaystyle\leq$ $\displaystyle\epsilon_{n,d}+C\sqrt{\frac{d}{n}}.$

The last inequality follows from Lemma 17 and the constant $C$ is described
therein. Note that $\widehat{\beta}_{\rm bc}^{(1)},\ldots,\widehat{\beta}_{\rm
bc}^{(B)}$ are independent bias-corrected estimators based on $\lfloor
n/B\rfloor$. For any $1\leq b\leq B$, we note from (39) that

$\displaystyle{\rm Medbias}\left(c^{\top}(\widehat{\beta}_{\rm
bc}^{(b)}-\beta)\right)\leq\epsilon_{\frac{n}{B},d}+C\sqrt{\frac{Bd}{n}}.$

An application of Theorem 2 of Kuchibhotla et al. [2021] results in

$\displaystyle\mathbb{P}\left(\beta\notin{\rm CI}_{\alpha}^{\rm
HULC}\right)-\alpha$ $\displaystyle\leq$ $\displaystyle\frac{\alpha
B(B-1)}{2}\left(\epsilon_{\frac{n}{B},d}+C\sqrt{\frac{Bd}{n}}\right)^{2}$
$\displaystyle\leq$ $\displaystyle\alpha
B^{2}\epsilon_{\frac{n}{B},d}^{2}+C^{2}\frac{\alpha B^{3}d}{n}.$

This completes the proof.

###  A.7 Consistency of the Sandwich Variance Estimator

####  A.7.1 Proof of Theorem 4

In this section, we collect various bounds that are used in the proof of
Theorem 4 about the consistency of the sandwich variance estimator with
respect to a given query vector $c\in\mathbb{R}^{d}\setminus\\{0_{d}\\}$. Note
that $|\widehat{\sigma}_{c}^{2}/\sigma_{c}^{2}-1|$ is invariant to the scaling
of $c$. Hence, we let $c$ be such that $\|c\|_{\Sigma}=1$, and write the
$L_{2}$-normalized version as $\tilde{c}$, so that $c=\Sigma^{1/2}\tilde{c}$
and $\|\tilde{c}\|_{2}=1$. Then, the sandwich variance estimate
$\widehat{\sigma}_{c}^{2}$ can be expressed as

$\widehat{\sigma}_{c}^{2}=\frac{1}{n}\sum_{i=1}^{n}(\tilde{c}^{\top}\Sigma^{1/2}\widehat{\Sigma}^{-1}X_{i})^{2}(Y_{i}-X_{i}^{\top}\widehat{\beta})^{2}.$

We commence by introducing an an intermediary quantity
$\tilde{\sigma}_{c}^{2}$, defined as

$\tilde{\sigma}_{c}^{2}=\frac{1}{n}\sum_{i=1}^{n}(\tilde{c}^{\top}\Sigma^{-1/2}X_{i})^{2}(Y_{i}-X_{i}^{\top}\beta)^{2}.$
(40)

It is immediate from its definition that
$\mathbb{E}[\tilde{\sigma}_{c}^{2}]=\sigma_{c}^{2}$. Moreover, we have from
Assumption 3 that

$\left|\frac{\widehat{\sigma}_{c}^{2}}{\sigma_{c}^{2}}-1\right|\leq\underline{\lambda}^{-2}\left|\widehat{\sigma}_{c}^{2}-\sigma_{c}^{2}\right|\leq\underline{\lambda}^{-2}\left(\left|\widehat{\sigma}_{c}^{2}-\tilde{\sigma}_{c}^{2}\right|+\left|\tilde{\sigma}_{c}^{2}-\sigma_{c}^{2}\right|\right).$

Since $(1/q_{x}+1/q)^{-1}\geq 4$, the quantity
$\left|\tilde{\sigma}_{c}^{2}-\sigma_{c}^{2}\right|$ can be controlled via
Chebyshev’s inequality as:

$\displaystyle\mathbb{P}\left(\left|\tilde{\sigma}_{c}^{2}-\sigma_{c}^{2}\right|\geq
t\right)$ $\displaystyle\leq$ $\displaystyle t^{-2}{\rm
Var}\left(\frac{1}{n}\sum_{i=1}^{n}(\tilde{c}^{\top}\Sigma^{-1/2}X_{i})^{2}(Y_{i}-X_{i}^{\top}\beta)^{2}\right)=(nt^{2})^{-1}{\rm
Var}\left((\tilde{c}^{\top}\Sigma^{-1/2}X)^{2}(Y-X^{\top}\beta)^{2}\right)$
$\displaystyle\overset{(i)}{\leq}$
$\displaystyle(nt^{2})^{-1}\mathbb{E}\left[(\tilde{c}^{\top}\Sigma^{-1/2}X)^{4}(Y-X^{\top}\beta)^{4}\right]$
$\displaystyle\overset{(ii)}{\leq}$
$\displaystyle(nt^{2})^{-1}\left(\mathbb{E}\left[(\tilde{c}^{\top}\Sigma^{-1/2}X)^{q_{x}}\right]\right)^{4/q_{x}}\left(\mathbb{E}\left[(Y-X^{\top}\beta)^{q}\right]\right)^{4/q}\overset{(iii)}{\leq}\frac{K_{x}^{4}K_{y}^{4}}{nt^{2}}.$

Here, $(i)$ is Jensen’s inequality and the inequality, and $(ii)$ follows from
$(1/q_{x}+1/q)^{-1}\geq 4$ and Hölder’s inequality. The last inequality
$(iii)$ is due to Assumption 1 and 2.a. The choice of
$t=K_{x}^{2}K_{y}^{2}n^{-1/3}$ yields

$\mathbb{P}\left(\left|\tilde{\sigma}_{c}^{2}-\sigma_{c}^{2}\right|\geq
K_{x}^{2}K_{y}^{2}n^{-1/3}\right)\leq n^{-1/3}.$

Our next lemma presents a deterministic bound for
$|\widehat{\sigma}_{c}^{2}-\tilde{\sigma}_{c}^{2}|$, implying that the rate of
convergence depends on several quantities, which will be described therein.

######  Lemma 7.

[Deterministic Bound for the Sandwich Estimator] For a given
$\tilde{c}\in\mathbb{S}^{d-1}$, define

$\displaystyle\mathcal{M}$ $\displaystyle=$
$\displaystyle\sup_{\theta\in\mathbb{S}^{d-1}}\frac{1}{n}\sum_{i=1}^{n}(\theta^{\top}\Sigma^{-1/2}X_{i})^{2}(\tilde{c}^{\top}\Sigma^{-1/2}X_{i})^{2},\quad\mathcal{N}=\sup_{\theta\in\mathbb{S}^{d-1}}\frac{1}{n}\sum_{i=1}^{n}(\theta^{\top}\Sigma^{-1/2}X_{i})^{2}(Y_{i}-X_{i}^{\top}\beta)^{2},$
$\displaystyle\mathcal{L}$ $\displaystyle=$
$\displaystyle\max\left\\{\left\|\frac{1}{n}\sum_{i=1}^{n}\Sigma^{-1/2}X_{i}(\tilde{c}^{\top}\Sigma^{-1/2}X_{i})^{2}(Y_{i}-X_{i}^{\top}\beta)\right\|_{2},\left\|\frac{1}{n}\sum_{i=1}^{n}\Sigma^{-1/2}X_{i}(\tilde{c}^{\top}\Sigma^{-1/2}X_{i})(Y_{i}-X_{i}^{\top}\beta)^{2}\right\|_{2}\right\\}.$

If $\\{\lambda_{\rm min}(\Sigma^{-1/2}\widehat{\Sigma}\Sigma^{-1/2})\geq
1/2\\}$ holds true, then

$\displaystyle\left|\widehat{\sigma}_{c}^{2}-\tilde{\sigma}_{c}^{2}\right|$
$\displaystyle\leq$ $\displaystyle
2\mathcal{L}\|\widehat{\beta}-\beta\|_{\Sigma}+3\mathcal{M}\|\widehat{\beta}-\beta\|_{\Sigma}^{2}+\left(8\mathcal{L}+12\max_{1\leq
i\leq
n}|X_{i}^{\top}(\widehat{\beta}-\beta)|^{2}\right)\mathcal{D}_{\Sigma}+8\mathcal{N}\mathcal{D}_{\Sigma}^{2}.$

###### proof of Lemma 7.

We begin by writing

$\displaystyle\widehat{\Sigma}^{-1}\Sigma^{1/2}\tilde{c}=\Sigma^{-1/2}(\tilde{c}+r_{1}),$
$\displaystyle
r_{1}=(\Sigma^{1/2}\widehat{\Sigma}^{-1}\Sigma^{1/2}-I_{d})\tilde{c},$
$\displaystyle
Y_{i}-X_{i}^{\top}\widehat{\beta}=Y_{i}-X_{i}^{\top}\beta+r_{2,i},$
$\displaystyle r_{2,i}=X_{i}^{\top}(\widehat{\beta}-\beta)\quad\mbox{for}\quad
i=1,\ldots,n.$

From their definitions, the difference
$|\widehat{\sigma}_{c}^{2}-\tilde{\sigma}_{c}^{2}|$ can be bounded by the sum
of eight quantities, denoted as $\mathcal{Q}_{k}$ for $k=1,\ldots,8$ in the
presented order.

$\displaystyle|\widehat{\sigma}_{c}^{2}-\tilde{\sigma}_{c}^{2}|\leq\left|\frac{2}{n}\sum_{i=1}^{n}|\tilde{c}^{\top}\Sigma^{-1/2}X_{i}|^{2}(Y_{i}-X_{i}^{\top}\beta)r_{2,i}\right|+\left|\frac{1}{n}\sum_{i=1}^{n}|\tilde{c}^{\top}\Sigma^{-1/2}X_{i}|^{2}r_{2,i}^{2}\right|$
$\displaystyle+\left|\frac{2}{n}\sum_{i=1}^{n}(\tilde{c}^{\top}\Sigma^{-1/2}X_{i})(r_{1}^{\top}\Sigma^{-1/2}X_{i})(Y_{i}-X_{i}^{\top}\beta)^{2}\right|+\left|\frac{4}{n}\sum_{i=1}^{n}(\tilde{c}^{\top}\Sigma^{-1/2}X_{i})(r_{1}^{\top}\Sigma^{-1/2}X_{i})(Y_{i}-X_{i}^{\top}\beta)r_{2,i}\right|$
$\displaystyle+\left|\frac{2}{n}\sum_{i=1}^{n}(\tilde{c}^{\top}\Sigma^{-1/2}X_{i})(r_{1}^{\top}\Sigma^{-1/2}X_{i})r_{2,i}^{2}\right|+\left|\frac{1}{n}\sum_{i=1}^{n}(r_{1}^{\top}\Sigma^{-1/2}X_{i})^{2}(Y_{i}-X_{i}^{\top}\beta)^{2}\right|$
$\displaystyle+\left|\frac{1}{n}\sum_{i=1}^{n}(r_{1}^{\top}\Sigma^{-1/2}X_{i})^{2}(Y_{i}-X_{i}^{\top}\beta)r_{2,i}\right|+\left|\frac{1}{n}\sum_{i=1}^{n}(r_{1}^{\top}\Sigma^{-1/2}X_{i})^{2}r_{2,i}^{2}\right|.$

We observe there exists following deterministic inequalities between the
quantities. They are merely inequalities of arithmetic and geometric means.

$\mathcal{Q}_{4}\leq\mathcal{Q}_{3}+\mathcal{Q}_{5},\quad\mathcal{Q}_{5}\leq\mathcal{Q}_{2}+\mathcal{Q}_{8},\quad\mathcal{Q}_{7}\leq\mathcal{Q}_{6}+\mathcal{Q}_{8}.$

Combining these lead to

$|\widehat{\sigma}_{c}^{2}-\tilde{\sigma}_{c}^{2}|\leq\mathcal{Q}_{1}+3\mathcal{Q}_{2}+2\mathcal{Q}_{3}+2\mathcal{Q}_{6}+4\mathcal{Q}_{8}$
(41)

Before we analyze each individual quantity, we observe that

$\displaystyle\|r_{1}\|_{2}=\|(\Sigma^{1/2}\widehat{\Sigma}^{-1}\Sigma^{1/2}-I_{d})\tilde{c}\|_{2}\leq\|\Sigma^{1/2}\widehat{\Sigma}^{-1}\Sigma^{1/2}-I_{d}\|_{\rm
op}\leq\frac{\mathcal{D}_{\Sigma}}{\lambda_{\rm
min}(\Sigma^{-1/2}\widehat{\Sigma}\Sigma^{-1/2})}\leq 2\mathcal{D}_{\Sigma}.$

The last inequality holds on the event $\\{\lambda_{\rm
min}(\Sigma^{-1/2}\widehat{\Sigma}\Sigma^{-1/2})\geq 1/2\\}$. Furthermore, the
following is useful in analyzing $\mathcal{Q}_{8}$:

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}(r_{1}^{\top}\Sigma^{-1/2}X_{i})^{2}$
$\displaystyle=$ $\displaystyle
r_{1}^{\top}\Sigma^{-1/2}\left(\frac{1}{n}\sum_{i=1}^{n}X_{i}X_{i}^{\top}\right)\Sigma^{-1/2}r_{1}=r_{1}^{\top}\left(\Sigma^{-1/2}\widehat{\Sigma}\Sigma^{-1/2}\right)r_{1}$
$\displaystyle=$
$\displaystyle\tilde{c}\left(\Sigma^{1/2}\widehat{\Sigma}^{-1}\Sigma^{1/2}-I_{d}\right)\left(\Sigma^{-1/2}\widehat{\Sigma}\Sigma^{-1/2}\right)\left(\Sigma^{1/2}\widehat{\Sigma}^{-1}\Sigma^{1/2}-I_{d}\right)\tilde{c}$
$\displaystyle=$
$\displaystyle\tilde{c}\left(\Sigma^{1/2}\widehat{\Sigma}^{-1}\Sigma^{1/2}-I_{d}+\Sigma^{-1/2}\widehat{\Sigma}\Sigma^{-1/2}-I_{d}\right)\tilde{c}$
$\displaystyle\leq$
$\displaystyle\mathcal{D}_{\Sigma}+\|\Sigma^{1/2}\widehat{\Sigma}^{-1}\Sigma^{1/2}-I_{d}\|_{\rm
op}\leq 3\mathcal{D}_{\Sigma}.$

We now present deterministic inequalities for each quantity, requiring no
further explanation.

$\displaystyle\mathcal{Q}_{1}$ $\displaystyle=$
$\displaystyle\left|\frac{2}{n}\sum_{i=1}^{n}|\tilde{c}^{\top}\Sigma^{-1/2}X_{i}|^{2}(Y_{i}-X_{i}^{\top}\beta)X_{i}^{\top}(\widehat{\beta}-\beta)\right|\leq
2\mathcal{L}\|\widehat{\beta}-\beta\|_{\Sigma},$
$\displaystyle\mathcal{Q}_{2}$ $\displaystyle=$
$\displaystyle\left|\frac{1}{n}\sum_{i=1}^{n}|\tilde{c}^{\top}\Sigma^{-1/2}X_{i}|^{2}(\widehat{\beta}-\beta)^{\top}X_{i}X_{i}^{\top}(\widehat{\beta}-\beta)\right|\leq\mathcal{M}\|\widehat{\beta}-\beta\|_{\Sigma}^{2},$
$\displaystyle\mathcal{Q}_{3}$ $\displaystyle=$
$\displaystyle\left|\frac{2}{n}\sum_{i=1}^{n}(\tilde{c}^{\top}\Sigma^{-1/2}X_{i})(r_{1}^{\top}\Sigma^{-1/2}X_{i})(Y_{i}-X_{i}^{\top}\beta)^{2}\right|\leq
2\mathcal{L}\|r_{1}\|_{2}\leq 4\mathcal{L}\mathcal{D}_{\Sigma}.$
$\displaystyle\mathcal{Q}_{6}$ $\displaystyle=$
$\displaystyle\left|\frac{1}{n}\sum_{i=1}^{n}(r_{1}^{\top}\Sigma^{-1/2}X_{i})^{2}(Y_{i}-X_{i}^{\top}\beta)^{2}\right|\leq\mathcal{N}\|r_{1}\|_{2}^{2}\leq
4\mathcal{N}\mathcal{D}_{\Sigma}^{2}.$ $\displaystyle\mathcal{Q}_{8}$
$\displaystyle=$
$\displaystyle\left|\frac{1}{n}\sum_{i=1}^{n}(r_{1}^{\top}\Sigma^{-1/2}X_{i})^{2}r_{2,i}^{2}\right|\leq\max_{1\leq
i\leq
n}r_{2,i}^{2}\frac{1}{n}\sum_{i=1}^{n}(r_{1}^{\top}\Sigma^{-1/2}X_{i})^{2}\leq
3\mathcal{D}_{\Sigma}\max_{1\leq i\leq n}r_{2,i}^{2}.$

Combining all in (41) gives the result. ∎

Lemma 7 underscores the intricacies involved in the convergence of the
sandwich variance, which hinge upon several quantities:
$\mathcal{D}_{\Sigma}$, $\|\widehat{\beta}-\beta\|_{\Sigma}$, $\mathcal{M}$,
$\mathcal{N}$, $\mathcal{L}$, and $\max_{1\leq i\leq
n}|X_{i}^{\top}(\widehat{\beta}-\beta)|$. Each of these elements is
individually addressed and examined in Proposition 22, Proposition 27, Lemma
8, Lemma 9, Lemma 10, and Theorem 11.

######  Lemma 8.

For $\tilde{c}\in\mathbb{S}^{d-1}$, define

$\displaystyle\overline{\mathcal{M}}$ $\displaystyle=$
$\displaystyle\sup_{\theta\in\mathbb{S}^{d-1}}\left|\frac{1}{n}\sum_{i=1}^{n}(\theta^{\top}\Sigma^{-1/2}X_{i})^{2}(\tilde{c}^{\top}\Sigma^{-1/2}X_{i})^{2}-\mathbb{E}(\theta^{\top}\Sigma^{-1/2}X)^{2}(\tilde{c}^{\top}\Sigma^{-1/2}X)^{2}\right|$

If Assumption 2.b holds with $q_{x}\geq 4$, then there exists a universal
constant $C$ such that

$\mathbb{E}\left[\overline{\mathcal{M}}\right]\leq
CK_{x}^{4}\left[\sqrt{\frac{\log(2n)d}{n^{1-4/q_{x}}}}+\frac{\log(2n)d}{n^{1-4/q_{x}}}\right].$

Consequently, we have

$\mathbb{E}\left[\mathcal{M}\right]\leq
CK_{x}^{4}\left[1+\frac{\log(2n)d}{n^{1-4/q_{x}}}\right],$

for some (possibly) different universal constant $C$.

###### proof of Lemma 8.

Let $Z_{i}=(\tilde{c}^{\top}\Sigma^{-1/2}X_{i})\Sigma^{-1/2}X_{i}$ for
$i=1,\ldots,n$. The unit ball $\mathbb{S}^{d-1}$ is a symmetric and convex
body of radius 1, and it has a modulus of convexity†††The definition of the
modulus of convexity of geometric objects can be found in Section 2 of Guédon
and Rudelson [2007]. of power type 2. Hence, Theorem 3 of Guédon and Rudelson
[2007] applies and yields,

$\displaystyle\mathbb{E}\sup_{\theta\in\mathbb{S}^{d-1}}\left|\frac{1}{n}\sum_{i=1}^{n}|\theta^{\top}Z_{i}|^{2}-\mathbb{E}|\theta^{\top}Z_{i}|^{2}\right|\leq$
$\displaystyle C\frac{\log(2n)}{n}\left[\mathbb{E}\max_{1\leq i\leq
n}\|Z_{i}\|_{2}^{2}\right]$ (42) $\displaystyle+$ $\displaystyle
C\sqrt{\frac{\log(2n)}{n}}\left[\mathbb{E}\max_{1\leq j\leq
n}\|Z_{i}\|_{2}^{2}\right]^{1/2}\sup_{\theta\in\mathbb{S}^{d-1}}\left(\mathbb{E}|\theta^{\top}Z_{1}|^{2}\right)^{1/2},$

where $C$ is a universal constant. We note that

$\displaystyle\sup_{\theta\in\mathbb{S}^{d-1}}\mathbb{E}[|\theta^{\top}Z_{i}|^{2}]\leq\sup_{\theta\in\mathbb{S}^{d-1}}\mathbb{E}\left[|\theta^{\top}\Sigma^{-1/2}X_{i}|^{4}\right]\leq
K_{x}^{4}.$ (43)

Furthermore, Jensen’s inequality yields that

$\displaystyle\mathbb{E}\max_{1\leq j\leq n}\|Z_{i}\|_{2}^{2}$
$\displaystyle\leq$ $\displaystyle\left(\mathbb{E}\max_{1\leq j\leq
n}\|Z_{i}\|_{2}^{q_{x}/2}\right)^{2/q_{x}}\leq\left(\sum_{j=1}^{n}\mathbb{E}\|Z_{i}\|_{2}^{q_{x}/2}\right)^{2/q_{x}}$
(44) $\displaystyle\leq$ $\displaystyle
n^{4/q_{x}}\left(\mathbb{E}\|Z_{i}\|_{2}^{q_{x}/2}\right)^{4/q_{x}}\leq
n^{4/q_{x}}K_{x}^{4}d.$

Combining inequalities (43) and (44) in (42) gives the desired result. ∎

######  Lemma 9.

For $\tilde{c}\in\mathbb{S}^{d-1}$, define

$\displaystyle\overline{\mathcal{N}}$ $\displaystyle=$
$\displaystyle\sup_{\theta\in\mathbb{S}^{d-1}}\left|\frac{1}{n}\sum_{i=1}^{n}(\theta^{\top}\Sigma^{-1/2}X_{i})^{2}(Y_{i}-X_{i}^{\top}\beta)^{2}-\mathbb{E}(\theta^{\top}\Sigma^{-1/2}X)^{2}(Y_{i}-X_{i}^{\top}\beta)^{2}\right|,$

If $s=(1/q_{x}+1/q)^{-1}\geq 4$, then there exists a universal constant $C$
such that

$\displaystyle\mathbb{P}\left(\overline{\mathcal{N}}\leq
C\frac{\overline{\lambda}K_{x}^{2}K_{y}^{2}}{\underline{\lambda}^{2}}\left[\sqrt{\frac{d\log(2d/\delta)}{n}}+\frac{d\log^{3/2}(2d/\delta)}{\delta^{2/s}n^{1-2/s}}\right]\right)\geq
1-\delta,$

for all $\delta\in(0,1)$.

###### proof of Lemma 9.

Let $W_{i}=\Sigma^{-1/2}X_{i}(Y_{i}-X_{i}^{\top}\beta)$, so that ${\rm
Var}(W_{i})=V$. We use Proposition 22 which proved the concentration of
$\|V^{-1/2}(\sum_{i=1}^{n}W_{i}W_{i}^{\top}/n)V^{-1/2}-I_{d}\|_{\rm op}$. This
implies the concentration of $\overline{\mathcal{N}}$ as Assumption 3
indicates that

$\overline{\mathcal{N}}\leq\|V\|_{\rm
op}\left\|V^{-1/2}\left(\frac{1}{n}\sum_{i=1}^{n}W_{i}W_{i}^{\top}\right)V^{-1/2}-I_{d}\right\|_{\rm
op}\leq\overline{\lambda}\left\|V^{-1/2}\left(\frac{1}{n}\sum_{i=1}^{n}W_{i}W_{i}^{\top}\right)V^{-1/2}-I_{d}\right\|_{\rm
op}.$

To ensure the assumptions made in Proposition 22, we observe from Assumption 3
and Hölder’s inequality that

$\displaystyle\sup_{\theta\in\mathbb{S}^{d-1}}\mathbb{E}\left[|\theta^{\top}V^{-1}W_{i}|^{s}\right]$
$\displaystyle\leq$ $\displaystyle\lambda_{\rm
min}(V)^{-s}\sup_{\theta\in\mathbb{S}^{d-1}}\mathbb{E}\left[|\theta^{\top}W_{i}|^{s}\right]\leq\underline{\lambda}^{-s}\sup_{\theta\in\mathbb{S}^{d-1}}\mathbb{E}\left[\left|(\theta^{\top}\Sigma^{-1/2}X_{i})(Y_{i}-X_{i}^{\top}\beta)\right|^{s}\right]$
$\displaystyle\leq$
$\displaystyle\underline{\lambda}^{-s}\sup_{\theta\in\mathbb{S}^{d-1}}\left[\mathbb{E}(\theta^{\top}\Sigma^{-1/2}X_{i})^{q_{x}}\right]^{s/q_{x}}\left[\mathbb{E}(Y_{i}-X_{i}^{\top}\beta)^{q}\right]^{s/q}\leq\left(\frac{K_{x}K_{y}}{\underline{\lambda}}\right)^{s}.$

Hence, an application of Proposition 22 yields the desired result. ∎

######  Lemma 10.

For $\tilde{c}\in\mathbb{S}^{d-1}$, define

$\displaystyle\mathcal{L}_{1}$ $\displaystyle:=$
$\displaystyle\left\|\frac{1}{n}\sum_{i=1}^{n}\Sigma^{-1/2}X_{i}(\tilde{c}^{\top}\Sigma^{-1/2}X_{i})^{2}(Y_{i}-X_{i}^{\top}\beta)\right\|_{2},$
$\displaystyle\mathcal{L}_{2}$ $\displaystyle:=$
$\displaystyle\left\|\frac{1}{n}\sum_{i=1}^{n}\Sigma^{-1/2}X_{i}(\tilde{c}^{\top}\Sigma^{-1/2}X_{i})(Y_{i}-X_{i}^{\top}\beta)^{2}\right\|_{2}.$

Suppose that Assumption 1 and Assumption 2.b hold.

  1. 1.

If $(3/q_{x}+1/q)^{-1}\geq 2$, then $\mathbb{P}\left(\mathcal{L}_{1}\leq
K_{x}^{3}K_{y}\\{1+\sqrt{d/(n\delta)}\\}\right)\geq 1-\delta$ for all
$\delta\in(0,1)$.

  2. 2.

If $s=(1/q_{x}+1/q)^{-1}\geq 4$, then $\mathbb{P}\left(\mathcal{L}_{2}\leq
K_{x}^{2}K_{y}^{2}\\{1+\sqrt{d/(n\delta)}\\}\right)\geq 1-\delta$ for all
$\delta\in(0,1)$.

Taking $\delta=d/(2n)$ in both inequalities shows that
$\mathbb{P}(\mathcal{L}\leq(1+\sqrt{2})(K_{x}^{3}K_{y}+K_{x}^{2}K_{y}^{2}))\geq
1-d/n$.

###### proof of Lemma 10.

To show the first claim, we let
$S_{i}=|\tilde{c}^{\top}\Sigma^{-1/2}X_{i}|^{2}(Y_{i}-X_{i}^{\top}\beta)\Sigma^{-1/2}X_{i}$
for $i=1,\ldots,n$, and observe that

$\displaystyle\left\|\frac{1}{n}\sum_{i=1}^{n}S_{i}\right\|_{2}\leq\left\|\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}S_{i}\right\|_{2}+\left\|\frac{1}{n}\sum_{i=1}^{n}S_{i}-\mathbb{E}S_{i}\right\|_{2}.$

The expectation of $S_{i}$ can be controlled as

$\displaystyle\|\mathbb{E}S_{i}\|_{2}$ $\displaystyle=$
$\displaystyle\sup_{u\in\mathbb{S}^{d-1}}\mathbb{E}\left[|\tilde{c}^{\top}\Sigma^{-1/2}X_{i}|^{2}(Y_{i}-X_{i}^{\top}\beta)(u^{\top}\Sigma^{-1/2}X_{i})\right]\leq
K_{x}^{3}K_{y},$

provided that $3/q_{x}+1/q\leq 1$. Meanwhile,

$\displaystyle\mathbb{E}\left\|\frac{1}{n}\sum_{i=1}^{n}S_{i}-\mathbb{E}S_{i}\right\|_{2}^{2}$
$\displaystyle=$ $\displaystyle{\rm tr}\left({\rm
Var}\left[\frac{1}{n}\sum_{i=1}^{n}S_{i}\right]\right)=\frac{1}{n}{\rm
tr}\left({\rm Var}\left[W_{1}\right]\right)$ $\displaystyle\leq$
$\displaystyle\frac{1}{n}\mathbb{E}\left[|\tilde{c}^{\top}\Sigma^{-1/2}X_{1}|^{4}(Y_{1}-X_{1}^{\top}\beta)^{2}\|X_{1}\|_{\Sigma^{-1}}^{2}\right]\leq
K_{x}^{6}K_{y}^{2}\frac{d}{n},$

provided that $3/q_{x}+1/q\leq 1/2$. An application of Chebyshev’s inequality
proves the claim.

The second claim follows similarly. Denote
$T_{i}=|\tilde{c}^{\top}\Sigma^{-1/2}X_{i}|(Y_{i}-X_{i}^{\top}\beta)^{2}\Sigma^{-1/2}X_{i}$,
and observe that

$\displaystyle\left\|\frac{1}{n}\sum_{i=1}^{n}T_{i}\right\|_{2}\leq\left\|\frac{1}{n}\sum_{i=1}^{n}\mathbb{E}T_{i}\right\|_{2}+\left\|\frac{1}{n}\sum_{i=1}^{n}T_{i}-\mathbb{E}T_{i}\right\|_{2}.$

The leading term on the right-hand side can be bounded as

$\displaystyle\|\mathbb{E}T_{i}\|_{2}$ $\displaystyle=$
$\displaystyle\sup_{u\in\mathbb{S}^{d-1}}\mathbb{E}\left[|\tilde{c}^{\top}\Sigma^{-1/2}X_{i}|(Y_{i}-X_{i}^{\top}\beta)^{2}(u^{\top}\Sigma^{-1/2}X_{i})\right]\leq
K_{x}^{2}K_{y}^{2},$

provided that $1/q_{x}+1/q\leq 1/2$. Furthermore, we observe that

$\displaystyle\mathbb{E}\left\|\frac{1}{n}\sum_{i=1}^{n}T_{i}-\mathbb{E}T_{i}\right\|_{2}^{2}$
$\displaystyle=$ $\displaystyle{\rm tr}\left({\rm
Var}\left[\frac{1}{n}\sum_{i=1}^{n}T_{i}\right]\right)=\frac{1}{n}{\rm
tr}\left({\rm Var}\left[W_{1}\right]\right)$ $\displaystyle\leq$
$\displaystyle\frac{1}{n}\mathbb{E}\left[|\tilde{c}^{\top}\Sigma^{-1/2}X_{1}|^{2}(Y_{1}-X_{1}^{\top}\beta)^{4}\|X_{1}\|_{\Sigma^{-1}}^{2}\right]\leq
K_{x}^{4}K_{y}^{4}\frac{d}{n},$

provided that $1/q_{x}+1/q\leq 1/4$. The Chebyshev’s inequality completes the
proof. ∎

####  A.7.2 Maximal Concentration of the Fitted Values

In this section, we aim to prove the following theorem. For the convenience of
notation and clarity, we write the weighted harmonic sum of $q_{x}$ and $q$ as
$q_{m,n}=(m/q_{x}+n/q)^{-1}$ for $m,n>0$.

######  Theorem 11.

Suppose that Assumption 1, 2.a, 3 hold for $q_{x}\geq 4$ and
$(1/q_{x}+1/q)^{-1}>2$. There exists constants