Robustness Auditing for Linear Regression: To Singularity and Beyond

## Abstract

Abstract It has recently been discovered that the conclusions of many highly influential econometrics studies can be overturned by removing a very small fraction of their samples (often less than 0.5 % percent 0.5 0.5\% 0.5 % ).
These conclusions are typically based on the results of one or more Ordinary Least Squares (OLS) regressions, raising the question: given a dataset, can we certify the robustness of an OLS fit on this dataset to the removal of a given number of samples?
Brute-force techniques quickly break down even on small datasets.
Existing approaches which go beyond brute force either can only find candidate small subsets to remove (but cannot certify their non-existence) Broderick et al. ( 2020 ); Kuschnig et al. ( 2021 ) , are computationally intractable beyond low dimensional settings Moitra and Rohatgi ( 2022 ) , or require very strong assumptions on the data distribution and too many samples to give reasonable bounds in practice Bakshi and Prasad ( 2021 ); Freund and Hopkins ( 2023 ) .
We present an efficient algorithm for certifying the robustness of linear regressions to removals of samples.
We implement our algorithm and run it on several landmark econometrics datasets with hundreds of dimensions and tens of thousands of samples, giving the first non-trivial certificates of robustness to sample removal for datasets of dimension 4 4 4 4 or greater.
We prove that under distributional assumptions on a dataset, the bounds produced by our algorithm are tight up to a 1 + o ⁢ ( 1 ) 1 𝑜 1 1+o(1) 1 + italic_o ( 1 ) multiplicative factor.

## 1 Introduction

Consider a supervised learning problem with feature vectors $X_{1},\ldots,X_{n}\in\mathbb{R}^{d}$ and labels $Y_{1},\ldots,Y_{n}\in\mathbb{R}$, to which we fit a model $f:\mathbb{R}^{d}\rightarrow\mathbb{R}$.
Robustness auditing addresses the question:

> *How would $f$ have differed if we had been missing a small fraction of the data?*

We study this question in the context of ordinary least squares (OLS) linear regression, where $f(X)=\langle\beta,X\rangle$ is the linear function minimizing the mean squared error $\tfrac{1}{n}\sum_{i\leq n}(f(X_{i})-Y_{i})^{2}$.
We focus on OLS for two reasons.
First, OLS is a statistics workhorse, with widespread use across economics, social science, finance, machine learning, and beyond.
Second, its relative simplicity affords us the opportunity to design algorithms with provable guarantees, and offers a stepping stone to more complex models (logistic regression, kernel methods, neural networks).

###### Problem 1 (Robustness Auditing for OLS Regression) .

Given a linear regression instance $(X_{1},Y_{1}),\ldots,(X_{n},Y_{n})\in\mathbb{R}^{d+1}$, a direction $e\in\mathbb{R}^{d}$, and an integer $k\leq n$, what is

$$ $\Delta_{k}(e)=\max_{\begin{subarray}{c}S\subseteq[n]\\ |S|=n-k\end{subarray}}\langle\beta_{[n]}-\beta_{S},e\rangle\,$ (1) $$

where for $T\subseteq[n]$, $\beta_{T}\in\mathbb{R}^{d}$ denotes the vector of OLS coefficients for the dataset $\{(X_{i},Y_{i})\}_{i\in T}$?

In particular, for a threshold $\theta\in\mathbb{R}$ what is the minimal number of removals $k_{\theta}(e)$ for which $\Delta_{k}(e)>\theta$?

###### Problem 1 (Robustness Auditing for OLS Regression) .

##### Context and Applications for Robustness Auditing

Problem 1 was introduced in this form by Broderick, Giordano, and Meager Broderick et al. (), who use a heuristic algorithm, AMIP, to identify very small subsets of landmark datasets from econometrics which can be removed to overturn important conclusions of the respective studies Finkelstein et al. (); Angelucci and De Giorgi (); often this can be achieved by removing less than $0.5\%$ of a dataset.
Researchers have subsequently used AMIP to audit a wide range of recent studies in economics Martinez (); Di and Xu (); Davies et al. (); Zachmann et al. (); Burton and Roach (); Beuermann et al. (); Bondy et al. ().
Subsequent algorithmic works Kuschnig et al. (); Moitra and Rohatgi (); Freund and Hopkins () develop additional algorithms for auditing a similar notion of robustness; we discuss prior work in detail below.

It is important that, similar to these prior works, we focus on robustness to a shift of $\beta$ in only a single user-specified direction $e$.
This is because the main conclusion of a regression is often determined by the projection of its result on a particular axis.
For instance:

###### Robustness of Parameter Estimate

A researcher may want to estimate the correlation $\langle\beta,e\rangle$ between a specific explanatory variable and a target variable, controlled for additional factors, where $e$ is the indicator vector corresponding to the explanatory variable.
Moreover, the sign and statistical significance of $\langle\beta,e\rangle$ is often of greatest interest.

This correlation can have a causal interpretation.
For instance, in a randomized control trial, where $e$ is the indicator for the treatment variable, the projection $\langle\beta,e\rangle$ can be used to estimate the “average treatment effect” (ATE) of a new drug or policy on the outcome $Y$, while including the control variables in the regression can help reduce the variance of this estimate.
Even more complex causal inferences (e.g., instrumental variables regression) can often be decomposed into a small number of OLS regressions, where the result of the causal inference depends on a single coefficient from the result of each regression.

Conclusions from a study where this shift $\Delta_{k}(e)$ is large when $k\ll n$ are therefore driven by a small number of data points, meriting at minimum reinspection of a dataset, and perhaps casting doubt on generalizability.
In many real world regressions, the sign of $\langle\beta,e\rangle$ is not robust to a small number of removals, even though it is statistically significant Broderick et al. ().
*Non-existence* of a small set of highly influential samples indicates robustness of the measured effect to an interesting class of distribution shifts – any removal of a small fraction of the population.

###### Data Attribution

Suppose that instead of looking for the effect of a particular feature on the label $Y$, instead we use the linear model $f$ to predict the label of a fresh point $X_{\text{new}}$, and we want to identify what part of the training data led to the prediction that $Y_{\text{new}}\approx f(X_{\text{new}})$.
Following the counterfactual formulation of this data attribution problem from Koh and Liang (); Ilyas et al. (), we arrive again at robustness auditing:
using $e=X_{\text{new}}$, we can find the smallest set of whose removal would significantly shift $f(X_{\text{new}})$.
We can evaluate the *brittleness* of the prediction by measuring the size of the smallest set of samples we could remove to cause $f(X_{\text{new}})$ to cross a decision boundary.

##### Intractability, Heuristics and Upper Bounds

As soon as $k$ exceeds single digits, robustness auditing by brute-force search over all $|S|=n-k$ takes times $\binom{n}{k}$, which is computationally intractable.
In fact, under standard computational complexity assumptions, some degree of intractability is inherent Moitra and Rohatgi (), so we need either some assumptions on $\{(X_{i},Y_{i})\}_{i\in[n]}$, a relaxed version of Problem 1, or both.

The works Broderick et al. (); Kuschnig et al. () relax the goal to finding upper bounds on $k_{\theta}(e)$, for which they use greedy/local search algorithms.
This approach leaves open the risk that $k_{\theta}(e)$ might be much smaller than the upper bound suggests.
Indeed, later experiments by us (see Figure 2), Moitra and Rohatgi (), and Freund and Hopkins () uncover numerous real-world examples where local search techniques give loose upper bounds.
Figure 1 illustrates a simple synthetic example where a small number of datapoint removals can shift $\langle\beta,e\rangle$, but greedy/local search algorithms only find much larger sets.

Following Moitra and Rohatgi (), we aim for algorithms which provide *unconditionally valid upper and lower bounds* on $k_{\theta}(e)$ for every dataset $\{(X_{i},Y_{i})\}_{i\in[n]}$, and which return *high-quality* upper/lower bounds (as close to matching each other as possible) under reasonable assumptions on $\{(X_{i},Y_{i})\}_{i\in[n]}$.
Prior approaches to go beyond greedy algorithms and provide lower bounds on $k_{\theta}(e)$ Klivans et al. (); Bakshi and Prasad (); Moitra and Rohatgi (); Freund and Hopkins () so far don’t yield results in practice for datasets of dimension $4$ or greater, due to running times which scale exponentially in $d$ and/or prohibitively strong assumptions on $X_{1},\ldots,X_{n}$.

Figure: (a) A robust dataset whose robustness cannot be certified by continuous algorithms such as Bakshi and Prasad (); Freund and Hopkins ().
Refer to caption: extracted/5914814/figures/original_dataset.png

### 1.1 Our Contributions

We present and analyze two new algorithms, ACRE and OHARE , which provide *lower bounds* on $k_{\theta}(e)$, *certifying robustness* of OLS regression.
Our algorithms provide the first nontrival bounds on the number of samples which must be removed to flip the signs of important parameter estimates for benchmark datasets studied in prior work with hundreds of dimensions and tens of thousands of samples.
We evaluate our algorithms experimentally and in theory.

ACRE (Algorithm for Certifying Robustness Efficiently) takes as input a dataset $X\in\mathbb{R}^{n\times d},Y\in\mathbb{R}^{n}$ and a vector $e\in\mathbb{R}^{d}$, runs in time $O\left(n^{2}d+n^{2}\log(n)\right)$, and outputs a set of upper and lower bounds $U,L\in\mathbb{R}^{n}$ on the removal effects such that $U_{k}\geq\Delta_{k}(e)\geq L_{k}$.
In particular, this runtime avoids exponential dependence on $k$, $d$, and $n$.
The upper and lower bounds ACRE provides are valid even making *no assumptions whatsoever* on $X$ and $Y$.
The bounds are *good*, meaning that the upper and lower bounds are close to matching, when $X$ and $Y$ are drawn from a sufficiently “nice” distribution (such as a linear model with subgaussian features and labels) – we formalize this below.

However, ACRE provides poor bounds – with a wide gap between $U_{k}$ and $L_{k}$ – on a very important class of datasets: those using one-hot encodings (also known as indicator or dummy variables) to express categorical features.
Even though one-hot encoded datasets can yield robust regressions, certifying this is challenging because of singularities which emerge in the covariance when samples are removed (see Figure 1).

OHARE (One-Hot aware Algorithm for certifying Robustness Efficiently) extends ACRE to certify robustness of datasets with a mix of continuous and categorical features.
It uses dynamic programming in conjunction with a fine-grained linear-algebraic analysis of the contribution of categorical features.
OHARE takes as input a dataset $X\in\mathbb{R}^{n\times(d+m)}$ (where the last $m$ features represent a one-hot encoding) and a direction $e\in\mathbb{R}^{d}$, and runs in time $O\left(n^{2}(d+m)+n^{2}m\log(n)\right)$.

The most important property of both ACRE and OHARE is that the bounds they produce are *valid* regardless of any assumptions on the dataset:

###### Theorem 1.1 (Correctness) .

Given $e,X$, and $Y$, ACRE and OHARE output lists of upper/lower bounds $U,L\in\mathbb{R}^{n}$ s.t.

$$ $\forall k\in[n]\;\;\;\;\;L_{k}\leq\Delta_{k}(e)=\max_{\begin{subarray}{c}S \subseteq[n]\\ {\left\lvert{S}\right\rvert}=n-k\end{subarray}}\left\langle{\beta-\beta_{S}},{ e}\right\rangle\leq U_{k}$ $$

The proof of Theorem 1.1 is given in Section 3.3 for ACRE and Section 4 for OHARE .
On its own, Theorem 1.1 says little about the usefulness of the upper and lower bounds $L_{k},U_{k}$.
We provide two types of evidence that the bounds produced by ACRE and OHARE are interesting.
First, we demonstrate on real-world econometric datasets studied in prior work on robustness auditing that ACRE and OHARE produce significantly better lower bounds bounds than were previously known (see Table
1 and Figure 2).
Second, we prove that ACRE and OHARE both produce nearly-matching upper and lower bounds under relatively mild distributional assumptions on $X$ and $Y$, for the interesting range of $k$.

###### Theorem 1.1 (Correctness) .

###### Interesting Values of k 𝑘 k italic_k

For a direction $e$, we emphasize two values of $k$: $k_{2\sigma}(e)$, the number of removals needed to shift $\langle\beta,e\rangle$ outside its $95\%$ confidence interval, and $k_{\textup{sign}}(e)=k_{\langle\beta,e\rangle}$, the number to flip the sign of $\langle\beta,e\rangle$.
In the parameter estimation setting, $\langle\beta,e\rangle$ and $2\sigma$ are often of similar magnitude, because rejecting a null hypothesis often involves placing the estimator $\langle\beta,e\rangle$ in a confidence interval which does not contain 0 0.

##### Experimental Results

For comparison with prior works, we focus our experiments on $k_{\textup{sign}}(e)$.
We provide lower bounds on $k_{\textup{sign}}(e)$ for benchmark datasets drawn from important studies in economics and social sciences, first investigated in the context of robustness auditing by Broderick et al. ().

We study real-world datasets corresponding to each of the parameter estimation use-cases listed above: *Nightlights* Martinez () (correlation controlled for additional features), *Cash Transfer* Angelucci and De Giorgi () (randomized control trial), and the *Oregon Health Insurance Experiment (OHIE)* Finkelstein et al. () (IV regression), with 14 distinct linear or instrumental-variables regressions drawn from the corresponding papers, all of which appeared in top econometrics journals.
In many cases, our lower bounds match known upper bounds up to a factor of $2$ or $3$, where no nontrivial lower bounds were previously known; see Table 1 and Figure 2.

To illustrate, consider Nightlights: Martinez () studies whether democractic countries publish more accurate economic growth estimates than dictatorships, after controlling for variables like regional stability and wealth in natural resources.
Martinez formulates this as a linear regression with dimension $d=209$, over $200$ of which correspond to one-hot encoded categorical variables.
The sign and statistical significance of a single coordinate of $\beta$ govern the conclusion of the study.
Algorithms from prior work find a subset of $2.8\%$ of the samples which can be removed to reverse Martinez ()’s main conclusion, but prior algorithms could not rule out the existence of much smaller subsets.
Our algorithm OHARE provides a certificate that no subset of $\leq 0.7\%$ of the samples would reverse the study’s conclusion.
See Table 1 for our results on Nightlights and the other 13 regressions we audit, and Section 6 for detailed discussion of our experiments.

Figure: (a) A comparison of OHARE and known upper bounds on benchmark datasets.
Refer to caption: extracted/5914814/figures/bar_plot.png

**Table 1: All experiments were run on a single core of an AMD EPYC 9654 96-Core Processor. For detailed notes on the data analysis and experimental procedures, see Section 6.**
| Paper | Regression | n | d | AMIP | KZC21 | OHARE | Runtime | Memory |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Nightlights |  | 3895 | 209 | 136 | 110 | 29 | 25 s | 3.81 GiB |
| Cash<br>Transfer | (np, 8) | 4543 | 18 | 5 | 5 | 5 | 6 s | 1.78 GiB |
|  | (np, 9) | 3769 | 18 | 21 | 21 | 17 | 4 s | 2.80 GiB |
|  | (np, 10) | 4191 | 18 | 26 | 26 | 20 | 5 s | 1.73 GiB |
|  | (p, 8) | 10781 | 18 | 225 | 224 | 119 | 50 s | 11.20 GiB |
|  | (p, 9) | 9489 | 18 | 321 | 314 | 126 | 24 s | 5.37 GiB |
|  | (p, 10) | 10368 | 18 | 570 | 555 | 178 | 29 s | 6.59 GiB |
| OHIE | Health genflip | 23361 | 18 | 257 | 257 | 77 | 9 m 5 s | 34.19 GiB |
|  | Health notpoor | 23361 | 18 | 149 | 149 | 40 | 9 m 10 s | 43.50 GiB |
|  | Health change flip | 23407 | 18 | 184 | 184 | 52 | 9 m 10 s | 43.55 GiB |
|  | Not bad days total | 21881 | 18 | 72 | 72 | 21 | 7 m 39 s | 38.64 GiB |
|  | Not bad days physical | 21384 | 18 | 84 | 84 | 25 | 7 m 17 s | 45.77 GiB |
|  | Not bad days mental | 21601 | 18 | 118 | 118 | 31 | 7 m 20 s | 46.21 GiB |
|  | Nodep Screen | 23147 | 18 | 116 | 116 | 32 | 8 m 57 s | 51.53 GiB |

An implementation of our algorithms is available via [Github](https://github.com/ittai-rubinstein/ols_robustness).
Our implementation is efficient enough to run our algorithms with $n$ up to $3\times 10^{4}$ and $d$ up to $10^{3}$ on a single CPU core with $<64GB$ of RAM.
The main bottleneck in practice is storing $3$ floating-point matrices of size $n\times n$ each.

##### Theory for ACRE

Without some assumptions on $X,Y$, finding matching upper and lower bounds on $\Delta_{k}(e)$ is computationally intractable, under standard complexity assumptions Moitra and Rohatgi ().
So, we analyze tightness of the bounds produced by ACRE and OHARE under some relatively mild distributional assumptions on $X$ and $Y$.
Our main assumption for ACRE will be that the samples $X_{1},\ldots,X_{n}$ are drawn iid from a well-behaved distribution:

###### Definition 1 (Well-Behaved Distribution) .

We say that a mean-zero distribution $\mathcal{X}$ on $\mathbb{R}^{d}$ with covariance $\Sigma={\mathop{\mathbb{E}}_{X\sim\mathcal{X}}\left[{XX^{\intercal}}\right]}$ is *well-behaved* if it has exponentially decaying tails in the sense that

$$ $\exists C>0\;\;\forall v\in\mathbb{S}^{d-1},t>0\;\;\;{\Pr_{X\sim\mathcal{X}} \left[{{\left\lvert{\left\langle{v},{\Sigma^{-1/2}X}\right\rangle}\right\rvert }>t}\right]}\leq\exp\left(-\Omega\left(t^{C}\right)\right)\,.$ $$

Any subgaussian or subexponential distribution is well-behaved (with $C=2,1$, respectively), and distributions with even heavier tails than subexponential can still be well-behaved.

We will assume that the labels $Y_{1},\ldots,Y_{n}$ are drawn according to a ground-truth linear model specified by an unknown vector $\beta_{\textup{gt}}\in\mathbb{R}^{d}$.
Concretely, we assume $Y_{i}=\langle\beta_{\textup{gt}},X_{i}\rangle+\epsilon_{i}$, where $\epsilon_{1},\ldots,\epsilon_{n}$ are iid from $\mathcal{N}(0,\sigma^{2})$ for an unknown variance parameter $\sigma^{2}>0$.

###### Theorem 1.2 (Tightness of ACRE , Proof in Section 7 ) .

Let $n,d\in\mathbb{N},\;e,\beta_{\textup{gt}}\in\mathbb{R}^{d},\;X,Y$, and $\sigma>0$ be as above.
There exists $k_{\textup{threshold}}=\widetilde{\Theta}\left(\min\left\{\frac{n}{\sqrt{d}},
\frac{n^{2}}{d^{2}}\right\}\right)$ such that for all $k\leq k_{\textup{threshold}}$, the bounds $L_{k},U_{k}$ produced by ACRE satisfy

$$ $\frac{U_{k}}{L_{k}}=1+\widetilde{O}\left(\frac{d+k\sqrt{d}}{n}\right)\,.$ (2) $$

Note that this theorem holds even when the covariance $\Sigma$ of $\mathcal{X}$ is unknown *a priori*.

###### Definition 1 (Well-Behaved Distribution) .

###### Theorem 1.2 (Tightness of ACRE , Proof in Section 7 ) .

###### Relationship of k 𝑘 k italic_k , d 𝑑 d italic_d , and n 𝑛 n italic_n

Theorem 1.2 guarantees that ACRE gives nearly optimal bounds so long as $n\gg d$ and the number of removals is at most $k\leq\widetilde{\Theta}(\min\{n/\sqrt{d},n^{2}/d^{2}\})$; note that in this regime the RHS of (2) is $1+o(1)$.
First of all, OLS is only appropriate when $n\geq d$, and all of the datasets on which we perform experiments have $n$ well in excess of $d$.
The range of interesting values of $k$ is more subtle – to illustrate, consider $k_{2\sigma}(e)$.
Assuming reasonably well-behaved samples, we expect $k_{2\sigma}(e)=O(\sqrt{n})$, so as long as we also have $n\gg d^{4/3}$, ACRE gives nearly tight bounds on $k_{2\sigma}(e)$.
This represents a mild multiplicative overhead of $d^{1/3}$ samples compared to the $n\geq d$ required anyway to use OLS regression.

##### Theory for OHARE

Even though they are mild, the assumptions for ACRE are not satisfied by the real-world datasets used in our experiments, because of the presence of one-hot encoded categorical features, which is the reason we design OHARE in the first place.
We also prove a tightness theorem for OHARE – we now turn to the assumptions on $X,Y$ which underlie it.

We study regression in $d+m$ dimensions – that is, $X\in\mathbb{R}^{n\times(d+m)}$.
The block of $m$ coordinates will be a one-hot encoding of a categorical variable, while the block of $d$ coordinates will act as in the ACRE setting.
Formally, let $n\in\mathbb{N}$ and let $B_{1},\ldots,B_{m}$ partition $[n]$ into $m$ *buckets*.
Let $e,\beta_{\textup{gt}}\in\mathbb{R}^{d+m}$, and $\sigma>0$.
For $i\in[n]$, let $X_{1}^{\prime},\ldots,X_{n}^{\prime}$ be iid draws from a $d$-dimensional well-behaved distribution and let $X_{i}$ be $X_{i}^{\prime}$ concatenated with the $j$-th indicator vector, where $j=j(i)$ is such that $i\in B_{j}$.
Finally, let $Y_{i}=\langle\beta_{\textup{gt}},X_{i}\rangle+\epsilon_{i}$ where $\epsilon_{1},\ldots,\epsilon_{n}\sim\mathcal{N}(0,\sigma^{2})$.

###### Theorem 1.3 (Tightness of OHARE , informal, see Section 8 ) .

Let $X,Y$ be as described above.
For any arbitrarily small $\varepsilon>0$, if all buckets have sizes $n^{\varepsilon}\sqrt{d}<|B_{j}|<0.49n$, and $n\geq d^{5/4+o(1)}$, then with high probability, for all $k<k_{\textup{threshold}}$, where
$k_{\textup{threshold}}=\widetilde{\Theta}\left(\min\left\{\frac{n}{\sqrt{d}},
\frac{n^{2}}{d^{2}},n^{1-\varepsilon}\right\}\right)$, the upper and lower bounds produced by OHARE satisfy

$$ $\frac{U_{k}}{L_{k}}=1+O\left(\frac{1}{\sqrt{\log n}}\right)\,.$ $$

Careful analysis is needed to prove Theorem 1.3 under the relatively weak assumptions $|B_{j}|>n^{\varepsilon}\sqrt{d}$ and $n\geq d^{5/4+o(1)}$.
The stronger assumptions $|B_{j}|>d$ or $n\geq d^{2}$ would simplify the analysis.
But neither assumption would be valid for all of our real-world datasets.
Getting away with such weak assumptions ultimately requires us to put together a number of technical tools, including novel matrix concentration arguments (e.g. Lemma 8.12).

Additionally, the error term $1/\sqrt{\log n}$, which we believe is nearly tight, goes to zero slowly compared to the error term in Theorem 1.2 – to capture this fine-grained behavior of $U_{k}/L_{k}$, our proof carefully exploits Gaussian anticoncentration.

###### Theorem 1.3 (Tightness of OHARE , informal, see Section 8 ) .

### 1.2 Prior Work and Future Directions

##### Prior Work

Robustness in linear regression is too vast to survey here, so we restrict attention to recent works in robustness auditing.
As previously discussed, our work is preceded by Broderick et al. (); Kuschnig et al. (); Moitra and Rohatgi (); Freund and Hopkins (), all studying robustness auditing for least-squares regression.
Unlike any of these prior works, our algorithms provide nontrivial lower bounds on $k_{\theta}$ in practice for datasets with tens/hundreds of dimensions.
All these algorithms and ours fit into the broader tradition of influence functions as a measure of robustness in regression, dating at least to Cook and Weisberg ().

Our algorithms are partly inspired by recent developments in algorithmic robust statistics – see Diakonikolas and Kane () and references therein.
While the main goal in robust statistics differs somewhat from robustness auditing, modern algorithms for robust regression typically contain subroutines for tasks very similar to robustness auditing.
But the subroutines in recent breakthroughs in robust regression Klivans et al. (); Bakshi and Prasad () are not practically implementable due to reliance on semidefinite programming.

##### Future Directions

We propose several directions for future work:

###### Beyond OLS

Given the prevalence of regression beyond ordinary least squares in machine learning, an important next step is to design algorithms to certify robustness of other regression methods which arise by minimizing a convex loss – e.g. logistic regression and LASSO Tibshirani ().

###### Applications to Differential Privacy (DP)

DP Dwork et al. () is the gold-standard mathematically rigorous approach to protecting privacy of individuals represented in a dataset.
Tighter certificates of robustness for regression have great potential to improve *privacy-accuracy tradeoffs* in private regression Dwork and Lei (); poor privacy-accuracy tradeoffs are a major roadblock to widespread adoption of private data analysis techniques.

###### Interpretations of Robustness Certificates

It is an appealing intuition that statistical conclusions which are robust to removing many samples should generalize better – formalizations of the relationship between stability and generalization have been highly influential, e.g. Bousquet and Elisseeff ().
Can rigorous interpretations of robustness certificates be formalized?
Can robustness certificates yield e.g. tighter empirical confidence intervals, or out-of-distribution generalization guarantees?

###### Beyond Sample Removal

Removing a small set of datapoints is just one of many potential ways to perturb a dataset.
Can we certify robustness of OLS or other regression algorithms to other types of dataset perturbation, e.g. $\ell_{2}$ or $\ell_{\infty}$-bounded perturbations of the feature vectors?

### 1.3 Organization of this Paper

In Section 2 we will define some of the notation used throughout the paper.

Sections 3, 4 and 5 will cover the main algorithmic contributions of our paper.
In Section 3 we present the ACRE algorithm for certifying robustness of OLS from continuous features and a proof of its correctness, as well as the main “MSN-bounding algorithm” we use throughout the paper.
In Section 4 we construct the OHARE algorithm by extending the ACRE framework to the case where some of the dimensions of the regression represent a one-hot encoding of a categorical feature, and in Section 5 we present additional MSN-bounding algorithms.

In Section 6, we give an overview of the econometrics datasets used for our evaluation of the OHARE algorithm, describe our experimental setup, and present the results of running our OHARE algorithm on these regressions (Table 1).

Sections 7 and 8 are devoted to a theoretical analysis of the ACRE and OHARE algorithms, and contain the proofs of Theorem 1.2 and Theorem 1.3 respectively.
We note that the proof of Theorem 1.3 is somewhat challenging and requires some non-trivial concentration bounds which may be of independent interest (see Section 8.3).

Finally, we defer some of the simpler proof to Appendix A and Appendix B.

## 2 Preliminaries

Let $(X_{1},Y_{1}),\ldots,(X_{n},Y_{n})\in\mathbb{R}^{(d+1)}$ represent features/covariates and labels/target variables of a linear regression instance.
We always assume $n\geq d$.
Let $\Sigma=X^{\intercal}X\in\mathbb{R}^{d\times d}$ denote the (un-normalized) empirical second moment of $X$, let $\beta$ denote the OLS fit $\beta=\Sigma^{-1}X^{\intercal}Y\in\mathbb{R}^{d}$, and
let $R=Y-\beta X\in\mathbb{R}^{n}$ denote the residuals on the complete regression.
Finally, let $e\in\mathbb{R}^{d}$ be some “direction of interest” along which we wish to certify the robustness of $\beta$.
Using standard normalization techniques, we may ensure that $e\in\mathbb{S}^{d-1}$ has norm $1$ and $\Sigma=I$.

For any $S\subseteq[n]$ representing some subset of the samples, let $X_{S},Y_{S}$ represent the samples limited to only those whose indices lie in $S$.
Similarly, let $\Sigma_{S},\beta_{S}$ denote the empirical second moment and regression when using only the samples in $S=\overline{T}$.
Note that while we could set $\Sigma=I$, removing some of the samples may change the covariance matrix $\Sigma_{S}=\Sigma-\Sigma_{T}\neq\Sigma$.

Finally, we use standard asymptotic notation $O(\cdot),\Theta(\cdot)$, and we write $f(n)=\widetilde{O}(g(n))$ if there is a constant $C$ such that $f(n)=O(g(n)\cdot\log^{C}n)$, and similarly for $\widetilde{\Theta}$.

## 3 ACRE: Certifying Robustness without Categorical Features

In this section we present ACRE (Algorithm for Certifying Robustness Efficiently), our algorithm for certifying robustness of regressions without categorical features.

### 3.1 Separating First Order and Higher Order Effects on the OLS fit β 𝛽 \beta italic_β

We split the effect of data removal on $\beta$ into a first order term and a higher order correction.
The first order term is linear in the datapoints $X_{i}$, allowing us to analyze it exactly.
For well-behaved datasets, the higher order term is smaller in magnitude, so even very loose bounds will suffice.

More concretely, let $T\subseteq[n]$ be a set of $k={\left\lvert{T}\right\rvert}\ll n$ samples we might remove from the regression data, and set $S=[n]\setminus T$.
A simple analysis yields the identity

$$ $\beta-\beta_{S}=\Sigma_{S}^{-1}\sum_{i\in T}R_{i}X_{i}$ (3) $$

where $R$ denotes the residuals of the original regression and $\Sigma_{S}$ denotes the empirical 2nd moment of the retained samples $\Sigma_{S}=X_{S}^{\intercal}X_{S}=\sum_{i\in S}X_{i}X_{i}^{\intercal}$.

The difficulty in analyzing equation (3) is the effect induced by the non-linear matrix inversion operation.
Recall that we normalized our datasets so that the empirical second moment over the entire dataset $\Sigma$ is the identity matrix.
$\Sigma_{S}$ is generated by removing some of the samples, so it is no longer normalized.

Because only a very small number of samples were removed, we might hope $\Sigma_{S}=I-\Sigma_{T}$ is still close to the identity matrix.
Therefore, it makes sense to try to develop equation (3) in orders of $\Sigma_{T}$.

More concretely, we use the identity $(I-\Sigma_{T})^{-1}=I+(I-\Sigma_{T})^{-1}\Sigma_{T}$ to derive:

$$ $\displaystyle\beta-\beta_{S}$ $\displaystyle=\Sigma_{S}^{-1}\sum_{i\in T}R_{i}X_{i}=\underbrace{\sum_{i\in T} R_{i}X_{i}}_{\text{first order term}}+\underbrace{\Sigma_{S}^{-1}\Sigma_{T} \sum_{i\in T}R_{i}X_{i}}_{\text{higher order correction}}$ (4) $$

Projecting the first order term onto some axis $e\in\mathbb{S}^{d-1}$ yields the gradients / influence scores used by the AMIP algorithm of Broderick et al. ().
Our analysis will focus on bounding the higher order term.

### 3.2 Maximal Subset Sum Norm (MSN) – The Backend of ACRE

Under the hood of ACRE is a simple algorithm which places upper/lower bounds on the following optimization problem.

###### Problem 2 ( M aximal S ubset sum N orm (MSN) 1 1 1 When the vectors have 0 0 mean ( $\sum_{i\in[n]}v_{i}=0$ ), MSN is equivalent to resilience Steinhardt et al. ( 2017 ) . ) .

Given a set of vectors $v_{1},\ldots,v_{n}\in\mathbb{R}^{d}$ with Gram matrix $G=\left(\left\langle{v_{i}},{v_{j}}\right\rangle\right)_{i,j\in[n]}$,
we define

$$ $\forall k\in[n]\;\;\;\;\textup{MSN}_{k}\left(G\right)=\max_{\begin{subarray}{c }T\subseteq[n]\\ {\left\lvert{T}\right\rvert}=k\end{subarray}}{\left\|{\sum_{i\in T}v_{i}} \right\|}=\max_{\begin{subarray}{c}T\subseteq[n]\\ {\left\lvert{T}\right\rvert}=k\end{subarray}}{\sqrt{1_{T}^{\intercal}G1_{T}}}$ $$

A constant-factor approximation to MSN would refute the small-set expansion hypothesis Hopkins and Li (), so we aim for *MSN-bounding* algorithms which place upper and lower bounds on the optimal value, with the aim that these bounds are close to tight on well-behaved $v_{i}$s.
For our purposes, a simple algorithm in Section 3.5 gives useful bounds, but ACRE can use any MSN-bounding algorithm as a subroutine, and improved MSN-bounding algorithms will lead to improved performance for ACRE .

For now, we treat MSN-bounding as a black-box, and show how ACRE produces its upper and lower bounds $U_{k},L_{k}$ by making a few calls to an MSN-bounding algorithm.

###### Problem 2 ( M aximal S ubset sum N orm (MSN) 1 1 1 When the vectors have 0 0 mean ( ∑ i ∈ [ n ] v i = 0 subscript 𝑖 delimited-[] 𝑛 subscript 𝑣 𝑖 0 \sum_{i\in[n]}v_{i}=0 ∑ start_POSTSUBSCRIPT italic_i ∈ [ italic_n ] end_POSTSUBSCRIPT italic_v start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT = 0 ), MSN is equivalent to resilience Steinhardt et al. ( 2017 ) . ) .

### 3.3 Reducing Robustness Certification to MSN-bounding

Recall equation (4) for the effect of removals on an OLS regression.
Projecting onto $e$, we have

$$ $\displaystyle\left\langle{e},{\beta-\beta_{S}}\right\rangle$ $\displaystyle=\left\langle{e},{\Sigma_{S}^{-1}\sum_{i\in T}R_{i}X_{i}}\right \rangle=\underbrace{\sum_{i\in T}R_{i}\left\langle{X_{i}},{e}\right\rangle}_{ \text{first order term}}+\underbrace{\left\langle{\Sigma_{S}^{-1}\Sigma_{T}e}, {\sum_{i\in T}R_{i}X_{i}}\right\rangle}_{\text{high order term}}$ (5) $$

We can find the maximizing $T$ for the first order term of equation (5) via a greedy algorithm, allowing us to focus on bounding the maximum value of the high order term for any choice of $T$.
We start with the following bound:

$$ $\displaystyle{\left\lvert{\text{high order term}}\right\rvert}$ $\displaystyle\leq\max(\lambda(\Sigma_{S}^{-1}))\left\|{(\Sigma_{S}-I)e}\right \|\left\|{\sum_{i\in T}X_{i}R_{i}}\right\|$ (6) $$

where $\max(\lambda(\Sigma_{S}^{-1}))$ is the largest eigenvalue of $\Sigma_{S}^{-1}$.
We show that each of the three terms in the RHS of equation (6) can be upper-bounded by the value of an MSN problem.
The last term, $\|\sum_{i\in T}X_{i}R_{i}\|$, is already corresponds to an MSN problem with the Gram matrix $G_{RX}=\mathop{\text{diag}}\left({R}\right)G_{X}\mathop{\text{diag}}\left({R}\right)$, so we focus on the other two terms.

##### MSN bound for max ⁡ ( λ ⁢ ( Σ S − 1 ) ) 𝜆 superscript subscript Σ 𝑆 1 \max(\lambda(\Sigma_{S}^{-1})) roman_max ( italic_λ ( roman_Σ start_POSTSUBSCRIPT italic_S end_POSTSUBSCRIPT start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT ) )

We start by simplifying the $\Sigma_{S}^{-1}$ term.
Recall that $\Sigma_{S}=I-\Sigma_{T}$, allowing us to use

$$ $\max\left\{{\lambda\left(\Sigma_{S}^{-1}\right)}\right\}\leq\frac{1}{1-\max \left\{{\lambda(\Sigma_{T})}\right\}}$ (7) $$

so long as $\max\left\{{\lambda\left(\Sigma_{T}\right)}\right\}\leq 1$.
We apply the inequality

$$ $\max\left\{{\lambda\left(\Sigma_{T}\right)}\right\}\leq\sqrt{\sum_{\lambda_{i} \in\lambda\left(\Sigma_{T}\right)}\lambda_{i}^{2}}=\left\|{\Sigma_{T}}\right\| _{F}=\left\|{\sum_{i\in T}X_{i}\otimes X_{i}}\right\|$ $$

where $\|M\|_{F}$ is the Frobenius norm of a matrix $M$ ($\ell_{2}$ norm of the vector of eigenvalues) and $\otimes$ denotes tensor product.
So we have an MSN problem with the vectors $X_{i}\otimes X_{i}$.

Even though this MSN is represented by $n$ vectors of dimension $d^{2}$, its Gram matrix representation $G_{X\otimes X}$ can still be computed in time $O(n^{2}d)$, by computing the Gram matrix $G_{X}=\left(\left\langle{X_{i}},{X_{j}}\right\rangle\right)_{i,j\in[n]}$, and squaring its entries to obtain

$$ $G_{X\otimes X}=\left(\left\langle{X_{i}\otimes X_{i}},{X_{j}\otimes X_{j}} \right\rangle\right)_{i,j\in[n]}=\left(\left\langle{X_{i}},{X_{j}}\right \rangle^{2}\right)_{i,j\in[n]}$ $$

##### MSN bound for Σ T ⁢ e subscript Σ 𝑇 𝑒 \Sigma_{T}e roman_Σ start_POSTSUBSCRIPT italic_T end_POSTSUBSCRIPT italic_e

Recall
$\Sigma_{T}=\sum_{i\in T}X_{i}X_{i}^{\intercal}$.
Let $G_{X\langle X,e\rangle}$ be the Gram matrix of $\{X_{i}\langle X_{i},e\rangle\}_{i\in[n]}$.
So,

$$ $\left\|{\Sigma_{T}e}\right\|=\left\|{\sum_{i\in T}X_{i}\langle X_{i},e\rangle} \right\|\leq\text{MSN}\left(G_{X\langle X,e\rangle}\right)\,,$ $$

### 3.4 Algorithm

Combining the results of the analysis above yields the following expression and the corresponding Algorithm 1:

$$ ${\left\lvert{\Delta_{k}(e)-\text{(first-order term)}_{k}}\right\rvert}\leq \frac{1}{1-\text{MSN}_{k}\left(G_{X\otimes X}\right)}\cdot\text{MSN}_{k}\left( G_{XR}\right)\cdot\text{MSN}_{k}\left(G_{XZ}\right)$ $$

Figure: Algorithm 1 ACRE (Algorithm for Certifying Robustness Efficiently)

### 3.5 Refined Triangle Inequality

In this section we discuss the main MSN-bounding algorithm we use as the backend of ACRE and OHARE for all our experiemnts on real-world data.
We call this algorithm Refined Triangle Inequality (RTI).
We implemented other MSN-bounding algorithms, in particular one based on eigenvalues and eigenvectors of the Gram matrix $G$, but found that they improved over RTI only on synthetic data, so we defer them to a later section.

RTI relies on the following inequality, evaluating the RHS via a greedy algorithm:

$$ $\displaystyle\max_{|T|=k}\left\|{\sum_{i\in T}X_{i}}\right\|^{2}$ $\displaystyle=\max_{|T|=k}\sum_{i,j\in T}\langle X_{i},X_{j}\rangle\leq\max_{| T|,|S_{1}|,\ldots,|S_{k}|=k}\sum_{i\in T,j\in S_{i}}\langle X_{i},X_{j}\rangle\,.$ (8) $$

Figure: Algorithm 2 Refined Triangle Inequality

## 4 Categorical-Aware Robustness Analysis (OHARE)

Suppose that $X_{1},\ldots,X_{n}\in\mathbb{R}^{d+m}$ consist of $d$ continuous-valued features and a single categorical feature with $m$ categories, one-hot encoded.

The bounds computed by ACRE are valid for such $X_{1},\ldots,X_{n}$.
But, consider removing the samples $T$ comprising any single category, encoded in coordinate $i$; let $S=[n]\setminus T$.
The matrix $\Sigma_{S}$ is singular, since all variance from category $i$ has been removed.
So our approach to bounding the high-order term from (6) by pulling out $\max(\lambda(\Sigma_{S}^{-1}))$ is doomed to failure once $k$ exceeds the size of the smallest category, since $\max(\lambda(\Sigma_{S}^{-1}))\rightarrow\infty$.
We will assume that the direct of interest $e$ lies orthogonal to the one-hot features, so $e\in\mathbb{R}^{d}$.

The key idea in OHARE is to rephrase the OLS algorithm as a two-phase process, first explicitly controlling for the categorical feature by computing a “controlled” dataset $\{(\widetilde{X}_{i},\widetilde{Y}_{i})\}_{i\in[n]}\subseteq\mathbb{R}^{d+1}$, then performing OLS on the controlled dataset to arrive at $\beta\in\mathbb{R}^{d}$.
We show that this process yields the same $\beta$ which would have been produced by running OLS on the original dataset $\{(X_{i},Y_{i})\}_{i\in[n]}$ (Claim 4.1) and restricting to the span orthogonal to the one-hot encoding.

We derive explicit formulae for $\widetilde{X}_{i},\widetilde{Y}_{i}$ by analyzing the Gram-Schmidt orthogonalization process which we call “reaveraging”.
The upshot is that we replace each term in (6) as well as the first order effect, with two terms: one corresponding to the direct effect of sample removals on the $\widetilde{X},\widetilde{Y}$ regression, and one corresponding to effects of removing $X_{T},Y_{T}$ on the remaining controlled samples $\widetilde{X}_{S},\widetilde{Y_{S}}$ through reaveraging.
Crucially, this process allows us to certify that the matrix $\widetilde{\Sigma}_{S}$ is nonsingular even in cases where $\Sigma_{S}$ can be singular.

To bound the new correction term coming from the influence of categorical features on the sample-removal effect, we use a knapsack-style dynamic program to combine bounds on the influence of data removals from each category into a single bound by searching over partitions $k=k_{1}+\ldots+k_{m}$.

### 4.1 Regressing over a Categorical Feature

Consider a linear regression over a set of feature vectors, $F_{i}\in\mathbb{R}^{d}$ for $i\in n$, and an additional $m$ dimensions corresponding to a one-hot encoding with buckets $B_{1}\sqcup\cdots\sqcup B_{m}=[n]$, where we want to fit our labels $L_{i}$ to a model of the form

$$ $L_{i}\approx\sum_{j\in[d]}\beta_{j}F_{i,j}+\sum_{j\in[m]}t_{j}1_{i\in B_{j}}$ $$

This linear regression would be represented by a matrix $X\in\mathbb{R}^{n\times(m+d)}$ whose rows are the samples $X_{i}=\left(F_{i},1_{i\in B_{1}},\ldots,1_{i\in B_{m}}\right)$.
Then, if we are interested in the correlation between one of the continuous features $F_{:,j}$ and the labels $Y$ (controlled for the rest of our categorical / continuous features), we may compute this correlation by running an OLS regression $\beta=(X^{\intercal}X)^{-1}X^{\intercal}Y$ and output the relevant entry $\beta_{j}$ of the fit.

Analyzing the robustness of this process is challenging, and we proceed by performing the Gram Schmidt orthogonalization between the dummy variables and the continuous features explicitly.
Before delving into the details of this analysis, we note that OHARE can be used to certify robustness for a slightly more general class of regressions.
Using this more general notation will help motivate our orthogonalization process and is crucial for certifying the robustness of weighted regressions with categorical features (such as the ones in the OHIE study - see Section 6.3).

Indeed, let $u_{1},\ldots,u_{m}$ be the columns representing the dummy variables.
So long as for any sample $i\in[n]$ there is a unique $j=b(i)$ (representing the bucket to which the $i$th sample belongs) such that $u_{j,i}\neq 0$ (and for all $j\neq j^{\prime}$, we have $u_{j^{\prime},i}=0$).
This property ensures that for any $S\subseteq[n]$, these columns are still perpendicular to one another $u_{j,S}\perp u_{j^{\prime},S}$.

- 1.
For each bucket $B_{j}\subseteq[n]$, we compute the weighted averages over the features $f_{j,i}=\frac{uu^{\intercal}}{\left\|{u}\right\|^{2}}F=\frac{\sum_{i\in B_{j}}
u_{j,i}f_{i,:}}{\left\|{u_{j}}\right\|^{2}}u_{j,i}$ and over the labels $\ell_{j}=\frac{\sum_{i\in B_{j}}u_{j,i}\ell_{i}}{\left\|{u_{j}}\right\|^{2}}u_
{j,i}$ for samples from this bucket.
When the samples are unweighted $u_{j,i}=1_{i\in B_{j}}$, these are equal to the averages over the features / labels $f_{j}={\mathop{\mathbb{E}}_{i\in B_{j}}\left[{F_{i,:}}\right]}$ and $\ell_{j}={\mathop{\mathbb{E}}_{i\in B_{j}}\left[{L_{i}}\right]}$, and when the samples are weighted, these are the projections of the continuous features / labels onto the space spanned by the indicator columns.
- 2.
Compute the normalized features $X\in\mathbb{R}^{d\times n}$ obtained by subtracting the feature averages of each bucket from its samples, and the normalized labels $Y\in\mathbb{R}^{n}$ obtained by subtracting the average label from each bucket:
$\displaystyle X_{i}=F_{i,:}-f_{b(i),i}$
$\displaystyle Y_{i}=L_{i}-\ell_{b(i),i}$
- 3.
Perform a linear regression on $Y$ as a function of $X$ and output the fit $\beta=(X^{\intercal}X)^{-1}X^{\intercal}Y$.

###### Claim 4.1 .

The output of the process described above $\beta$ is equal to the coefficients for the continuous features on a full OLS with the dummy variables.

The proof of Claim 4.1 follows directly by performing the Gram-Schmidt orthogonalization (to compute $\Sigma^{-1}$) explicitly while taking the one-hot encoding columns into account (see Section A).
This approach of explicitly writing out a Gram-Schmidt orthogonalization to remove the effects of some of the controls could potentially be used whenever we have two sets of features, one of which is easier to analyze than the other.

Our focus on one-hot encodings is due to their prevalence in econometrics datasets, along with the fact that, as we will see over the next few pages, the perpendicular structure of the indicator variables makes them particularly amenable to such a divide and conquer strategy (and much harder to deal with using previous techniques).

For ease of notation, we will normalize the dummy variables so that $\left\|{u_{j}}\right\|=1$ for all $j\in[m]$.

###### Claim 4.1 .

### 4.2 First Order Effects

Claim 4.1 allows us to reduce the problem of performing a regression on $d$ continuous feature and a categorical feature with $m$ potential values into a regression on just a reaveraged version of the $d$ continuous features.
Using this reduction, we can split the effects of removing samples from the smaller linear regression into direct removal effects that change the regression fit directly by removing samples from the $d$-dimensional regression $X,Y$, and reaveraging effects that shift the regression fit by causing the expectation of the data within each bucket to change (effectively shifting the values of the retained rows of $X,Y$).

In particular, we may again write

$$ $\displaystyle\left\langle{e},{\beta-\beta_{S}}\right\rangle$ $\displaystyle=-\left\langle{e},{\hat{\Sigma}_{S}^{-1}\sum_{i\in S}\hat{X}_{i} \hat{R}_{i}}\right\rangle$ (9) $$

where $\hat{\cdot}$ is the value of $\cdot$ after the new reaveraging due to the removal of the elements in $T$.
More concretely, we have

$$ $\displaystyle\hat{X}_{i}=X_{i}-x_{b(i)}u_{b(i),i}\text{ where }x_{j}=\sum_{i^{ \prime}\in B_{j}\cap S}X_{i}u_{j,i^{\prime}}$ (10) $\displaystyle\hat{R}_{i}=R_{i}-r_{b(i)}u_{b(i),i}\text{ where }r_{j}=\sum_{i^{ \prime}\in B_{j}\cap S}R_{i}u_{j,i^{\prime}}$ $$

From here, we can derive closed form formula for the new versions of all the terms in the previous analysis.
For instance, the “first order effect” (i.e., the effect we would have seen on the regression, had $\Sigma_{S}$ been equal to $I$), which was previously given by $\text{First Order}_{\text{continuous}}=\sum_{i\in T}R_{i}\left\langle{e},{X_{i
}}\right\rangle$ is now given by

$$ $\displaystyle\text{First Order}_{\text{one-hot}}$ $\displaystyle=-\left\langle{e},{\sum_{i\in S}\hat{X}_{i}\hat{R}_{i}}\right \rangle=-\left\langle{e},{\sum_{i\in S}(X_{i}-x_{b(i)}u_{b(i),i})\hat{R}_{i}} \right\rangle=$ (11) $\displaystyle=-\left\langle{e},{\sum_{j\in[m]}\sum_{i\in B_{j}\cap S}(X_{i}-x_ {j}u_{j,i})\hat{R}_{i}}\right\rangle=$ $\displaystyle=-\left\langle{e},{\sum_{j\in[m]}\sum_{i\in B_{j}\cap S}X_{i}\hat {R}_{i}}\right\rangle-\left\langle{e},{\sum_{j\in[m]}x_{j}\underbrace{\sum_{i \in B_{j}\cap S}u_{j,i}\hat{R}_{i}}_{=0}}\right\rangle=$ $\displaystyle=-\left\langle{e},{\sum_{j\in[m]}\sum_{i\in B_{j}\cap S}X_{i}R_{i }}\right\rangle+\left\langle{e},{\sum_{j\in[m]}r_{j}\underbrace{\sum_{i\in B_{ j}\cap S}u_{j,i}X_{i}}_{=-\sum_{i\in B_{j}\cap T}X_{i}u_{j,i}}}\right\rangle=$ $\displaystyle=\underbrace{\left\langle{e},{\sum_{i\in T}X_{i}R_{i}}\right \rangle}_{\text{First Order}_{\text{continuous}}}+\underbrace{\sum_{j\in[m]} \frac{1}{\left\|{u_{j,S}}\right\|^{2}}\left(\sum_{i\in T\cap B_{j}}R_{i}u_{j,i }\right)\left(\sum_{i\in T\cap B_{j}}\left\langle{e},{X_{i}}\right\rangle u_{j ,i}\right)}_{\text{Correction Term}}$ $$

In order to bound the right-hand-side of equation (11), we consider each of the $3$ terms $\sum_{i\in T\cap B_{j}}R_{i}\left\langle{e},{X_{i}}\right\rangle$, $\sum_{i\in T\cap B_{j}}R_{i}u_{j,i}$ and $\sum_{i\in T\cap B_{j}}\left\langle{e},{X_{i}}\right\rangle u_{j,i}$ separately.
Using a greedy algorithm, we can maximize / minimize each of these terms separately as a function of $k_{j}={\left\lvert{T\cap B_{j}}\right\rvert}$.
We can then combine these into a bound of the form

$$ $\displaystyle\text{bound}^{\pm}_{j}(k_{j})=d_{j}(k_{j})+\frac{c^{\pm}_{j}(k_{j })}{\min_{S_{j}\in\binom{B_{j}}{k_{j}}}\left\{{\left\|{u_{j,S_{j}}}\right\|^{2 }}\right\}}$ $$

where

$$ $\displaystyle d_{j}(k_{j})=\max_{\begin{subarray}{c}T_{j}\subseteq B_{j}\\ {\left\lvert{T_{j}}\right\rvert}=k_{j}\end{subarray}}\sum_{i\in T_{j}}R_{i} \left\langle{e},{X_{i}}\right\rangle$ $\displaystyle c^{+}_{j}(k_{j})=\max\left\{{\left(\max_{\begin{subarray}{c}T_{j }\subseteq B_{j}\\ {\left\lvert{T_{j}}\right\rvert}=k_{j}\end{subarray}}\sum_{i\in T_{j}}R_{i}u_{ j,i}\right)\left(\max_{\begin{subarray}{c}T_{j}\subseteq B_{j}\\ {\left\lvert{T_{j}}\right\rvert}=k_{j}\end{subarray}}\sum_{i\in T_{j}}\left \langle{e},{X_{i}}\right\rangle u_{j,i}\right),\left(\min_{\begin{subarray}{c} T_{j}\subseteq B_{j}\\ {\left\lvert{T_{j}}\right\rvert}=k_{j}\end{subarray}}\sum_{i\in T_{j}}R_{i}u_{ j,i}\right)\left(\min_{\begin{subarray}{c}T_{j}\subseteq B_{j}\\ {\left\lvert{T_{j}}\right\rvert}=k_{j}\end{subarray}}\sum_{i\in T_{j}}\left \langle{e},{X_{i}}\right\rangle u_{j,i}\right)}\right\}$ $\displaystyle c^{-}_{j}(k_{j})=\max\left\{{\left(\min_{\begin{subarray}{c}T_{j }\subseteq B_{j}\\ {\left\lvert{T_{j}}\right\rvert}=k_{j}\end{subarray}}\sum_{i\in T_{j}}R_{i}u_{ j,i}\right)\left(\max_{\begin{subarray}{c}T_{j}\subseteq B_{j}\\ {\left\lvert{T_{j}}\right\rvert}=k_{j}\end{subarray}}\sum_{i\in T_{j}}\left \langle{e},{X_{i}}\right\rangle u_{j,i}\right),\left(\max_{\begin{subarray}{c} T_{j}\subseteq B_{j}\\ {\left\lvert{T_{j}}\right\rvert}=k_{j}\end{subarray}}\sum_{i\in T_{j}}R_{i}u_{ j,i}\right)\left(\min_{\begin{subarray}{c}T_{j}\subseteq B_{j}\\ {\left\lvert{T_{j}}\right\rvert}=k_{j}\end{subarray}}\sum_{i\in T_{j}}\left \langle{e},{X_{i}}\right\rangle u_{j,i}\right)}\right\}$ $$

For $k_{j}=n_{j}$, there is no reaveraging effect and $\text{bound}^{\pm}_{j}(n_{j})=d_{j}(n_{j})$.

We can then combine these bounds on the individual buckets by using a dynamic programming algorithm to solve the integer knapsack problem of

$$ $\max_{k_{1}+\cdots+k_{m}=k}\sum_{j\in[m]}\text{bound}^{\pm}_{j}(k_{j})\,.$ $$

From equation (11), it is clear that these maximizations yield upper and lower bounds on the maximal first order removal effect

$$ $\max_{k_{1}+\cdots+k_{m}=k}\sum_{j\in[m]}\text{bound}^{-}_{j}(k_{j})\leq\max_{ S\in\binom{[n]}{n-k}}\left\{{-\left\langle{e},{\sum_{i\in S}\hat{X}_{i}\hat{R} _{i}}\right\rangle}\right\}\leq\max_{k_{1}+\cdots+k_{m}=k}\sum_{j\in[m]}\text{ bound}^{+}_{j}(k_{j})\,.$ $$

The upper bound on the overall first order effect tends to be larger than the continuous first order effect but only by a $1+o(1)$ factor.
This is because $d_{j}$ tends to be the dominant effect, while $c_{j}$ are typically much smaller.

For instance, consider an unweighted regression on normally distributed samples $(X_{i},Y_{i})$, and let $\sigma_{R}=\sqrt{\frac{1}{n}\sum_{i\in[n]}R_{i}^{2}}$ denote the root-mean-square (RMS) / scale of the residuals.
With high probability, we can find removal sets that would have $d_{j}\approx\frac{k_{j}}{\sqrt{n}}\sigma_{R}\log(n/k)$ (by taking only samples which are on the $\varepsilon$ tail end of the distribution of having both large residual and large inner product with the direction of interest).
On the other hand, for sufficiently large buckets ${\left\lvert{B_{j}}\right\rvert}\gg\log(n)$, we expect to have $c_{j}=O\left(\frac{k_{j}^{2}}{{\left\lvert{B_{j}}\right\rvert}\sqrt{n}}\sigma_
{R}\right)<O\left(\frac{k_{j}}{\sqrt{n}}\sigma_{R}\right)$, so they tend to be somewhat smaller than $d_{j}$.

Moreover, by using a dynamic programming algorithm, we can enforce the constraint that the number of samples removed from each bucket $k_{j}$ has to be the same for the direct effects and for the reaveraging effect.
This constraint limits causes the contribution of the reaveraging effects to be even smaller, as they require taking many samples from the same bucket to enjoy the quadratic scaling of $k_{j}$, while the more dominant direct effects are typically optimized by selecting samples evenly from the buckets.

Figure: Algorithm 3 Dynamic Programming Algorithm for Integer Knapsack

### 4.3 High Order Terms

We continue our analysis of equation (9).
As in the continuous analysis, we have

$$ $\displaystyle\left\langle{e},{\beta-\beta_{S}}\right\rangle$ $\displaystyle=-\left\langle{e},{\hat{\Sigma}_{S}^{-1}\sum_{i\in S}\hat{X}_{i} \hat{R}_{i}}\right\rangle=$ (12) $\displaystyle=-\left\langle{e},{I\sum_{i\in S}\hat{X}_{i}\hat{R}_{i}}\right \rangle-\left\langle{e},{\left(\hat{\Sigma}_{S}^{-1}-I\right)\sum_{i\in S}\hat {X}_{i}\hat{R}_{i}}\right\rangle=$ $\displaystyle=\text{First Order}_{\text{one-hot}}+\left\langle{\hat{\Sigma}_{S }^{-1}\left(I-\hat{\Sigma}_{S}\right)e},{\sum_{i\in S}\hat{X}_{i}\hat{R}_{i}} \right\rangle\leq$ $\displaystyle\leq\text{First Order}_{\text{one-hot}}+\max\left\{{\lambda\left( \hat{\Sigma}_{S}^{-1}\right)}\right\}\times\left\|{\left(I-\hat{\Sigma}_{S} \right)e}\right\|\times\left\|{\sum_{i\in S}\hat{X}_{i}\hat{R}_{i}}\right\|$ $$

To analyze the first 2 higher order terms, we begin with an analysis of $\hat{\Sigma}_{S}$.
This analysis shows that $\hat{\Sigma}_{S}$ essentially behaves like $\Sigma_{S}$ with a minor correction for each bucket:

$$ $\displaystyle\hat{\Sigma}_{S}$ $\displaystyle=\sum_{i\in S}\hat{X}_{i}\hat{X}_{i}^{\intercal}=\sum_{j\in[m]} \sum_{i\in B_{j}\cap S}(X_{i}-x_{j}u_{j,i})\hat{X}_{i}^{\intercal}=$ (13) $\displaystyle=\sum_{j\in[m]}\sum_{i\in B_{j}\cap S}X_{i}(X_{i}-x_{j}u_{j,i})^{ \intercal}-\sum_{j\in[m]}x_{j}\underbrace{\sum_{i\in B_{j}\cap S}\hat{X}_{i}^{ \intercal}u_{j,i}}_{=0}=$ $\displaystyle=\sum_{j\in[m]}\sum_{i\in B_{j}\cap S}X_{i}X_{i}^{\intercal}-\sum _{j\in[m]}\left(\sum_{i\in B_{j}\cap T}X_{i}u_{j,i}\right)\left(\sum_{i\in B_{ j}\cap T}X_{i}\right)^{\intercal}=$ $\displaystyle=I-\Sigma_{T}-\sum_{j\in[m]}\frac{1}{\left\|{u_{j,S}}\right\|^{2} }\left(\sum_{i\in B_{j}\cap T}X_{i}u_{j,i}\right)\left(\sum_{i\in B_{j}\cap T} X_{i}u_{j,i}\right)^{\intercal}$ $$

In order to bound $\max\left\{{\lambda\left(\hat{\Sigma}_{S}^{-1}\right)}\right\}=\min\left\{{
\lambda\left(\hat{\Sigma}_{S}\right)}\right\}^{-1}$ from above, we begin by bounding $\max\left\{{\lambda\left(\Sigma_{T}\right)}\right\}$ from above (using the same MSN-bounding reductions from Section 3).
We then use a MSN-bounding algorithm and the same dynamic programming as above to bound the term in equation (14) (MSN-bounding to bound the individual terms and dynamic programming to combine them).

$$ $\displaystyle\max\left\{{\lambda\left(\sum_{j\in[m]}\frac{1}{\left\|{u_{j,S}} \right\|^{2}}\left(\sum_{i\in B_{j}\cap T}X_{i}u_{j,i}\right)\left(\sum_{i\in B _{j}\cap T}X_{i}u_{j,i}\right)^{\intercal}\right)}\right\}=$ (14) $\displaystyle\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;=\max_{k_ {1}+\cdots+k_{m}=k}\sum_{j\in[m]}\frac{1}{\min_{S_{j}}\left\{{\left\|{u_{j,S_{ j}}}\right\|^{2}}\right\}}\max_{\begin{subarray}{c}T_{j}\subseteq B_{j}\\ {\left\lvert{T_{j}}\right\rvert}=k_{j}\end{subarray}}\left\{{\left\|{\sum_{i \in T_{j}}X_{i}u_{j,i}}\right\|^{2}}\right\}$ $$

For each $j$, we can bound $\max_{\begin{subarray}{c}T_{j}\subseteq B_{j}\\
{\left\lvert{T_{j}}\right\rvert}=k_{j}\end{subarray}}\left\{{\left\|{\sum_{i
\in T_{j}}X_{i}u_{j,i}}\right\|}\right\}$ from above by performing a call to an MSN bounding algorithm.
However, we can refine the results of this MSN bounding call (especially with regards to larger values of $k_{j}\geq\frac{n_{j}}{2}$), by utilizing the fact that $\sum_{i\in B_{j}}X_{i}=0$, which implies that

$$ $\max_{\begin{subarray}{c}T_{j}\subseteq B_{j}\\ {\left\lvert{T_{j}}\right\rvert}=k_{j}\end{subarray}}\left\{{\left\|{\sum_{i \in T_{j}}X_{i}u_{j,i}}\right\|}\right\}=\max_{\begin{subarray}{c}T_{j} \subseteq B_{j}\\ {\left\lvert{T_{j}}\right\rvert}=n_{j}-k_{j}\end{subarray}}\left\{{\left\|{ \sum_{i\in T_{j}}X_{i}u_{j,i}}\right\|}\right\}\,.$ $$

Therefore, if $M_{j}$ is the result of an MSN bounding algorithm (such as RTI) on the $n_{j}\times n_{j}$ Gram matrix of the samples in the $j$th bucket, we have the bound

$$ $\max_{\begin{subarray}{c}T_{j}\subseteq B_{j}\\ {\left\lvert{T_{j}}\right\rvert}=k_{j}\end{subarray}}\left\{{\left\|{\sum_{i \in T_{j}}X_{i}u_{j,i}}\right\|}\right\}\leq\min\left\{{M_{j}(k_{j}),M_{j}(n_{ j}-k_{j})}\right\}\,.$ $$

Similarly, we may bound the other terms.
Using equation (13), we have

$$ $\displaystyle\left(I-\hat{\Sigma}_{S}\right)e$ $\displaystyle=\left(\Sigma_{T}+\sum_{j\in[m]}\left(\sum_{i\in B_{j}\cap T}X_{i }u_{j,i}\right)\left(\sum_{i\in B_{j}\cap T}X_{i}u_{j,i}\right)^{\intercal} \right)e=$ (15) $\displaystyle=\sum_{i\in T}X_{i}\left\langle{X_{i}},{e}\right\rangle+\sum_{j \in[m]}\frac{1}{\left\|{u_{j,S}}\right\|^{2}}\left(\sum_{i\in T\cap B_{j}}X_{i }u_{j,i}\right)\left(\sum_{i\in T\cap B_{j}}\left\langle{X_{i}},{e}\right \rangle u_{j,i}\right)$ $$

As usual, we bound the norm of the first term $\sum_{i\in T}X_{i}\left\langle{X_{i}},{e}\right\rangle$ using an MSN-bounding algorithm, and for each bucket, we bound $\left\|{\sum_{i\in T\cap B_{j}}X_{i}u_{j,i}}\right\|$ and ${\left\lvert{\sum_{i\in T\cap B_{j}}\left\langle{X_{i}},{e}\right\rangle u_{j,
i}}\right\rvert}$ as a function of $k_{j}={\left\lvert{T\cap B_{j}}\right\rvert}$.
We then combine these bounds using the triangle inequality, and the same symmetry and dynamic programming algorithm as above.

Finally, to bound $\left\|{\sum_{i\in S}\hat{X}_{i}\hat{R}_{i}}\right\|$, we use the same analytic techniques to write

$$ $\displaystyle\sum_{i\in S}\hat{X}_{i}\hat{R}_{i}$ $\displaystyle=\sum_{j\in[m]}\sum_{i\in B_{j}\cap S}(X_{i}-x_{j}u_{j,i})\hat{R} _{i}=\cdots$ (16) $\displaystyle\cdots=\sum_{i\in T}X_{i}R_{i}+\sum_{j\in[m]}\frac{1}{\left\|{u_{ j,S}}\right\|^{2}}\left(\sum_{i\in B_{j}\cap T}X_{i}u_{j,i}\right)\left(\sum_{ i\in B_{j}\cap T}R_{i}u_{j,i}\right)$ $$

which can be bounded by the same MSN + symmetry + dynamic programming above.

### 4.4 The OHARE Algorithm

Figure: Algorithm 4 The OHARE Algorithm

### 4.5 Supplementary Algorithms

We begin by an extension of the RTI Algorithm 2 that bounds the triangle inequality terms as a function of $k$ (number of removals) and $u$ (number of unique buckets from which samples were removed).
The basic idea is the same as with the original RTI algorithm and requires only minor adaptations to track which row-column pairs are largest in their respective buckets.

Figure: Algorithm 5 KU Triangle Inequality

We also adapt Algorithm 3 to the case where we wish to keep track of both the total number of removals and the number of unique buckets from which we remove samples.

Figure: Algorithm 6 Dynamic Programming 2D

## 5 Additional MSN-Bounding Algorithms

### 5.1 Spectral Decomposition

At a high-level, this algorithm similarly tries to bound $\left\|{\sum_{i\in T}Z_{i}}\right\|^{2}=1_{T}^{\intercal}G1_{T}$ for the same Gram matrix $G$, but here we use a spectral decomposition. We refine on the standard spectral algorithms which might output $\left\|{1_{T}}\right\|\max\lambda(G)=\sqrt{k}\max\lambda(G)$, by computing naive bounds on the inner products between $1_{G}$ and each of the top eigenvalues using a standard greedy algorithm.

Figure: Algorithm 7 Spectral Bound Algorithm

### 5.2 Greedy Lower Bound

In order to get some additional diagnostic capabilities and improve the interpretability of our results, we also implemented a simple greedy algorithm that helps us generate a *lower-bound* on an MSN problem.
Algorithm 8 initializes candidate set $T$ to contain the longest vector $v_{i}$ in our MSN instance, and greedily adds vectors to this set in an attempt to increase $\left\|{\sum_{i\in T}v_{i}}\right\|$ as fast as possible.

Figure: Algorithm 8 Greedy Algorithm for Lower Bound and Candidate Set

In real-world datasets, we found that the upper bounds obtained by Algorithms 2 and 7 tend to be very close to the greedy lower bounds on their respective MSNs.

Moreover, when analyzing the Cash Transfer study (see Section 6.2), we noticed that our analysis behaved poorly in some regimes.
Using this greedy lower bound technique, we were able to identify that this behaviour was caused by a small number of households ($<0.5\%$) dominating the variance in land-ownership ($>80\%$), which was included in the study in linear scale.
Converting the land-ownership to a logarithmic scale resolved this issue (allowing us to certify robustness to removing many more samples).

## 6 Applied Experiments

##### AMIP

The AMIP algorithm by Broderick et al. Broderick et al. () can be used to estimate the robustness of a linear regression in one of two ways.
The fastest option is to compute the AMIP gradients / influence scores, and estimate that the number of samples one must remove in order to change the fit by a certain amount by taking a cummulative sum over the sorted gradients.

However, this method both typically tends to overestimate the robustness of a dataset, and despite this does not even provide a formal upper bound on $k_{\textup{sign}}$ or $k_{2\sigma}$.
Instead, for our experiments we sort the samples from most influential to least influential and remove them one at a time until the fit crosses the threshold we are considering.

### 6.1 How Much Should We Trust the Dictator’s GDP Growth Estimates?

#### 6.1.1 Paper Overview

Reliable estimates of GDP figures are crucial for analysts to assess the performance of an economy.
However, leaders often have a variety of political and financial incentives manipulate GDP figures in order to improve their perception.
Therefore, economists often use proxies to obtain independent estimates of these figures that may be harder to manipulate.
One such method which has gained a lot of attention in recent years is to simply measure the amount of light emitted from a region at night (nightlights or NTL), as observed by satellite imaging Henderson et al. ().

Martinez Martinez () uses this proxy in conjunction with several well-known measures for the democracy of a country (such as the freedom-in-the-world or FiW metric), in an effort to find evidence of GDP figure manipulation in autocratic regimes.
One methods used by Martinez is to measure the “autocracy gradient in the NTL elasticity of GDP”.
In practice, this translates to running a regression of the form

$$ $\ln(\text{GDP})_{i,t}=\mu_{i}+\delta_{t}+\phi_{0}\ln(\text{NTL})_{i,t}+\phi_{1 }\text{FiW}_{i,t}+\phi_{2}\text{FiW}^{2}_{i,t}+\phi_{3}[\ln(\text{NTL})\times \text{FiW}]_{i,t}+\varepsilon_{i,t}$ (17) $$

where $i$ represents the index of a country, $t$ the year from which this sample was taken, $\mu_{i}$ and $\delta_{t}$ control for country-specific or time-specific effects, $\phi_{0},\phi_{1},\phi_{2}$ control for the direct correlation between nightlights and GDP, as well as for any $\leq 2$nd order dependence between GDP and democracy.
$\phi_{3}$ is the main effect we wish to observe, and a positive value of $\phi_{3}$ could be explained by autocratic regimes being more prone to embellishing their GDP figures.

Martinez reports a statistically significant positive value of $\phi_{3}$, as well as additional evidence for GDP figure manipulation in autocratic regimes (such as a larger difference between NTL-based estimates and reported GDP in years leading up to IMF evaluations of autocratic regimes).
Martinez hypothesises that the separation of powers and cross-examination of figures by opposition contribute to making GDP figures more reliable in democratic countries.

#### 6.1.2 Robustness

Martinez uses the AMIP tool Broderick et al. () to assess the robustness of the regression (17).
AMIP finds a set of 136 samples whose removal flips the sign of the OLS fit on the $[\ln(\text{NTL})\times\text{FiW}]_{i,t}$ term.
Kusching et al. Kuschnig et al. () improve this upper bound on the stability by finding a set of 110 samples that flip the sign of the fit parameter.

Running our continuous regression toolkit on the Martinez dataset can only certify robustness to the removal of at most $7$ samples, as the dataset contains a one-hot encoding of the country and it contains only $8$ samples from Monetenegro.
To overcome this, we apply our one-hot aware algorithms to the same data, and are able to certify the sign of this parameter in Martinez’ regression requires at least $k_{\textup{sign}}\geq 29$ to overturn.

### 6.2 Indirect Effects of an Aid Program: How Do Cash Transfers Affect Ineligibles’ Consumption?

#### 6.2.1 Paper Overview

Angelucci and De Giorgi Angelucci and De Giorgi () study the indirect effect of the Progresa welfare program in Mexico.
This program gave financial aid to eligible households within “treated” villages, and did not give aid to ineligible households or households from “untreated” villages.

Angelucci and De Giorgi then track spending patterns of eligible and ineligible households from both treated and untreated villages.
They then use linear regressions to estimate the treatment effects on both eligible and ineligible households when controlling on various other features.

#### 6.2.2 Minor Difficulties in Using the Data

We found two different regression formula that have been attributed to this paper.
The formula use slightly different control, with one option controlling for: head of household age, sex, literacy and education level, household poverty index and amount of land owned, local poverty index, number of households in the village and “average shock” (to the best of our knowledge, this was the regression used in Angelucci and De Giorgi’s original paper).
The regression used in the later paper by Broderick et al. Broderick et al. () also controls for head of household marital status and which region the samples came from (using a one-hot encoding).
Ultimately, these additional controls had very little effect on the fit parameters or their error bars, and we chose to focus on the latter to maintain consistency with previous benchmarks.

A more significant issue is that the dataset also contained many clear outliers.
For instance, almost all the entries of the column for head of household sex were either “hombre” or “mujer” with a very small fraction (roughly 16 out of 59455 samples) having the value $9.0$, and many columns had entries of “nr” (presumably no response).

We removed these “bad samples” for the sake of our analysis (both to avoid regressing over a clearly problematic dataset, and also because not doing so would cause the regression to include nearly empty categories in several different one-hot encodings which is beyond the scope of the OHARE algorithm).
This had only a minor effect on the regression results.

**Table 2: Number and values of outliers we removed from each of the covariates in Angelucci and De Giorgi’s study Angelucci and De Giorgi ().**
| Column | Values Removed | Samples Removed |
| --- | --- | --- |
| hhhsex | 9.0, nr | 16 |
| hhhalpha | nr | 21 |
| hhhspouse | nr, 2.0 | 546 |
| p16 | nr | 146 |
| hhhage | 97 y más, no sabe, nr | 205 |
| Total |  | 921 |

#### 6.2.3 Robustness

We ran our tools on the $6$ regressions from Angelucci and De Giorgi’s paper (treatment effects on eligible / ineligible over $3$ periods).
Our robustness lower bounds nearly match the AMIP upper bound for the ineligible studies, but left room for improvement on the eligible regressions.

To find the root of this gap, we first noted that the main component of our bound that failed was that after removing $\approx 50$ samples, our bound on $\max\left\{{\lambda(\Sigma_{S}^{-1})}\right\}$ becomes very large.
Using Algorithm 8, we were able to find the samples responsible for these loose bounds: $\approx 50$ households (out of $>10000$) that account for more than $80\%$ of the variance in the amount of land owned by each household.

When one of the columns of a linear regression is heavy-tailed, it is not uncommon to replace this column with its logarithm.
This presented a slight challenge in this case, as many of the household in the eligible studies owned no land at all.
To overcome this, we replace the hectacres (amount of land owned) column with $\log\left(\text{hectacres}+\text{median hectacres}\right)$.
We then reran the regression and robustness analysis to see the effects of this change to produce the results in Table 1 and Figure 2.

### 6.3 The Oregon Health Insurance Experiment: Evidence from the First Year

#### 6.3.1 Paper Overview

The Oregon Health Insurance Experiment provides a unique opportunity to assess the effects of expanding access to public health insurance on low-income adults using a randomized controlled design. In 2008, Oregon implemented a lottery to select uninsured low-income adults to apply for Medicaid. This random assignment allows researchers to compare the outcomes of the treatment group (those selected by the lottery) with the control group (those not selected).

Finkelstein et al. Finkelstein et al. () analyze data from the first year after the lottery to evaluate the impacts on health care utilization, financial strain, and health outcomes. The study finds that individuals in the treatment group were about 25 percentage points more likely to have health insurance compared to the control group. The results show that the treatment group experienced higher health care utilization, including increased primary and preventive care visits, and hospitalizations. They also faced lower out-of-pocket medical expenditures and medical debt, evidenced by fewer bills sent to collection agencies.

Moreover, the treatment group reported better physical and mental health than the control group. The authors suggest that the increase in health care utilization due to insurance coverage led to improved health outcomes and reduced financial strain. These findings provide significant evidence on the benefits of expanding Medicaid coverage to low-income populations.

#### 6.3.2 Regression Analysis

The regression analysis conducted by Finkelstein et al. utilized instrumental variable (IV) regression to estimate the impact of Medicaid coverage on various health and financial outcomes. They used a treatment variable as the instrument and several dummy variable controls corresponding to different medical metrics (e.g., “Not bad days physical”).

IV regression is particularly useful in this context as it helps address potential endogeneity issues, where the treatment (Medicaid enrollment) might be correlated with unobserved factors affecting the outcomes. By using the lottery selection as an instrument for Medicaid enrollment, they ensured a more accurate estimation of the causal effect.

The IV regression can be computed as the ratio of two ordinary least squares (OLS) regressions: one estimating the relationship between the instrument (lottery selection) and the endogenous variable (Medicaid enrollment), and the other estimating the relationship between the instrument and the outcome (health or financial metrics). This implies that if both OLS regressions are robust to sign changes, so is their ratio. In this study, the correlation between the instrument and the endogenous variable was very strong and robust, making the primary robustness concern the OLS regression between the endogenous variable and the outcome.

Finkelstein et al. used weighted OLS regressions with several dummy variables, some of which had very sparse categories. This sparsity caused singularities after 13-14 removals, making it difficult to certify robustness with the ACRE algorithm. Therefore, the OHARE algorithm was used to certify the robustness of these regressions.

| Outcome | Type | ACRE | OHARE | AMIP | KZC21 |
| --- | --- | --- | --- | --- | --- |
| Health genflip | Instrument vs Endogenous | 14 | 752 | $\geq 10\%$ | $\geq 10\%$ |
| Instrument vs Outcome | 14 | 77 | 257 | 257 |  |
| Health notpoor | Instrument vs Endogenous | 14 | 752 | $\geq 10\%$ | $\geq 10\%$ |
| Instrument vs Outcome | 14 | 40 | 149 | 149 |  |
| Health change flip | Instrument vs Endogenous | 14 | 755 | $\geq 10\%$ | $\geq 10\%$ |
| Instrument vs Outcome | 14 | 52 | 184 | 184 |  |
| Not bad days total | Instrument vs Endogenous | 13 | 685 | $\geq 10\%$ | $\geq 10\%$ |
| Instrument vs Outcome | 13 | 21 | 72 | 72 |  |
| Not bad days physical | Instrument vs Endogenous | 13 | 669 | $\geq 10\%$ | $\geq 10\%$ |
| Instrument vs Outcome | 13 | 25 | 84 | 84 |  |
| Not bad days mental | Instrument vs Endogenous | 13 | 676 | $\geq 10\%$ | $\geq 10\%$ |
| Instrument vs Outcome | 13 | 31 | 118 | 118 |  |
| Nodep Screen | Instrument vs Endogenous | 14 | 742 | $\geq 10\%$ | $\geq 10\%$ |
| Instrument vs Outcome | 14 | 32 | 116 | 116 |  |

## 7 Tightness of ACRE

In this section, we will prove that the ACRE algorithm produces tight bounds on “well-behaved” distributions, proving Theorem 1.2.

### 7.1 Preliminaries

Throughout this section, we use $\widetilde{O}/\widetilde{\Theta}/\widetilde{\Omega}$ to denote big-O statements that hold up to a factor of $\textup{polylog}(n)$.
We will make no attempt to optimize the $\textup{polylog}(n)$ factors in this analysis.

Moreover, we say that a term $\eta$ is *negligible* if $\eta^{-1}=n^{\omega(1)}$ is superpolynomial in $n$.
Finally, we say an event happens with very high probability if it happens with probability $1-\eta$ for negligible $\eta$.

### 7.2 Main Result

Our main goal for this section will be to prove that ACRE produces good bounds with high probability when the regression data is drawn from a well behaved distribution:

###### Theorem 7.1 ( ACRE Bounds are Tight on Well-behaved Data) .

Let $X\in\mathbb{R}^{n\times d},Y\in\mathbb{R}^{n}$ be a linear regression problem such that the covariates (i.e., the rows of $X$) are drawn iid from a well-behaved distribution $X_{i}\sim\mathcal{X}$ and the outcomes $Y$ are drawn iid from $Y\sim\mathcal{N}\left(X\beta_{\textup{gt}},I_{n}\right)$.

Then, for any axis $e\in\mathbb{S}^{d-1}$, with very high probability, the upper and lower bounds produced by ACRE on this regression are close to tight

$$ $\frac{U_{k}}{L_{k}}=1+\widetilde{O}\left(\frac{d+k\sqrt{d}}{n}\right)$ $$

for all $k<k_{\textup{threshold}}$, where $k_{\textup{threshold}}=\widetilde{\Theta}\left(\min\left\{\frac{n}{\sqrt{d}},
\frac{n^{2}}{d^{2}}\right\}\right)$

In particular, if the samples $X_{i},Y_{i}$ are drawn iid from some normal distribution $\mathcal{N}(0,\Sigma^{\prime})$ (for some covariance $\Sigma^{\prime}\in\mathbb{R}^{(d+1)\times(d+1)}$), or if the covariates $X_{i}$ are drawn iid from the hypercube or unit sphere and the target variable $Y_{i}\sim\left\langle{\beta_{\textup{gt}}},{X_{i}}\right\rangle+\mathcal{N}(0,
\sigma_{R})$ are drawn iid from a normal distribution around some ground truth model, then Theorem 7.1 holds for them.
Therefore, Theorem 1.2 follows from Theorem 7.1.

Our goal for the rest of this section will be to prove Theorem 7.1.

###### Theorem 7.1 ( ACRE Bounds are Tight on Well-behaved Data) .

### 7.3 Proof Sketch

First, we note that if $n=\widetilde{O}(d)$, Theorem 7.1 holds vacuously, as $k_{\textup{threshold}}<1$.
Therefore, we limit our analysis to the cases where $n>n_{\textup{threshold}}$ for some $n_{\textup{threshold}}=d\times\textup{polylog}(n)$ (the exact power of this polylogarithmic factor will depend on the specific constant $C$ in the exponential decay assumption in Definition 1.

To prove Theorem 7.1, we first define a set of condition under which we can prove that ACRE will produce good bounds:

###### Definition 2 (ACRE-friendly) .

Let $X,Y$ be the covariates and target variable of a regression. Let $\Sigma=X^{\intercal}X$ denote the unnormalized empirical covariance of the covariates, and let $R=Y-X\Sigma^{-1}X^{T}Y$ denote the residuals.

We say that this regression is *ACRE-friendly* for direction $e\in\mathbb{R}^{d}$ with $k$ removals and parameters $P_{1},P_{2},P_{3},P_{4},P_{5}\geq 0$ if

- 1.
The covariates are bounded in Mahalanobis distance $\max_{i\in[n]}X_{i}^{\intercal}\Sigma^{-1}X_{i}\leq P_{1}\frac{d}{n}$.
- 2.
The inner products of samples are bounded $\max_{i\neq j\in[n]}X_{i}^{\intercal}\Sigma^{-1}X_{j}\leq P_{2}\frac{\sqrt{d}}
{n}$.
- 3.
The residuals are bounded $\max_{i\in[n]}R_{i}^{2}\leq P_{3}\sigma_{R}^{2}$, where $\sigma_{R}=\sqrt{\frac{1}{n}\sum_{i\in[n]}R_{i}^{2}}$.
- 4.
The inner products between the covariates and the axis of interest are bounded $\max_{i\in[n]}{\left(e^{\intercal}\Sigma^{-1}X_{i}\right)^{2}}\leq\frac{P_{4}}
{n}\sum_{i\in[n]}{\left(e^{\intercal}\Sigma^{-1}X_{i}\right)^{2}}=\frac{P_{4}}
{n}e^{\intercal}\Sigma^{-1}e$.
- 5.
Let $\alpha_{i}=e^{\intercal}\Sigma^{-1}X_{i}R_{i}$ be the AMIP influence scores. We require that the sum over the $k$ largest influence scores is at least
$a_{k}=\max_{T\in\binom{[n]}{k}}\left\{{\sum_{i\in T}\alpha_{i}}\right\}\geq
\frac{1}{P_{5}}\sigma_{Z}\sigma_{R}k$
where $\sigma_{Z}=\sqrt{\frac{1}{n}\sum_{i\in[n]}Z_{i}^{2}}$, for $Z_{i}=e^{\intercal}\Sigma^{-1}X_{i}$.

When $P_{1},\ldots,P_{5}$ are at most polylogarithmic in $n$, we say that the regression is ACRE-friendly.

Our proof of Theorem 7.1 will have two main components.
The first and smaller portion of the proof will be to show that the bounds produced by ACRE are close to tight when a regression is ACRE-friendly.
It should not be surprising that Definition 2 is sufficient condition for producing good bounds with ACRE, since ACRE essentially checks a slightly more robust version of these conditions.

###### Claim 7.2 (ACRE bounds are nearly tight on ACRE-friendly regressions) .

Let $X,Y,e$ be an ACRE-friendly regression with parameters $P_{1},\ldots,P_{5}$ for all $k\leq k_{0}$.
Then, there exists $k_{\textup{threshold}}=\Theta\left(\nicefrac{{\min\left\{\frac{n}{\sqrt{d}},
\frac{n^{2}}{d^{2}},k_{0}\right\}}}{{\textnormal{poly}\left(P_{1},\ldots,P_{5}
\right)}}\right)$ such that for all $k\leq k_{\textup{threshold}}$, the bounds $L_{k},U_{k}$ produced by ACRE satisfy

$$ $\frac{U_{k}}{L_{k}}\leq 1+O\left(\frac{d+k\sqrt{d}}{n}\times\textnormal{poly} \left(P_{1},\ldots,P_{5}\right)\right)$ $$

The second and longer portion of our proof will be devoted to showing that well-behaved distributions yield ACRE-friendly regressions with high probability.

###### Claim 7.3 (Well-behaved distributions yield ACRE-friendly regressions with high probability) .

Let $n,d$ be as above, and let $X,Y$ be as in Theorem 7.1.

Then, for any axis $e\in\mathbb{S}^{d-1}$, with very high probability, the regression $X,Y$ is ACRE-friendly for all $k\leq k_{0}$, for $k_{0}=\widetilde{\Omega}(n)$.

Combined, Claims 7.2 and 7.3 yield Theorem 7.1.

###### Definition 2 (ACRE-friendly) .

###### Claim 7.2 (ACRE bounds are nearly tight on ACRE-friendly regressions) .

###### Claim 7.3 (Well-behaved distributions yield ACRE-friendly regressions with high probability) .

### 7.4 Symmetries

Before proceeding to the proofs of Claims 7.2 and 7.3, we note that our definitions of well-behaved distributions and of ACRE-friendly regressions are normalized in a way that permits some symmetries.
Indeed, it is easy to see that for any invertible matrix $L\in\mathbb{R}^{d\times d}$ and positive scalars $\alpha_{e},\alpha_{y}\in\mathbb{R}^{+}$, a regression $X,Y,e$ is well-behaved / ACRE-friendly if and only if the regression $\widetilde{X}=XL,\widetilde{Y}=\alpha_{y}Y,\widetilde{e}=\alpha_{e}Le$ is well-behaved / ACRE-friendly.
Moreover, such a reparametrizations also has no effect on the robustness of the regression, so we may apply this symmetry on the input and output of Theorem 7.1 or Claims 7.2 and 7.3.

For Theorem 7.1 or Claim 7.3, we may apply this symmetry with $L=\Sigma^{-1/2}$, where $\Sigma$ is the *ground truth* covariance of the distribution $\mathcal{X}$, resulting in a distribution with covariance identity.
Note that this renormalization is *not* the same as the renormalization process in the ACRE algorithm.
In the ACRE algorithm, we renormalize by the empirical covariance, whereas the renormalization above is by the ground-truth covariance.
One of the steps of our analysis will be to show that when the regression is drawn from a well-behaved distribution, the two are close to one another with high probability, but this is not immediate.

For the proof of Claim 7.2, we renormalize the samples according to their empirical covariance.
This would result in a set of samples $X$ such that $\Sigma^{*}=X^{\intercal}X=I$.
Moreover, we use the $\alpha_{Y}$ symmetry to ensure that the residuals to have standard deviation $\sigma_{R}=1$ and the $\alpha_{e}$ symmetry to ensure that $e$ has norm $1$.

### 7.5 Proof of Claim 7.2

After this renormalization step above, the conditions on the normalized regression simplify to:

- 1.
The covariates are bounded in $\ell_{2}$ norm $\max_{i\in[n]}\left\|{X_{i}}\right\|^{2}\leq P_{1}\frac{d}{n}$
- 2.
The inner products of samples are bounded $\max_{i\neq j\in[n]}\left\langle{X_{i}},{X_{j}}\right\rangle\leq P_{2}\frac{
\sqrt{d}}{n}$
- 3.
The residuals are bounded $\max_{i\in[n]}R_{i}^{2}\leq P_{3}$.
- 4.
The inner products between the covariates and the axis of interest are bounded $\max_{i\in[n]}{\left\langle{e},{X_{i}}\right\rangle^{2}}\leq\frac{P_{4}}{n}$
- 5.
Let $\alpha_{i}=\left\langle{e},{X_{i}}\right\rangle R_{i}$ be the AMIP influence scores. We require that the sum over the $k$ largest influence scores is at least
$a_{k}=\max_{T\in\binom{[n]}{k}}\left\{{\sum_{i\in T}\alpha_{i}}\right\}\geq
\frac{1}{P_{5}}\times\frac{k}{\sqrt{n}}$

We prove Claim 7.2, working with the normalized regression.
Recall that ACRE produces its bounds by combining RTI bounds on the following Gram matrices:

- 1.
$G_{X\otimes X}$ whose entries are the squared entries of the Gram matrix $G_{X}$ of the original covariates. By our assumptions on the maximal norm and inner products of covariates, we know that the diagonal entries of this matrix are bounded by $\frac{d^{2}}{n^{2}}P_{1}^{2}$ and its off-diagonal entries are bounded by $\frac{d}{n^{2}}P_{2}^{2}$.
- 2.
$G_{XR}$ whose entries are inner products between covariates multiplied by the product of two residuals. Therefore, the diagonal entries of this matrix are bounded (in absolute value) by $\frac{d}{n}P_{1}P_{3}$ and its off-diagonal entries by $\frac{\sqrt{d}}{n}P_{2}P_{3}$.
- 3.
and $G_{XZ}$ whose entries are inner products between covariates, rescaled by the product of their weights on the axis of interest $e$. Therefore, its diagonal entries are bounded (in absolute value) by $\frac{d}{n^{2}}P_{1}P_{4}$ and its off-diagonal entries by $\frac{\sqrt{d}}{n^{2}}P_{2}P_{4}$.

Note that the output of the RTI algorithm is an upper bound on the resulting MSN problem.
The output of the RTI algorithm squared is always equal to the sum of $k$ diagonal entries of the Gram matrix plus fewer than $k^{2}$ off-diagonal entries.
Therefore, we have:

$$ $\text{RTI Bound}\leq\sqrt{k\times\text{Largest Diagonal Entry}+k^{2}\times \text{Largest Off-Diagonal Entry}}$ (18) $$

##### The X ⊗ X tensor-product 𝑋 𝑋 X\otimes X italic_X ⊗ italic_X Term

Let $M_{X\otimes X}$ denote the MSN-bound obtained by running the RTI algorithm on the $G_{X\otimes X}$.
Combining equation (18) with our knowledge of the diagonal and off-diagonal entries of $G_{X\otimes X}$, we clearly have:

$$ $M_{X\otimes X}\leq\sqrt{k\frac{d^{2}}{n^{2}}P_{1}^{2}+k^{2}\frac{d}{n^{2}}P_{2 }^{2}}$ $$

Therefore, for all $k\leq k_{\text{threshold}}=\min\left\{{\frac{n}{\sqrt{d}P_{2}2},\frac{n^{2}}{d
^{2}P_{1}^{2}4}}\right\}$, $M_{X\otimes X}$ is at most

$$ $M_{X\otimes X}\leq\sqrt{\frac{1}{2}}<1$ $$

##### The X ⁢ R 𝑋 𝑅 XR italic_X italic_R and X ⁢ Z 𝑋 𝑍 XZ italic_X italic_Z Terms

Similarly, let $M_{XR}$ and $M_{XZ}$ denote the MSN bounds obtained by the RTI algorithm on $G_{XR}$ and $G_{XZ}$ respectively.
As before, we combine equation (18) with our bounds on the diagonal and off-diagonal terms of $G_{XR}$ and $G_{XZ}$ to show that

$$ $M_{XR}\leq\sqrt{k\frac{d}{n}P_{1}P_{3}+k^{2}\frac{\sqrt{d}}{n}P_{2}P_{3}}\leq \sqrt{k\frac{d}{n}P_{1}P_{3}}+\sqrt{k^{2}\frac{\sqrt{d}}{n}P_{2}P_{3}}$ $$

and

$$ $M_{XZ}\leq\sqrt{k\frac{d}{n^{2}}P_{1}P_{4}+k^{2}\frac{\sqrt{d}}{n^{2}}P_{2}P_{ 4}}\leq\sqrt{k\frac{d}{n^{2}}P_{1}P_{4}}+\sqrt{k^{2}\frac{\sqrt{d}}{n^{2}}P_{2 }P_{4}}$ $$

##### Wrapping Up

Therefore, for all $k\leq k_{\textup{threshold}}$, we have

$$ $b_{k}=\frac{1}{1-M_{X\otimes X}}M_{XR}M_{XZ}\leq\frac{\sqrt{2P_{3}P_{4}}}{ \sqrt{2}-1}\times\left(\frac{kd}{n^{3/2}}P_{1}+2\frac{k^{3/2}d^{3/4}}{n^{3/2}} \sqrt{P_{1}P_{2}}+\frac{k^{2}\sqrt{d}}{n^{3/2}}P_{2}\right)=O\left(\frac{kd+k^ {2}\sqrt{d}}{n^{3/2}}\times\textnormal{poly}\left(P_{1},\ldots,P_{4}\right)\right)$ $$

Finally, applying the 5th condition, we have

$$ $\frac{U_{k}}{L_{k}}=\frac{a_{k}+b_{k}}{a_{k}-b_{k}}=1+O\left(\frac{b_{k}}{a_{k }}\right)=1+O\left(\frac{d+k\sqrt{d}}{n}\times\textnormal{poly}\left(P_{1}, \ldots,P_{5}\right)\right)$ $$

completing the proof of Claim 7.2.

### 7.6 Proof of Claim 7.3

We now move on to the main portion of the proof which will be devoted to showing that when the covariates $X_{i}\sim\mathcal{X}$ are drawn from a well-behaved distribution and the target variable $Y_{i}\sim N(X_{i}^{\intercal}{\beta}_{\textup{gt}},1)$ are drawn from a normal distribution around some ground truth linear model, the resulting regression is ACRE-friendly with very high probability (i.e., to prove Claim 7.3).
Recall that as we showed in Section 7.4, it suffices to prove this result for the case where the ground-truth covariance of the distribution $\mathcal{X}$ is equal to identity.

##### Main Challenge and Proof Strategy

The main challenge will be that the requirements of Definition 2 (ACRE-friendly) require statements like $v^{\intercal}\Sigma^{-1}w=\textup{small}$, where $v$ and $w$ are related to the samples of the regression (e.g., $v,w$ might be two samples $X_{i},X_{j}$) and from here on out, $\Sigma=X^{\intercal}X$ denotes the unnormalized empirical covariance (recall that we normalized the ground truth covariance to be the identity).

Because of our assumption that $n>d\times\textup{polylog}(n)$, we can use matrix Bernstein to prove that the empirical covariance is close to its expectation $\Sigma\approx nI$, but even for moderate dimensions $d$, our bounds on the overall error of this approximation (i.e., $\left\|{\Sigma-nI}\right\|$) would not be strong enough to prove Claim 7.3.
For any fixed $v,w\in\mathbb{R}^{d}$ that do not depend on the sample distributions, we can easily prove sufficiently strong concentration bounds on $v^{\intercal}\Sigma^{-1}w$.
However, the $v,w$ pairs for which we will need to prove our concentration bounds *do* depend on the samples, creating the potential for an alignment between the large eigenvectors of $\Sigma-nI$ and this pair.

Our proof strategy will be to expand $\Sigma^{-1}$ into terms that depend on $v,w$ and terms that do not.
The terms that do not depend on $v,w$ can be tightly bounded using simple concentration bounds, and the terms that do depend on $v,w$ will be so small to begin with that we can bound them very loosely using the matrix Bernstein inequality (see Lemma 7.4).

#### 7.6.1 Matrix Bernstein Inequality

Throughout our analysis we will often make use of the matrix Bernstein inequality:

###### Lemma 7.4 (Matrix Bernstein Tropp et al. ( 2015 ) ) .

Let $Z_{1},\ldots,Z_{n}$ be independently-distributed random symmetric $d\times d$ matrices with mean ${\mathop{\mathbb{E}}{\left[{Z_{i}}\right]}}=0$.
Moreover, assume that $\left\|{Z_{i}}\right\|\leq L$ and $\left\|{\sum_{i\in[n]}{\mathop{\mathbb{E}}{\left[{Z_{i}^{2}}\right]}}}\right\|
\leq\sigma^{2}$, then

$$ ${\Pr\left[{\left\|{\sum_{i\in[n]}Z_{i}}\right\|\geq t}\right]}\leq 2d\cdot\exp \left(\frac{-t^{2}/2}{\sigma^{2}+Lt/3}\right)$ $$

In many cases however, our input will not directly fit the assumptions of this inequality and instead of a statement of the form $\left\|{Z_{i}}\right\|\leq L$, we will only have a statement of the form ${\Pr\left[{\left\|{Z_{i}}\right\|\geq L}\right]}\leq\delta$ for some negligible probability $\delta$.
Therefore it will be useful for us to have an extension of this theorem to cases when the assumption fails, but with a very small probability:

###### Lemma 7.5 (Approximate Matrix Bernstein) .

Let $Z_{1},\ldots,Z_{n}$ be independently-distributed random symmetric $d\times d$ matrices with mean ${\mathop{\mathbb{E}}{\left[{Z_{i}}\right]}}=0$.
Moreover, assume that for all $i$, ${\Pr\left[{\left\|{Z_{i}}\right\|\geq L}\right]}\leq\delta$ and $\left\|{{\mathop{\mathbb{E}}{\left[{Z_{i}^{2}}\right]}}}\right\|\leq\tau^{2}$, and that $\left\|{\sum_{i\in[n]}{\mathop{\mathbb{E}}{\left[{Z_{i}^{2}}\right]}}}\right\|
\leq\sigma^{2}$, then

$$ ${\Pr\left[{\left\|{\sum_{i\in[n]}Z_{i}}\right\|\geq t+n\tau\sqrt{\delta}} \right]}\leq 2d\cdot\exp\left(\frac{-t^{2}/2}{\sigma^{2}+Lt/3}\right)+n\delta$ $$

###### Proof of Lemma 7.5 .

Let $\mathcal{S}^{d}=\left\{Z\in\mathbb{R}^{d\times d}\middle|Z=Z^{T}\right\}$ denote the space of symmetric real valued $d\times d$ matrices.

Consider the random variable $\zeta_{i}=Z_{i}\times 1_{\left\|{Z_{i}}\right\|\leq L}$.
At first it might seem like we can directly apply the matrix Bernstein inequality to these new variables, and then use the fact that $Z_{i}=\zeta_{i}$ for all $i$ w.p. $\geq 1-n\delta$.

The issue is that we no longer know that these $\zeta_{i}$ maintain the other assumptions of the matrix Bernstein inequality.
In particular, we will bound $\left\|{{\mathop{\mathbb{E}}{\left[{\zeta_{i}}\right]}}}\right\|$ and $\left\|{\sum_{i\in[n]}{\mathop{\mathbb{E}}{\left[{{\zeta_{i}}^{2}}\right]}}}\right\|$ from above.

Let $f_{i}$ denote the probability density function of $Z_{i}$, and consider

$$ ${\mathop{\mathbb{E}}{\left[{\zeta_{i}}\right]}}={\mathop{\mathbb{E}}{\left[{ \zeta_{i}}\right]}}-{\mathop{\mathbb{E}}{\left[{Z_{i}}\right]}}=\int_{\left\|{ Z}\right\|>L}Zf_{i}(Z)$ $$

Let $v,w\in\mathbb{S}^{d-1}$ be any two points on the unit sphere.
Using the CS inequality, we have

$$ $v^{\intercal}{\mathop{\mathbb{E}}{\left[{\zeta_{i}}\right]}}w=\int_{\left\|{Z} \right\|>L}v^{\intercal}Zwf_{i}(Z)=\int_{Z\in\mathcal{S}^{d}}1_{\left\|{Z} \right\|>L}v^{\intercal}Zwf_{i}(Z)\leq\sqrt{\int_{Z\in\mathcal{S}^{d}}1_{\left \|{Z}\right\|>L}f_{i}(Z)}\times\sqrt{\int_{Z\in\mathcal{S}^{d}}\left(v^{ \intercal}Zw\right)^{2}f_{i}(Z)}$ $$

The first term in the RHS is bounded by

$$ $\sqrt{\int_{Z\in\mathcal{S}^{d}}1_{\left\|{Z}\right\|>L}f_{i}(Z)}=\sqrt{{\Pr \left[{\left\|{Z_{i}}\right\|\geq L}\right]}}\leq\sqrt{\delta}$ $$

For the latter term, we use the CS inequality again:

$$ $\sqrt{\int_{Z\in\mathcal{S}^{d}}\left(v^{\intercal}Zw\right)^{2}f_{i}(Z)}\leq \sqrt{\int_{Z\in\mathcal{S}^{d}}\left\|{Zv}\right\|^{2}\left\|{w}\right\|^{2}f _{i}(Z)}=\sqrt{\int_{Z\in\mathcal{S}^{d}}v^{\intercal}Z^{2}vf_{i}(Z)}=\sqrt{v^ {\intercal}{\mathop{\mathbb{E}}{\left[{Z_{i}^{2}}\right]}}v}\leq\tau$ $$

Therefore

$$ $\left\|{{\mathop{\mathbb{E}}{\left[{\zeta_{i}}\right]}}}\right\|\leq\tau\sqrt{\delta}$ $$

Finally, we bound the change in the variance of these variables.

$$ $\textup{Cov}\left(\zeta_{i}\right)={\mathop{\mathbb{E}}{\left[{\zeta_{i}^{2}} \right]}}-{\mathop{\mathbb{E}}{\left[{\zeta_{i}}\right]}}^{2}\preceq{\mathop{ \mathbb{E}}{\left[{\zeta_{i}^{2}}\right]}}\preceq{\mathop{\mathbb{E}}{\left[{Z _{i}^{2}}\right]}}\Rightarrow\left\|{\sum_{i\in[n]}\textup{Cov}\left(\zeta_{i} \right)}\right\|\leq\left\|{\sum_{i\in[n]}{\mathop{\mathbb{E}}{\left[{Z_{i}^{2 }}\right]}}}\right\|\leq\sigma^{2}$ $$

Therefore, we may apply the matrix Bernstein inequality (Lemma 7.4) on the matrices $z_{i}=\zeta_{i}-{\mathop{\mathbb{E}}{\left[{\zeta_{i}}\right]}}$ to obtain the bound

$$ ${\Pr\left[{\left\|{\sum_{i\in[n]}\zeta_{i}}\right\|\geq t+n\tau\sqrt{\delta}} \right]}\leq{\Pr\left[{\left\|{\sum_{i\in[n]}\zeta_{i}-\sum_{i\in[n]}{\mathop{ \mathbb{E}}{\left[{\zeta_{i}}\right]}}}\right\|\geq t}\right]}\leq 2d\cdot\exp \left(\frac{-t^{2}/2}{\sigma^{2}+Lt/3}\right)\,,$ $$

yielding the main claim.

∎

###### Lemma 7.4 (Matrix Bernstein Tropp et al. ( 2015 ) ) .

###### Lemma 7.5 (Approximate Matrix Bernstein) .

###### Proof of Lemma 7.5 .

#### 7.6.2 A Useful Corollary of Lemma 7.5

Recall that we defined a “well-behaved distribution” to be a distribution whose tails decay rapidly when projected on any direction.
It will be beneficial to our use-case to work with a more relaxed condition on the decay of these tails, by allowing an additional $\textnormal{poly}(n)$ factor.
Note that any well-behaved distribution is clearly also almost well-behaved.

###### Definition 3 .

We say that a mean-zero distribution $\mathcal{X}$ on $\mathbb{R}^{d}$ is *almost well-behaved* with respect to the scaling parameter $n$, if it has exponentially decaying tails in the sense that

$$ $\exists C>0\;\;\forall v\in\mathbb{S}^{d-1},t>0\;\;\;{\Pr_{X\sim\mathcal{X}} \left[{{\left\lvert{\left\langle{v},{\Sigma^{-1/2}X}\right\rangle}\right\rvert }>t}\right]}\leq\textnormal{poly}(n)\times\exp\left(-\Omega\left(t^{C}\right)\right)$ $$

where

$$ $\Sigma={\mathop{\mathbb{E}}_{X\sim\mathcal{X}}\left[{XX^{\intercal}}\right]}= \textup{Covariance}(\mathcal{X})$ $$

In the previous section, we proved a generalisation of the matrix Bernstein inequality that can deal with a small probability that the norm bound assumption of the Bernstein inequality is violated.
A corollary of this result is that almost well-behaved distributions are closed to summation over polynomially many iid samples.
More concretely:

###### Lemma 7.6 .

Let $\mathcal{X}$ be a distribution that is almost well-behaved with respect to the scaling parameter $n$.
Let $k=\textnormal{poly}(n)$ and let $X_{1},\ldots,X_{k}\sim\mathcal{X}$ be iid random variables in $\mathbb{R}^{d}$, drawn from $\mathcal{X}$.

Then, the distribution of their empirical mean (i.e., of the variable $\hat{X}=\frac{1}{k}\sum_{i\in[k]}X_{i}$) is also almost well-behaved.

###### Proof of Lemma 7.6 .

Let $\hat{\mathcal{X}}$ denote the distribution of $\hat{X}$.

First, we note that from linearity of expectation ${\mathop{\mathbb{E}}{\left[{\hat{X}}\right]}}=\frac{1}{k}\sum_{i\in[k]}{
\mathop{\mathbb{E}}{\left[{X_{i}}\right]}}=0$, so $\hat{\mathcal{X}}$ is indeed a mean zero distribution.
Let $\Sigma={\mathop{\mathbb{E}}_{X\sim\mathcal{X}}\left[{XX^{\intercal}}\right]}$.
Again, using linearity of expectation, we have

$$ $\hat{\Sigma}={\mathop{\mathbb{E}}_{X\sim\hat{\mathcal{X}}}\left[{XX^{\intercal }}\right]}=\frac{1}{k^{2}}\sum_{i\in[k]}\Sigma=\frac{1}{k}\Sigma$ $$

Finally, we need to show that the tails of $\hat{\mathcal{X}}$ are bounded.
Let $v\in\mathbb{S}^{d-1}$ be any vector on the unit sphere.
From our assumption that $\mathcal{X}$ is well-behaved, for any $t>0$, we have

$$ ${\Pr_{X\sim\mathcal{X}}\left[{{\left\lvert{\left\langle{v},{\Sigma^{-1/2}X} \right\rangle}\right\rvert}\geq t}\right]}\leq\textnormal{poly}(n)\times\exp \left(-\Omega\left(t^{C}\right)\right)$ $$

Now, consider

$$ $Z=\left\langle{v},{\hat{\Sigma}^{-1/2}\hat{X}}\right\rangle=\frac{1}{\sqrt{k}} \sum_{i\in[k]}\left\langle{v},{\Sigma^{-1/2}X_{i}}\right\rangle\,.$ $$

We would like to bound the probability that ${\left\lvert{Z}\right\rvert}$ is greater than some threshold $T$.

Let $Z_{i}=\left\langle{v},{\Sigma^{-1/2}X_{i}}\right\rangle$.
These are iid random variables and our goal is to prove a concentration bound on their sum, so we may try to use Bernstein-type inequalities to do so.

Since ${\mathop{\mathbb{E}}{\left[{X_{i}}\right]}}=0$ and ${\mathop{\mathbb{E}}{\left[{X_{i}X_{i}^{\intercal}}\right]}}=\Sigma$, we have that ${\mathop{\mathbb{E}}{\left[{Z_{i}}\right]}}=0$ and ${\mathop{\mathbb{E}}{\left[{Z_{i}}\right]}}^{2}=1$.
Moreover, from the assumption that $\mathcal{X}$ is well behaved, we have concentration bounds on these individual variables.
Indeed, for any $t$, we have that

$$ $\delta_{t}\stackrel{{\scriptstyle\textup{def}}}{{=}}{\Pr\left[{{\left\lvert{X_ {i}}\right\rvert}\geq t}\right]}=\textnormal{poly}(n)\times\exp\left(-\Omega \left(t^{C}\right)\right)$ $$

Therefore, from the $d=1$ dimensional case of the approximate matrix Bernstein inequality (Lemma 7.5), we have that

$$ ${\Pr\left[{{\left\lvert{\left\langle{v},{\hat{\Sigma}^{-1/2}\hat{\mathcal{X}}} \right\rangle}\right\rvert}\geq\tau+k\sqrt{\delta_{t}}}\right]}\leq\exp\left(- \frac{\frac{1}{2}\tau^{2}k}{k+\frac{1}{3}t\tau\sqrt{k}}\right)+k\times \textnormal{poly}(n)\times\exp\left(-\Omega\left(t^{C}\right)\right)$ $$

Setting $t=\sqrt{\frac{T}{2}}$, we consider two regimes.
When $t<t_{\textup{threshold}}=\Theta\left(\log(n)^{1/C}\right)$, the bound

$$ ${\Pr\left[{{\left\lvert{\left\langle{v},{\hat{\Sigma}^{-1/2}\hat{X}}\right \rangle}\right\rvert}\geq T}\right]}\leq\textnormal{poly}(n)\times\exp\left(- \Omega\left(t^{C}\right)\right)$ $$

holds vacuously, since the RHS is greater than $1$ and the LHS is a probability.

When $t\geq t_{\textup{threshold}}$, we can ensure that $\delta_{t}<1/k^{2}=1/\textnormal{poly}(n)$.
Therefore in this case for $\tau=t^{2}$, we have

$$ ${\Pr\left[{{\left\lvert{\left\langle{v},{\hat{\Sigma}^{-1/2}\hat{X}}\right \rangle}\right\rvert}\geq T}\right]}\leq{\Pr\left[{{\left\lvert{\left\langle{v },{\hat{\Sigma}^{-1/2}\hat{X}}\right\rangle}\right\rvert}\geq\tau+k\sqrt{ \delta_{t}}}\right]}\leq\exp\left(-\Omega\left(\tau^{1/2}\right)\right)+k \times\textnormal{poly}(n)\times\exp\left(-\Omega\left(\tau^{C/2}\right)\right)$ $$

completing our proof of Lemma 7.6.

∎

###### Definition 3 .

###### Lemma 7.6 .

###### Proof of Lemma 7.6 .

#### 7.6.3 Renormalization

Our goal in this portion of the proof will be to show that the empirical covariance matrix $\Sigma=X^{\intercal}X$ is not far from the expectation ${\Sigma}_{\textup{gt}}=\mathbb{E}_{X\sim\mathcal{X}^{n}}[X^{\intercal}X]=nI$.
In particular, we prove the following claim:

###### Claim 7.7 .

Let $X_{1},\ldots,X_{n}\sim\mathcal{X}$ be samples drawn iid from an almost well-behaved distribution $\mathcal{X}$ with covariance identity.
Let $\Sigma=\sum_{i\in[n]}X_{i}X_{i}^{\intercal}$ denote the unnormalized empirical covariance, and set ${\Sigma}_{\textup{gt}}\stackrel{{\scriptstyle\textup{def}}}{{=}}{\mathop{
\mathbb{E}}{\left[{\Sigma}\right]}}=nI$.
Then,

With very high probability, $\left\|{\Sigma-{\Sigma}_{\textup{gt}}}\right\|<\widetilde{O}\left(\sqrt{nd}+d
\right)=o(n)$.

Moreover, for all $i$, we have $\left\|{{\mathop{\mathbb{E}}{\left[{\left\|{X_{i}}\right\|^{2}X_{i}X_{i}^{
\intercal}}\right]}}}\right\|=O\left({d}\right)$ and with very high probability $\left\|{X_{i}}\right\|^{2}=\widetilde{O}\left({d}\right)$.

###### Proof of Claim 7.7 .

We prove Claim 7.7 using our adaptation of the matrix Bernstein inequality (see Lemma 7.5).
We will apply this inequality for the $Z_{i}=X_{i}X_{i}^{\intercal}-I_{d}$, allowing us to obtain probabilistic bounds on

$$ $\Sigma-{\Sigma}_{\textup{gt}}=\sum_{i=1}^{n}Z_{i}$ $$

To use the approximate matrix Bernstein, we begin by proving a bound on the norm of $\max_{i}\left\|{Z_{i}}\right\|$ that holds with probability $1-\frac{1}{\textup{superpoly}(n)}$.
Note that

$$ $\left\|{Z_{i}}\right\|\leq\left\|{X_{i}X_{i}^{\intercal}}\right\|+\left\|{I} \right\|=\left\|{X_{i}}\right\|^{2}+1=\left(\sum_{j\in[d]}\left\langle{e_{j}}, {X_{i}}\right\rangle^{2}\right)+1$ $$

Next, we use the fact that $X_{i}$ is drawn from a well-behaved distribution $\mathcal{X}$ to prove strong tail bounds on the distribution of its norm.
In particular for any primary axis $e_{j}$ (for $j\in[d]$), we have $\left\langle{e_{j}},{X_{i}}\right\rangle\leq\textup{polylog}(n)$ with very high probability.
Therefore, using the union bound over all $i,j$, we also have that with very high probability

$$ $\max_{\begin{subarray}{c}i\in[n]\\ j\in[d]\end{subarray}}\left\{{\left\langle{e_{j}},{X_{i}}\right\rangle^{2}} \right\}\leq\textup{polylog}\left(n\right)\,.$ $$

Therefore, with very high probability

$$ $\left\|{Z_{i}}\right\|\leq\widetilde{O}\left(d\right)$ (19) $$

Next, we consider the second moment of $Z_{i}$.

$$ ${\mathop{\mathbb{E}}{\left[{Z_{i}^{2}}\right]}}={\mathop{\mathbb{E}}{\left[{ \left\|{X_{i}}\right\|^{2}X_{i}X_{i}^{\intercal}-2X_{i}X_{i}^{\intercal}+I} \right]}}={\mathop{\mathbb{E}}{\left[{\left\|{X_{i}}\right\|^{2}X_{i}X_{i}^{ \intercal}}\right]}}-I$ $$

Fix some pair of primary axis $e_{s}$ and unit vector $v\in\mathbb{S}^{d-1}$.
Recall that by our definition of $\mathcal{X}$ being a well-behaved distribution, we have an exponentially decaying concentration bound on the projection of our samples onto either of these axes

$$ ${\Pr_{X\sim\mathcal{X}}\left[{\max\left\{{{\left\lvert{\left\langle{X},{e_{s}} \right\rangle}\right\rvert},{\left\lvert{\left\langle{X},{v}\right\rangle} \right\rvert}}\right\}\geq t}\right]}=\exp\left(-\Omega\left(t^{C}\right)\right)$ $$

Therefore, a similar exponential tail bound also holds on the square of the product of these projections

$$ ${\Pr_{X\sim\mathcal{X}}\left[{{\left\lvert{\left\langle{X},{e_{s}}\right \rangle}\right\rvert}^{2}{\left\lvert{\left\langle{X},{v}\right\rangle}\right \rvert}^{2}\geq t}\right]}=\exp\left(-\Omega\left(t^{C/4}\right)\right)$ $$

In particular, we can conclude the far milder bound that the expectation of the squared product of these projections has bounded mean:

$$ ${\mathop{\mathbb{E}}_{X\sim\mathcal{X}}\left[{{\left\lvert{\left\langle{X},{e_ {s}}\right\rangle}\right\rvert}^{2}{\left\lvert{\left\langle{X},{v}\right \rangle}\right\rvert}^{2}}\right]}=\int_{t\in[0,\infty)}{\Pr_{X\sim\mathcal{X} }\left[{{\left\lvert{\left\langle{X},{e_{s}}\right\rangle}\right\rvert}^{2}{ \left\lvert{\left\langle{X},{v}\right\rangle}\right\rvert}^{2}\geq t}\right]} \leq\int_{t\in[0,\infty)}\exp\left(-\Omega\left(t^{C/4}\right)\right)=O(1)$ (20) $$

Note that equation (20) is no longer a concentration bound that holds with high probability, but a bound on the expectation of a random variable, that holds for *any* such pair $v,e_{s}$.
In particular, we have

$$ $\max_{v,e_{s}\in\mathbb{S}^{d-1}}\left\{{{\mathop{\mathbb{E}}_{X\sim\mathcal{X }}\left[{{\left\lvert{\left\langle{X},{e_{s}}\right\rangle}\right\rvert}^{2}{ \left\lvert{\left\langle{X},{v}\right\rangle}\right\rvert}^{2}}\right]}}\right \}=O(1)$ $$

Therefore, for any $v\in\mathbb{S}^{d-1}$,

$$ $v^{\intercal}{\mathop{\mathbb{E}}{\left[{\left\|{X_{i}}\right\|^{2}X_{i}X_{i}^ {\intercal}}\right]}}v={\mathop{\mathbb{E}}{\left[{\left\langle{v},{X_{i}} \right\rangle^{2}\left\|{X_{i}}\right\|^{2}}\right]}}={\mathop{\mathbb{E}}{ \left[{\sum_{s\in[d]}\left\langle{v},{X_{i}}\right\rangle^{2}\left\langle{e_{s }},{X_{i}}\right\rangle^{2}}\right]}}=\sum_{s\in[d]}{\mathop{\mathbb{E}}{\left [{\left\langle{v},{X_{i}}\right\rangle^{2}\left\langle{e_{s}},{X_{i}}\right \rangle^{2}}\right]}}=O(d)$ $$

and since this holds for all $v$, it is also true when maximizing over the unit sphere

$$ $\left\|{{\mathop{\mathbb{E}}{\left[{\left\|{X_{i}}\right\|^{2}X_{i}X_{i}^{ \intercal}}\right]}}}\right\|=\max_{v\in\mathbb{S}^{d-1}}\left\{{v^{\intercal} {\mathop{\mathbb{E}}{\left[{\left\|{X_{i}}\right\|^{2}X_{i}X_{i}^{\intercal}} \right]}}v}\right\}=O(d)$ $$

Therefore

$$ $\left\|{\sum_{i\in[n]}{\mathop{\mathbb{E}}{\left[{Z_{i}^{2}}\right]}}}\right\| \leq\sum_{i\in[n]}\left\|{{\mathop{\mathbb{E}}{\left[{Z_{i}^{2}}\right]}}} \right\|=O(nd)$ (21) $$

Combining equations (19) and (21) with Lemma 7.5 yields Claim 7.7
∎

###### Claim 7.7 .

###### Proof of Claim 7.7 .

#### 7.6.4 Maximal Norm

We now proceed to prove that each of the conditions required for the regression to be well-behaved occurs with very high probability.
The first (and easiest to prove) is the condition that $\max_{i\in[n]}X_{i}^{\intercal}\Sigma^{-1}X_{i}$ is bounded.

###### Claim 7.8 .

Let $X_{1},\ldots,X_{n}\sim\mathcal{X}$ be samples drawn iid from an almost well-behaved distribution $\mathcal{X}$ with covariance identity.
Let $\Sigma=\sum_{i\in[n]}X_{i}X_{i}^{\intercal}$ denote the unnormalized empirical covariance, and set ${\Sigma}_{\textup{gt}}\stackrel{{\scriptstyle\textup{def}}}{{=}}{\mathop{
\mathbb{E}}{\left[{\Sigma}\right]}}=nI$.

Then, with very high probability

$$ $\max_{i\in[n]}\left\{{X_{i}^{\intercal}\Sigma^{-1}X_{i}}\right\}=\widetilde{O} \left(\frac{d}{n}\right)$ $$

###### Proof of Claim 7.8 .

Let $\lambda$ denote the spectrum of a matrix.
In the proof of Claim 7.7, we already showed that with very high probability $\left\|{\Sigma-{\Sigma}_{\textup{gt}}}\right\|=o(n)=o\left(\min\lambda\left({
\Sigma}_{\textup{gt}}\right)\right)$.
When this holds, we also have

$$ $\lambda\left(\Sigma^{-1}\right)\subseteq\left(1\pm o(1)\right)\times\frac{1}{n }\,.$ $$

Moreover, the second part of Claim 7.7 states that with very high probability $\max_{i\in[n]}\left\|{X_{i}}\right\|^{2}=\widetilde{O}(d)$.

Combining these results, we have that with very high probability

$$ $\max_{i\in[n]}\left\{{X_{i}^{\intercal}\Sigma^{-1}X_{i}}\right\}\leq\max \lambda\left(\Sigma^{-1}\right)\times\max_{i\in[n]}\left\{{\left\|{X_{i}} \right\|^{2}}\right\}=\widetilde{O}\left(\frac{d}{n}\right)$ $$

∎

###### Claim 7.8 .

###### Proof of Claim 7.8 .

#### 7.6.5 Bounded Inner Products

For the next step of our proof, we show that the second condition of ACRE-friendliness of the regression holds with very high probability.
In particular, we will show that

###### Claim 7.9 .

Let $X_{1},\ldots,X_{n}\sim\mathcal{X}$ be samples drawn iid from an almost well-behaved distribution $\mathcal{X}$ with covariance identity.
Let $\Sigma=\sum_{i\in[n]}X_{i}X_{i}^{\intercal}$ denote the unnormalized empirical covariance, and set ${\Sigma}_{\textup{gt}}\stackrel{{\scriptstyle\textup{def}}}{{=}}{\mathop{
\mathbb{E}}{\left[{\Sigma}\right]}}=nI$.

Then, with very high probability

$$ $\max_{i\neq j\in[n]}\left\{{{\left\lvert{X_{i}^{\intercal}\Sigma^{-1}X_{j}} \right\rvert}}\right\}=\widetilde{O}\left(\frac{\sqrt{d}}{n}\right)\,.$ $$

###### Proof of Claim 7.9 .

Denote $A=\Sigma-X_{i}X_{i}^{\intercal}=\Sigma_{[n]\setminus\left\{{i}\right\}}$, and $v=X_{i}$.
Note that from Claim 7.7, we know that with very high probability $v^{\intercal}A^{-1}v\leq\left\|{X_{i}}\right\|^{2}\left\|{A^{-1}}\right\|=
\widetilde{O}\left(\frac{d}{n}\right)=o(1)$.

Therefore, in this regime we may apply the Sherman-Morrison formula to show that

$$ $\Sigma^{-1}=\left(A+vv^{\intercal}\right)^{-1}=A^{-1}-\frac{A^{-1}vv^{ \intercal}A^{-1}}{1+v^{\intercal}A^{-1}v^{\intercal}}\,.$ $$

Note that neither $A$ nor $X_{j}$ depend on $X_{i}$, so from our assumption that $\mathcal{X}$ is well behaved, with very high probability

$$ ${\left\lvert{X_{i}^{\intercal}A^{-1}X_{j}}\right\rvert}=\widetilde{O}\left( \left\|{A^{-1}X_{j}}\right\|\right)=\widetilde{O}\left(\frac{\sqrt{d}}{n} \right)\,.$ $$

Similarly, for our target expression, with very high probability,

$$ $\displaystyle{\left\lvert{X_{i}^{\intercal}\Sigma^{-1}X_{j}}\right\rvert}$ $\displaystyle={\left\lvert{X_{i}^{\intercal}\left(A^{-1}-\frac{A^{-1}X_{i}X_{i }^{\intercal}A^{-1}}{1+X_{i}^{\intercal}A^{-1}X_{i}^{\intercal}}\right)X_{j}} \right\rvert}=$ $\displaystyle=\left(1-\frac{X_{i}^{\intercal}A^{-1}X_{i}}{1-X_{i}^{\intercal}A ^{-1}X_{i}}\right){\left\lvert{X_{i}^{\intercal}A^{-1}X_{j}}\right\rvert}= \widetilde{O}\left(\frac{\sqrt{d}}{n}\right)$ $$

∎

###### Claim 7.9 .

###### Proof of Claim 7.9 .

#### 7.6.6 Projection on the e 𝑒 e italic_e Axis

For the next property of well-behaved regressions, we will want to show that with very high probability $e^{\intercal}\Sigma^{-1}X_{i}$ is bounded for all $i$.
Note that from the exponential decay assumption due to $\mathcal{X}$ being well-behaved would suffice to give a good bound on $e^{\intercal}\hat{\Sigma}^{-1}X_{i}$, but as before, the challenge will be to show that $\Sigma^{-1}$ doesn’t rotate $X_{i}$ onto $e$.

###### Claim 7.10 .

Let $X_{1},\ldots,X_{n}\sim\mathcal{X}$ be samples drawn iid from an almost well-behaved distribution $\mathcal{X}$ with covariance identity.
Let $\Sigma=\sum_{i\in[n]}X_{i}X_{i}^{\intercal}$ denote the unnormalized empirical covariance, and set ${\Sigma}_{\textup{gt}}\stackrel{{\scriptstyle\textup{def}}}{{=}}{\mathop{
\mathbb{E}}{\left[{\Sigma}\right]}}=nI$.
Let $e\in\mathbb{S}^{d-1}$ be any fixed vector independent of the $X_{i}$.

Then, with very high probability

$$ $\forall i\in[n]\;\;\;\;e^{\intercal}\Sigma^{-1}X_{i}=\frac{e^{\intercal}X_{i}} {n}\pm o\left(\frac{1}{n}\right)\,.$ $$

In particular, with very high probability

$$ $\max_{i\in[n]}{\left\lvert{\left\langle{\Sigma^{-1}X_{i}},{e}\right\rangle} \right\rvert}=\widetilde{O}\left(\frac{1}{n}\right)=\sqrt{\widetilde{O}\left( \frac{1}{n}\right)\times e^{\intercal}\Sigma^{-1}e}$ $$

###### Proof of Claim 7.10 .

As in the proof of Claim 7.9, let $v=X_{i}$ and $A=\Sigma-vv^{\intercal}$.
Moreover, because with very high probability $v^{\intercal}A^{-1}v=o(1)<1$, we may apply the Sherman-Morrison formula

$$ $\Sigma^{-1}=\left(A+vv^{\intercal}\right)^{-1}=A^{-1}-\frac{A^{-1}vv^{ \intercal}A^{-1}}{1+v^{\intercal}A^{-1}v^{\intercal}}\,.$ $$

Therefore, with high very high probability,

$$ $\displaystyle e^{\intercal}\Sigma^{-1}X_{i}$ $\displaystyle=e^{\intercal}A^{-1}X_{i}-\frac{e^{\intercal}A^{-1}X_{i}X_{i}^{ \intercal}A^{-1}X_{i}^{\intercal}}{1+X_{i}^{\intercal}A^{-1}X_{i}}=$ $\displaystyle=e^{\intercal}{\Sigma}_{\textup{gt}}^{-1}X_{i}+e^{\intercal}\left (A^{-1}-{\Sigma}_{\textup{gt}}^{-1}\right)X_{i}-\frac{e^{\intercal}A^{-1}X_{i} X_{i}^{\intercal}A^{-1}X_{i}^{\intercal}}{1+X_{i}^{\intercal}A^{-1}X_{i}}$ $$

From Claim 7.7, we know that with very high probability

$$ $\displaystyle\left\|{A^{-1}-{\Sigma}_{\textup{gt}}^{-1}}\right\|\leq\frac{1}{2 n^{2}}\left\|{A-{\Sigma}_{\textup{gt}}}\right\|=\widetilde{O}\left({\frac{d+ \sqrt{nd}}{n^{2}}}\right)$ $\displaystyle X_{i}^{\intercal}A^{-1}X_{i}\leq\left\|{X_{i}}\right\|^{2}\left \|{A^{-1}}\right\|=\widetilde{O}\left({\frac{d}{n}}\right)\,.$ $$

Therefore, using the fact that $X_{i}$ is well-behaved and independent of $e$ and $A$, we have that with very high probability

$$ ${\left\lvert{e^{\intercal}\Sigma^{-1}X_{i}-e^{\intercal}{\Sigma}_{\textup{gt}} ^{-1}X_{i}}\right\rvert}\leq{\left\lvert{e^{\intercal}\left(A^{-1}-{\Sigma}_{ \textup{gt}}^{-1}\right)X_{i}}\right\rvert}+{\left\lvert{\frac{e^{\intercal}A^ {-1}X_{i}X_{i}^{\intercal}A^{-1}X_{i}^{\intercal}}{1+X_{i}^{\intercal}A^{-1}X_ {i}}}\right\rvert}=\widetilde{O}\left({\frac{d+\sqrt{nd}}{n^{2}}}\right)=o \left(\frac{1}{n}\right)\,.$ $$

∎

###### Claim 7.10 .

###### Proof of Claim 7.10 .

#### 7.6.7 Bounded Residuals

For the next step of our analysis we will show that the residuals are bounded with very high probability.
Recall that under the assumptions of Claim 7.3, we assume that the labels are drawn from the distribution

$$ $Y=X{\beta}_{\textup{gt}}+\zeta\sim X{\beta}_{\textup{gt}}+\mathcal{N}\left( \vec{0},I_{n}\right)$ $$

Therefore, we clearly have that with very high probability $\max_{i\in[n]}{\left\lvert{\zeta_{i}}\right\rvert}\leq\log(n)=\widetilde{O}(1)$ as required.
The issue is that the residuals of the regression are not necessarily equal to $\zeta$.

Recall that the residuals are equal to

$$ $R\stackrel{{\scriptstyle\textup{def}}}{{=}}Y-X\Sigma^{-1}X^{\intercal}Y=\zeta- X\Sigma^{-1}X^{\intercal}\zeta$ $$

###### Claim 7.11 .

With very high probability

$$ $\forall i\;\;\;\;{\left\lvert{R_{i}-\zeta_{i}}\right\rvert}=o(1)$ $$

###### Proof of Claim 7.11 .

Fix some index $i\in[n]$.

$$ $R_{i}=\zeta_{i}-X_{i}^{\intercal}\Sigma^{-1}X^{\intercal}\zeta=\zeta_{i}-\sum_ {j}X_{i}^{\intercal}\Sigma^{-1}X_{j}\zeta_{j}=\zeta_{i}-\zeta_{i}X_{i}^{ \intercal}\Sigma^{-1}X_{i}-\sum_{j\neq i}\zeta_{j}X_{i}^{\intercal}\Sigma^{-1} X_{j}$ (22) $$

Therefore, from Claim 7.8, we have

$$ ${\left\lvert{R_{i}-\zeta_{i}}\right\rvert}={\left\lvert{\sum_{j\neq i}\zeta_{j }X_{i}^{\intercal}\Sigma^{-1}X_{j}}\right\rvert}+\widetilde{O}\left({\frac{d}{ n}}\right)={\left\lvert{\sum_{j\neq i}\zeta_{j}X_{i}^{\intercal}\Sigma^{-1}X_{ j}}\right\rvert}+o(1)$ $$

This leaves us with the task of analyzing the term

$$ $X_{i}^{\intercal}\Sigma^{-1}\left(\sum_{j\neq i}\zeta_{j}X_{j}\right)=\sum_{j \neq i}Z_{j}$ $$

where $Z_{j}=X_{i}^{\intercal}\Sigma^{-1}X_{j}\zeta_{j}$.
To bound this term, we view the process of generating the samples as first generating the covariates $X_{i}$, and then after fixing some values for the $X_{i}$, it generates the errors $\zeta_{j}$.

In other words, we will show that with very high probability over the $X_{i}$

$$ ${\Pr_{\zeta_{j}\sim\mathcal{N}(0,1)}\left[{{\left\lvert{\sum_{j\neq i}Z_{j}} \right\rvert}>o(1)}\right]}<\frac{1}{\textup{super-poly}(n)}\,.$ $$

In particular, with very high probability over the covariates, we have

$$ $\forall j\neq i\;\;\;\;\left(X_{i}^{\intercal}\Sigma^{-1}X_{j}\right)^{2}= \widetilde{O}\left({\frac{d}{n^{2}}}\right)\Rightarrow\sum_{j\in[n]\setminus \left\{{i}\right\}}\left(X_{i}^{\intercal}\Sigma^{-1}X_{j}\right)^{2}= \widetilde{O}\left({\frac{d}{n}}\right)\,.$ $$

Therefore, fixing the covariates $X_{i}$, and viewing $\sum_{j\neq i}Z_{j}$ as a random variable dependent on the randomness of the errors $\zeta_{j}$, we have

$$ $\sum_{j\neq i}Z_{j}\sim\mathcal{N}\left(0,\sum_{j\in[n]\setminus\left\{{i} \right\}}\left(X_{i}^{\intercal}\Sigma^{-1}X_{j}\right)^{2}\right)=\mathcal{N} \left(0,\widetilde{O}\left({\frac{d}{n}}\right)\right)\,,$ $$

yielding the claim.

∎

###### Claim 7.11 .

###### Proof of Claim 7.11 .

#### 7.6.8 Large Influence Scores

For the final step of our proof of Claim 7.3, we will show that with very high probability, there are many samples in the regression that have relatively high AMIP influence scores.

###### Claim 7.12 .

Let $\alpha_{i}=e^{\intercal}\Sigma^{-1}X_{i}R_{i}$ denote the AMIP influence score of the $i$th sample.
Then with very high probability, there are is a set $T_{0}\subseteq[n]$ of least $k_{0}=\widetilde{\Omega}(n)$ “influential samples” – i.e., such that $\forall i\in T\;\;\alpha_{i}\geq\frac{1}{10n}$.

Proving Claim 7.12 will also conclude our proof of Claim 7.3, as this will show that all conditions required for a regression to be well-behaved are fulfilled with very high probability.

###### Proof of Claim 7.12 .

To prove Claim 7.12 we first show that a very large number of samples must have a relatively high inner product with the axis of interest.
In other words, we will show that with very high probability

$$ ${\left\lvert{\left\{i\middle|{\left\lvert{e^{\intercal}\Sigma^{-1}X_{i}}\right \rvert}\geq\frac{1}{2n}\right\}}\right\rvert}=\widetilde{\Omega}(n)$ (23) $$

To prove equation (23), note that:

- •
$e^{\intercal}\Sigma^{-1}e=\frac{1\pm o(1)}{n}$ (this follows immediately from Claim 7.7).
- •
$\sum_{i\in[n]}{\left\lvert{e^{\intercal}\Sigma^{-1}X_{i}}\right\rvert}^{2}=
\sum_{i\in[n]}e^{\intercal}\Sigma^{-1}X_{i}X_{i}^{\intercal}\Sigma^{-1}e=e^{
\intercal}\Sigma^{-1}e$.
- •
In Claim 7.10, we showed that wvhp $\forall i\;\;\left(e^{\intercal}\Sigma^{-1}X_{i}\right)^{2}\leq\widetilde{O}
\left(\frac{1}{n}\right)\times e^{\intercal}\Sigma^{-1}e$.

Therefore, we must have at least $\widetilde{\Omega}(n)$ samples with ${\left\lvert{e^{\intercal}\Sigma^{-1}X_{i}}\right\rvert}\geq\frac{1}{2n}$.
Let $i$ be the index of such a sample.
If $\textup{sign}\left({e^{\intercal}\Sigma^{-1}X_{i}}\right)\times R_{i}\geq\frac
{1}{5}$, then we will also have $\alpha_{i}\geq\frac{1}{10n}$.

In the proof of Claim 7.11, we showed that wvhp ${\left\lvert{R_{i}-\zeta_{i}}\right\rvert}=o(1)$ for all $i$ (where $\zeta_{i}=Y_{i}-X_{i}^{\intercal}{\beta}_{\textup{gt}}$ are the “ground truth residuals”, and are drawn iid from a normal distribution).
In particular, wvhp $\forall i\;\;{\left\lvert{R_{i}-\zeta_{i}}\right\rvert}<\frac{1}{4}-\frac{1}{5}$, so as long as $\textup{sign}\left({e^{\intercal}\Sigma^{-1}X_{i}}\right)\zeta_{i}\geq\frac{1}
{4}$, we have $\alpha_{i}\geq\frac{1}{10n}$.

But $\zeta_{i}$ is drawn iid from a normal distribution, so $\textup{sign}\left({e^{\intercal}\Sigma^{-1}X_{i}}\right)\zeta_{i}\geq\frac{1}
{4}$ has constant probability and is independent of $X_{i}$.
Therefore, applying the Hoeffding-Chernoff bound, we can easily see that wvhp at least a constant fraction of the $\widetilde{\Omega}(n)$ samples for which ${\left\lvert{e^{\intercal}\Sigma^{-1}X_{i}}\right\rvert}\geq\frac{1}{2n}$ also have $\alpha_{i}\geq\frac{1}{10n}$, thus concluding our proof of Claims 7.12 and 7.3.

∎

###### Claim 7.12 .

###### Proof of Claim 7.12 .

## 8 Tightness of OHARE

In the previous section, we proved Theorem 7.1 which says that for “well-behaved” data, the ACRE algorithm outputs nearly tight bounds on the removal effects for a range of removal set sizes $k$.
In this section, we will extend those results to the one-hot aware version of the algorithm – OHARE .

###### Theorem 8.1 ( OHARE Bounds are Tight on Well-behaved Data) .

Consider a linear regression from a set of continuous features $X\in\mathbb{R}^{n\times d}$ and a set of $m$ dummy variables, representing a categorical feature $B_{1}\sqcup\cdots\sqcup B_{m}=[n]$, to a target variable $Y$.

For any fixed $\varepsilon>0$, there exists $\nu\in\textup{polylog}(n)$ such that:

If $n_{j}={\left\lvert{B_{j}}\right\rvert}$ denote the number of samples that take the value $j$ in the categorical feature, and for all $j\in[m]$, we have

$$ $n^{\varepsilon}+\nu\sqrt{d}<n_{j}<0.49n\,,$ $$

that the dimension of the continuous features $d$ is at most $d\leq n^{4/5}/\nu$, and that the continuous features are then drawn iid from a well-behaved distribution $X_{i}\sim\mathcal{X}$ independently of their value on the categorical feature.

And if the outcomes $Y$ are drawn iid from a normal distribution around a linear model of the features

$$ $Y_{i}\sim\underbrace{\mu_{j(i)}}_{\textup{categorical contribution}}+ \underbrace{\left\langle{X_{i}},{{\beta}_{\textup{gt}}}\right\rangle}_{\textup {continuos contribution}}+\underbrace{\mathcal{N}\left(0,1\right)}_{\textup{ error}}\,,$ $$

for some unknown ground truth linear model $(\mu,{\beta}_{\textup{gt}})\in\mathbb{R}^{m+d}$.

Then, for any axis $e\in\mathbb{S}^{d-1}$, with very high probability, the upper and lower bounds produced by OHARE on this regression are close to tight

$$ $\frac{U_{k}}{L_{k}}=1+O\left(\frac{\textup{polyloglog}(n)}{\sqrt{\log(n)}}\right)$ $$

for all $k<k_{\textup{threshold}}$, where

$$ $k_{\textup{threshold}}=\widetilde{\Theta}\left(\min\left\{\frac{n}{\sqrt{d}}, \frac{n^{2}}{d^{2}},n^{1-\varepsilon}\right\}\right)\,.$ $$

Our goal for the rest of this section will be to prove Theorem 8.1.

###### Theorem 8.1 ( OHARE Bounds are Tight on Well-behaved Data) .

### 8.1 Main Challenges and Proof Structure

Recall from Section 4 that the key idea of the OHARE algorithm is to analyse a process that is equivalent to the regression with a one-hot encoding.
In this alternative formulation of one-hot controlled regression, we first split our samples into buckets $B_{j}\subseteq[n]$ corresponding to each of the potential values of the categorical feature, reaverage the samples in each bucket

$$ $\displaystyle\widetilde{X}_{i}=X_{i}-{\mathop{\mathbb{E}}_{i^{\prime}\in B_{j( i)}}\left[{X_{i^{\prime}}}\right]}\in\mathbb{R}^{d}$ $\displaystyle\widetilde{Y}_{i}=Y_{i}-{\mathop{\mathbb{E}}_{i^{\prime}\in B_{j( i)}}\left[{Y_{i^{\prime}}}\right]}\in\mathbb{R}$ $$

and perform a regression with just the reaveraged continuous features.
The OHARE algorithm then computes the same MSN bounds as the ACRE algorithm would, but on these reaveraged continous features and combines them with terms corresponding to the effect a removal might have on the reaveraging process.

Our proof of Theorem 8.1 will follow a similar path.
We will first prove a claim very similar to Claim 7.3 adapted to the OHARE case:

###### Claim 8.2 (Well-behaved distributions yield well-behaved regressions with high probability after reaveraging) .

Let $n,d,m,X,Y$ be as in Theorem 8.1.

Then, for any axis $e\in\mathbb{S}^{d-1}$, with very high probability, the regression $\widetilde{X},\widetilde{Y}$ is ACRE-friendly for all $k\leq k_{0}$, for $k_{0}=\widetilde{\Omega}(n^{1-\varepsilon})$.

Note that Claim 8.2 does not follow immediately from the corresponding Claim 7.3 for ACRE, since the continuous features $\widetilde{X}$ are no longer drawn iid from a well-behaved distribution as the reaveraging step could have changed them and similarly the reaveraged labels $\widetilde{Y}$ are not drawn iid from a normal distribution around a linear combination of the continous features.
The proof of Claim 8.2 will follow a very similar path to the proof of Claim 7.3, but will also have to account for these additional corrections.

Finally, we will prove that the additional corrections taken into account by OHARE will not change the upper and lower bounds too much, yielding Theorem 8.1.

###### Claim 8.2 (Well-behaved distributions yield well-behaved regressions with high probability after reaveraging) .

### 8.2 Proof of Claim 8.2

We begin by adapting the analysis from the continuous features (Claim 7.3) to the reaveraged samples.
Throughout this section, let $X\in\mathbb{R}^{n\times d}$ denote just the continuous covariates and $\widetilde{X}\in\mathbb{R}^{n\times d}$ denote the reaveraged continuous covariates.

As in the proof of Claim 7.3, we normalize our samples so that the ground truth covariance of the continuous features is equal to the identity, and denote by $\Sigma=X^{\intercal}X\in\mathbb{R}^{d\times d}$ the unnormalized empirical covariance of $X$ and by $\widetilde{\Sigma}=\widetilde{X}^{\intercal}\widetilde{X}\in\mathbb{R}^{d
\times d}$ denote the unnormalized empirical covariance of the reaveraged samples covariates.

#### 8.2.1 Reaveraging

The first step of our analysis will be to show that with very high probability the reaveraging step makes only a small change to the covariates as well as the target variable.
Let $j\in[m]$ be the index of any bucket of samples and let

$$ $\displaystyle\xi_{j}=\frac{1}{n_{j}}\sum_{i\in B_{j}}X_{i}={\mathop{\mathbb{E} }_{i\in B_{j}}\left[{X_{i}}\right]}$ $\displaystyle y_{j}=\frac{1}{n_{j}}\sum_{i\in B_{j}}Y_{i}={\mathop{\mathbb{E}} _{i\in B_{j}}\left[{Y_{i}}\right]}$ $$

denote the averaging effect for this bucket.

###### Claim 8.3 .

The mean and covariance of $\xi_{j}$ are

$$ $\displaystyle{\mathop{\mathbb{E}}{\left[{\xi_{j}}\right]}}=0$ $\displaystyle{\mathop{\mathbb{E}}{\left[{\xi_{j}\xi_{j}^{\intercal}}\right]}}= \frac{1}{n_{j}}I$ $$

Moreover, with very high probability

$$ $\displaystyle\left\|{\xi_{j}}\right\|\leq\widetilde{O}\left(\sqrt{d/n_{j}} \right)=o(\sqrt{d})$ $\displaystyle{\left\lvert{y_{j}-\mu_{j}-\left\langle{\xi_{j}},{{\beta}_{ \textup{gt}}}\right\rangle}\right\rvert}\leq\widetilde{O}\left(1/\sqrt{n_{j}} \right)=o(1)$ $$

Finally, the fourth moments of $\xi_{j}$ are also bounded

$$ $\left\|{{\mathop{\mathbb{E}}{\left[{\left\|{\xi_{j}}\right\|^{2}\xi_{j}\xi_{j} ^{\intercal}}\right]}}}\right\|=O\left(\frac{d}{n_{j}^{2}}\right)$ $$

###### Proof of Claim 8.3 .

First, recall that our covariate distribution $\mathcal{X}$ was normalized to have mean 0 0 and covariance identity, yielding the first part of Claim 8.3 immediately from the fact that $\xi_{j}$ is the empirical average of $n_{j}$ iid samples drawn from $\mathcal{X}$.

Recall that we assumed the target variable was drawn from a normal distribution

$$ $Y_{i}\sim\mu_{j}+\left\langle{X_{i}},{{\beta}_{\textup{gt}}}\right\rangle+ \mathcal{N}(0,1)$ $$

In particular,

$$ ${\mathop{\mathbb{E}}_{i\in B_{j}}\left[{Y_{i}}\right]}-\mu_{j}-{\mathop{ \mathbb{E}}_{i\in B_{j}}\left[{\left\langle{X_{i}},{{\beta}_{\textup{gt}}} \right\rangle}\right]}$ $$

is the empirical average over $n_{j}$ samples of this normal distribution, yielding the claim that with very high probability

$$ ${\left\lvert{y_{j}-\mu_{j}-\left\langle{\xi_{j}},{{\beta}_{\textup{gt}}}\right \rangle}\right\rvert}\leq\widetilde{O}\left(1/\sqrt{n_{j}}\right)=o(1)\,.$ $$

The rest of Claim 8.3 will follow from a combination of Claim 7.7 (which bounds the norms and higher moments of well-behaved distributions), and Lemma 7.6 (which states that the sum over iid samples from a well-behaved distribution is also well-behaved).

Indeed $\xi_{j}$ is the empirical average over $n_{j}=\textnormal{poly}(n)$ samples from a well-behaved distribution $\mathcal{X}$, so Lemma 7.6 ensures that for all $j$, the variable $\sqrt{n_{j}}\xi_{j}$ is well behaved, and, as noted above, it has covariance identity.
Therefore, it follows from Claim 7.7 that

$$ $\left\|{{\mathop{\mathbb{E}}{\left[{\left\|{\xi_{j}}\right\|^{2}\xi_{j}\xi_{j} ^{\intercal}}\right]}}}\right\|\leq O\left(\frac{d}{n_{j}^{2}}\right)$ $$

and that with very high probability

$$ $\left\|{\xi_{j}}\right\|\leq\widetilde{O}\left(\sqrt{d/n_{j}}\right)$ $$

∎

###### Claim 8.3 .

###### Proof of Claim 8.3 .

#### 8.2.2 Renormalization

Let $\Sigma=X^{\intercal}X=\sum_{i\in[n]}X_{i}X_{i}^{\intercal}$ be the unnormalized empirical covariance of the continuous features, and let ${\Sigma}_{\textup{gt}}=nI={\mathop{\mathbb{E}}{\left[{\Sigma}\right]}}$ be the ground truth mean of these features.
We continue along the same lines as the proof of Claim 7.3 by adapting Claim 7.7 to the reaveraged setting:

###### Claim 8.4 .

With very high probability, $\left\|{\widetilde{\Sigma}-\Sigma}\right\|=\widetilde{O}\left(d+m\right)$.

In particular, due to Claim 7.7 and the triangle inequality

$$ $\left\|{\widetilde{\Sigma}-{\Sigma}_{\textup{gt}}}\right\|\leq\left\|{ \widetilde{\Sigma}-{\Sigma}}\right\|+\left\|{{\Sigma}-{\Sigma}_{\textup{gt}}} \right\|=\widetilde{O}\left(\sqrt{nd}+d+m\right)=\widetilde{O}\left(\sqrt{nd}+ m\right)=o(n)$ $$

###### Proof of Claim 8.4 .

We have

$$ $\widetilde{\Sigma}=\sum_{i\in[n]}\widetilde{X}_{i}\widetilde{X}_{i}^{\intercal }=\sum_{i\in[n]}\left(X_{i}-\xi_{j(i)}\right)\left(X_{i}-\xi_{j(i)}\right)^{ \intercal}=\sum_{i\in[n]}X_{i}X_{i}^{\intercal}-\sum_{i\in[n]}\xi_{j(i)}\xi_{j (i)}^{\intercal}=\Sigma-\sum_{j\in[m]}n_{j}\xi_{j}\xi_{j}^{\intercal}$ $$

Therefore, it only remains to bound

$$ $\left\|{\sum_{i\in[n]}\xi_{j(i)}\xi_{j(i)}^{\intercal}}\right\|=\left\|{\sum_{ j\in[m]}n_{j}\xi_{j}\xi_{j}^{\intercal}}\right\|$ $$

We do this using our approximate matrix Bernstein inequality – Lemma 7.5.
Let $Z_{j}=n_{j}\xi_{j(i)}\xi_{j(i)}^{\intercal}$.

In Claim 8.3, we showed that ${\mathop{\mathbb{E}}{\left[{Z_{j}}\right]}}=I$, that ${\mathop{\mathbb{E}}{\left[{Z_{j}^{2}}\right]}}={\mathop{\mathbb{E}}{\left[{n_
{j}^{2}\left\|{\xi_{j}}\right\|^{2}\xi_{j}\xi_{j}^{\intercal}}\right]}}$ has bounded norm

$$ $\left\|{{\mathop{\mathbb{E}}{\left[{n_{j}^{2}\left\|{\xi_{j}}\right\|^{2}\xi_{ j}\xi_{j}^{\intercal}}\right]}}}\right\|=O\left(d\right)\,,$ $$

and that with very high probability

$$ $\left\|{Z_{j}}\right\|=n_{j}\left\|{\xi_{j}}\right\|^{2}=\widetilde{O}(d)\,.$ $$

Therefore, from Lemma 7.5, with very high probability

$$ $\left\|{\sum_{j\in[m]}Z_{j}}\right\|\leq\left\|{\sum_{j\in[m]}Z_{j}-\sum_{j\in [m]}{\mathop{\mathbb{E}}{\left[{Z_{j}}\right]}}}\right\|+m=\widetilde{O}\left( d+\sqrt{md}+m\right)=\widetilde{O}\left(d+m\right)$ $$

Finally, recall that we assumed to have at least $n^{\varepsilon}$ samples in the smallest category.
Therefore, $m\leq n^{1-c}\ll n$, completing the proof.

∎

###### Claim 8.4 .

###### Proof of Claim 8.4 .

#### 8.2.3 Maximal Norm

We proceed to prove that each of the conditions required for the regression to be well-behaved occurs with very high probability.
The first (and easiest to prove) is the condition that $\max_{i\in[n]}X_{i}^{\intercal}\Sigma^{-1}X_{i}$ is bounded.

###### Claim 8.5 .

With very high probability

$$ $\max_{i\in[n]}\left\{{\widetilde{X}_{i}^{\intercal}\widetilde{\Sigma}^{-1} \widetilde{X}_{i}}\right\}=\widetilde{O}\left(\frac{d}{n}\right)$ $$

###### Proof of Claim 8.5 .

Let $\lambda$ denote the spectrum of a matrix.
In the proof of Claim 8.4, we already showed that with very high probability $\left\|{\widetilde{\Sigma}-{\Sigma}_{\textup{gt}}}\right\|=o(n)=o\left(\min
\lambda\left({\Sigma}_{\textup{gt}}\right)\right)$ (where ${\Sigma}_{\textup{gt}}=nI$).
When this holds, we also have

$$ $\lambda\left(\widetilde{\Sigma}^{-1}\right)\subseteq\left(1\pm o(1)\right) \times\frac{1}{n}\,.$ $$

Moreover, in Claim 7.7 we show that with very high probability $\max_{i\in[n]}\left\|{X_{i}}\right\|^{2}=\widetilde{O}(d)$.
Combined with Claim 8.3 and the triangle inequality, we can derive a bound on $\left\|{\widetilde{X}_{i}}\right\|<\left\|{X_{i}}\right\|+\left\|{\xi_{j(i)}}
\right\|=\widetilde{O}\left({\sqrt{d}}\right)$ with very high probability.

Therefore, with very high probability

$$ $\max_{i\in[n]}\left\{{\widetilde{X}_{i}^{\intercal}\widetilde{\Sigma}^{-1} \widetilde{X}_{i}}\right\}\leq\max\left\{{\lambda\left(\widetilde{\Sigma}^{-1} \right)}\right\}\times\max_{i\in[n]}\left\{{\left\|{\widetilde{X}_{i}}\right\| ^{2}}\right\}=\widetilde{O}\left(\frac{d}{n}\right)$ $$

∎

###### Claim 8.5 .

###### Proof of Claim 8.5 .

#### 8.2.4 Bounded Inner Products

For the next step of our proof, we show that the second condition of ACRE-friendliness of the regression holds with very high probability.
In particular, we will show that

###### Claim 8.6 .

For $M\in\left\{{\Sigma^{-1},\widetilde{\Sigma}^{-1}}\right\}$, the following inequalities hold with very high probability

- •
$\max_{j_{1}\neq j_{2}\in[m]}\left\{{{\left\lvert{\xi_{j_{1}}^{\intercal}M\xi_{
j_{2}}}\right\rvert}}\right\}=\widetilde{O}\left(\frac{\sqrt{d}}{n\sqrt{n_{j_{
1}}n_{j_{2}}}}\times\left(1+\frac{d\sqrt{n_{j_{1}}}}{n}\right)\right)$
- •
$\forall i\in[n],j\in[m]\;\;\;\;{\left\lvert{X_{i}^{\intercal}M\xi_{j}}\right
\rvert}=\widetilde{O}\left(\frac{\sqrt{d}}{n\sqrt{n_{j}}}+\frac{d}{nn_{j}}1_{i
\in B_{j}}\right)$
- •
$\max_{i_{1}\neq i_{2}\in[n]}\left\{{{\left\lvert{\widetilde{X}_{i_{1}}^{
\intercal}M^{-1}\widetilde{X}_{i_{2}}}\right\rvert}}\right\}=\widetilde{O}
\left(\frac{\sqrt{d}}{n}\right)$

Claim 7.9 proves the third inequality of Claim 8.6 for the case where $M=\Sigma^{-1}$.
The following Lemma 8.7 will prove the first inequality for this same case.

###### Lemma 8.7 .

Let $j_{1}\neq j_{2}\in[m]$ be the indices of two distinct buckets, and let $\Sigma=X^{\intercal}X=\sum_{i\in[n]}X_{i}X_{i}^{\intercal}$ denote the unnormalized empirical covariance of the unaveraged samples.

With very high probability

$$ ${\left\lvert{\xi_{j_{1}}^{\intercal}\Sigma^{-1}\xi_{j_{2}}}\right\rvert}= \widetilde{O}\left(\frac{\sqrt{d}}{\sqrt{n_{j_{1}}n_{j_{2}}}\times n}\times \left(1+\frac{d\sqrt{n_{j_{1}}}}{n}\right)\right)\,.$ $$

The proof of Lemma 8.7 is very long and technical, and we devote Section 8.3 to it.
For now, let us continue with our proof of the rest of the inequalities in Claim 8.6 assuming Lemma 8.7.

###### Proof of Claim 8.6 .

We begin by bounding the inner product between the bucket averages through the covariance $\widetilde{\Sigma}^{-1}$.
Define

$$ $\eta\stackrel{{\scriptstyle\textup{def}}}{{=}}\max_{j_{1},j_{2}}\left\{{\frac{ {\left\lvert{\xi_{j_{1}}^{\intercal}\widetilde{\Sigma}^{-1}\xi_{j_{2}}}\right \rvert}}{\frac{\sqrt{d}}{\sqrt{n_{j_{1}}n_{j_{2}}}\times n}\times\left(1+\frac {d\sqrt{n_{j_{1}}}}{n}\right)}}\right\}\,.$ $$

Our first goal will be to show that with very high probability $\eta=\widetilde{O}(1)$.

Consider the following identity (where $A$ and $B$ are $d\times d$ matrices such that $A$ and $A-B$ are invertible):

$$ $A=\left(A-B\right)+B\,.$ $$

Multiplying this equation by $\left(A-B\right)^{-1}$ from the left and $A^{-1}$ yields the following identity that we will use in our analysis.

$$ $\left(A-B\right)^{-1}=A^{-1}+\left(A-B\right)^{-1}BA^{-1}$ (24) $$

We apply equation (24) for $A=\Sigma,B=\Sigma-\widetilde{\Sigma}=\sum_{j\in[m]}n_{j}\xi_{j}\xi_{j}^{\intercal}$, giving us the identity

$$ $\widetilde{\Sigma}^{-1}=\Sigma^{-1}+\sum_{j\in[m]}\widetilde{\Sigma}^{-1}\xi_{ j}n_{j}\xi_{j}^{\intercal}\Sigma^{-1}\,.$ $$

Therefore, from the triangle inequality,

$$ $\eta\times\frac{\sqrt{d}}{\sqrt{n_{j_{1}}n_{j_{2}}}n}\times\left(1+\frac{d \sqrt{n_{j_{1}}}}{n}\right)={\left\lvert{\xi_{j_{1}}^{\intercal}\widetilde{ \Sigma}^{-1}\xi_{j_{2}}}\right\rvert}\leq{\left\lvert{\xi_{j_{1}}^{\intercal}{ \Sigma}^{-1}\xi_{j_{2}}}\right\rvert}+\sum_{j\in[m]}{\left\lvert{\xi_{j_{1}}^{ \intercal}\widetilde{\Sigma}^{-1}\xi_{j}n_{j}\xi_{j}^{\intercal}\Sigma^{-1}\xi _{j_{2}}}\right\rvert}\,.$ $$

For all $j\neq j_{1},j_{2}$, from Lemma 8.7, with very high probability

$$ ${\left\lvert{\xi_{j_{1}}^{\intercal}\widetilde{\Sigma}^{-1}\xi_{j}n_{j}\xi_{j} ^{\intercal}\Sigma^{-1}\xi_{j_{2}}}\right\rvert}=\widetilde{O}\left({\eta\frac {d}{n^{2}\sqrt{n_{j_{1}}n_{j_{2}}}}\left(1+\frac{d^{2}n_{j}}{n^{2}}\right)} \right)\,.$ $$

Therefore, with very high probability

$$ $\sum_{j\neq j_{1},j_{2}}{\left\lvert{\xi_{j_{1}}^{\intercal}\widetilde{\Sigma} ^{-1}\xi_{j}n_{j}\xi_{j}^{\intercal}\Sigma^{-1}\xi_{j_{2}}}\right\rvert}= \widetilde{O}\left({\frac{md}{n^{2}\sqrt{n_{j_{1}}n_{j_{2}}}}+\frac{d^{3}}{n^{ 3}\sqrt{n_{j_{1}}n_{j_{2}}}}}\right)\times\eta=\widetilde{O}\left({\frac{\sqrt {d}}{\sqrt{n_{j_{1}}n_{j_{2}}}n}\eta}\right)\,,$ $$

where the last step utilized our assumptions that $m\leq\frac{n}{n_{\textup{min}}}\ll\frac{n}{\sqrt{d}}$ and that $n^{4}\gg d^{5}$.

When $j=j_{1}$, we have

$$ $\displaystyle{\left\lvert{\xi_{j_{1}}^{\intercal}\widetilde{\Sigma}^{-1}\xi_{j }n_{j}\xi_{j}^{\intercal}\Sigma^{-1}\xi_{j_{2}}}\right\rvert}$ $\displaystyle={\left\lvert{\xi_{j_{1}}^{\intercal}\widetilde{\Sigma}^{-1}\xi_{ j_{1}}n_{j_{1}}\xi_{j_{1}}^{\intercal}\Sigma^{-1}\xi_{j_{2}}}\right\rvert}= \widetilde{O}\left({\frac{d}{n\times n_{j_{1}}}\times n_{j_{1}}\times\frac{ \sqrt{d}}{\sqrt{n_{j_{1}}n_{j_{2}}}}\left(1+\frac{d\sqrt{n_{j_{1}}}}{n}\right) }\right)=$ $\displaystyle=\widetilde{O}\left({\frac{\sqrt{d}}{\sqrt{n_{j_{1}}n_{j_{2}}}} \left(1+\frac{d\sqrt{n_{j_{1}}}}{n}\right)}\right)\,.$ $$

Similarly, for $j=j_{2}$, we have

$$ $\displaystyle{\left\lvert{\xi_{j_{1}}^{\intercal}\widetilde{\Sigma}^{-1}\xi_{j }n_{j}\xi_{j}^{\intercal}\Sigma^{-1}\xi_{j_{2}}}\right\rvert}$ $\displaystyle={\left\lvert{\xi_{j_{1}}^{\intercal}\widetilde{\Sigma}^{-1}\xi_{ j_{2}}n_{j_{2}}\xi_{j_{2}}^{\intercal}\Sigma^{-1}\xi_{j_{2}}}\right\rvert}= \widetilde{O}\left({\frac{d}{n\times n_{j_{2}}}\times n_{j_{2}}\times\frac{ \sqrt{d}\eta}{\sqrt{n_{j_{1}}n_{j_{2}}}}\left(1+\frac{d\sqrt{n_{j_{1}}}}{n} \right)}\right)=$ $\displaystyle=o\left(\frac{\sqrt{d}\eta}{\sqrt{n_{j_{1}}n_{j_{2}}}}\left(1+ \frac{d\sqrt{n_{j_{1}}}}{n}\right)\right)\,.$ $$

Therefore, with very high probability

$$ $\eta=\frac{{\left\lvert{\xi_{j_{1}}^{\intercal}\widetilde{\Sigma}^{-1}\xi_{j_{ 2}}}\right\rvert}}{\frac{\sqrt{d}}{\sqrt{n_{j_{1}}n_{j_{2}}}n}\times\left(1+ \frac{d\sqrt{n_{j_{1}}}}{n}\right)}\leq\frac{{\left\lvert{\xi_{j_{1}}^{ \intercal}{\Sigma}^{-1}\xi_{j_{2}}}\right\rvert}+\sum_{j\in[m]}{\left\lvert{ \xi_{j_{1}}^{\intercal}\widetilde{\Sigma}^{-1}\xi_{j}n_{j}\xi_{j}^{\intercal} \Sigma^{-1}\xi_{j_{2}}}\right\rvert}}{\frac{\sqrt{d}}{\sqrt{n_{j_{1}}n_{j_{2}} }n}\times\left(1+\frac{d\sqrt{n_{j_{1}}}}{n}\right)}=\widetilde{O}\left({1} \right)+o(\eta)\,.$ $$

Therefore, with very high probability $\eta=\widetilde{O}\left({1}\right)$, proving the first portion of our claim.

Next, consider a term of the form $X_{i}^{\intercal}\widetilde{\Sigma}^{-1}\xi_{j}$.
To bound this term, we open $\widetilde{\Sigma}^{-1}$ again using equation (24).
Indeed, we have

$$ ${\left\lvert{X_{i}^{\intercal}\widetilde{\Sigma}^{-1}\xi_{j}}\right\rvert}\leq {\left\lvert{X_{i}^{\intercal}{\Sigma}^{-1}\xi_{j}}\right\rvert}+\sum_{j^{ \prime}}{\left\lvert{X_{i}^{\intercal}\Sigma^{-1}\xi_{j^{\prime}}n_{j^{\prime} }\xi_{j^{\prime}}^{\intercal}\widetilde{\Sigma}^{-1}\xi_{j}}\right\rvert}\,.$ $$

Define

$$ $\zeta_{i,j}=\xi_{j}-\frac{1}{n_{j}}X_{i}1_{i\in B_{j}}=\begin{cases}\xi_{j}- \frac{1}{n_{j}}X_{i}&i\in B_{j}\\ \xi_{j}&i\notin B_{j}\end{cases}$ $$

$\zeta_{i,j}$ over $B_{j}$ bucket had the $X_{i}$ sample been replaced with the 0 0 vector.
Note that with very high probability

$$ ${\left\lvert{X_{i}^{\intercal}{\Sigma}^{-1}\left(\xi_{j}-\zeta_{i,j}\right)} \right\rvert}\leq\frac{1_{i\in B_{j}}}{n_{j}}\left\|{X_{i}}\right\|^{2}\left\| {{\Sigma}^{-1}}\right\|=\widetilde{O}\left(\frac{d1_{i\in B_{j}}}{nn_{j}} \right)\,.$ $$

To proceed, we apply the Sherman-Morrison identity with $v=X_{i}$ and $A=\Sigma_{[n]\setminus\left\{{i}\right\}}=\Sigma-vv^{\intercal}$, to show that

$$ $\Sigma^{-1}=A^{-1}-\frac{A^{-1}vv^{\intercal}A^{-1}}{1-v^{\intercal}A^{-1}v}\,.$ $$

Moreover, because both $A$ and $\zeta_{i,j}$ are independent of the $i$th sample and $\mathcal{X}$ is well-behaved, with very high probability

$$ ${\left\lvert{X_{i}^{\intercal}A^{-1}\zeta_{i,j}}\right\rvert}=\widetilde{O} \left(\left\|{A^{-1}\zeta_{i,j}}\right\|\right)=\widetilde{O}\left(\frac{\sqrt {d}}{\sqrt{n_{j}}n}\right)$ $$

Combining this with the Sherman-Morrison formula, we have

$$ ${\left\lvert{X_{i}^{\intercal}\Sigma^{-1}\zeta_{i,j}}\right\rvert}=\left(1- \frac{X_{i}^{\intercal}A^{-1}X_{i}}{1-X_{i}^{\intercal}A^{-1}X_{i}}\right){ \left\lvert{X_{i}^{\intercal}A^{-1}\zeta_{i,j}}\right\rvert}=\widetilde{O} \left(\frac{\sqrt{d}}{\sqrt{n_{j}}n}\right)\,.$ $$

Therefore, for all $j\in[m]$, with very high probability

$$ ${\left\lvert{X_{i}^{\intercal}\Sigma^{-1}\xi_{j}}\right\rvert}=\widetilde{O} \left(\frac{\sqrt{d}}{\sqrt{n_{j}}n}+\frac{1_{i\in B_{j}}d}{n_{j}n}\right)\,.$ $$

Therefore, with very high probability

$$ $\displaystyle\sum_{j^{\prime}\in[m]}{\left\lvert{X_{i}^{\intercal}\Sigma^{-1} \xi_{j^{\prime}}n_{j^{\prime}}\xi_{j^{\prime}}^{\intercal}\widetilde{\Sigma}^{ -1}\xi_{j}}\right\rvert}\leq$ $\displaystyle\leq\sum_{j^{\prime}\in[m]}\widetilde{O}\left(\left(\frac{\sqrt{d }}{n\sqrt{n_{j^{\prime}}}}+\frac{d1_{i\in B_{j^{\prime}}}}{nn_{j^{\prime}}} \right)\times n_{j^{\prime}}\times\left(\frac{\sqrt{d}}{\sqrt{n_{j}n_{j^{ \prime}}}n}\left(1+\frac{d\sqrt{n_{j^{\prime}}}}{n}\right)+\frac{1_{j=j^{ \prime}}d}{\sqrt{n_{j}n_{j^{\prime}}}n}\right)\right)=$ $\displaystyle=\sum_{j^{\prime}\in[m]}\widetilde{O}\left(\frac{d}{n^{2}\sqrt{n_ {j}}}+\frac{d^{2}\sqrt{n_{j^{\prime}}}}{n^{3}\sqrt{n_{j}}}+\frac{d^{3/2}1_{j=j ^{\prime}}}{\sqrt{n_{j}}n^{2}}+\frac{d^{3/2}1_{i\in B_{j^{\prime}}}}{\sqrt{n_{ j}n_{j^{\prime}}}n^{2}}\times\left(1+\frac{d\sqrt{n_{j^{\prime}}}}{n}\right)+ \frac{d^{2}1_{i\in B_{j}}1_{j=j^{\prime}}}{n_{j}n^{2}}\right)=$ $\displaystyle=\widetilde{O}\left({\frac{\sqrt{d}}{\sqrt{n_{j}}n}+\frac{1_{i\in B _{j}}d}{n_{j}n}}\right)+\widetilde{O}\left({\frac{md}{n^{2}\sqrt{n_{j}}}+\frac {\sqrt{mn}d^{2}}{n^{3}\sqrt{n_{j}}}+\frac{d^{3/2}}{\sqrt{n_{j^{\prime}}}n^{2}} +\frac{d1_{i\in B_{j}}}{nn_{j}}}\right)=\widetilde{O}\left(\frac{\sqrt{d}}{ \sqrt{n_{j}}n}+\frac{1_{i\in B_{j}}d}{nn_{j}}\right)\,,$ $$

where the last step utilizes our assumptions that $m=\widetilde{O}\left(\frac{n}{\sqrt{d}}\right)$ and that $d^{5}=\widetilde{O}\left({n^{4}}\right)$.

Finally, let $i_{1}\neq i_{2}$ be the indices of two samples.
From Claim 7.9, we have that with very high probability ${\left\lvert{X_{i_{1}}^{\intercal}\Sigma^{-1}X_{i_{2}}}\right\rvert}=
\widetilde{O}\left(\frac{\sqrt{d}}{n}\right)$.
Opening $\widetilde{\Sigma}^{-1}$ again using equation (24), we see that with very high probability

$$ $\displaystyle{\left\lvert{X_{i_{1}}^{\intercal}\widetilde{\Sigma}^{-1}X_{i_{2} }}\right\rvert}$ $\displaystyle\leq{\left\lvert{X_{i_{1}}^{\intercal}\Sigma^{-1}X_{i_{2}}}\right \rvert}+\sum_{j\in[m]}{\left\lvert{X_{i_{1}}^{\intercal}\Sigma^{-1}\xi_{j}n_{j }\xi_{j}^{\intercal}\widetilde{\Sigma}^{-1}X_{i_{2}}}\right\rvert}=$ $\displaystyle=\widetilde{O}\left(\frac{\sqrt{d}}{n}\right)+\sum_{j\in[m]} \widetilde{O}\left(\left(\frac{\sqrt{d}}{\sqrt{n_{j}}n}+\frac{1_{i_{1}\in B_{j }}d}{n_{j}n}\right)\times\left(\frac{\sqrt{d}}{\sqrt{n_{j}}n}+\frac{1_{i_{2} \in B_{j}}d}{n_{j}n}\right)\times n_{j}\right)=$ $\displaystyle=\widetilde{O}\left(\frac{\sqrt{d}}{n}\right)+\widetilde{O}\left( \frac{md}{n^{2}}+\frac{d^{2}}{n^{2}n_{\textup{min}}}\right)=\widetilde{O}\left (\frac{\sqrt{d}}{n}\right)\,.$ $$

where $n_{\textup{min}}\stackrel{{\scriptstyle\textup{def}}}{{=}}\min_{j\in[m]}\left
\{{n_{j}}\right\}=\widetilde{\Omega}\left(\sqrt{d}\right)$.

∎

###### Claim 8.6 .

###### Lemma 8.7 .

###### Proof of Claim 8.6 .

#### 8.2.5 Projection on the e 𝑒 e italic_e Axis

For the next property of well-behaved regressions, we will want to show that with very high probability $e^{\intercal}\Sigma^{-1}X_{i}$ is bounded for all $i$.
Note that from the exponential decay assumption due to $\mathcal{X}$ being well-behaved would suffice to give a good bound on $e^{\intercal}{\Sigma}_{\textup{gt}}^{-1}X_{i}$ (where ${\Sigma}_{\textup{gt}}=nI$), but as before, the challenge will be to show that $\Sigma^{-1}$ doesn’t rotate $X_{i}$ onto $e$.

###### Claim 8.8 .

With very high probability

$$ $\forall i\in[n]\;\;\;\;\left\langle{\widetilde{\Sigma}^{-1}\widetilde{X}_{i}}, {e}\right\rangle=\frac{e^{\intercal}X_{i}}{n}\pm o\left(\frac{1}{n}\right)\,.$ $$

In particular, with very high probability

$$ $\max_{i\in[n]}{\left\lvert{\left\langle{\widetilde{\Sigma}^{-1}\widetilde{X}_{ i}},{e}\right\rangle}\right\rvert}=\widetilde{O}\left(\frac{1}{n}\right)=\sqrt {\widetilde{O}\left(\frac{1}{n}\right)\times e^{\intercal}\widetilde{\Sigma}^{ -1}e}$ $$

Our proof of Claim 8.8, will make use of the following Lemma 8.9 that will also be proved in Section 8.3.

###### Lemma 8.9 .

For any $j\in[m]$, with very high probability,

$$ ${\left\lvert{e^{\intercal}\Sigma^{-1}\xi_{j}}\right\rvert}=\widetilde{O}\left( \frac{1}{n\sqrt{n_{j}}}\times\left(1+\frac{d\sqrt{n_{j}}}{n}\right)\right)\,.$ $$

###### Proof of Claim 8.8 .

As in the proof of Claim 8.6, we will use the matrix identity in equation (24) to show that

$$ $\widetilde{\Sigma}^{-1}=\Sigma^{-1}+\Sigma^{-1}C\widetilde{\Sigma}^{-1}\,,$ $$

where $C=\Sigma-\widetilde{\Sigma}=\sum_{j\in[m]}n_{j}\xi_{j}\xi_{j}^{\intercal}$.

Define

$$ $\eta\stackrel{{\scriptstyle\textup{def}}}{{=}}\max_{j\in[m]}\left\{{{\left \lvert{e^{\intercal}\widetilde{\Sigma}^{-1}\xi_{j}}\right\rvert}}\right\}$ $$

Using Lemma 8.9, we see that with very high probability

$$ $\eta={\left\lvert{e^{\intercal}\widetilde{\Sigma}^{-1}\xi_{j}}\right\rvert} \leq{\left\lvert{e^{\intercal}{\Sigma}^{-1}\xi_{j}}\right\rvert}+\sum_{j^{ \prime}}{\left\lvert{e^{\intercal}\Sigma^{-1}\xi_{j^{\prime}}n_{j^{\prime}}\xi _{j^{\prime}}^{\intercal}\widetilde{\Sigma}^{-1}\xi_{j}}\right\rvert}= \widetilde{O}\left(\frac{1}{n}\times\left(\frac{m\sqrt{d}+d}{n\sqrt{n_{j}}} \right)\right)=\widetilde{O}\left(\frac{1}{n\sqrt{n_{j}}}\right)\,.$ $$

Opening $\widetilde{\Sigma}^{-1}$ again, we have

$$ ${\left\lvert{e^{\intercal}\widetilde{\Sigma}^{-1}X_{i}-e^{\intercal}\Sigma^{-1 }X_{i}}\right\rvert}=\sum_{j\in[m]}{\left\lvert{e^{\intercal}\widetilde{\Sigma }^{-1}\xi_{j}n_{j}\xi_{j}^{\intercal}\Sigma^{-1}X_{i}}\right\rvert}$ $$

From Claim 7.10, we know that with very high probability

$$ ${\left\lvert{e^{\intercal}\Sigma^{-1}X_{i}-\frac{e^{\intercal}X_{i}}{n}}\right \rvert}=o\left(\frac{1}{n}\right)\,.$ $$

Therefore, from Claim 8.6, with very high probability

$$ $\displaystyle{\left\lvert{e^{\intercal}\widetilde{\Sigma}^{-1}X_{i}-e^{ \intercal}\Sigma^{-1}X_{i}}\right\rvert}$ $\displaystyle\leq\sum_{j\in[m]}{\left\lvert{e^{\intercal}\widetilde{\Sigma}^{- 1}\xi_{j}n_{j}\xi_{j}^{\intercal}\Sigma^{-1}X_{i}}\right\rvert}=$ $\displaystyle=\sum_{j\in[m]}\widetilde{O}\left(\frac{1}{n\sqrt{n_{j}}}\times n _{j}\times\left(\frac{\sqrt{d}}{n\sqrt{n_{j}}}+\frac{d1_{i\in B_{j}}}{n{n_{j}} }\right)\right)=$ $\displaystyle=\widetilde{O}\left(\frac{\sqrt{d}m}{n^{2}}+\frac{d}{n^{2}\sqrt{n _{\textup{min}}}}\right)=o\left(\frac{1}{n}\right)\,.$ $$

Altogether, we have that with very high probability

$$ ${\left\lvert{e^{\intercal}\widetilde{\Sigma}^{-1}\widetilde{X}_{i}-\frac{e^{ \intercal}X_{i}}{n}}\right\rvert}\leq{\left\lvert{e^{\intercal}\widetilde{ \Sigma}^{-1}\xi_{j(i)}}\right\rvert}+{\left\lvert{e^{\intercal}\widetilde{ \Sigma}^{-1}X_{i}-e^{\intercal}\Sigma^{-1}X_{i}}\right\rvert}+{\left\lvert{e^{ \intercal}\Sigma^{-1}X_{i}-\frac{e^{\intercal}X_{i}}{n}}\right\rvert}=o\left( \frac{1}{n}\right)\,.$ $$

∎

###### Claim 8.8 .

###### Lemma 8.9 .

###### Proof of Claim 8.8 .

#### 8.2.6 Bounded Residuals

Recall that we assumed that the samples of our regression were drawn from a ground truth linear model plus an iid normally distributed error.
Therefore, in order to bound the empirical residuals, it suffices to show that they are close to the ground truth residuals with very high probability.

In particular, in the setting of Theorem 8.1, we assumed that

$$ $Y_{i}\sim\underbrace{\mu_{j(i)}}_{\textup{categorical contribution}}+ \underbrace{\left\langle{X_{i}},{{\beta}_{\textup{gt}}}\right\rangle}_{\textup {continuous contribution}}+\underbrace{\mathcal{N}\left(0,1\right)}_{\textup{ error}}$ $$

Define the ground truth residuals to be

$$ ${R}_{\textup{gt}}=\left(Y_{i}-\mu_{j(i)}-\left\langle{X_{i}},{{\beta}_{\textup {gt}}}\right\rangle\right)_{i\in[n]}\,,$ $$

and the empirical residuals to be

$$ $R=\widetilde{Y}-\widetilde{X}\beta=\widetilde{Y}-\widetilde{X}\widetilde{ \Sigma}^{-1}\widetilde{X}^{\intercal}\widetilde{Y}\,.$ $$

###### Claim 8.10 .

With very high probability,

$$ $\left\|{R-{R}_{\textup{gt}}}\right\|_{\infty}=o(1)\,.$ $$

###### Proof of Claim 8.10 .

Denote

$$ $\mu^{*}_{j}\stackrel{{\scriptstyle\textup{def}}}{{=}}\frac{1}{n_{j}}\sum_{i\in B _{j}}Y_{i}=\mu_{j}+\xi_{j}^{\intercal}{\beta}_{\textup{gt}}+\frac{1}{n_{j}} \sum_{i\in B_{j}}\left({R}_{\textup{gt}}\right)_{i}$ $$

From the definitions of the residuals and the ground truth residuals, we have

$$ $\displaystyle R_{i}-\left({R}_{\textup{gt}}\right)_{i}$ $\displaystyle=\left(Y_{i}-\mu^{*}_{j(i)}\right)-\left(X_{i}-\xi_{j(i)}\right)^ {\intercal}\beta-\left(Y_{i}-\mu_{j(i)}-X_{i}^{\intercal}{\beta}_{\textup{gt}} \right)=$ $\displaystyle=-\xi_{j(i)}^{\intercal}{\beta}_{\textup{gt}}-\frac{1}{n_{j(i)}} \sum_{i^{\prime}\in B_{j(i)}}\left({R}_{\textup{gt}}\right)_{i^{\prime}}-X_{i} ^{\intercal}\left(\beta-{\beta}_{\textup{gt}}\right)+\xi_{j(i)}^{\intercal}\beta=$ $\displaystyle=-\frac{1}{n_{j(i)}}\sum_{i^{\prime}\in B_{j(i)}}\left({R}_{ \textup{gt}}\right)_{i^{\prime}}-\widetilde{X}_{i}^{\intercal}\left(\beta-{ \beta}_{\textup{gt}}\right)\,.$ $$

Expanding on the difference between the results of the reaveraged OLS and the ground truth linear model, we have

$$ $\displaystyle\beta-{\beta}_{\textup{gt}}$ $\displaystyle=\widetilde{\Sigma}^{-1}\widetilde{X}^{\intercal}\left(\widetilde {Y}-\widetilde{X}{\beta}_{\textup{gt}}\right)=\widetilde{\Sigma}^{-1} \widetilde{X}^{\intercal}\left(Y_{i}-\mu^{*}_{j(i)}-\left(Y_{i}-\mu_{j(i)}- \left({R}_{\textup{gt}}\right)_{i}\right)+\xi_{j(i)}^{\intercal}{\beta}_{ \textup{gt}}\right)_{i\in[n]}=$ $\displaystyle=\widetilde{\Sigma}^{-1}\widetilde{X}^{\intercal}\left(\left({R}_ {\textup{gt}}\right)_{i}-\xi_{j(i)}^{\intercal}{\beta}_{\textup{gt}}-\frac{1}{ n_{j(i)}}\sum_{i^{\prime}\in B_{j(i)}}\left({R}_{\textup{gt}}\right)_{i^{ \prime}}+\xi_{j(i)}^{\intercal}{\beta}_{\textup{gt}}\right)_{i\in[n]}=$ $\displaystyle=\widetilde{\Sigma}^{-1}\widetilde{X}^{\intercal}\left(\left({R}_ {\textup{gt}}\right)_{i}-\frac{1}{n_{j(i)}}\sum_{i^{\prime}\in B_{j(i)}}\left( {R}_{\textup{gt}}\right)_{i^{\prime}}\right)_{i\in[n]}\,.$ $$

Therefore,

$$ $R_{i}-\left({R}_{\textup{gt}}\right)_{i}=-\frac{1}{n_{j(i)}}\sum_{i^{\prime} \in B_{j(i)}}\left({R}_{\textup{gt}}\right)_{i^{\prime}}-\widetilde{X}_{i}^{ \intercal}\widetilde{\Sigma}^{-1}\widetilde{X}^{\intercal}\left(\left({R}_{ \textup{gt}}\right)_{i^{*}}-\frac{1}{n_{j({i^{*}})}}\sum_{i^{\prime}\in B_{j({ i^{*}})}}\left({R}_{\textup{gt}}\right)_{i^{\prime}}\right)_{{i^{*}}\in[n]}\,.$ $$

Summing over the contributions of the second term, we have

$$ $\displaystyle\sum_{{i^{*}}\in[n]}\widetilde{X}_{i}^{\intercal}\widetilde{ \Sigma}^{-1}\widetilde{X}_{i^{*}}\left(\left({R}_{\textup{gt}}\right)_{i^{*}}- \frac{1}{n_{j({i^{*}})}}\sum_{i^{\prime}\in B_{j({i^{*}})}}\left({R}_{\textup{ gt}}\right)_{i^{\prime}}\right)=\sum_{{i^{*}}\in[n]}\widetilde{X}_{i}^{ \intercal}\widetilde{\Sigma}^{-1}\left(\widetilde{X}_{i^{*}}-\xi_{j(i^{*})} \right)\left({R}_{\textup{gt}}\right)_{i^{*}}\sim$ $\displaystyle\;\;\;\;\;\;\;\;\;\;\sim\mathcal{N}\left(0,\sum_{i^{*}\in[n]} \left(\widetilde{X}_{i}^{\intercal}\widetilde{\Sigma}^{-1}\left(\widetilde{X}_ {i^{*}}-\xi_{j(i^{*})}\right)\right)^{2}\right)=\mathcal{N}\left(0,\widetilde{ O}\left({\frac{d}{n}+\frac{d}{nn_{\textup{min}}}+\frac{d^{2}}{n^{2}}}\right) \right)\,,$ $$

so with very high probability this contribution is bounded in absolute value by $o(1)$.

Similarly,

$$ $\frac{1}{n_{j(i)}}\sum_{i^{\prime}\in B_{j(i)}}\left({R}_{\textup{gt}}\right)_ {i^{\prime}}\sim\mathcal{N}\left(0,\frac{1}{n_{j(i)}}\right)=\mathcal{N}\left( 0,O\left(\frac{1}{n^{\varepsilon}}\right)\right)\,,$ $$

so with very high probability this term is also $o(1)$, concluding our proof of Claim 8.10.

∎

###### Claim 8.10 .

###### Proof of Claim 8.10 .

#### 8.2.7 Large Influence Scores

For the final step of our proof of Claim 8.2, we will show that with very high probability, there are many samples in the regression that have relatively high AMIP influence scores.

###### Claim 8.11 .

Let $\alpha_{i}=e^{\intercal}\widetilde{\Sigma}^{-1}\widetilde{X}_{i}R_{i}$ denote the AMIP influence score of the $i$th sample.
Then with very high probability, there are is a set $T_{0}\subseteq[n]$ of least $k_{0}=\widetilde{\Omega}(n^{1-\varepsilon})$ “influential samples” – i.e., such that

$$ $\forall i\in T\;\;\alpha_{i}\geq\left(\frac{\sqrt{\log(n)}}{n}\right)=\omega \left(\frac{1}{n}\right)$ $$

Proving Claim 8.11 will also conclude our proof of Claim 8.2, as this will show that all conditions required for a regression to be well-behaved are fulfilled with very high probability.

###### Proof of Claim 8.11 .

The proof of Claim 8.11 will follow the same approach as our proof of its ACRE counterpart – Claim 7.12.
Indeed, we first note that from the definition of $\widetilde{\Sigma}$ and Claim 8.4, with very high probability,

$$ $\sum_{i\in[n]}\left(e^{\intercal}\widetilde{\Sigma}^{-1}\widetilde{X}_{i} \right)^{2}=\sum_{i\in[n]}e^{\intercal}\widetilde{\Sigma}^{-1}\widetilde{X}_{i }\widetilde{X_{i}}^{\intercal}\widetilde{\Sigma}^{-1}e=e^{\intercal}\widetilde {\Sigma}e=\frac{1\pm o(1)}{n}\,.$ $$

Moreover, from Claim 8.8, with very high probability, the contribution of each individual sample to this sum is at most $\widetilde{O}\left({\frac{1}{n}}\right)$-fraction of this total.
Therefore, with very high probability

$$ ${\left\lvert{\left\{i\middle|{\left\lvert{e^{\intercal}\widetilde{\Sigma}^{-1} \widetilde{X}_{i}}\right\rvert}\geq\frac{1}{2n}\right\}}\right\rvert}= \widetilde{\Omega}\left(n\right)\,.$ $$

As in the proof of Claim 7.12, we note that Claim 8.10 guarantees that with very high probability $\forall i\;\;{\left\lvert{R_{i}-\left({R}_{\textup{gt}}\right)_{i}}\right
\rvert}=o(1)<\frac{1}{20}$.
Moreover, because the ground truth residuals $\left({R}_{\textup{gt}}\right)_{i}$ are normally distributed independently of anything else, it follows that so are $\rho_{i}\stackrel{{\scriptstyle\textup{def}}}{{=}}\textup{sign}\left({e^{
\intercal}\widetilde{\Sigma}^{-1}\widetilde{X}_{i}}\right)\times\left({R}_{
\textup{gt}}\right)_{i}$.

Therefore, with very high probability, at least $\widetilde{\Omega}\left(n^{1-\varepsilon}\right)$ of the samples such that ${\left\lvert{e^{\intercal}\widetilde{\Sigma}^{-1}\widetilde{X}_{i}}\right
\rvert}\geq\frac{1}{2n}$ have $\rho_{i}\geq\Omega(\sqrt{\log(n)})$ (this is because we can set the constants in the $\Omega$ to be such that the probability of each of $\rho_{i}$ being above this threshold is $\gg n^{-\varepsilon}$, allowing us to apply Hoeffding on the $\widetilde{\Omega}\left(n\right)$ iid $\rho_{i}$).

Therefore, for this set of samples it holds that

$$ $\alpha_{i}=e^{\intercal}\widetilde{\Sigma}^{-1}\widetilde{X}_{i}R_{i}=\left(1 \pm o(1)\right){\left\lvert{e^{\intercal}\widetilde{\Sigma}^{-1}\widetilde{X}_ {i}}\right\rvert}\rho_{i}=\Omega\left(\frac{\sqrt{\log(n)}}{n}\right)\,.$ $$

∎

###### Claim 8.11 .

###### Proof of Claim 8.11 .

### 8.3 Proof of Lemmas 8.7 and 8.9

Recall the lemma we wish to prove:

###### Lemma 8.12 .

Let $X_{1},\ldots,X_{n}\sim\mathcal{X}$ be $n$ iid samples of a well-behaved distribution $\mathcal{X}$ with covariance $I_{d\times d}$.
Let $\xi=\frac{1}{k}\sum_{i\in[k]}X_{i}$ be the empirical average over the first $k<0.49n$ of these samples, and let $v\in\mathbb{R}^{d}$ be any vector that is independent of these first $k$ samples (but may depend on the other samples).
Finally, denote by $\Sigma\stackrel{{\scriptstyle\textup{def}}}{{=}}\sum_{i\in[n]}X_{i}X_{i}^{\intercal}$ the unnormalized empirical second moment of these samples.

Then, with very high probability,

$$ ${\left\lvert{\xi^{\intercal}\Sigma^{-1}v}\right\rvert}=\widetilde{O}\left( \frac{\left\|{v}\right\|}{\sqrt{k}n}\times\left(1+\frac{d\sqrt{k}}{n}\right) \right)\,.$ $$

Clearly, Lemma 8.12 implies Lemma 8.7 (set $\xi=\xi_{j_{1}}$ and $v=\xi_{j_{2}}$) and Lemma 8.9 (set $\xi=\xi_{j}$ and $v=e$).

###### Lemma 8.12 .

#### 8.3.1 Proof Sketch

We will split the proof of Lemma 8.7 into 3 main steps.
To do this, let $S=\sum_{i=k+1}^{n}X_{i}X_{i}^{\intercal}=\Sigma_{[n]\setminus[k]}$ be the contributions to the empirical covariance due to samples not amongst the first $k$, and let $C=\Sigma-S$ be the contributions from within the bucket.
From a standard analysis, so long as $C\prec S$, we have

$$ $\Sigma^{-1}=S^{-1}-\Sigma^{-1}CS^{-1}=S^{-1}-S^{-1}CS^{-1}-\left(\Sigma^{-1}-S ^{-1}\right)CS^{-1}\,,$ $$

and these components will correspond to the main components of our analysis.
We will show that the inequality $C\prec S$ holds with very high probability, before proving the following claims over the rest of this section:

###### Claim 8.13 .

With very high probability

$$ ${\left\lvert{\xi^{\intercal}S^{-1}v}\right\rvert}=\widetilde{O}\left(\frac{ \left\|{v}\right\|}{\sqrt{k}\times n}\right)$ $$

###### Claim 8.14 .

Let $w\in\mathbb{R}^{d}$ be any vector that does not depend on the first $k$ samples.

Then, with very high probability,

$$ ${\left\lvert{\xi^{\intercal}S^{-1}Cw}\right\rvert}=\widetilde{O}\left(\frac{ \left\|{w}\right\|}{\sqrt{k}}\left(1+\frac{d\sqrt{k}}{n}\right)\right)$ $$

###### Claim 8.15 .

Let $w\in\mathbb{R}^{d}$ be any vector that does not depend on the first $k$ samples.
Then, with very high probability,

$$ ${\left\lvert{\xi^{\intercal}\left(S^{-1}-\Sigma^{-1}\right)Cw}\right\rvert}= \widetilde{O}\left(\frac{\left\|{w}\right\|}{\sqrt{k}}\left(1+\frac{d\sqrt{k}} {n}\right)\right)$ $$

The proofs of Claims 8.13, 8.14 and 8.15 will grow progressively more complex and each claim will build on ideas from the previous one.
Throughout the latter two, the key challenge will be to deal with cases where both the $\zeta\stackrel{{\scriptstyle\textup{def}}}{{=}}k\xi$ term and the $C$ or $\Sigma$ multiplicand may depend on the same samples.

In these cases, we will proceed by applying a sort of divide-and-conquer approach by splitting the bucket into two subsets.
For instance, instead of analyzing $\zeta^{\intercal}S^{-1}Cw$ directly, we will split the samples in the bucket into two subsets and track their contributions to both $\zeta=\zeta_{0}+\zeta_{1}$ and $C=C_{0}+C_{1}$.

$$ $\displaystyle\zeta^{\intercal}S^{-1}Cw$ $\displaystyle=\left(\zeta_{0}+\zeta_{1}\right)S^{-1}\left(C_{0}+C_{1}\right)w=$ (25) $\displaystyle=\underbrace{\zeta_{0}^{\intercal}S^{-1}C_{0}w+\zeta_{1}^{ \intercal}S^{-1}C_{1}w}_{\textup{diagonal terms}}+\underbrace{\zeta_{1}^{ \intercal}S^{-1}C_{0}w+\zeta_{0}^{\intercal}S^{-1}C_{1}w}_{\textup{off- diagonal terms}}$ $$

The “off-diagonal” subsets will be relatively simple to bound as they contain inner products of independent vectors in $\mathbb{R}^{d}$ (and this independence will give us a $1/\sqrt{d}$ scaling to their inner product), and the diagonal elements will be split again recursively.
This will also leave us with a large number of “single sample” diagonal elements of the form

$$ $X_{i}^{\intercal}S^{-1}X_{i}X_{i}^{\intercal}w$ $$

These single-sample terms will no longer enjoy the same $1/\sqrt{d}$ scaling the other terms gain due to independence, but will instead gain a sort of $1/k$ scaling, because instead of having $k^{2}$ sample-times-sample contributions in

$$ $\zeta^{\intercal}S^{-1}C=\sum_{i,i^{\prime}\in[k]}X_{i}^{\intercal}AX_{i^{ \prime}}X_{i^{\prime}}^{\intercal}w\,,$ $$

we have only $k$ terms in the sum

$$ $\sum_{i\in[k]}X_{i}^{\intercal}AX_{i}X_{i}^{\intercal}w$ $$

Finally, in Claim 8.15, we will have 2 types of diagonal vs off-diagonal splits.
The first will be to track the cases where $\zeta$ and $C$ may depend on the same samples and will be very similar to our analysis of Claim 8.14.
The second and much more difficult of the two will be dealing with dependencies between $\Sigma$ and $\zeta$.

###### Claim 8.13 .

###### Claim 8.14 .

###### Claim 8.15 .

#### 8.3.2 Setup and Proof of Claim 8.13

Recall our assumption from Lemma 8.7 that $k\leq 0.49n$.
Therefore $n-k\geq 0.51n>k+\Omega(n)$.

Let $S=\sum_{i\in[n]\setminus[k]}X_{i}X_{i}^{\intercal}=\Sigma_{[n]\setminus[k]}$.
From Claim 7.7, we have that with very high probability

$$ $\left\|{S-\left(n-k\right)I}\right\|=o(n)\,.$ $$

Combined with Lemma 7.6, which shows that $\xi$ is well-behaved, we have that with high probability

$$ ${\left\lvert{\xi^{\intercal}S^{-1}v}\right\rvert}=\widetilde{O}\left(\frac{1}{ \sqrt{k}}\times\left\|{S^{-1}v}\right\|\right)=\widetilde{O}\left(\frac{1}{ \sqrt{k}}\times\max\left\{{\lambda\left(S^{-1}\right)}\right\}\times\left\|{v} \right\|\right)=\widetilde{O}\left(\frac{\left\|{v}\right\|}{\sqrt{k}\times n}\right)$ $$

Our goal for the rest of the proof will be to show a similar bound for $\xi^{\intercal}\Sigma^{-1}v$.
Let $S$ and $C=\Sigma-S$ be our “main” and “correction” terms.
Applying Claim 7.7 again, we have that with very high probability $C\preceq kI+o(n)\prec(n-k)I-o(n)\preceq S$, so both $S$ and $I\pm S^{-1}C$ are invertible.
We have

$$ $\Sigma^{-1}=\left(S+C\right)^{-1}=\left(I+S^{-1}C\right)^{-1}S^{-1}=S^{-1}- \left(I+S^{-1}C\right)^{-1}S^{-1}CS^{-1}=S^{-1}-\Sigma^{-1}CS^{-1}$ $$

Therefore, it remains to bound $\xi^{\intercal}\Sigma^{-1}CS^{-1}v=\xi^{\intercal}\Sigma^{-1}Cw$ (where $w=S^{-1}v$) in absolute value.
To do this, we once again use the intuition that in some sense $S\approx\Sigma$, and first bound $\xi^{\intercal}S^{-1}CS^{-1}v$.
We will then slowly break down the difference between $\xi^{\intercal}S^{-1}CS^{-1}v$ and $\xi_{j_{1}}^{\intercal}\Sigma^{-1}CS^{-1}v$ into a series of corrections and bound each of these corrections in absolute value.

#### 8.3.3 Proof of Claim 8.14

Before proving Claim 8.14, we prove a lemma that will help us in our analysis.

###### Lemma 8.16 .

With very high probability

$$ $\left\|{Cw}\right\|=\widetilde{O}\left(\left(k+\sqrt{kd}\right)\left\|{w} \right\|\right)\,.$ $$

###### Proof of Lemma 8.16 .

At first glance, it might seem like Lemma 8.16 should follow immediately from Claim 7.7, but this is only true when $k=\widetilde{\Omega}\left(d\right)$, and we will want to apply Lemma 8.16 even when $k\ll d$.

We prove Lemma 8.16 using the approximate matrix Bernstein inequality (Lemma 7.5).
Indeed,

$$ $Cw=\sum_{i\in[k]}X_{i}X_{i}^{\intercal}w=\sum_{i\in[k]}v_{i}\,,$ $$

can be written as the sum of $k$ iid vectors $v_{i}\stackrel{{\scriptstyle\textup{def}}}{{=}}X_{i}X_{i}^{\intercal}w$.

From the fact that the $X_{i}$s are well behaved and independent of $w$, we have that with very high probability,

$$ $\left\|{v_{i}}\right\|\leq\sqrt{d}\left\|{w}\right\|\,.$ $$

Moreover, from our assumption that $\mathcal{X}$ has covariance identity, we have that

$$ ${\mathop{\mathbb{E}}{\left[{v_{i}}\right]}}={\mathop{\mathbb{E}}{\left[{X_{i}X _{i}^{\intercal}}\right]}}w=Iw=w\,.$ $$

Finally, from Claim 7.7, we have that

$$ ${\mathop{\mathbb{E}}{\left[{\left\|{v_{i}}\right\|^{2}}\right]}}=w^{\intercal} {\mathop{\mathbb{E}}{\left[{\left\|{X_{i}}\right\|^{2}X_{i}X_{i}^{\intercal}} \right]}}w=\widetilde{O}\left(d\left\|{w}\right\|^{2}\right)\,.$ $$

Therefore, applying the approximate matrix Bernstein inequality on the standard embedding matrix

$$ $V_{i}=\begin{pmatrix}0&v_{i}^{\intercal}\\ v_{i}&0\end{pmatrix}\in\mathbb{R}^{(d+1)\times(d+1)}\,,$ $$

yields the desired result.
∎

We now return to the proof of of Claim 8.14.

###### Proof of Claim 8.14 .

Let $\zeta\stackrel{{\scriptstyle\textup{def}}}{{=}}k\xi=\sum_{i\in[k]}X_{i}$.
Using this notation, we have

$$ $\xi^{\intercal}S^{-1}Cw=\frac{1}{k}\zeta^{\intercal}S^{-1}Cw$ (26) $$

We will bound the RHS of equation (26) by breaking the contributions to $\zeta$ and to $C$ into smaller and smaller subsets of the bucket $[k]$.

We will split the contributions from the first $k$ samples into “clusters” based on their index modulus some number, and the assumption above is just to ensure that the subsets are of roughly equal size.

Indeed, for any string of bits $a=(a_{0},\ldots,a_{t})\in\left\{{0,1}\right\}^{*}$, define the $a$th *cluster* of samples to be the set of indices from $[k]$ whose bitwise representation ends with the string $a$:

$$ $\mathcal{C}_{a}\stackrel{{\scriptstyle\textup{def}}}{{=}}\left\{i\in[k]\middle |i\mod 2^{t+1}=a_{0}+a_{1}2^{1}+\cdots a_{t}2^{t}\right\}$ $$

If $\epsilon$ is the empty string, then $\mathcal{C}_{\epsilon}=[k]$, and for all $a$, we have

$$ $\mathcal{C}_{a}=\mathcal{C}_{a0}\sqcup\mathcal{C}_{a1}$ $$

Similarly, we may split the contributions of samples in $[k]$ to $C$ and $\zeta$ based on their cluster

$$ $\displaystyle C_{a}\stackrel{{\scriptstyle\textup{def}}}{{=}}\sum_{i\in \mathcal{C}_{a}}X_{i}X_{i}^{\intercal}$ (27) $\displaystyle\zeta_{a}\stackrel{{\scriptstyle\textup{def}}}{{=}}\sum_{i\in \mathcal{C}_{a}}X_{i}$ $$

Moreover, we have the property that $C_{a}=C_{a0}+C_{a1}$ and $\zeta_{a}=\zeta_{a0}+\zeta_{a1}$.
Using this property, we may begin to split the RHS of equation (26) to smaller components

$$ $\displaystyle\zeta^{\intercal}S^{-1}Cw$ $\displaystyle=\left(\zeta_{0}+\zeta_{1}\right)A\left(C_{0}+C_{1}\right)w=$ (28) $\displaystyle=\underbrace{\zeta_{0}^{\intercal}S^{-1}C_{0}w+\zeta_{1}^{ \intercal}S^{-1}C_{1}w}_{\textup{diagonal terms}}+\underbrace{\zeta_{1}^{ \intercal}S^{-1}C_{0}w+\zeta_{0}^{\intercal}S^{-1}C_{1}w}_{\textup{off- diagonal terms}}$ $$

We split the terms in the RHS of equation (28) into “diagonal” terms which correspond to the contributions where the $C$ term and the $\zeta$ term correspond to the same samples and “off-diagonal” terms where $\zeta$ and $C$ depend on disjoint sets of samples.

To bound the contribution of the off-diagonal terms, we note that for any bitstring $a=(a_{0},\ldots,a_{t})$, it holds that $\zeta_{a}$ is well-behaved (as it is the sum over ${\left\lvert{\mathcal{C}_{a}}\right\rvert}$ iid samples from $\mathcal{X}$) and has covariance ${\left\lvert{\mathcal{C}_{a}}\right\rvert}I=\Theta\left({\frac{k}{2^{t}}}
\right)I$ (this is because ${\left\lvert{\mathcal{C}_{a}}\right\rvert}\approx\frac{k}{2^{t+1}}$).
Moreover, the other terms in the product do not depend on the samples in $\mathcal{C}_{a}$, so with very high probability

$$ ${\left\lvert{\zeta_{0}^{\intercal}S^{-1}C_{1}w}\right\rvert}=\widetilde{O} \left(\sqrt{k}\times\left\|{S^{-1}C_{1}w}\right\|\right)\,.$ $$

Moreover, applying Claim 7.7, we have that with very high probability

$$ $\left\|{S-(n-k)I}\right\|=o(n)\Rightarrow\left\|{S^{-1}}\right\|=O\left(\frac{ 1}{n}\right)\,,$ $$

and applying Lemma 8.16, we have that with very high probability,

$$ $\left\|{C_{1}w}\right\|=\widetilde{O}\left(\left(k+\sqrt{kd}\right)\left\|{w} \right\|\right)\,.$ $$

Therefore, with very high probability

$$ ${\left\lvert{\zeta_{0}^{\intercal}S^{-1}C_{1}w}\right\rvert}=\widetilde{O} \left(\sqrt{k}\times\left\|{S^{-1}C_{1}w}\right\|\right)=\widetilde{O}\left( \sqrt{k}\frac{k+\sqrt{kd}}{n}\left\|{w}\right\|\right)\,,$ $$

and similarly for the other off-diagonal term, with very high probability

$$ ${\left\lvert{\zeta_{1}^{\intercal}S^{-1}C_{0}w}\right\rvert}=\widetilde{O} \left(\sqrt{k}\frac{k+\sqrt{kd}}{n}\left\|{w}\right\|\right)\,.$ $$

It now remains to bound the diagonal terms.
Consider $\zeta_{0}^{\intercal}S^{-1}C_{0}w$.
We can open up the next bit of the indices of the samples to obtain

$$ $\displaystyle\zeta_{0}^{\intercal}S^{-1}C_{0}w$ $\displaystyle=\left(\zeta_{00}+\zeta_{01}\right)^{\intercal}A\left(C_{00}+C_{0 1}\right)w=$ (29) $\displaystyle=\zeta_{00}^{\intercal}S^{-1}C_{00}w+\zeta_{01}^{\intercal}S^{-1} C_{01}w+\zeta_{00}^{\intercal}S^{-1}C_{01}w+\zeta_{01}^{\intercal}S^{-1}C_{00}w$ $$

We split the RHS of equation (29) again into diagonal and off-diagonal terms.
The off-diagonal terms can be bounded again in exactly the same manner and the off diagonal can again be split by specifying another bit of the sample indices.
Applying this logic recursively, we have

$$ $\displaystyle\textup{Diagonal}_{\epsilon}$ $\displaystyle=\zeta^{\intercal}S^{-1}Cw=\underbrace{\zeta_{0}^{\intercal}S^{-1 }C_{0}w}_{\textup{Diagonal}_{0}}+\underbrace{\zeta_{1}^{\intercal}S^{-1}C_{1}w }_{\textup{Diagonal}_{1}}+\underbrace{\zeta_{1}^{\intercal}S^{-1}C_{0}w}_{ \textup{Off-Diagonal}_{1,0}}+\underbrace{\zeta_{0}^{\intercal}S^{-1}C_{1}w}_{ \textup{Off-Diagonal}_{0,1}}$ (30) $\displaystyle=\textup{Diagonal}_{00}+\textup{Diagonal}_{01}+\textup{Diagonal}_ {10}+\textup{Diagonal}_{11}+\textup{Off-Diagonal}_{1,0}+\textup{Off-Diagonal}_ {0,1}+$ $\displaystyle\;\;\;\;\;\;\;\;+\textup{Off-Diagonal}_{00,01}+\textup{Off- Diagonal}_{01,00}+\textup{Off-Diagonal}_{10,11}+\textup{Off-Diagonal}_{11,10}=\cdots$ $\displaystyle\cdots=\sum_{i\in[k]}\textup{Diagonal}_{i}+\sum_{t=1}^{t=\left \lceil\log_{2}(k)\right\rceil}\sum_{a\in\left\{{0,1}\right\}^{t}}\textup{Off- Diagonal}_{a\mid 0,a\mid 1}+\textup{Off-Diagonal}_{a\mid 1,a\mid 0}$ $$

From the same analysis as the one above, we see that for any bitstring $a\in\left\{{0,1}\right\}^{t}$ and bit $b\in\left\{{0,1}\right\}$, it holds that with very high probability

$$ ${\left\lvert{\textup{Off-Diagonal}_{a\mid b,a\mid\overline{b}}}\right\rvert}= \widetilde{O}\left(\sqrt{{\left\lvert{\mathcal{C}_{a}}\right\rvert}}\frac{{ \left\lvert{\mathcal{C}_{a}}\right\rvert}+\sqrt{{\left\lvert{\mathcal{C}_{a}} \right\rvert}d}}{n}\left\|{w}\right\|\right)=\widetilde{O}\left(2^{-t}\sqrt{k} \frac{k+\sqrt{kd}}{n}\left\|{w}\right\|\right)$ $$

Union bounding over the $O\left(2^{t}\right)=\textnormal{poly}(n)$ off-diagonal combinations, we see that with very high probability they are all bounded.
Summing over these off-diagonal terms gives us

$$ $\displaystyle{\left\lvert{\sum_{t=1}^{t=\left\lceil\log_{2}(k)\right\rceil} \sum_{a\in\left\{{0,1}\right\}^{t}}\textup{Off-Diagonal}_{a\mid 0,a\mid 1}+ \textup{Off-Diagonal}_{a\mid 1,a\mid 0}}\right\rvert}\leq$ $\displaystyle\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\leq\sum_{t=1}^{t=\left \lceil\log_{2}(k)\right\rceil}2^{t}\widetilde{O}\left(2^{-t}\sqrt{k}\frac{k+ \sqrt{kd}}{n}\left\|{w}\right\|\right)=\widetilde{O}\left(\sqrt{k}\frac{k+ \sqrt{kd}}{n}\left\|{w}\right\|\right)\,.$ $$

Now, consider a $\text{Diagonal}_{i}$ term $X_{i}^{\intercal}S^{-1}X_{i}X_{i}^{\intercal}w$.
Because $w$ is independent of the $i$th sample $X_{i}$, with very high probability,

$$ ${\left\lvert{X_{i}^{\intercal}w}\right\rvert}=\widetilde{O}\left(\left\|{w} \right\|\right)\,,$ $$

and from Claim 7.7, we have that with very high probability,

$$ ${\left\lvert{X_{i}^{\intercal}S^{-1}X_{i}}\right\rvert}\leq\left\|{X_{i}} \right\|^{2}\left\|{S^{-1}}\right\|=\widetilde{O}\left(\frac{d}{n}\right)\,.$ $$

Therefore, from the triangle inequality, with very high probability

$$ ${\left\lvert{\sum_{i\in[k]}\textup{Diagonal}_{i}}\right\rvert}\leq\sum_{i\in[k ]}{\left\lvert{\textup{Diagonal}_{i}}\right\rvert}=\widetilde{O}\left(\frac{kd }{n}\left\|{w}\right\|\right)\,.$ $$

Altogether, we have

$$ ${\left\lvert{\xi^{\intercal}S^{-1}Cw}\right\rvert}=\widetilde{O}\left({\frac{ \left\|{w}\right\|}{\sqrt{k}}\times\left(1+\frac{d\sqrt{k}}{n}\right)}\right)\,,$ $$

concluding the proof of Claim 8.14.
∎

###### Lemma 8.16 .

###### Proof of Lemma 8.16 .

###### Proof of Claim 8.14 .

#### 8.3.4 Proof of Claim 8.15

In the previous portion of the proof, we bounded $\xi^{\intercal}S^{-1}Cw$, when $w$ is independent of $\xi$.
The rest of our analysis will be devoted to bounding the effect that replacing $S^{-1}$ with $\Sigma^{-1}$ will not make this inner product much larger.

As in the previous portion of the proof, we set $\zeta=k\xi$, and separate samples into clusters based on the least significant bits of the bitwise representations of their indices.
In particular, for any bitstring $a$, let $\mathcal{C}_{a}$, $\zeta_{a}$ and $C_{a}$ be as defined above, and define

$$ $S_{a}=\Sigma-C_{a}=\Sigma_{[n]\setminus\mathcal{C}_{a}}=\sum_{i\in[n]\setminus \mathcal{C}_{a}}X_{i}X_{i}^{\intercal}$ $$

to be the unnormalized empirical covariance of the samples *not* in the $\mathcal{C}_{a}$ cluster.

As in the proof of Claim 8.14, we will separate the contributions to $\zeta^{\intercal}\left(S^{-1}-\Sigma^{-1}\right)Cw$ based on the cluster of the samples and label these contributions as diagonal or off-diagonal based on whether or not the same samples were used in $\zeta$ and $C$.

$$ $\displaystyle\zeta^{\intercal}\left(S^{-1}-\Sigma^{-1}\right)Cw$ $\displaystyle=\left(\zeta_{0}+\zeta_{1}\right)^{\intercal}\left(S^{-1}-\Sigma^ {-1}\right)\left(C_{0}+C_{1}\right)w=$ (31) $\displaystyle=\underbrace{\zeta_{0}^{\intercal}\left(S^{-1}-\Sigma^{-1}\right) C_{0}w+\zeta_{1}^{\intercal}\left(S^{-1}-\Sigma^{-1}\right)C_{1}w}_{\textup{ diagonal terms}}+$ $\displaystyle\;\;\;\;\;\;\;\;\;\;\;\;+\underbrace{\zeta_{0}^{\intercal}\left(S ^{-1}-\Sigma^{-1}\right)C_{1}w+\zeta_{1}^{\intercal}\left(S^{-1}-\Sigma^{-1} \right)C_{0}w}_{\textup{off-diagonal terms}}$ $$

In other words, we have the recursive formula that for all $a\in\left\{{0,1}\right\}^{*}$,

$$ $\displaystyle\textup{Diagonal Term}_{a}$ $\displaystyle=\textup{Diagonal Term}_{a\mid 0}+\textup{Diagonal Term}_{a\mid 1}+$ (32) $\displaystyle+\textup{Off-Diagonal Term}_{a\mid 0,a\mid 1}+\textup{Off- Diagonal Term}_{a\mid 1,a\mid 0}$ $$

Applying equation (32) recursively, we have

$$ $\displaystyle\zeta^{\intercal}\left(S^{-1}-\Sigma^{-1}\right)Cw$ $\displaystyle=\textup{Diagonal Term}_{\epsilon}=\cdots=\sum_{i\in[k]}\textup{ Diagonal Term}_{i}+$ $\displaystyle\;\;\;\;\;\;\;\;\;\;\;+\sum_{t=0}^{\left\lceil\log(k)\right\rceil }\sum_{a\in\left\{{0,1}\right\}^{t}}\left(\textup{Off-Diagonal Term}_{a\mid 0, a\mid 1}+\textup{Off-Diagonal Term}_{a\mid 1,a\mid 0}\right)$ $$

In Claim 8.17, we will bound the off-diagonal terms.
We will show that with very high probability,

$$ $\displaystyle{\left\lvert{\textup{Off-Diagonal Term}_{a,b}}\right\rvert}$ $\displaystyle=\widetilde{O}\left(\frac{\sqrt{{\left\lvert{\mathcal{C}_{a}} \right\rvert}}}{n}\times\left({\left\lvert{\mathcal{C}_{b}}\right\rvert}+\sqrt {{\left\lvert{\mathcal{C}_{b}}\right\rvert}d}\right)\times\left\|{w}\right\| \times\left(1+\frac{d\sqrt{k}}{n}\right)\right)=$ $\displaystyle=\widetilde{O}\left(2^{-\left({\left\lvert{a}\right\rvert}+{\left \lvert{b}\right\rvert}\right)/2}\frac{\sqrt{k}}{n^{2}}\times\left(k+\sqrt{kd} \right)\times\left\|{w}\right\|\times\left(1+\frac{d\sqrt{k}}{n}\right)\right)\,.$ $$

Since in our case ${\left\lvert{a}\right\rvert}={\left\lvert{b}\right\rvert}=t$, we may conclude that with very high probability

$$ ${\left\lvert{\textup{Off-Diagonal Term}_{a,b}}\right\rvert}=\widetilde{O}\left (2^{-t}\frac{\sqrt{k}}{n^{2}}\times\left(k+\sqrt{kd}\right)\times\left\|{w} \right\|\times\left(1+\frac{d\sqrt{k}}{n}\right)\right)\,.$ $$

Applying the triangle inequality, with very high probability the total contribution of all the off-diagonal terms is of order

$$ $\displaystyle{\left\lvert{\sum_{t=0}^{\left\lceil\log(k)\right\rceil}\sum_{a \in\left\{{0,1}\right\}^{t}}\left(\textup{Off-Diagonal Term}_{a\mid 0,a\mid 1} +\textup{Off-Diagonal Term}_{a\mid 1,a\mid 0}\right)}\right\rvert}=$ $\displaystyle\;\;\;\;\;\;\;\;\;=\sum_{t=0}^{\left\lceil\log(k)\right\rceil} \sum_{a\in\left\{{0,1}\right\}^{t}}{\left\lvert{\textup{Off-Diagonal Term}_{a \mid 0,a\mid 1}}\right\rvert}+{\left\lvert{\textup{Off-Diagonal Term}_{a\mid 1 ,a\mid 0}}\right\rvert}=$ $\displaystyle\;\;\;\;\;\;\;\;\;=\sum_{t=0}^{\left\lceil\log(k)\right\rceil} \sum_{a\in\left\{{0,1}\right\}^{t}}\widetilde{O}\left(2^{-t}\frac{\sqrt{k}}{n^ {2}}\times\left(k+\sqrt{kd}\right)\times\left\|{w}\right\|\times\left(1+\frac{ d\sqrt{k}}{n}\right)\right)=$ $\displaystyle\;\;\;\;\;\;\;\;\;=\sum_{t=0}^{\left\lceil\log(k)\right\rceil} \widetilde{O}\left(\frac{\sqrt{k}}{n^{2}}\times\left(k+\sqrt{kd}\right)\times \left\|{w}\right\|\times\left(1+\frac{d\sqrt{k}}{n}\right)\right)=$ $\displaystyle\;\;\;\;\;\;\;\;\;=\widetilde{O}\left(\frac{\sqrt{k}}{n^{2}} \times\left(k+\sqrt{kd}\right)\times\left\|{w}\right\|\times\left(1+\frac{d \sqrt{k}}{n}\right)\right)$ $$

This leaves us with only the “single-sample” diagonal terms

$$ $\textup{Diagonal Term}_{i}=X_{i}^{\intercal}\left(S^{-1}-\Sigma^{-1}\right)X_{ i}X_{i}^{\intercal}w\,.$ $$

To analyse this term, simply note that from Claim 7.7 and the CS inequality, with very high probability,

$$ ${\left\lvert{X_{i}^{\intercal}\left(S^{-1}-\Sigma^{-1}\right)X_{i}}\right \rvert}\leq\left\|{X_{i}}\right\|^{2}\left(\left\|{S^{-1}}\right\|+\left\|{ \Sigma^{-1}}\right\|\right)=\widetilde{O}\left({\frac{d}{n}}\right)\,.$ $$

Moreover, from our assumption that $X_{i}\sim\mathcal{X}$ is well-behaved and that $w$ is independent of $X_{i}$, we have that with very high probability

$$ ${\left\lvert{X_{i}^{\intercal}w}\right\rvert}=\widetilde{O}\left({\left\|{w} \right\|}\right)\,.$ $$

Therefore, with very high probability

$$ $\frac{1}{k}\sum_{i\in[k]}{\left\lvert{\textup{Diagonal Term}_{i}}\right\rvert} =\widetilde{O}\left({\frac{\left\|{w}\right\|}{\sqrt{k}}\times\frac{d\sqrt{k}} {n}}\right)\,.$ $$

##### Bounding the Off-Diagonal Terms

###### Claim 8.17 .

Consider the off-diagonal term

$$ $\zeta_{a}^{\intercal}\left(S^{-1}-\Sigma^{-1}\right)C_{b}w\,,$ $$

where $a\neq b$ are bitstrings representing disjoint clusters $\mathcal{C}_{a}\cap\mathcal{C}_{b}=\emptyset$.

With very high probability,

$$ $\zeta_{a}^{\intercal}\left(S^{-1}-\Sigma^{-1}\right)C_{b}w=\widetilde{O}\left( \frac{\sqrt{{\left\lvert{\mathcal{C}_{a}}\right\rvert}}}{n}\times\left({\left \lvert{\mathcal{C}_{b}}\right\rvert}+\sqrt{{\left\lvert{\mathcal{C}_{b}}\right \rvert}d}\right)\times\left\|{w}\right\|\times\left(1+\frac{d\sqrt{k}}{n} \right)\right)\,.$ $$

###### Proof of Claim 8.17 .

Our goal is to bound the following term in absolute value

$$ $\zeta_{a}^{\intercal}\left(S^{-1}-\Sigma^{-1}\right)C_{b}w\,.$ $$

The difficulty in analysing this term is that both $\zeta_{a}$ and $\left(S^{-1}-\Sigma^{-1}\right)$ may depend on the same samples.
To circumvent this, we begin by splitting $\left(S^{-1}-\Sigma^{-1}\right)$ into a term that is easy to deal with and a small correction:

$$ $\left(S^{-1}-\Sigma^{-1}\right)=\left(S^{-1}-S_{a}^{-1}\right)+\left(S_{a}^{-1 }-\Sigma^{-1}\right)\,,$ $$

where $S_{a}=\Sigma-C_{a}=\sum_{i\notin\mathcal{C}_{a}}X_{i}X_{i}^{\intercal}$.

The first component has the property that it does not depend on the samples in $\mathcal{C}_{a}$, while $\zeta_{a}$ depends only on the samples in $a$.
Therefore, because $\zeta_{a}$ is well-behaved with covariance ${\left\lvert{\mathcal{C}_{a}}\right\rvert}I$, with very high probability

$$ ${\left\lvert{\zeta_{a}^{\intercal}\left(S^{-1}-S_{a}^{-1}\right)C_{b}w}\right \rvert}=\widetilde{O}\left(\sqrt{{\left\lvert{\mathcal{C}_{a}}\right\rvert}} \left\|{\left(S^{-1}-S_{a}^{-1}\right)C_{b}w}\right\|\right)$ $$

We begin by bounding the norm of $\left(S^{-1}-S_{a}^{-1}\right)C_{b}w$.
First, we note that Lemma 8.16, with very high probability

$$ $\left\|{C_{b}w}\right\|=\widetilde{O}\left({\left({\left\lvert{\mathcal{C}_{b} }\right\rvert}+\sqrt{{\left\lvert{\mathcal{C}_{b}}\right\rvert}d}\right)\left \|{w}\right\|}\right)\,.$ $$

To continue, we bound the norm of $\left(S^{-1}-S_{a}^{-1}\right)$.
Let $A=S_{a}$ and $B=S_{a}-S=\sum_{i\in[k]\setminus\mathcal{C}_{a}}X_{i}X_{i}^{\intercal}$.
We utilize the identity

$$ $\left(A-B\right)^{-1}-A^{-1}=\left(A-B\right)^{-1}BA^{-1}\,,$ $$

as well as Claim 7.7 which states that with very high probability $\left\|{\left(A-B\right)^{-1}}\right\|,\left\|{A^{-1}}\right\|=O\left(\frac{1}
{n}\right)$ and that $\left\|{B}\right\|=\widetilde{O}\left(k+d\right)$ to show that with very high probability

$$ $\left\|{S^{-1}-S_{a}^{-1}}\right\|=\widetilde{O}\left(\frac{k+d}{n^{2}}\right)$ $$

Therefore, with very high probability

$$ $\displaystyle{\left\lvert{\zeta_{a}^{\intercal}\left(S^{-1}-S_{a}^{-1}\right)C _{b}w}\right\rvert}$ $\displaystyle=\widetilde{O}\left(\sqrt{{\left\lvert{\mathcal{C}_{a}}\right \rvert}}\times\frac{k+d}{n^{2}}\times\left({\left\lvert{\mathcal{C}_{b}}\right \rvert}+\sqrt{{\left\lvert{\mathcal{C}_{b}}\right\rvert}d}\right)\times\left\| {w}\right\|\right)=$ $\displaystyle=\widetilde{O}\left(\frac{\sqrt{{\left\lvert{\mathcal{C}_{a}} \right\rvert}}}{n}\times\left({\left\lvert{\mathcal{C}_{b}}\right\rvert}+\sqrt {{\left\lvert{\mathcal{C}_{b}}\right\rvert}d}\right)\times\left\|{w}\right\| \right)\,.$ $$

It remains to bound the contribution of the $S_{a}^{-1}-\Sigma^{-1}$ term, which brings with it the added difficulty that it may depend on the samples in $\mathcal{C}_{a}$.
To bound the effect of this term, we split the contributions to $\zeta$ once more

$$ $\displaystyle\textup{Hard Term}_{a,b}$ $\displaystyle\stackrel{{\scriptstyle\textup{def}}}{{=}}\zeta_{a}^{\intercal} \left(S_{a}^{-1}-\Sigma^{-1}\right)C_{b}w=$ (33) $\displaystyle=\left(\zeta_{a0}+\zeta_{a1}\right)^{\intercal}\left(S_{a}^{-1}- \Sigma^{-1}\right)C_{b}w=$ $\displaystyle=\underbrace{\zeta_{a0}^{\intercal}\left(S_{a}^{-1}-S_{a0}^{-1} \right)C_{b}w+\zeta_{a1}^{\intercal}\left(S_{a}^{-1}-S_{a1}^{-1}\right)C_{b}w} _{\textup{Easy Terms}}+$ $\displaystyle\;\;\;\;\;\;\;\;\;\;\;\;+\underbrace{\zeta_{a0}^{\intercal}\left( S_{a0}^{-1}-\Sigma^{-1}\right)C_{b}w+\zeta_{a1}^{\intercal}\left(S_{a1}^{-1}- \Sigma^{-1}\right)C_{b}w}_{\textup{Hard Terms}}=$ $\displaystyle=\textup{Easy Term}_{a0,b}+\textup{Easy Term}_{a1,b}+\textup{Hard Term}_{a0,b}+\textup{Hard Term}_{a1,b}$ $$

We split the right hand side of equation (33) into “easy” terms which can be dealt with using the same logic above and hard terms which can again be split into smaller easy and hard terms.
Applying equation (33) recursively, we have that:

$$ $\displaystyle\textup{Hard Term}_{a,b}$ $\displaystyle=\textup{Easy Term}_{a,0,b}+\textup{Easy Term}_{a,1,b}+\textup{ Hard Term}_{a0,b}+\textup{Hard Term}_{a1,b}=$ (34) $\displaystyle\;\;\;\;\;\vdots$ $\displaystyle=\sum_{i\in\mathcal{C}_{a}}\textup{Hard Term}_{i,b}+\sum_{t=1}^{ \log\left({\left\lvert{\mathcal{C}_{a}}\right\rvert}\right)+1}\sum_{\begin{ subarray}{c}a^{\prime}\in\left\{{0,1}\right\}^{t-1}\\ z\in\left\{{0,1}\right\}\end{subarray}}\textup{Easy Term}_{a\mid a^{\prime},z,b}$ $$

To bound the easy terms, we note that as in the analysis above, the samples included in their $\zeta$ are disjoint from the samples included in their other factors.
Therefore, with very high probability

$$ ${\left\lvert{\textup{Easy Term}_{a,z,b}}\right\rvert}\stackrel{{\scriptstyle \textup{def}}}{{=}}{\left\lvert{\zeta_{az}^{\intercal}\left(S_{a}^{-1}-S_{az}^ {-1}\right)C_{b}w}\right\rvert}=\widetilde{O}\left(\sqrt{{\left\lvert{\mathcal {C}_{az}}\right\rvert}}\times\left\|{\left(S_{a}^{-1}-S_{az}^{-1}\right)C_{b}w }\right\|\right)$ $$

Recall that we showed above that with very high probability

$$ $\left\|{C_{b}w}\right\|=\widetilde{O}\left(\left({\left\lvert{\mathcal{C}_{b}} \right\rvert}+\sqrt{{\left\lvert{\mathcal{C}_{b}}\right\rvert}d}\right)\times \left\|{w}\right\|\right)\,.$ $$

From here, we set $A=S_{az}$ and $B=C_{a\overline{z}}=S_{az}-S_{a}$, and recall the identity

$$ $\left(A-B\right)^{-1}-A^{-1}=A^{-1}B\left(A-B\right)^{-1}$ $$

to obtain the equation

$$ $S_{a}^{-1}-S_{az}^{-1}=S_{az}^{-1}C_{a\overline{z}}S_{a}^{-1}\,.$ $$

From here, it might be tempting to simply bound the norm of this product, but that would result in too loose of a bound.
Instead, we perform the somewhat finer relaxation

$$ $\left\|{\left(S_{a}^{-1}-S_{az}^{-1}\right)C_{b}w}\right\|=\left\|{S_{az}^{-1} C_{a\overline{z}}S_{a}^{-1}C_{b}w}\right\|\leq\left\|{S_{az}^{-1}}\right\| \times\left\|{C_{a\overline{z}}S_{a}^{-1}C_{b}w}\right\|\,.$ $$

The advantage of this finer analysis is that we can now use Lemma 8.16 again, but this time on the matrix $C_{a\overline{z}}$.
Indeed, $S_{a}$ depends only on the samples outside $\mathcal{C}_{a}\supseteq\mathcal{C}_{a\overline{z}}$, $C_{b}$ depends only on the samples in $\mathcal{C}_{b}$ which is disjoint from $\mathcal{C}_{a}$, and $S$ and $w$ do not depend on the first $k$ samples.
Therefore, $C_{a\overline{z}}$ is independent of them all, so applying Lemma 8.16, we have that with very high probability

$$ $\displaystyle\left\|{C_{a\overline{z}}S_{a}^{-1}C_{b}w}\right\|$ $\displaystyle=\widetilde{O}\left({\left\lvert{\mathcal{C}_{a\overline{z}}} \right\rvert}+\sqrt{{\left\lvert{\mathcal{C}_{a\overline{z}}}\right\rvert}d} \right)\times\left\|{S_{a}^{-1}C_{b}w}\right\|\leq\widetilde{O}\left({\left \lvert{\mathcal{C}_{a\overline{z}}}\right\rvert}+\sqrt{{\left\lvert{\mathcal{C }_{a\overline{z}}}\right\rvert}d}\right)\times\left\|{S_{a}^{-1}}\right\|\left \|{C_{b}w}\right\|=$ $\displaystyle=\widetilde{O}\left(\frac{{\left\lvert{\mathcal{C}_{a\overline{z} }}\right\rvert}+\sqrt{{\left\lvert{\mathcal{C}_{a\overline{z}}}\right\rvert}d} }{n}\times\left({\left\lvert{\mathcal{C}_{b}}\right\rvert}+\sqrt{{\left\lvert{ \mathcal{C}_{b}}\right\rvert}d}\right)\times\left\|{w}\right\|\right)$ $$

Altogether, we have that with very high probability

$$ ${\left\lvert{\textup{Easy Term}_{a,z,b}}\right\rvert}=\widetilde{O}\left(\sqrt {{\left\lvert{\mathcal{C}_{a}}\right\rvert}}\times\frac{{\left\lvert{\mathcal{ C}_{a}}\right\rvert}+\sqrt{{\left\lvert{\mathcal{C}_{a}}\right\rvert}d}}{n^{2} }\times\left({\left\lvert{\mathcal{C}_{b}}\right\rvert}+\sqrt{{\left\lvert{ \mathcal{C}_{b}}\right\rvert}d}\right)\times\left\|{w}\right\|\right)\,.$ $$

Therefore, for $a^{\prime}$ of length $t-1$, we have ${\left\lvert{\mathcal{C}_{a\mid a^{\prime}}}\right\rvert}\approx 2^{-t}{\left
\lvert{\mathcal{C}_{a}}\right\rvert}$, so with very high probability

$$ $\displaystyle{\left\lvert{\textup{Easy Term}_{a\mid a^{\prime},z,b}}\right\rvert}$ $\displaystyle=\widetilde{O}\left(\sqrt{{\left\lvert{\mathcal{C}_{a\mid a^{ \prime}}}\right\rvert}}\times\frac{{\left\lvert{\mathcal{C}_{a\mid a^{\prime}} }\right\rvert}+\sqrt{{\left\lvert{\mathcal{C}_{a\mid a^{\prime}}}\right\rvert} d}}{n^{2}}\times\left({\left\lvert{\mathcal{C}_{b}}\right\rvert}+\sqrt{{\left \lvert{\mathcal{C}_{b}}\right\rvert}d}\right)\times\left\|{w}\right\|\right)=$ $\displaystyle=\widetilde{O}\left(2^{-t}\sqrt{{\left\lvert{\mathcal{C}_{a}} \right\rvert}}\times\frac{2^{-t/2}{\left\lvert{\mathcal{C}_{a}}\right\rvert}+ \sqrt{{\left\lvert{\mathcal{C}_{a}}\right\rvert}d}}{n^{2}}\times\left({\left \lvert{\mathcal{C}_{b}}\right\rvert}+\sqrt{{\left\lvert{\mathcal{C}_{b}}\right \rvert}d}\right)\times\left\|{w}\right\|\right)$ $$

Therefore, summing over all of the easy terms and applying the triangle inequality, we have that with very high probability

$$ $\displaystyle{\left\lvert{\sum_{t=1}^{\log\left({\left\lvert{\mathcal{C}_{a}} \right\rvert}\right)+1}\sum_{\begin{subarray}{c}a^{\prime}\in\left\{{0,1} \right\}^{t-1}\\ z\in\left\{{0,1}\right\}\end{subarray}}\textup{Easy Term}_{a\mid a^{\prime},z, b}}\right\rvert}$ $\displaystyle\leq\sum_{t=1}^{\log\left({\left\lvert{\mathcal{C}_{a}}\right \rvert}\right)+1}\sum_{\begin{subarray}{c}a^{\prime}\in\left\{{0,1}\right\}^{t -1}\\ z\in\left\{{0,1}\right\}\end{subarray}}{\left\lvert{\textup{Easy Term}_{a\mid a ^{\prime},z,b}}\right\rvert}=$ $\displaystyle=\widetilde{O}\left(\sqrt{{\left\lvert{\mathcal{C}_{a}}\right \rvert}}\frac{{\left\lvert{\mathcal{C}_{a}}\right\rvert}+\sqrt{d{\left\lvert{ \mathcal{C}_{a}}\right\rvert}}}{n^{2}}\times\left({\left\lvert{\mathcal{C}_{b} }\right\rvert}+\sqrt{{\left\lvert{\mathcal{C}_{b}}\right\rvert}d}\right)\times \left\|{w}\right\|\right)=$ $\displaystyle=\widetilde{O}\left(\frac{\sqrt{{\left\lvert{\mathcal{C}_{a}} \right\rvert}}}{n}\times\left({\left\lvert{\mathcal{C}_{b}}\right\rvert}+\sqrt {{\left\lvert{\mathcal{C}_{b}}\right\rvert}d}\right)\times\left\|{w}\right\|\right)$ $$

It remains to bound the single-sample hard terms.
That is, we want to bound

$$ ${\left\lvert{X_{i}^{\intercal}\left(\Sigma_{[n]\setminus\left\{{i}\right\}}^{- 1}-\Sigma^{-1}\right)C_{b}w}\right\rvert}\,.$ $$

Using the same matrix identity as in the previous terms, we have that

$$ $\Sigma_{[n]\setminus\left\{{i}\right\}}^{-1}-\Sigma^{-1}=\Sigma^{-1}X_{i}X_{i} ^{\intercal}\Sigma_{[n]\setminus\left\{{i}\right\}}^{-1}$ $$

Therefore, with very high probability, we have

$$ $\displaystyle{\left\lvert{X_{i}^{\intercal}\left(\Sigma_{[n]\setminus\left\{{i }\right\}}^{-1}-\Sigma^{-1}\right)C_{b}w}\right\rvert}$ $\displaystyle\leq\left\|{X_{i}}\right\|\left\|{\left(\Sigma_{[n]\setminus\left \{{i}\right\}}^{-1}-\Sigma^{-1}\right)C_{b}w}\right\|=$ (35) $\displaystyle=\left\|{X_{i}}\right\|\left\|{\Sigma^{-1}X_{i}X_{i}^{\intercal} \Sigma_{[n]\setminus\left\{{i}\right\}}^{-1}C_{b}w}\right\|\leq$ $\displaystyle\leq\left\|{X_{i}}\right\|\left\|{\Sigma^{-1}}\right\|\left\|{X_{ i}}\right\|{\left\lvert{X_{i}^{\intercal}\Sigma_{[n]\setminus\left\{{i}\right \}}^{-1}C_{b}w}\right\rvert}\leq$ $\displaystyle\leq\left\|{X_{i}}\right\|\left\|{\Sigma^{-1}}\right\|\left\|{X_{ i}}\right\|\times\widetilde{O}\left(\left\|{\Sigma_{[n]\setminus\left\{{i} \right\}}^{-1}C_{b}w}\right\|\right)=$ $\displaystyle=\widetilde{O}\left(\left\|{X_{i}}\right\|^{2}\left\|{\Sigma^{-1} }\right\|\left\|{\Sigma_{[n]\setminus\left\{{i}\right\}}^{-1}}\right\|\left\|{ C_{b}w}\right\|\right)=\widetilde{O}\left(\frac{d\left({\left\lvert{\mathcal{C }_{b}}\right\rvert}+\sqrt{{\left\lvert{\mathcal{C}_{b}}\right\rvert}d}\right)} {n^{2}}\left\|{w}\right\|\right)$ $$

Therefore, from the triangle inequality, with very high probability

$$ ${\left\lvert{\sum_{i\in\mathcal{C}_{a}}\textup{Hard Term}_{i,b}}\right\rvert} \leq\sum_{i\in\mathcal{C}_{a}}{\left\lvert{\textup{Hard Term}_{i,b}}\right \rvert}=\widetilde{O}\left({\left\lvert{\mathcal{C}_{a}}\right\rvert}\frac{d \left({\left\lvert{\mathcal{C}_{b}+\sqrt{\mathcal{C}_{b}d}}\right\rvert}\right )}{n^{2}}\left\|{w}\right\|\right)\,,$ $$

completing our proof of Claim 8.17.
∎

###### Claim 8.17 .

###### Proof of Claim 8.17 .

### 8.4 Proof of Theorem 8.1

Recall that the OHARE algorithm works by computing each of the MSN style bounds produced by the ACRE algorithm and then adds a correction term to each one, where these correction terms correspond to the indirect removal effects due to the change in the reaveraging step.

In the previous subsection, we proved Claim 8.2 which shows that with very high probability the ACRE components of the OHARE algorithm will produce good bounds on well-behaved regressions with categorical features.
In order to conclude Theorem 8.1, we would also need to bound the higher order corrections that OHARE takes into account.

Finally, note that for the case of an unweighted one-hot encoding, we have $u_{i,j}=1_{i\in B_{j}}$ (i.e., the columns corresponding to the dummy variables are indicators of their respective buckets).

#### 8.4.1 Indirect Contributions to the First Order Term

Analysing the indirect contributions to the first order term in OHARE will require significantly more care than our analysis of the higher order terms.
This is because the first order term is the dominant one to begin with and the indirect contributions to it are smaller than the main effect by only a $\sqrt{\log(n)}$ factor, forcing us to track $\textup{polylog}(n)$ factors much more carefully.

Recall from Section 4.2, that the OHARE algorithm bounds the first order effect of removals on the regression result from above/below by

$$ $\displaystyle\text{bound}^{\pm}_{j}(k_{j})=d_{j}(k_{j})+\frac{c^{\pm}_{j}(k_{j })}{n_{j}-k_{j}}$ $$

where $d_{j}$ represents AMIP gradients on the reaveraged samples $\widetilde{X}_{i}$:

$$ $\displaystyle d_{j}(k_{j})=\max_{\begin{subarray}{c}T_{j}\subseteq B_{j}\\ {\left\lvert{T_{j}}\right\rvert}=k_{j}\end{subarray}}\sum_{i\in T_{j}}R_{i}e^{ \intercal}\widetilde{\Sigma}^{-1}\widetilde{X}_{i}\,,$ $$

and $c^{\pm}_{j}(k_{j})$ represents the effect that reaveraging after the removals can have on the AMIP gradients of the retained samples:

$$ $\displaystyle c^{+}_{j}(k_{j})=\max\left\{{\left(\max_{\begin{subarray}{c}T_{j }\subseteq B_{j}\\ {\left\lvert{T_{j}}\right\rvert}=k_{j}\end{subarray}}\sum_{i\in T_{j}}R_{i} \right)\left(\max_{\begin{subarray}{c}T_{j}\subseteq B_{j}\\ {\left\lvert{T_{j}}\right\rvert}=k_{j}\end{subarray}}\sum_{i\in T_{j}}e^{ \intercal}\widetilde{\Sigma}^{-1}\widetilde{X}_{i}\right),\left(\min_{\begin{ subarray}{c}T_{j}\subseteq B_{j}\\ {\left\lvert{T_{j}}\right\rvert}=k_{j}\end{subarray}}\sum_{i\in T_{j}}R_{i} \right)\left(\min_{\begin{subarray}{c}T_{j}\subseteq B_{j}\\ {\left\lvert{T_{j}}\right\rvert}=k_{j}\end{subarray}}\sum_{i\in T_{j}}e^{ \intercal}\widetilde{\Sigma}^{-1}\widetilde{X}_{i}\right)}\right\}$ $\displaystyle c^{-}_{j}(k_{j})=\max\left\{{\left(\min_{\begin{subarray}{c}T_{j }\subseteq B_{j}\\ {\left\lvert{T_{j}}\right\rvert}=k_{j}\end{subarray}}\sum_{i\in T_{j}}R_{i} \right)\left(\max_{\begin{subarray}{c}T_{j}\subseteq B_{j}\\ {\left\lvert{T_{j}}\right\rvert}=k_{j}\end{subarray}}\sum_{i\in T_{j}}e^{ \intercal}\widetilde{\Sigma}^{-1}\widetilde{X}_{i}\right),\left(\max_{\begin{ subarray}{c}T_{j}\subseteq B_{j}\\ {\left\lvert{T_{j}}\right\rvert}=k_{j}\end{subarray}}\sum_{i\in T_{j}}R_{i} \right)\left(\min_{\begin{subarray}{c}T_{j}\subseteq B_{j}\\ {\left\lvert{T_{j}}\right\rvert}=k_{j}\end{subarray}}\sum_{i\in T_{j}}e^{ \intercal}\widetilde{\Sigma}^{-1}\widetilde{X}_{i}\right)}\right\}$ $$

For $k_{j}=n_{j}$, there is no reaveraging effect and $\text{bound}^{\pm}_{j}(n_{j})=d_{j}(n_{j})$.

Our goal will be to show that $d_{j}(k_{j})$, which is the contribution of the AMIP gradients to this first order effect, is the dominant effect.
In particular, the main claim we will prove in this subsection is

###### Claim 8.18 .

Let $k_{\textup{threshold}}=\widetilde{\Omega}\left(n^{1-\varepsilon}\right)$ be as promised by Claim 8.11.

Then, with very high probability

$$ $\forall k\leq k_{\textup{threshold}},\;\;\;\;\max_{k_{1}+\cdots+k_{m}=k}\left \{{\sum_{j\in[m]}\textup{bound}_{j}^{\pm}(k_{j})}\right\}=\left(1\pm O\left( \frac{1}{\sqrt{\log(n)}}\right)\right)\textup{AMIP}(k)\,,$ $$

where

$$ $\textup{AMIP}(k)=\max_{T\in\binom{[n]}{k}}\left\{{\sum_{i\in T}\alpha_{i}}\right\}$ $$

is the sum over the $k$ largest AMIP influence scores

$$ $\alpha_{i}=e^{\intercal}\widetilde{\Sigma}^{-1}\widetilde{X}_{i}R_{i}\,.$ $$

###### Proof of Claim 8.18 .

We expect $d_{j}(k_{j})$ to grow roughly linearly with $k_{j}$, which motivates us to focus on the expression

$$ $\textup{Scaled Indirect Effect}=\eta_{j}(k_{j})\stackrel{{\scriptstyle\textup{ def}}}{{=}}\frac{\max\left\{{{\left\lvert{c^{+}_{j}(k_{j})}\right\rvert},{ \left\lvert{c^{-}_{j}(k_{j})}\right\rvert}}\right\}}{k_{j}(n_{j}-k_{j})}\,.$ $$

We will bound these $\eta_{j}$ in the following lemma:

###### Claim 8.19 .

There exists some $\nu=\textup{polylog}(n)$ such that for all $k_{j}\leq\frac{n_{j}}{\nu}$, with very high probability

$$ $\eta_{j}(k_{j})=O\left(\frac{1}{n}\right)\,.$ $$

Moreover, for any $\nu^{\prime}=\textup{polylog}(n)$ there exists some threshold $\tau=\textup{polyloglog}(n)$, such that with very high probability,

$$ $\displaystyle{\left\lvert{\left\{i\in B_{j}\mid{\left\lvert{R_{i}}\right\rvert }>\tau\right\}}\right\rvert}<\frac{n_{j}}{\nu^{\prime}}$ $\displaystyle{\left\lvert{\left\{i\in B_{j}\mid{\left\lvert{e^{\intercal} \widetilde{\Sigma}^{-1}\widetilde{X}_{i}}\right\rvert}>\frac{\tau}{n}\right\}} \right\rvert}<\frac{n_{j}}{\nu^{\prime}}$ $$

In particular, for all $k_{j}\geq\frac{n_{j}}{\nu}$, with very high probability

$$ $\displaystyle\max_{T_{j}\in\binom{B_{j}}{k_{j}}}\left\{{\sum_{i\in T_{j}} \alpha_{i}}\right\}=O\left(\frac{\textup{polyloglog}(n)}{n}\times k_{j}\right)$ $\displaystyle\eta_{j}(k_{j})=O\left(\frac{\textup{polyloglog}(n)}{n}\right)\,,$ $$

where $\alpha_{i}=e^{\intercal}\widetilde{\Sigma}^{-1}\widetilde{X}_{i}R_{i}$ are the AMIP influence scores.

The first part of Claim 8.19 promises that the $\eta_{j}$ components to the bound are very small in any bucket for which there are not too many removals $k_{j}<\frac{n_{j}}{\textup{polylog}(n)}$, while the second portion of the claim will help us bound the total contribution of buckets from which more than $\frac{n_{j}}{\textup{polylog}(n)}$ have been removed.

Let $T=\mathop{\textup{argmax}}_{T\in\binom{[n]}{k}}\left\{{\sum_{i\in T}\alpha_{i}
}\right\}$ denote the set of $k$ samples with largest AMIP influences, and let $\kappa_{j}={\left\lvert{B_{j}\cap T}\right\rvert}$ denote the distribution of these samples across the buckets.
By definition, we have

$$ $\sum_{i\in T}\alpha_{i}=\textup{AMIP}(k)\,.$ $$

From Claim 8.11, we know that with very high probability all the samples in $T$ must have influence at least $\min_{i\in T}\left\{{\alpha_{i}}\right\}=\Omega\left(\frac{\sqrt{\log(n)}}{n}\right)$.
But from the second part of Claim 8.19, we know that with very high probability, for all $j$, the $j$th bucket cannot have more than $\frac{n_{j}}{\nu}$ such samples, so with very high probability $\kappa_{j}<\frac{n_{j}}{\nu}$.

Consider our maximization problem

$$ $\textup{MaxScore}=\max_{k_{1}+\cdots+k_{m}=k}\left\{{\sum_{j\in[m]}\textup{ bound}_{j}^{\pm}(k_{j})}\right\}\,.$ $$

This maximization is lower bounded by every valid assignment to $k_{1},\ldots,k_{m}$, so in particular it is lower bounded by the score of $\kappa_{1},\ldots,\kappa_{m}$.
Utilising the first part of Claim 8.19 to bound $\eta_{j}(\kappa_{j})$, we have

$$ $\max_{k_{1}+\cdots+k_{m}=k}\left\{{\sum_{j\in[m]}\textup{bound}_{j}^{\pm}(k_{j })}\right\}\geq\sum_{j\in[m]}\textup{bound}_{j}^{\pm}(\kappa_{j})=\sum_{i\in T }\alpha_{i}-\sum_{j\in[m]}\kappa_{j}\eta_{j}(\kappa_{j})\in\textup{AMIP}(k)-O \left(\frac{k}{n}\right)\,.$ $$

Now, consider any other assignment $k_{1},\ldots,k_{m}$.
If we still have $k_{j}\leq\frac{n_{j}}{\nu}$ for all $j$, then $\sum_{j\in[m]}\textup{bound}_{j}^{\pm}(k_{j})$ will still be bounded by $\textup{AMIP}(k)\pm O\left(\frac{k}{n}\right)$ following the same logic as above.

Otherwise, let ${j_{1}},\ldots,{j_{\ell}}$ denote the set of buckets for whick $k_{j_{i}}>\frac{n_{j_{i}}}{\nu}$, and define $k^{\prime}\stackrel{{\scriptstyle\textup{def}}}{{=}}k-k_{j_{1}}-\cdots-k_{j_{
\ell}}$.
From the same analysis as above, we have

$$ $\sum_{j\in[m]}\textup{bound}_{j}^{\pm}(k_{j})\leq\textup{AMIP}(k^{\prime})+ \sum_{i\in[\ell]}\left\{{\textup{bound}_{j_{i}}^{\pm}(k_{j_{i}})}\right\}+O \left(\frac{k^{\prime}}{n}\right)\,.$ $$

Next, note that we know that with very high probability

$$ $\textup{AMIP}(k)\geq\textup{AMIP}(k^{\prime})+\Omega\left(\frac{\sqrt{\log(n)} \left(k-k^{\prime}\right)}{n}\right)\,,$ $$

because Claim 8.11 tells us that with very high probability, there are at least $k_{\textup{threshold}}\geq k$ samples, each of which has a sufficiently large contribution to the AMIP score.

Therefore, using the last part of Claim 8.19, which bounds with very high probability every term in the bounds of buckets with more than $\frac{n_{j}}{\nu}$ removals, we have

$$ $\displaystyle\sum_{j\in[m]}\textup{bound}_{j}^{\pm}(k_{j})$ $\displaystyle\leq\textup{AMIP}(k^{\prime})+\sum_{i\in[\ell]}\left\{{\textup{ bound}_{j_{i}}^{\pm}(k_{j_{i}})}\right\}+O\left(\frac{k^{\prime}}{n}\right)\leq$ $\displaystyle\leq\textup{AMIP}(k)+\underbrace{O\left(\frac{\textup{polyloglog} (n)\times\left(k-k^{\prime}\right)}{n}\right)-\Omega\left(\frac{\sqrt{\log(n)} \left(k-k^{\prime}\right)}{n}\right)}_{\leq 0}\pm O\left(\frac{k^{\prime}}{n} \right)\leq$ $\displaystyle\leq\textup{AMIP}(k)+O\left(\frac{k}{n}\right)\,.$ $$

Altogether, we have bounded our maximization target from above and from below by $\textup{AMIP}(k)+O\left(\frac{k}{n}\right)=\left(1\pm O\left(\frac{1}{\sqrt{
\log(n)}}\right)\right)\textup{AMIP}(k)$, completing our proof.

∎

###### Proof of Claim 8.19 .

If $k_{j}<\frac{n_{j}}{\nu}\ll n_{j}$, then we use the bound

$$ ${\left\lvert{c^{\pm}_{j}(k_{j})}\right\rvert}<\frac{k_{j}^{2}}{n_{j}-k_{j}} \max_{i\in[n]}\left\{{{\left\lvert{R_{i}}\right\rvert}}\right\}\max_{i\in[n]} \left\{{{\left\lvert{e^{\intercal}\widetilde{\Sigma}^{-1}\widetilde{X}_{i}} \right\rvert}}\right\}=\widetilde{O}\left({\frac{k_{j}}{\nu n}}\right)=o\left( \frac{k_{j}}{n}\right)\,.$ $$

Next, note that due to Claim 8.10 (which states that with very high probability the empirical residuals $R_{i}$ is close to the ground truth residual $\left({R}_{\textup{gt}}\right)_{i}$), we have that with very high probability for all $i$:

$$ ${\left\lvert{R_{i}-\left({R}_{\textup{gt}}\right)_{i}}\right\rvert}=o(1)<1\,.$ $$

From our assumption that the ground truth residuals were normally distributed, we have that

$$ ${\Pr_{\left({R}_{\textup{gt}}\right)_{i}\sim\mathcal{N}(0,1)}\left[{{\left \lvert{\left({R}_{\textup{gt}}\right)_{i}}\right\rvert}>\tau-1}\right]}<\frac{ 1}{2\nu^{\prime}}\,.$ $$

Combining the two, along with the Hoeffding bound which promises us that with very high probability

$$ ${\left\lvert{\left\{i\in B_{j}\mid{\left\lvert{{R}_{i}}\right\rvert}>\tau-1 \right\}}\right\rvert}\leq n_{j}{\Pr_{\left({R}_{\textup{gt}}\right)_{i}\sim \mathcal{N}(0,1)}\left[{{\left\lvert{\left({R}_{\textup{gt}}\right)_{i}}\right \rvert}>\tau-1}\right]}\pm\widetilde{O}\left({\sqrt{n_{j}}}\right)\,,$ $$

we have that with very high probability

$$ ${\left\lvert{\left\{i\in B_{j}\mid{\left\lvert{{R}_{i}}\right\rvert}>\tau \right\}}\right\rvert}<{\left\lvert{\left\{i\in B_{j}\mid{\left\lvert{ \widetilde{R}_{i}}\right\rvert}>\tau-1\right\}}\right\rvert}=\frac{n_{j}}{\nu^ {\prime}}\,.$ $$

We bound ${\left\lvert{\left\{i\in B_{j}\mid{\left\lvert{e^{\intercal}\widetilde{\Sigma}
^{-1}\widetilde{X}_{i}}\right\rvert}>\frac{\tau}{n}\right\}}\right\rvert}$ in much the same manner, by utilizing Claim 8.8 which states that with very high probability ${\left\lvert{e^{\intercal}\widetilde{\Sigma}^{-1}\widetilde{X}_{i}-\frac{e^{
\intercal}X_{i}}{n}}\right\rvert}=o\left(\frac{1}{n}\right)$, and our assumption that the $X_{i}$ are very well behaved, which shows that there can’t be too many samples in a given bucket for which ${\left\lvert{e^{\intercal}X_{i}}\right\rvert}>\textup{polyloglog}(n)$.

∎

###### Claim 8.18 .

###### Proof of Claim 8.18 .

###### Claim 8.19 .

###### Proof of Claim 8.19 .

#### 8.4.2 Indirect Contributions to the Higher Order Terms

##### Indirect Contributions to the Covariance Shift Term

Recall from Section 4.3 that we bound the covariance shift as

$$ $\max\left\{{\lambda\left(\hat{\Sigma}_{S}^{-1}\right)}\right\}\leq\frac{1}{1- \max\left\{{\lambda\left(\widetilde{\Sigma}_{T}\right)}\right\}-\max_{k_{1}+ \cdots+k_{m}=k}\left\{{\sum_{j\in[m]}\frac{1}{n_{j}-k_{j}}\left\|{\sum_{i\in T \cap B_{j}}\widetilde{\Sigma}^{-1/2}\widetilde{X}_{i}}\right\|}\right\}}\,.$ $$

Each component in the denominator is bounded separately by running an MSN-bounding algorithm.
The first MSN is run on the Gram matrix $G_{X\otimes X}$ whose $i,j$ entry is:

$$ $\left(\widetilde{X}_{i}^{\intercal}\widetilde{\Sigma}^{-1}\widetilde{X}_{j} \right)^{2}\,.$ $$

From Claim 8.2, the regression $\widetilde{X},\widetilde{Y}$ is well-behaved, allowing us to use the same analysis as in Claim 7.2 to bound the output of this MSN by

$$ $M_{X\otimes X}(k)=\widetilde{O}\left({\sqrt{k\frac{d^{2}}{n^{2}}+k^{2}\frac{d} {n^{2}}}}\right)\,.$ $$

Similarly, we define $M_{j}$ to be the MSN bound achieved by RTI on the Gram matrix

$$ $G_{j}[i_{1},i_{2}]=\widetilde{X}_{i_{1}}^{\intercal}\widetilde{\Sigma}^{-1} \widetilde{X}_{i_{2}}\,.$ $$

Claims 8.5 and 8.6 promise us that with very high probability the largest diagonal entry of this Gram matrix is at most $\widetilde{O}\left({\frac{d}{n}}\right)$ and its largest off-diagonal entry is at most $\widetilde{O}\left({\frac{\sqrt{d}}{n}}\right)$.
Therefore, the resulting MSN bounds are at most

$$ $M_{j}(k_{j})=\widetilde{O}\left({\sqrt{\frac{k_{j}d+k_{j}^{2}\sqrt{d}}{n}}} \right)\,.$ $$

Recall that for our actual OHARE bound, we also utilize the symmetry that allows us to replace $M_{j}(k_{j})$ with $\overline{M}_{j}(k_{j})=\min\left\{{M_{j}(k_{j}),M_{j}(n_{j}-k_{j})}\right\}$.

For all $k_{j}\leq\frac{n_{j}}{2}$, we have

$$ $\frac{\overline{M}_{j}(k_{j})^{2}}{(n_{j}-k_{j})k_{j}}\leq\frac{M_{j}(k_{j})^{ 2}}{(n_{j}-k_{j})k_{j}}=\widetilde{O}\left({\frac{d+k_{j}\sqrt{d}}{n(n_{j}-k_{ j})}}\right)=\widetilde{O}\left({\frac{\sqrt{d}}{n}}\right)\,,$ $$

and for all $k_{j}>\frac{n_{j}}{2}$, we have

$$ $\frac{\overline{M}_{j}(k_{j})^{2}}{(n_{j}-k_{j})k_{j}}\leq\frac{M_{j}(n_{j}-k_ {j})^{2}}{(n_{j}-k_{j})k_{j}}=\widetilde{O}\left({\frac{\sqrt{d}}{n}}\right)\,.$ $$

Therefore, for all $k<k_{\textup{threshold}}$

$$ $\max_{k_{1}+\cdots+k_{m}=k}\left\{{\sum_{j\in[m]}\frac{\overline{M}_{j}(k_{j}) ^{2}}{n_{j}-k_{j}}}\right\}=\widetilde{O}\left({\frac{\sqrt{d}k}{n}}\right)=o \left(M_{X\otimes X}(k)^{2}\right)=o\left(1\right)\,.$ $$

Therefore, in this regime the covariance shift term does not contribute a factor of more than $3$.

##### Indirect Contributions to the X ⁢ R 𝑋 𝑅 XR italic_X italic_R Term

Recall from Section 4.3 that the residual contributions to the indirect XR term are bounded by

$$ $\rho_{j}(k_{j})\stackrel{{\scriptstyle\textup{def}}}{{=}}\min\left\{{\max_{T_{ j}\in\binom{B_{j}}{k_{j}}}\left\{{\sum_{i\in T_{j}}{\left\lvert{R_{i}}\right \rvert}}\right\},\max_{T_{j}\in\binom{B_{j}}{k_{j}}}\left\{{\sum_{i\in B_{j} \setminus T_{j}}{\left\lvert{R_{i}}\right\rvert}}\right\}}\right\}\,,$ $$

where $R_{i}$ are the empirical residuals.

With an analysis similar to the one above, we can utilize Claim 8.10 which states that with very high probability $\max_{i\in[n]}\left\{{{\left\lvert{R_{i}}\right\rvert}}\right\}=\widetilde{O}
\left({1}\right)$ to show that with very high probability

$$ $\frac{\rho_{j}(k_{j})}{\sqrt{k_{j}(n_{j}-k_{j})}}=\widetilde{O}\left({1}\right )\,.$ $$

Combining this with our bound on $\frac{\overline{M}_{j}(k_{j})}{\sqrt{k_{j}(n_{j}-k_{j})}}$, we have

$$ $\max_{k_{1}+\cdots+k_{m}=k}\left\{{\sum_{j\in[m]}\frac{\overline{M}_{j}(k_{j}) \rho_{j}(k_{j})}{n_{j}-k_{j}}}\right\}=\widetilde{O}\left({\sqrt{\frac{\sqrt{d }}{n}}k}\right)\,.$ $$

Note that this is also smaller than the bound we proved for the direct contribution to the XR term in Section 7.5,

$$ $M_{XR}=\widetilde{O}\left({\sqrt{\frac{kd}{n}}+\sqrt{\frac{\sqrt{d}}{n}}k} \right)\,,$ $$

and which we can apply after reaveraging due to Claim 8.2 which states that $\widetilde{X},\widetilde{Y}$ is an ACRE-friendly with very high probability.

##### Indirect Contributions to the X ⁢ Z 𝑋 𝑍 XZ italic_X italic_Z Term

Finally, recall from Section 4.3 that the $Z$ component of the indirect effect on the $XZ$ term was bounded by

$$ $\zeta_{j}(k_{j})=\min\left\{{\max_{T_{j}\in\binom{B_{j}}{k_{j}}}\left\{{\sum_{ i\in T_{j}}{\left\lvert{Z_{i}}\right\rvert}}\right\},\max_{T_{j}\in\binom{B_{j }}{k_{j}}}\left\{{\sum_{i\in B_{j}\setminus T_{j}}{\left\lvert{Z_{i}}\right \rvert}}\right\}}\right\}\,,$ $$

where $Z_{i}=e^{\intercal}\widetilde{\Sigma}^{-1}\widetilde{X}_{i}$

Using the same analysis as above and Claim 8.8 which states that with very high probability ${\left\lvert{Z_{i}}\right\rvert}=\widetilde{O}\left({\frac{1}{n}}\right)$ for all $i\in[n]$, we have that with very high probability

$$ $\frac{\zeta_{j}(k_{j})}{\sqrt{k_{j}(n_{j}-k_{j})}}=\widetilde{O}\left({\frac{1 }{n}}\right)\,.$ $$

Combining this with our bound on $M_{j}(k_{j})$, we have

$$ $\max_{k_{1}+\cdots+k_{m}=k}\left\{{\sum_{j\in[m]}\frac{\overline{M}_{j}(k_{j}) \zeta_{j}(k_{j})}{n_{j}-k_{j}}}\right\}=\widetilde{O}\left({\sqrt{\frac{\sqrt{ d}}{n^{2}}}k}\right)\,.$ $$

As before, this is smaller than our bound on the direct effect

$$ $M_{XZ}=\widetilde{O}\left({\sqrt{\frac{kd}{n^{2}}+\frac{k^{2}\sqrt{d}}{n^{2}}} }\right)\,.$ $$

##### Putting it all together

In Claim 8.18 and over the last few paragraphs, we have bounded all the individual terms that go into generating the OHARE bounds.
In particular, we have shown that

$$ $\displaystyle U_{k},L_{k}$ $\displaystyle=\underbrace{\text{First Order}}_{\left(1\pm\left(\frac{1}{\sqrt{ \log(n)}}\right)\right)\times\textup{AMIP}(k)}\pm\frac{\lvert\text{Direct XR}+ \text{Indirect XR}\rvert\times\lvert\text{Direct XZ}+\text{Indirect XZ}\rvert} {1-\text{Direct CS}-\text{Indirect CS}}=$ $\displaystyle=\left(1\pm\left(\frac{1}{\sqrt{\log(n)}}\right)\right)\times \textup{AMIP}(k)\pm\widetilde{O}\left({\frac{kd}{n^{3/2}}+\frac{k^{2}\sqrt{d}} {n^{3/2}}}\right)=$ $\displaystyle=\left(1\pm\left(\frac{1}{\sqrt{\log(n)}}\right)\right)\times \textup{AMIP}(k)\,,$ $$

for all $k\leq k_{\textup{threshold}}$, concluding our proof of Theorem 8.1.

## Acknowledgments

Sam Hopkins was partially supported by NSF CAREER award no. 2238080.
Both authors were partially supported by MLA@CSAIL.

## Appendix A Proof of Claim 4.1

The main analytic idea we use in the OHARE algorithm is an alternative description of what happens when we perform a linear regression while controlling for a categorical feature.
We formalise this process in Claim 4.1, which we prove here:

###### Proof of Claim 4.1 .

Let $X=\left(F\mid U\right)$ be the covariates for the original regression.
Consider the Gram-Schmidt orthogonalization process, and define

$$ $\hat{X}=\left(F-D^{-1}U^{\intercal}F\mid U\right)$ $$

where $D=U^{\intercal}U$ is the diagonal matrix whose $j$th entry is $\left\|{u_{j}}\right\|$.

$\hat{X}$ now has a block-diagonal covariance matrix

$$ $\hat{\Sigma}=\hat{X}^{\intercal}\hat{X}=\begin{pmatrix}(F-D^{-1}U^{\intercal}F )^{\intercal}(F-D^{-1}U^{\intercal}F)&0\\ 0&D\end{pmatrix}$ $$

Similarly, let $\hat{Y}=\hat{Y}=Y-D^{-1}U^{\intercal}Y$.
These labels are now also perpendicular to the dummy variables $U$.
Therefore, for $\hat{X}=F-D^{-1}U^{\intercal}F$, $\hat{\Sigma}=\hat{X}^{\intercal}\hat{X}$, we have

$$ $\displaystyle\hat{\beta}$ $\displaystyle=\hat{\Sigma}^{-1}\hat{X}^{\intercal}\hat{Y}$ (36) $\displaystyle=\begin{pmatrix}(F-D^{-1}U^{\intercal}F)^{\intercal}(F-D^{-1}U^{ \intercal}F)&0\\ 0&D\end{pmatrix}^{-1}\begin{pmatrix}F-D^{-1}U^{\intercal}F\\ U\end{pmatrix}^{\intercal}\left(Y-D^{-1}U^{\intercal}Y\right)$ $\displaystyle=\begin{pmatrix}\widetilde{\Sigma}^{-1}\widetilde{X}^{\intercal} \widetilde{Y}\\ D^{-1}U^{\intercal}\left(Y-D^{-1}U^{\intercal}Y\right)\end{pmatrix}=\begin{ pmatrix}\hat{\beta}\\ \star\end{pmatrix}$ $$

∎

###### Proof of Claim 4.1 .

## Appendix B One-Hot Encodings are Almost Brittle

In this section, we will prove our claim from the introduction that datasets with one-hot encodings are arbitrarily close to being extremely brittle.
In particular, we will show that

###### Claim B.1 .

Let $X\in\mathbb{R}^{n\times d}$ be an array of features and let $Y\in\mathbb{R}^{n}$ be some labels, such that one of the features is 0 0 on all but $k_{\textup{bucket}}<n-d$ of the samples. In other words, for some set $S\subseteq[n]$ of size ${\left\lvert{S}\right\rvert}=n-k_{\textup{bucket}}$, we have

$$ $\forall i\in S\;\;\;\;\;X_{i,d}=0$ $$

Then, for all $\gamma\in\mathbb{R}^{d-1}$ and for all $c>0$, there exists a linear regression $X^{\prime},Y^{\prime}$ such that $\left\|{X^{\prime}-X}\right\|+\left\|{Y^{\prime}-Y}\right\|<c$, and $\text{OLS}\left(X^{\prime}_{S},Y^{\prime}_{S}\right)_{[d-1]}=\gamma$.

The reason regressions become so close to brittle is that with one of the features being always 0 0, we can make a very small change to its value, in a way that creates very strong correlations.

###### Proof of Claim B.1 .

The proof of Claim B.1 is relatively simple.

First, we want to ensure that the original regression problem has no other degeneracies.
We do this to ensure that the resulting OLS has a unique solution.

Let $V\subseteq\mathbb{R}^{n-k_{\textup{bucket}}}$ be the linear space spanned by the columns of $X_{S,[d-1]}$ and $Y_{S}$.
$V$ is spanned by $d$ vectors.
Therefore, $\dim V\leq d$.

If $\dim V=d$, we can skip this step, and if this inequality is strict we say that $X_{S,[d-1]},Y_{S}$ are degenerate.
If $\dim V<d$, it is easy to see that almost all $X_{S,[d-1]}^{\prime},Y_{S}^{\prime}$ in $B_{c/2}\left(X_{S,[d-1]},Y_{S}\right)$ (i.e. in the ball of radius $c/2$ around the original regression) are non-degenerate.
Therefore, let $X_{S,[d-1]}^{\prime},Y_{S}^{\prime}$ be such a non-degenerate pair.

Now, consider the vector $R=Y_{S}^{\prime}-X_{S,[d-1]}^{\prime}\gamma$.
This is the residual vector for the linear model $Y_{S}^{\prime}\approx X_{S,[d-1]}^{\prime}\gamma$, and by our assumption that $X_{S,[d-1]}^{\prime},Y_{S}^{\prime}$ are non-degenerate, $R\neq 0$ and is not in the span of the columns of $X_{S,[d-1]}^{\prime}$.

It is only left to decide the values of $X^{\prime}_{S,d}$.
Setting $X^{\prime}_{S,d}=\frac{c}{2\left\|{R}\right\|}R$, our regression has a perfect fit on the samples in $S$

$$ $Y^{\prime}_{S}=X^{\prime}_{S}\begin{pmatrix}\gamma\\ -\frac{2\left\|{R}\right\|}{c}\end{pmatrix}$ $$

By our construction, $\Sigma=\left(X_{S}^{\prime}\right)^{\intercal}X_{S}^{\prime}$ is full rank, making this OLS solution unique.
∎

###### Claim B.1 .

###### Proof of Claim B.1 .

## Peer Reviews

confidence: 4
contribution: 3
initial_rating: 8
presentation: 4
questions: Do you have a comment on how to expand your theoretical tightness results to broader types of assumptions on $X$ and $Y$. In particular, can you show a more general (but perhaps weaker) bound on $\frac{U_k}{L_k}$.
soundness: 4
strengths: This paper is well written and proposes a nice solution to a simple fundamental problem. Due to the challenging nature of this problem, their results still hinge on a reduction to an NP-complete problem. However, they are able to give performance bounds on their algorithm's tightness (both theoretically and experimentally), and crucially their algorithm improves over the naive approach (trying all subsets of size $k$) enough to provide meaningful audits of how robust the results from real-life linear regressions are.
summary: This paper studies data-based robustness for linear regression. Given an OLS solution that predicts $Y$ from $X$ (based on dataset $(X_1, Y_1), \dots, (X_n, Y_n)$), target direction $e$, and $1 \leq k \leq n$, they estimate $\Delta_k(e)$, the maximum possible change along the direction $e$ in $\beta$, the solution to the OLS, that can be achieved by removing $k$ points from the dataset. In essence, $\Delta_k(e)$ measures the robustness of a linear regression to potential training-time adversarial attacks.

Their main contribution is an algorithm, ACRE, which computes upper and lower bounds $U_k$ and $L_k$ for $\Delta_k(e)$. ACRE operates by bounding the difference between $\beta_S$, the value of $\beta$ after restricting to a subset, $S$ of the data, and $\beta$ by expressing it with two terms: first, a combination of the residuals of the original regression, and second, a higher order term, that factors in matrix terms based on the dataset $X$.

Next, ACRE bounds each of these terms by reducing them to Maximum Subset Norm (MSN) problems. The MSN problem consists of computing a subset of vectors from a set of vectors with the largest possible $L_2$-norm of their sum. This problem, while still posing significant combinatorial difficulty, is subsequently solved by a Refined Triangle Inequality (RTI)-based method which essentially uses a greedy approach.

Their first theoretical result shows that the upper and lower bounds absolutely hold for any data-set. The quality of the solution, however, is based on the difference between $U_k$ and $L_k$. They then show that for exponentially tailed distributions, the ratio $\frac{U_k}{L_k}$ is at most $1 + O \left( \frac{d + k \sqrt{d}}{n} \right)$. 

This paper also includes an adaptation of their algorithm for regressions that involve categorical data called OHARE. For this algorithm, they show that for datasets satisfying certain technical properties, they can obtain a bound of $\frac{U_k}{L_k}$ of at most $1 + O\left(\frac{1}{\sqrt{\log n}}\right)$. They then experimentally validate this algorithm by demonstrating its effectiveness in the Nightlights study, where they obtain a provable lower bound on the minimum number of examples that need to be changed in order to reverse the main findings of the study. 

They also include audits of 13 other linear regressions.
weaknesses: There is, of course, a gap between the theoretical bounds on $\frac{U_k}{L_k}$ and its behavior in practice (probably due to the nature of the assumptions needed to prove bounds on this quantity). Additionally, this algorithm is probably still not feasible for very large values of $k$. However, I do view this last weakness to be quite thoroughly addressed with the comprehensive empirical evaluation provided on real regressions.

---

confidence: 3
contribution: 2
initial_rating: 5
presentation: 2
questions: 1) Line 79: The authors talk about crossing a decision boundary. In terms of the regression problem addressed in the paper, a decision boundary is slightly inappropriate. It suits better for classification problem. 
2) Line 158: ACRE and OHARE are claimed to produce nearly matching upper and lower bounds. When do they (ACRE and OHARE) produce the same bounds, and when do they perform differently? 
3) Line 267: ACRE and OHARE are proposed for n>>d. What happens for the more interesting case when n is comparable to d or a high dimensional setting? 
4) Line 352: The norm is not specified in the equation. Is it Euclidean? 
5) Line 385: The authors assume the knowledge of the population covariance matrix.  What is the approach when it is unknown? 
6) Line 464: The authors use a greedy algorithm to obtain the bound. What is the computational complexity?
soundness: 2
strengths: 1) The paper addresses an interesting problem of auditing in linear regression. 
2) Rigorous theoretical analysis using interesting proof techniques like bounds via maximal subset sum norm (MSN) or knapsack-style dynamic programming.  
3) The authors have tried experiments on real datasets.
summary: The paper presents two algorithms, ACRE and OHARE, designed to audit the robustness of Ordinary Least Squares (OLS) regression models. This robustness is assessed by determining how removing a small subset of data points can alter the regression parameter vector. For detailed comments, please see strengths, weakness, questions section.
weaknesses: 1) Line 22: Instead of saying hundreds of dimensions, the authors should specify the exact value of d and n. Later in the paper, the authors have specified that they assume n>>d. Hence, saying hundreds of dimensions sounds slightly misleading. 
2) Line 135: The authors claim no assumptions whatsoever on X and Y. Aren’t you assuming well-behaved in Definition 1 and also iid? Also, see line 159 about assumptions. 
3) Line 253: Authors claim distributions with heavier tails than sub-exponential can still be well-behaved. 
While this is true, the authors should provide more information or intuition on which heavy-tailed distributions can specifically be well-behaved and which are not. For example, the counter-statement with some heavy-tailed distribution not being well-behaved can also be true.
4) Line 295: The authors believe that the error term 1/\sqrt(\log(n)) is nearly tight. Is there a formal proof for the tightness of the bound?
5) Line 322: The authors claim that loose bounds suffice. Why? Are the proposed bounds loose? Because the authors have also claimed the bounds to be nearly tight. This is slightly contradictory. What are the challenges in deriving tighter bounds? 
6) Line 374: The authors bound the first term using a greedy algorithm. What will the computational complexity be for that? Is this greedy step just for theoretical analysis? 
7) Line 390: The authors assume the max eigenvalue of a matrix which is a function of T (an unknown variable or variable of interest). This is slightly unconventional. 
8) Line 399: Can the computational complexity be improved from O(n^2 d) to something better? 
9) Line 440: Typo: direct -> direction 
10) Line 441: The assumption on $e$ should be specified in the main theorem. 
11) Line 497: The authors criticize existing works due to their reliance on semi-definite programming. Is the computational complexity of ACRE and OHARE better than existing works? 

Although the contributions are noteworthy, it is slightly challenging to follow the paper: 
1) It is tough to follow the story. The authors propose ACRE and then criticize it for proposing OHARE. 
2) Theorem 1.2 and 1.3 are presented even before describing the respective methods.  
3) This makes it tough to develop an intuitive understanding of theorems. For example, please explain the implication of bucket size from Theorem 1.3. What role does it play? What happens if all the bucket sizes are different or if all of them are the same? 
4) It feels like the authors have tried to squeeze in a lot of information, breaking the story's flow.

---

confidence: 3
contribution: 2
initial_rating: 5
presentation: 2
questions: * How sensitive are the bounds values $U_k$ and $L_k$ to the choice of MSN-bounding algorithm?
* How do the algorithms perform when the assumptions about data distribution are violated? In lines 252-253, the authors mention that the result is valid also without the exponential decaying tail assumption. Have the authors tried the algorithm on PowerLaw synthetic data?
soundness: 2
strengths: * The authors provide an implementation in their anonymous repository. The algorithm proposed by them is able to run on consumer grade machines for datasets of sizes of order $n = 10^{4}$.
* The main text provides an excellent introduction to ACRE, explaining the idea behind its definition. The connection to MSN (Maximal Subset-sum Norm) is clearly explained.
summary: The paper presents two algorithms, ACRE (Algorithm for Certifying Robustness Efficiently) and OHARE (One-Hot aware Algorithm for certifying Robustness Efficiently), for certifying the robustness of Ordinary Least Squares (OLS) regression to sample removal in the dataset with upper and lower bounds. The first algorithm is shown to work for continuous features and the second one with a mixture of categorical and one-hot encoded features.
Given a dataset with features $X_1,...,X_n \in \mathbb{R}^d$ and labels $Y_1,...,Y_n \in \mathbb{R}$, the algorithms provide bounds on how much the OLS coefficients can change when removing $k$ samples. The authors call the maximal discrepancy between the OLS with full dataset and the one where $k$ data points are removed $\Delta_k(e)$ along the direction $e$.

The algorithms provide upper and lower bounds $U_k, L_k$ such that $L_k \leq \Delta_k(e) \leq U_k$, where $\Delta_k(e)$ represents the maximum possible change in the projection of $\beta$ onto direction $e$ when removing $k$ samples. The bounds are proven to be (logarithmically) tight under mild distributional assumptions. One main characteristic of the algorithm is Minimal Subset-sum Norm bounding.

The authors provide some empirical analysis of the efficiency of these algorithms and an explanation of the correctness of ACRE.
The analysis is based on the decomposition of $\beta - \beta_S$ into linear order and higher-order terms where $\beta_S$ is the OLS estimator after the removal of samples.
weaknesses: * The clarity of the paper is a central issue. In Section 3, the intuition for the actual definition of the algorithm explains the setting and assumptions under which the result is obtained, and the authors believe the algorithm's work is not precise. 
The connection between theoretical guarantees and practical performance needs to be better explained. for example Theorem1.2 and 1.3 could be tested in varying the numbers $n,d$ used for Table 1

* The process of obtaining the values in Table 1 needs to be clearly explained, and the table needs to be clearer. For example, while the algorithms output bounds $(U_k, L_k)$ for each $k$, how these translate to the final numbers reported needs to be clarified. The comparison methodology with baseline methods needs more detail. Any form of cross-validation missing from the experimental results

* Additional experiments can be conducted to strengthen the results.
Performance could be assessed under different types of data distributions using synthetic datasets with known ground truth.
Comparison could also be made with simple baselines (e.g., random sampling).

Meta-review: The paper considers the following very natural and simple problem regarding robustness of linear regression: Given n points, and a budget k for points to remove and a direction e, what is the maximum change that you can make along this direction by removing k points from the training data. This is a very natural robustness question which has only recently been studied. The authors introduce an algorithm to compute upper and lower bounds on this measure. There is also substantial empirical evaluation.

Decision: Accept (Poster)