Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
44,464 characters · 11 sections · 45 citation commands
Provably Auditing Ordinary Least Squares in Low Dimensions
In applied statistics and machine learning, a common task is to estimate the effect of some real-valued treatment variable (e.g. years of education) on a real-valued response variable (e.g. income) in the presence of a small number of real-valued controls. In particular, it might be especially important to understand whether the effect is positive or not; this may inform e.g. governmental policy decisions. The most basic model for this task is linear regression, and the most common estimator for this model is Ordinary Least Squares. In their seminal 1980 book Regression Diagnostics belsley2005regression, Belsley, Kuh, and Welsch observe:
While the intervening decades have seen the development of increasingly sophisticated statistical estimators, and increasingly complex machine learning models, this observation still holds.
For samples $(X_i,y_i)_{i=1}^n$ where each sample consists of a treatment $X_{i1} \in \RR$, a response $y_i \in \RR$, and $d-1$ controls $X_{i2},\dots,X_{id} \in \RR$, the linear regression model posits that there is a linear relationship $y_i = \langle X_i,\beta^*\rangle + \text{noise}$ where $\beta^*$ is the ground truth regressor. Then, $\beta^*_1$ is the posited effect of the treatment. Of course, $\beta^*$ is unknown. The Ordinary Least Squares (OLS) estimator for $\beta^*$ is
and then $(\hat{\beta}_\text{OLS})_1$ is the estimated effect of the treatment.
Under strong statistical assumptions (e.g. that the model is well-specified with exogenous error term), the OLS estimator has compelling statistical guarantees (e.g. consistency, meaning that $\hat{\beta}_\text{OLS} \to \beta^*$ in probability), and it is possible to compute statistical confidence that the estimated effect is real rather than produced by chance (e.g. by $t$-statistics). But real-world data is (almost de facto) misspecified, with nonlinearities, endogeneity, heavy-tailed errors, or even gross outliers. When these irregularities cannot be easily modeled, practitioners may simply apply the OLS estimator anyways. Hence,
Ultimately, the goal of a scientific study is to extract conclusions from data that will generalize to broader populations. Perhaps the best way of validating studies is to collect new data and replicate the results, but this is labor-intensive and often inconclusive baker20161. The simpler alternative, pursued in this work, is to audit conclusions based on the data which has already been collected.\footnote{We are auditing for robustness, which should not be confused with the (philosophically related) notion of auditing for fairness kearns2018preventing.}
In recent work, Broderick, Giordano, and Meager propose assessing the sensitivity of a given statistical analysis by computing the minimum number of samples that need to be removed from the dataset to overturn the conclusion broderick2020automatic; we call this the finite-sample stability of the conclusion. Intuitively, if the stability is very small, this means the conclusion is driven by samples from a small subpopulation, and may not generalize well. Conversely, if the stability is large, then we may have more confidence that the conclusion is broadly true. Finite-sample stability is distinct from common statistical notions of confidence broderick2020automatic, and in fact the authors discovered several important economics studies where statistically significant conclusions can be overturned by dropping less than $1\%$ of the data. Their finding underscores the importance of auditing conclusions for stability. While the oft-cited scientific “replication crisis” ioannidis2005most will not be easily solved, at the least it may be partially alleviated by thorough audits.
Unfortunately, the obvious algorithm for determining whether a conclusion derived from $n$ samples can be overturned by removing e.g. $1\%$ of the data would require rerunning the analysis $\binom{n}{n/100}$ times, which is computationally intractable even for fairly small $n$. Instead, broderick2020automatic proposes a heuristic (applicable to any differentiable $Z$-estimator) for guessing which size-$n/100$ subset of the data will most heavily influence the conclusion: essentially, they choose the $n/100$ samples with the largest individual influences. Whenever removing this subset overturns the conclusion, then the heuristic has proven that the conclusion was not stable.
What the heuristic cannot do is provably certify that a particular conclusion is stable to dropping a few samples. Even for the OLS estimator, there are simple examples where the estimate can be changed dramatically by dropping a few samples, but the greedy approach fails to identify these samples (see Section (ref) for discussion). This motivates the following question: Are there efficient algorithms with provable guarantees for estimating the finite-sample stability of an OLS estimate (concretely, the number of samples that need to be dropped to change the sign of the first coordinate of the estimated regressor)?
This question is even interesting in the low-dimensional regime. In applications to e.g. econometrics, the number of covariates may be very small maddala1992introduction; nevertheless, as evidenced by broderick2020automatic, the conclusions of these studies may be very sensitive. Even high-dimensional machine learning algorithms frequently build upon low-dimensional OLS paired with feature selection methods hall2000correlation, tropp2007signal, heinze2018variable; thus, estimating stability of OLS is a prerequisite for estimating stability of these algorithms. Our work provides provable and efficient algorithms to not only better diagnose sensitive conclusions obtained by low-dimensional OLS regressions, but also validate robust conclusions.
We are given a deterministic and arbitrary set of $n$ samples $(X_i,y_i)_{i=1}^n$, where each $X_i$ is a vector of $d$ real-valued covariates, and each $y_i$ is a real-valued response. We are interested in a single coefficient of the OLS regressor (without loss of generality, the first coordinate): in an application, the first covariate may be the treatment and the rest may be controls. The sign of this coefficient is important because it estimates whether the treatment has a positive or negative effect. Thus, we want to determine if it can be changed by dropping a few samples from the regression. Formally, we consider the fractional relaxation, where we allow dropping fractions of samples:\footnote{This relaxation does not change the interpretation of the metric. Obviously, the fractional stability lower bounds the integral stability. Moreover, there is always an optimal weight vector $w \in [0,1]^n$ with at most $d+1$ non-integral weights.}
This is the minimum number of samples (in a fractional sense) which need to be removed to zero out the first coordinate of the OLS regressor. If the OLS solution set contains multiple regressors, then it suffices if any regressor $\beta$ in the solution set has $\beta_1=0$. Our algorithmic goal is to compute $\operatorname{Stability}(X,y)$, or at least to approximate $\operatorname{Stability}(X,y)$ up to an additive $\epsilon n$ error.
By brute-force search, the (integral) stability can be computed in time $2^n \cdot \poly(n)$. However, because the complexity is exponential in the number of samples, it is computationally infeasible even when the dimension $d$ of the data is low, which is a common situation in many scientific applications. Similarly, the fractional stability (Definition (ref)) is the solution to a non-convex optimization problem in more than $n$ variables, which seems no simpler. Can we still hope for a polynomial-time algorithm in constant dimensions? We show that the answer is yes.
We also show that the exponential dependence on dimension $d$ is necessary under standard complexity assumptions:
This lower bound in particular rules out fixed-parameter tractability, i.e. algorithms with time complexity $f(d) \cdot \poly(n)$. However, it only applies to exact algorithms. In practice, it is unlikely to matter whether $\operatorname{Stability}(X,y) = 0.01n$ or $\operatorname{Stability}(X,y) = 0.02n$; in both cases, the conclusion is sensitive to dropping a very small fraction of the data. This motivates our next two algorithmic results on $\epsilon n$-additive approximation of the stability (where we think of $\epsilon>0$ as a constant). First, we make a mild anti-concentration assumption, under which the stability can $\epsilon n$-approximated in time roughly $n^{d+O(1)}$. While still not fixed-parameter tractable, this algorithm can now be run on moderate sized problems in low dimensions, unlike the algorithm in Theorem (ref).
Under this assumption, we present an $O(\epsilon n)$-approximation algorithm:
In fact, \lpalg can detect failure of the anti-concentration assumption (see Theorem (ref) for the precise statement). Moreover, the required anti-concentration is very mild. If $\epsilon,\eta>0$ are constants, the algorithm has time complexity $n^{d+O(1)}$, so long as the samples satisfy $(\epsilon,\exp(-\Omega(n)))$-anti-concentration. This is true for arbitrary smoothed data (Appendix (ref)). Finally, unlike the exact algorithm, \lpalg avoids heavy algorithmic machinery; it only requires solving linear programs.
\paragraph{Fixed-parameter tractability?} Our final result is that $\epsilon n$-approximation of the stability is in fact fixed-parameter tractable, under a stronger anti-concentration assumption.
Although this assumption is stronger than the previous, it still holds with constant $\epsilon,\delta>0$ under certain distributional assumptions on $(X_i,y_i)_{i=1}^n$, e.g. centered Gaussian mixtures with uniformly bounded condition number (Appendix (ref)).
\paragraph{Extensions.} Another model, frequently used in causal inference and econometrics, is instrumental variables (IV) linear regression. When the noise $\eta$ in a hypothesized causal relationship $y = \langle X,\beta^*\rangle + \eta$ is believed to be endogenous (i.e. correlated with $X$), a common approach sargan1958estimation, angrist1996identification, card2001estimating is to find a $p$-dimensional variable $Z$ (the instrument) for which domain knowledge suggests that $\EE[\eta|Z] = 0$. Positing that $\beta^*$ is identified by the moment condition $\EE[Z(y - \langle X,\beta\rangle)] = 0$, the IV estimator set given samples $(X_i,y_i,Z_i)_{i=1}^n$ is then \[\textsf{IV}(X,y,Z) = \{\beta \in \RR^d: Z^T (w \star (X\beta - y) = 0\}\]
where $a \star b$ denotes elementwise product, and $Z: n \times p$ and $X: n \times d$ are the matrices of instruments and covariates respectively. Stability can be defined as in Definition (ref). Although for simplicity we state all of our results for OLS (i.e. the special case $Z = X$), it can be seen that Theorem (ref) and Theorem (ref) both extend directly to the IV regression setting. See Appendix (ref) for further discussion.
\paragraph{Experiments.} We implement modifications of \netalg and \lpalg which give unconditional, exact upper and lower bounds on stability, respectively. We use these algorithms to obtain tight data-dependent bounds on stability of isotropic Gaussian datasets for a broad range of signal-to-noise ratios, and we demonstrate heterogeneous synthetic datasets where our algorithms' upper bounds are an order of magnitude better than upper bounds obtained by the prior heuristic. On the Boston Housing dataset harrison1978hedonic, we regress house values against all pairs of features. For the majority of these regressions, we bound the stability within a factor of two. On the one hand, we detect many sensitive conclusions (including some which the greedy heuristic claims are stable); on the other hand, we certify that some conclusions are stable to dropping as much as half the dataset.
\paragraph{Broader Context.} Taking a step back, OLS is important not only in its own right but also as a key building block in more complex machine learning systems, ranging from regression trees loh2011classification to generative adversarial networks mao2017least and policy iteration in linear MDPs lagoudakis2003least. Our work on estimating stability of OLS is also a first step towards estimating stability for these systems.
In Section (ref) we review related work. In Section (ref) we collect notation and formulas that will be useful later. In Section (ref) we sketch the intuition behind our algorithmic results. Section (ref) covers our experiments. In Appendices (ref), (ref), (ref), and (ref) we prove Theorems (ref), (ref), (ref), and (ref) respectively.
\paragraph{Local and global sensitivity metrics.} Post-hoc evaluation of the sensitivity of a statistical inference to various types of model misspecification has long been recognized as an important research direction. Within this area, there is a distinction between local sensitivity metrics, which measure the sensitivity of the inference to infinitesimal misspecifications of the assumed model $M_0$ (e.g. polasek1984regression, castillo2004general, belsley2005regression), and global sensitivity metrics, which measure the set of possible inferences as the model ranges in some fixed set $\mathcal{M}$ around $M_0$ (e.g. leamer1984global, tanaka1989possibilistic, vcerny2013possibilistic). For OLS in particular, there is a well-established literature on the influences of individual data points cook1977detection, chatterjee1986influential, which falls under local sensitivity analysis, since deleting a single data point is an infinitesimal perturbation to a dataset of size $n$ as $n \to \infty$. In contrast, identifying jointly influential subsets of the data (the “global” analogue) has been a long-standing challenge due to computational issues (see e.g. page 274 of belsley2005regression). Existing approaches typically focus on identifying outliers in a generic sense rather than with respect to a specific inference hadi1993procedures, or study computationally tractable variations of deletion (e.g. constant-factor reweighting leamer1984global).
\paragraph{Finite-sample stability.} Finite-sample stability is a global sensitivity metric; it essentially seeks to identify jointly influential subsets. Most directly related to our work is the prior work on heuristics for the finite-sample stability broderick2020automatic, kuschnig2021hidden. The heuristic given by broderick2020automatic (to approximate the most-influential $k$ samples) is simply the local approximation: compute the local influence of each sample at $w = \mathbbm{1}$, sort the samples from largest to smallest influence, and output the top $k$ samples. Subsequent work kuschnig2021hidden proposed refining this heuristic by recomputing the influences after removing each sample, which alleviates issues such as masking chatterjee1986influential. But this is still just a greedy heuristic, and it may fail when samples are jointly influential but not individually influential. Except under the strong assumption that the sample covariance remains nearly constant when we remove any $\epsilon n$ samples (see Theorem 1 in broderick2020automatic, which relies on Condition 1 in giordano2019swiss), the local influence approach can upper bound the finite-sample stability but cannot provably lower bound it. In fact, in Section (ref) we provide examples where the greedy heuristic of kuschnig2021hidden is very inaccurate due to instability in the sample covariance.
\paragraph{The $s$-value.} Closely related to finite-sample stability, the $s$-value gupta2021r is the minimum Kullback-Leibler divergence $D(P||P_0)$ over all distributions $P$ for which the conclusion is null, where $P_0$ is the empirical distribution of the samples. Unfortunately, while the $s$-value is an interesting and well-motivated metric, computing the $s$-value for OLS estimation appears to be computationally intractable, and the algorithms given by gupta2021r lack provable guarantees.
\paragraph{Robustified estimators.} Ever since the work of Tukey and Huber, one of the central areas of statistics has been robustifying statistical estimators to be resilient to outliers (see, e.g. huber2004robust). While a valuable branch of research, we view robust statistics as incomparable if not orthogonal to post-hoc sensitivity evaluation, for three reasons. First, samples that drive the conclusion (in the sense that deleting them would nullify the conclusion) are not synonymous with outliers: removing an outlier that works against the conclusion only makes the conclusion stronger. Indeed, outlier-trimmed datasets are not necessarily finite-sample robust broderick2020automatic. Rather, finite-sample stability (along with the $s$-value gupta2021r), in the regime where a constant fraction of samples is removed, may be thought of as a measure of resilience to heterogeneity and distribution shift.
Second, it is unreasonable to argue that using robustified estimators obviates the need for sensitivity evaluation. Robust statistics has seen a recent algorithmic revival, with the development of computationally efficient estimators, for problems such as linear regression, that are robust in the strong contamination model (e.g. klivans2018efficient,diakonikolas2019sever, bakshi2021robust). However, even positing that the strong contamination model is correct, estimation guarantees for these algorithms require strong, unverifiable (and unavoidable klivans2018efficient) assumptions about the uncorrupted data, such as hypercontractivity. Sensitivity analyses should support modeling assumptions, not depend upon them.
Third and perhaps most salient, classical estimators such as OLS are ubiquitous in practice, despite the existence of robust estimators. This alone justifies sensitivity analysis of the resulting scientific conclusions.
\paragraph{Distributionally robust optimization.} Our motivations mirror a recent line of work in machine learning sinha2017certifying, duchi2018learning, cauchois2020robust, jeong2020robust which suggests that the lack of resilience of Empirical Risk Minimization to distribution shift can be mitigated by minimizing the supremum of risks with respect to distributions near the empirical training distribution (under e.g. Wasserstein distance or an $f$-divergence). Again, this approach of robustifying the estimator is valuable but incomparable to sensitivity analysis.
For vectors $u,v \in \RR^m$, we let $u \star v$ denote the elementwise product $(u\star v)_i = u_i v_i$. Throughout the paper, we will frequently use the closed-form expression for the (weighted) OLS solution set \[\OLS(X,y,w) = \{\beta \in \RR^d: X^T (w \star (X\beta - y)) = 0\}\] where $X: n \times d$ is the matrix with rows $X_1,\dots,X_n$. In particular, setting $\lambda = \beta_{2:d}$, this means that the finite-sample stability can be rewritten as
where (here and throughout the paper) $\tilde{X}: n \times (d-1)$ is the matrix with columns $(X^T)_2,\dots,(X^T)_d$.
\paragraph{An exact algorithm.} Our main tool for Theorem (ref) is the following special case of an important result due to renegar1992computational on solving quantified polynomial systems of inequalities:
Roughly, for a constant number of quantifier alternations, a quantified polynomial system can be decided in time exponential in the number of variables. Unfortunately, a naive formulation of the expression $\operatorname{Stability}(X,y) \leq k$, by direct evaluation of Equation (ref), has $n+d-1$ variables: \[\exists \lambda \in \RR^{d-1}, w \in [0,1]^n: \sum_{i=1}^n w_i \geq n-k \land X^T(w\star (\tilde{X}\lambda-y)) = 0.\] Intuitively, it may not be necessary to search over all $w \in [0,1]^n$; for fixed $\lambda$, the maximum-weight $w$ is described by a simple linear program. Formally, the linear program can be rewritten (Lemma (ref)) by the separating hyperplane theorem, so that the overall expression becomes:
Now, for fixed $\lambda$ and $u$, the maximum-weight $w$ has very simple description: it only depends on the relative ordering of the $n$ summands $(\langle \tilde{X},\lambda\rangle - y_i)\langle X_i,u\rangle$. By classical results on connected components of varieties, since the summands have only $2d-1$ variables, the number of achievable orderings is only $n^{\Omega(d)}$ rather than $n!$, and the orderings can be enumerated efficiently milnor1964betti, renegar1992computational. This allows the quantifier over $w \in [0,1]^n$ to be replaced by a quantifier over the $n^{\Omega(d)}$ achievable orderings, after which Theorem (ref) implies that the overall expression can be decided in time $n^{\Omega(d^3)}$. See Appendix (ref) for details.
\paragraph{Approximation via partitioning.} Next, we show how to avoid the heavy algorithmic machinery used in the previous result. For Theorem (ref), the strategy is to partition the OLS solution space $\RR^{d-1}$ into roughly $n^d$ regions, such that if we restrict $\lambda$ to any one region, the bilinear program which defines the stability can be approximated by a linear program.
Concretely, we start by writing the formulation ((ref)) as
This has a nonlinear (and nonconvex) constraint due to the pointwise product between $w$ and the residual vector $\tilde{X}\lambda - y$. Thus, we can introduce the change of variables $g_i = w_i(\langle \tilde{X}_i,\lambda\rangle - y_i)$ for $i \in [n]$. This causes two issues. First, the constraint $0 \leq w_i \leq 1$ becomes $0 \leq g_i/(\langle X_i,\lambda\rangle - y_i) \leq 1$, which is no longer linear. To fix this, suppose that instead of maximizing over all $\lambda \in \RR^{d-1}$, we maximize over a region $R \subseteq \RR^{d-1}$ where each residual $\langle \tilde{X},\lambda\rangle - y_i$ has constant sign $\sigma_i$. The constraint $0 \leq w_i \leq 1$ then becomes one of two linear constraints, depending on $\sigma_i$. Let $V_R$ denote the value of Program (ref) restricted to $\lambda \in R$. Then with the change of variables, we have that \[V_R = \sup_{g \in \RR^n, \lambda \in R} \left\{\sum_{i \in [n]} \frac{g_i}{\langle X_i,\lambda\rangle - y_i} \middle |
\right\},\] with the convention that $0/0 = 1$. Now the constraints are linear. Unfortunately, (and this is the second issue), the objective is no longer linear. The solution is to partition the region $R$ further: if the region were small enough that every residual $\langle X_i,\lambda\rangle - y_i$ had at most $(1\pm \epsilon)$-multiplicative variation, then the objective could be approximated to within $1\pm\epsilon$ by a linear objective.
How many regions do we need? Let $M = \norm{X\beta^{(0)}-y}_2$ be the unweighted OLS error. If all the residuals were bounded between $\delta M/\sqrt{n}$ and $M$ in magnitude, for all $\lambda \in \RR^{d-1}$, then the regions could be demarcated by $O(n\log_{1+\epsilon}(n/\delta))$ hyperplanes, for a total of $O(n\log_{1+\epsilon}(n/\delta))^d$ regions. Of course, for some $\lambda$, some residuals may be very small or very large. But $(\epsilon,\delta)$-anti-concentration implies that for every $\lambda$, at most $\epsilon n$ residuals are very small, and it can be shown that if $\lambda$ is a weighted OLS solution, the total weight on samples with large residuals is low. Thus, for any region, we can exclude from the objective function the samples with residuals that are not well-approximated within the region, and this only affects the objective by $O(\epsilon n)$.
This gives an algorithm with time complexity $(n\epsilon^{-1}\log(1/\delta))^{d+O(1)}$. To achieve the time complexity in Theorem (ref), where the $\log(n/\delta)$ is additive rather than multiplicative, we use subsampling. Every residual is still partitioned by sign, but we multiplicatively partition only a random $\tilde{O}(d/\epsilon)$-size subset of the residuals. Intuitively, most residuals will still be well-approximated in any given region. This can roughly be formalized via a VC dimension argument, albeit with some technical complications. See Section (ref) for details and Appendix (ref) for formal pseudocode of the algorithm \lpalg.
\paragraph{Net-based approximation.} The algorithm for Theorem (ref) is intuitively the simplest. For any fixed $\lambda \in \RR^{d-1}$, Program (ref) reduces to a linear program with value denoted $S(\lambda)$. Thus, an obvious approach is to construct a net $\mathcal{N} \subseteq \RR^{d-1}$ in some appropriate metric, and compute $\min_{\lambda \in \mathcal{N}} S(\lambda)$. This always upper bounds the stability, but to prove that it's an approximate lower bound, we need $S(\lambda)$ to be Lipschitz under the metric.
The right metric turns out to be \[d(\lambda,\lambda') = \norm{\frac{\tilde{X}\lambda - y}{\norm{\tilde{X}\lambda-y}_2} - \frac{\tilde{X}\lambda' - y}{\norm{\tilde{X}\lambda'-y}_2}}_2.\] Under this metric, $\RR^{d-1}$ essentially embeds into a $d$-dimensional subspace of the Euclidean sphere $\mathcal{S}^{n-1}$, and therefore has a $\gamma$-net of size $O(1/\gamma)^d$. Why is $S(\lambda)$ Lipschitz under $d$? First, if $\tilde{X}\lambda-y$ equals $\tilde{X}\lambda' - y$ up to rescaling, then it can be seen from Program (ref) that $S(\lambda) = S(\lambda')$. More generally, if the residuals are close up to rescaling, we apply the dual formulation of $S(\lambda)$ from expression ((ref)): \[n-S(\lambda) = \inf_{u \in \RR^d} \sup_{w \in [0,1]^n} \norm{w}_1: \sum_{i=1}^n w_i(\langle \tilde{X}_i,\lambda\rangle - y_i)\langle X_i,u\rangle \geq 0.\] For any $u$, the optimal $w$ for $\lambda$ and $u$ can be rounded to some feasible $w'$ for $\lambda'$ and $u$ without decreasing the $\ell_1$ norm too much, under strong anti-concentration. This shows that $S(\lambda)$ and $S(\lambda')$ are close. See Appendix (ref) for details and Appendix (ref) for formal pseudocode of the algorithm \netalg.
In this section, we apply (modifications of) \netalg and \lpalg to various datasets. The modifications are made for practical reasons, and while our theory will no longer apply, we emphasize that the modified \netalg will still provide an unconditional upper bound on stability (referred to henceforth as “net upper bound”), and the modified \lpalg now provides an unconditional lower bound (“LP lower bound”).
As a baseline upper bound, we implement the greedy heuristic of kuschnig2021hidden which refines broderick2020automatic. We are not aware of any prior work on lower bounding stability, so we implement a simplification of our full lower bound algorithm as a baseline. See Appendix (ref) for implementation details and hyperparameter choices of our algorithms and baselines. All error bars are at $25$th and $75$th percentiles over independent trials.
\paragraph{Heterogeneous data.} We start with a simple two-dimensional dataset with two disparate subpopulations, where the greedy baseline fails to estimate the stability but our algorithms give tight estimates. For parameters $n$, $k$, and $\sigma$, we generate $k$ independent samples $(X_i,y_i)$, where $X_i \in \RR^2$ has independent coordinates $X_{i1} \sim N(-1, 0.01)$ and $X_{i2} \sim N(0,1)$, and $y_i = X_{i1}$. Then, we generate $n-k$ independent samples $(X_i,y_i)$ where $X_{i1} = 0$ and $X_{i2} \sim N(0,1)$, and $y \sim N(0,1)$. It always suffices to remove the first subpopulation, so the stability is at most $k$. However, the first subpopulation has small individual influences, because the OLS regressor on the whole dataset can nearly interpolate the first subpopulation. Thus, we expect that the greedy algorithm will fail to notice the first subpopulation, and therefore remove far more than $k$ samples.
Indeed, this is what happens. For $n = 1000$ and $k$ varying from $10$ to $500$, we compare our net upper bound and LP lower bound with the baselines. As seen in Figure (ref), our methods are always better than the baselines, and certifiably approximate the stability within a small constant factor. In the regime where $k$ is small, our upper bound outperforms the greedy upper bound by a factor of $30$.
\paragraph{Covariance shift.} In the previous example, removing $k$ samples caused a pathological change in the sample covariance; it became singular. However, even modest, constant-factor instability in the sample covariance can cause the greedy algorithm to fail; see Appendix (ref) for details.
\paragraph{Isotropic Gaussian data.} Instability can arise even in homogeneous data, as a result of a low signal-to-noise ratio broderick2020automatic. But when the noise level is low, can we certify stability? For a broad range of noise levels, we experimentally show that this is the case. Specifically, for $d \in \{2,3\}$ and noise parameter $\sigma$ ranging from $0.1$ to $10$, we generate $n$ independent samples $(X_i,y_i)_{i=1}^n$ where $X_i \sim N(0,I_d)$ and $y_i = \langle X_i,\mathbbm{1}\rangle + N(0,\sigma^2)$. For $d = 2$ and $n = 1000$ (Figure (ref)), our LP lower bound is nearly tight with the upper bounds, particularly as the noise level increases (in comparison, the baseline lower bound quickly degenerates towards zero). For $d=3$ and $n = 500$ (Figure (ref)), the bounds are looser for small noise levels but still always within a small constant factor.
The Boston Housing dataset harrison1978hedonic,gilley1996harrison consists of data from $506$ census tracts of Greater Boston in the 1970 Census. There are $14$ real-valued features, one of which---the median house value in USD 1000s---we designate as the response. Unfortunately the entire set of features is too large for our algorithms, so for our experiments we pick various subsets of two or three features to use as covariates.
\paragraph{A Tale of Two Datasets.} We exemplify our results with two particular feature subsets. First, we investigate the effect of zn (percentage of residential land zoned for large lots) on house values, controlling for rm (average number of rooms per home) and rad (highway accessibility index) but no bias term. On the entire dataset, we find a modest positive effect: the estimated coefficient of zn is roughly $0.06$. Both the greedy heuristic and our net algorithm find subsets of just $8\%$ of the data ($38$-$40$ samples) which, if removed, would nullify the effect. But is this tight, or could there be a much smaller subset with the same effect? Our LP lower bound certifies that removing at least $22.4$ samples is necessary.
Second, we investigate the effect of zn on house values, this time controlling only for crim (per capita crime rate). Our net algorithm finds a subset of just $27\%$ of the data which was driving the effect, and the LP lower bound certifies that the stability is at least $8\%$. But this time, the greedy algorithm removes $90\%$ of the samples, a clear failure. What happened? Plotting zn against crim reveals a striking heterogeneity in the data: $73\%$ of the samples have $\textsf{zn} = 0$, and the remaining $27\%$ of the samples (precisely those removed by the net algorithm) have $\textsf{crim} < 0.83$, i.e. very low crime rates. As in the synthetic example, this heterogeneity explains the greedy algorithm's failure. But heterogeneity is very common in real data: in this case, it's between the city proper and the suburbs, and in fact the OLS regressors of these two subpopulations on all $13$ features are markedly different (Appendix (ref)). Thus, it's important to have algorithms with provable guarantees for detecting when heterogeneity causes (or doesn't cause) unstable conclusions.
\paragraph{All-feature-pairs analysis.} To be thorough, we also apply our algorithms to all $156$ ordered pairs of features. For each pair, we regress the response (i.e. median house value) against the two features by Ordinary Least Squares, and we use our algorithms on this $2$-dimensional dataset to estimate how many samples need to be removed to nullify the effect of the first feature on the response. We also compare to the greedy upper bound. See Figure (ref) for a perspective on the results. In each figure, each point corresponds to the results of one dataset. The left figure plots the net upper bound against the greedy upper bound: we can see that our net algorithm substantially outperforms the greedy heuristic on some datasets (i.e. finds a much smaller upper bound) and never performs much worse. The right figure plots the LP lower bound against the net upper bound (along with the line $y=x$). For a majority of the datasets, the upper bound and lower bound are close. Concretely, for $116$ of the $156$ datasets, we certifiably estimate the stability up to a factor of two \--- some are sensitive to removing less than $10$ samples, and some are stable to removing even a majority of the samples.
\paragraph{Acknowledgments.} We thank Allen Liu, Morris Yau, Fred Koehler, Jon Kelner, Sam Hopkins, and Vasilis Syrgkanis for valuable conversations on this and related topics.