Provably Auditing Ordinary Least Squares in Low Dimensions
Authors: Ankur Moitra
Dhruv Rohatgi
Venue: ICLR 2023 poster (2023)
Keywords: stability
linear regression
ordinary least squares
robustness

Abstract: Auditing the stability of a machine learning model to small changes in the training procedure is critical for engendering trust in practical applications. For example, a model should not be overly sensitive to removing a small fraction of its training data. However, algorithmically validating this property seems computationally challenging, even for the simplest of models: Ordinary Least Squares (OLS) linear regression. Concretely, recent work defines the stability of a regression as the minimum number of samples that need to be removed so that rerunning the analysis overturns the conclusion (Broderick et al., 2020), specifically meaning that the sign of a particular coefficient of the OLS regressor changes. But the only known approach for estimating this metric, besides the obvious exponential-time algorithm, is a greedy heuristic that may produce severe overestimates and therefore cannot certify stability. We show that stability can be efficiently certified in the low-dimensional regime: when the number of covariates is a constant but the number of samples is large, there are polynomial-time algorithms for estimating (a fractional version of) stability, with provable approximation guarantees. Applying our algorithms to the Boston Housing dataset, we exhibit regression analyses where our estimator outperforms the greedy heuristic, and can successfully certify stability even in the regime where a constant fraction of the samples are dropped.

Published as a conference paper at ICLR 2023
PROVABLY AUDITING ORDINARY LEAST SQUARES IN LOW DIMENSIONS
Ankur Moitra & Dhruv Rohatgi Massachusetts Institute of Technology (moitra, drohatg! @mit edu
ABSTRACT
Auditing the stability of a machine learning model to small changes in the training procedure is critical for engendering trust in practical applications. For example, a model should not be overly sensitive removing small fraction of its training data. However, algorithmically validating this property seems computationally challenging, even for the simplest of models: Ordinary Least Squares (OLS) linear regression. Concretely, recent work defines the stability of a regression as the minimum number of samples that need to be removed so that rerunning the analysis overturns the conclusion (Broderick et al., 2020), specifically meaning that the sign of particular coefficient of the OLS regressor changes. But the only known approach for estimating this metric, besides the obvious exponential- time algorithm, is a greedy heuristic that may produce severe overestimates and therefore cannot certify stability. We show that stability can be efficiently certified in the low-dimensional regime: when the number of covariates is constant but the number of samples is large, there are polynomial-time algorithms for estimating (a fractional version of) stability, with provable approximation guarantees. Applying our algorithmss the Boston Housing dataset, we exhibit regression analyses where our estimator outperforms the greedy heuristic, and can successfully certify stability even in the regime where constant fraction of the samples are dropped.
INTRODUCTION
key facet of interpretability of machine learning models is understanding how different subsets of the training data influence the Iearned model and its predictions. Computing the influences of individual training points has been shown to a useful tool for enhancing trust in the model (Zhou et al., 2019), for tracing the origins of model bias (Brunet et al., 2019), and for identifying mislabelled training data and other model debugging (Koh & Liang, 2017). Modelling the influence of groups of training points has applications to measuring fairness (Chen et al., 2018), vulnerability to contamination of multi- sourti- training data (Hayes & Ohrimenko, 2018). and (most relevant to this paper) identification of unstable predictions (Ilyas et al., 2022) and models (Broderick et al., 2020). In high-stakes machine learning application, it would likely be alarming if some data points were so influential that the removal of, say, 1% of the training data dramatically changed the model. An ideal, trustworthy machine learning pipeline therefore should include validation that this does not happen.
But the obvious algorithm for checking if model trained n data points exhibits this instability would require computing the group influences of (n//1o)) different subsets of the data, which is computationally infeasible even for fairly small n. Instead, current methods for estimating the stability of a model simply use the first-order approximation of group influence: namely, the sum of individual influences of data points in the group. With this approximation, vulnerability of model to dropping an data points is heuristically estimated by dropping the an most individually influential data points (Broderick et al., 2020; llyas et al., 2022).
This heuristic can be thought of as using "local" stability as a proxy for "global" stability, and has found substantial anecdotal success in diagnosing unstable models. Unfortunately, for correlated groups of data points, the first-order approximation the group influence is often an underestimate (Koh et al., 2019), so large local stability does not actually certify that model is provably stable to

---

Published as a conference paper at ICLR 2023
removing small subsets of data. In fact, stability certification is challenging and open problem even in the simplest of models: linear regression via Ordinary Least Squares (OLS).
Concretely, given regression dataset, natural metric for the stability of the OLS regressor is the minimum number of data points that need to be removed from the dataset to flip the sign of a particular coefficient of the regressor (e.g., in causal inference, the coefficient measuring the treatment effect). Recent work has used the local stability heuristic to diagnose unstable OLS regressions in several prominent economics studies (Broderick et al., 2020), identifying examples where even statistically significant conclusion can be overturned by removing less than 1% of the data points. However, the converse question of validating stable conclusions remains unaddressed:
Given a regression dataset, can we efficiently certify non-trivial lower bounds on the stability of the OLS regressor?
Our work takes steps towards addressing this question, via the following contributions:
We introduce natural fractional relaxation of the above notion of OLS stability, where we allow removing fractions of data points, and seek to minimize the total removed weight. We call this finite finite-sample stability, and henceforth refer to the prior notion as "integral" stability,
We develop approximation algorithms for the finite-sample stability, with (a) provable guarantees under reasonable anti-concentration assumptions on the dataset, and (b) running time polynomial in the size of the dataset, so long the dimension of the data is constant (in contrast, the naive algorithm is exponential in the size of the dataset). Moreover, we prove that (at least for exact algorithms) exponential dependence of the running time on the dimension is unavoidable under standard complexity assumptions.
We use modifications of our algorithms compute assumption-free upper and lower bounds on the finite-sample stability of several simple synthetic and real datasets, achieving tighter upper bounds than prior work and the first non-trivial lower bounds, i.e. certifications that the OLS regressor is stable.
Why define stability this way? The definition of integral stability was introduced in (Broderick al., 2020), along with several variants (e.g. smallest perturbation which causes the first coordinate to lose significance). We choose the definition based on the sign of the first coordinate, because has clear practical interpretation does the first covariate positively negatively affect the response? which does not depend on choice of additional parameters such as significance level.
We study the fractional relaxation so that the stability is defined by a continuous optimization problem. Note that certifying a lower bound on fractional stability immediately certifies a lower bound on the integral stability; we will see later (Remark 3.1) that near-converse holds in low dimensions.
Why is low-dimensional regression important? Given that much of machine learning happens in high-dimensional settings, where the number of covariates can even be larger than the number of datapoints, it is natural to wonder why low-dimensional settings are still important. First, in application areas such as econometrics, linear regressions with as few as two to four covariates are very common (Britto al., 2022; Bianchi & Bigio, 2022: Hopenhayn et al., 2022), often serving as proofs-of-concep for more complex models. Second, even settings where the number of covariates larger, it is often the expectation that few covariates relevant. In such applications, analysis often consists of a variable selection step followed by regression on much-reduced set of covariates (Cai & Wang, 2011). In all these settings, understanding the stability of an estimator is important, and our work gives some of the first provable guarantees that avoid making strong distributional assumptions. Moreover our lower bounds show that certifying stability of truly high- nigh-dimensional models, even linear ones, is intractable.
1.1 FORMAL PROBLEM STATEMENT
We are given a deterministic and arbitrary set of n samples (X,, y.) where each X, is a vector of d rellevalued covariates, and each y/ is real-valued response. We are interested in 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
2

---

Published as a conference paper at ICLR 2023
is important because it estimates whether the treatment has positive negative effect. Thus, 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:
Definition 1.1. Fix (x),3) with XI...., Xx E Fd and yI..... R. For any w E [0, 1]", the weight-w OLS solution set (X,; (x,;,))==1 is
OLS(X, w)
arg min BER 1=1
(., y.)
The finite-sample stability of (X,, y:) is Stability(X,2) y) = "E(o.I]".BBRR BERR fn- |wll : 31 = and 3 E OLS(X, y, w))
inf
This is the minimum number of samples (in fractional sense) which need to be removed zero out the first coordinate of the OLS regressor. If the OLS solution set contains multiple regressors, then it suffices if any regressor 3 in the solution set has 31 0. Our algorithmic goal is to compute Stability (X, v), or at least to approximate Stability (X, y) up to an additive en error.
1.2 RESULTS
By brute-force search, the (integral) stability can be computed in time 2" poly((). 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 1.1) is the solution 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.
Theorem 1.2. There is an n ((ddd) -time algorithm which, given n arbitrary samples (Xi, W)) with X,, XI,...,,, E R and ,y....,yn E IR, and given k 2 0, decides whether Stability (X, v) S k.
We also show that the exponential dependence on dimension d is necessary under standard complexity assumptions:
Theorem 1.3. Under the Exponential Time Hypothesis, there is no no -time algorithm which, given (1) and k 2 0, decides whether Stability(X, y) $ k.
This theorem in particular rules out seeed-parameter tractability, i.c. algorithms with time complexity f(d) (n)y(n). However, it only applies to exact algorithms. In practice, it is unlikely to matter whether Stability(X, (X,,itty(x,)) y) 0.0ln or Stability (X, y) 0.02n; in both cases, the conclusion is sensitive to dropping very small fraction of the data. This motivates our next two algorithmic results on en-additive approximation of the stability (where we think of c > 0 as a constant). First, we make mild anti-concentratiin assumption, under which the stability can en-approximated in time roughly n d+1011) While still not fixed- -xedmmerrmeter tractable, this algorithm can now be run on moderate sized problems low dimensions, unlike the algorithm in Theorem 1.2.
Assumption A. Let 8 > 0. We say that samples (x)y satisfy (c, anti-coonnnnnninooonnnn if for every 3 E [dd it holds that
[iel)::X,,. - -yll < Vn X3 (880)~111 -112] S en,
where X X d is the matrix with rows X1.... Xn, and ((0) E OLS(X, y, 1) is any unweighted OLS regressor of y against X.
See Appendix 1.1 for discussion of when Assumption holds. Under this assumption, we present an aaproxxapooiinaioonn algorithm:
Theorem 1.4, For any e, 5,1) > O, there is an algorithm PARTITIONANDAPROOX with time complex- ity
1 d+0(1)
Cd 2
n
log
log en
3

---

Published as a conference paper at ICLR 2023
which, given e, 6, 7,, and samples (X,, M1 satisfying (c.0))aatticconceentin returns estimate S such that with probability least 7).
|$ Stability(X, S 12cn + 1.
In fact, PARTITIONANDAPPRO, also can detect when Assumption A fails (see Theorem D.6 for precise statement), so can be used to compute unconditional lower bounds on stability with high probability (where the lower bound is provably tight if the data satisfies anti-concentration, Moreover, as discussed in Appendix F.I, the required anti concentration is very mild. If c, > are constants, the algorithm has time complexity n nd)0(1)) so long as the samples satisfy (e, -)()))))))mm concentration. This is true for arbitrary smoothed data. Finally, unlike the exact algorithm, PARTI- TIONANDAPPROX avoids heavy algorithmic machinery; it only requires solving linear programs.
Fixed-pare sieed-parameter tractability? Our final result is that en-approximation of the stability is in fact fixed-pe skeddparaamter tractable, under a stronger anti-concentration assumption. Assumption B. Eet > O. We say that samples (X,, y:) satisfy (e, --strong anti-concenttattion for every 3 E Rd+1,, it holds that
{ i [
< Vn
X 1x3112/
S Sen
where X (d+ 1) is the matrix with columns (X ) (x))
This assumption holds with constant e, under certain distributional assumptions on (X), e.g. centered Gaussian mixtures with uniformly bounded condition number (Appendix F.2).
Theorem 1.5. For any 8 > O, there is a (Vd/ (c(2)))) poly(()
(ply()) (timme algoritim NETAPPROX which, given e, and samples (X,, y/ ) satisfying (e, 8) -sstrong anti-concentration returns estimate $
satisfying Stability(X, < $ < Stability(X, y) 3en 1. Moreover, Stability(X, y) $ holds for arbitrary (i,y:)]==
Extensions. Another model, frequently used in causal inference and econometrics, is instrumental variables (IV) linear regression. When the noise in a hypothesized causal relationship believed to be endogenous (i.e. correlated with X), common approach (Sargan, 1958; Angrist et al., 1996; Card, 2001) 1) is to find a p-dimensional variable Z (the instrument) for which domain knowledge suggests that Emn) Z] O. Positing that 5* is identified by the moment condition EyZ(y- - (X, B] = 0, the IV estimator set given samples (X,, Vi. 2,) is then
(X, (x,8*)+ 1
+
IV(X, 2) {8 e æd _T (w (X3 y) 0}
where a denotes elementwise product, and Z : n x and X n x d are the matrices of instruments and covariates respectively. Stability can be defined in Definition 1.1. Although for simplicity state all of our results for OLS (i.e. the special case Z = X), it can be seen that Theorem 1.2 and Theorem 1.5 both extend directly the IV regression setting. See Appendix G for further discussion.
Experiments. We implement modifications NETAPPROX and PARTITIONA NDAPROO 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-nowl 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 (Harrison & Rubinfeld, 1978). we regress house values against all pairs of features. For the majority of these regressions, we bound the stability within a factor of two. the one hand, we detect many sensitive conclusions (including some which the greedy heuristic claims stable); the other hand, certify that some conclusions stable to dropping much as half the dataset.
1.3 ORGANIZATION
In Section 2 review related work. In Section 3 collect notation and formulas that will be useful later. In Section 4 we sketch the intuition behind our algorithmic results. Section covers our experiments. In Appendices B, C, D, and E we prove Theorems 1.2, 1.3, 1.4, and 1.5 respectively.
4

---

Published as a conference paper at ICLR 2023
2 RELATED WORK
There a rich literature on topics related to finite-sample stability, including sensitivity analysis and robustness to distribution shift and contamination. Due space constraints, here only discuss the works most relevant to ours, and we posipone broader discussion to Appendix A.
Most directly related is the prior work heuristics for the (integral) stability (Broderick al., 2020; Kuschnig et al., 2021). The heuristic given by Broderick et al. (2020) (to approximate the most-influentia& & samples) is simply the local approximation: compute the local influence of each sample at u = 1, sort the samples from largest to smallest influence, and output the top h samples. Subsequent work (Kuschnig al., 2021) refines this heuristic by recomputing the influences after removing each sample, which alleviates issues such as masking (Chatterjee & Hadi, 1986). But this still just a greedy heuristic, and may fail when samples are jointly but not individually influential. Except under the strong assumption that the sample covariance remains nearly constant when we remove any en samples (see Theorem Giordano al. (2019)), the local influence approach can upper bound the maite- stability but cannot provably lower bound it. In fact, in Section we provide examples where the greedy heuristic of Kuschnig et al. (2021) is very inaccurate due instability in the sample covariance.
Broderick et al. (2020), which relies on Condition in
Closely related to finite-sample stability, the s-value (Gupta & Rothenhäusler, 2021) is the minimum Kullback-Leibe divergence ((P over all distributions P for which the conclusion is null, where Po 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 Gupta & Rothenhäusler (2021) lack provable guarantees.
3 PRELIMINARIES
For vectors u, v E IR",, we let u * o denote the elementwise product (u * v); = u, Throughout the paper, we will frequently the closed-form expression for the (weighted) OLS solution set
OLS(X, y, w) {B E Rd X (w (X3 - v)) 0}
where : n d is the matrix with rows XI, Xn. particular, setting = 33.ds this means that the finite-sample stability can be rewritten as
inf Stability(X,v) "EI
(n - ||||| x2 tux (XA- 4)) 0}
(I)
where (here and throughout the paper) X x (d-1) 1) is the matrix with columns (X ....
(X")a
Remark 3.1. As previously mentioned, the finite-sample stability always lower bounds the integral stability (the minimum number of samples that need to be removed to make the first coordinate of the regressor change sign), by continuity of the OLS solution in w. Additionally, it can be seen from Equation 1 that partial converse holds in low dimensions. For any feasible (w, A), the set of w' such that (a) (w', X) is feasible, and (b) || = ||w0 I| has the form [0, Td V for some subspace V S T of codimension at most d + 1. Thus, there is some w' [0, 1) n V with at most d + 1 non-integral weights. If Stability v) = an, then w' witnesses that the first coordinate of the OLS regressor can be zeroed out by downweighting at most an + d + samples.
4 OVERVIEW OF ALGORITHMS
An exact algorithm. Our main tool for Theorem 1.2 is the following special case of an important result due to Renegar (1992) on solving quantified polynomial systems of inequalities: Theorem 4.1 (Renegar (1992)). Given an expression
Vx E I : By E R P(x, y).
where P(z, is system of polynomial inequalities with maximum degree d, the truth value of the expression can be decided in time ((1m)))))))
"Thi in the real number model; similar statement be made in the bit complexity model.
5

---

Published as a conference paper at ICLR 2023
Roughly, for constant number of quantifier alternations, quantified polynomial system can be decided in time exponential in the number of variables. Unfortunately, a naive formulation of the expression Stability(XX, y) k. by direct evaluation of Equation has + d - 1 variables:
YE R [0, 1]"::"""" Σ 2 and X" (w (XA -y) 0.
Intuitively, it may not be necessary to search over all w E [0, 1]"; for fixed X, the maximun-weight w described by simple linear program. Formally, the linear program be rewritten (Lemma B.1) by the separating hyperplane theorem, so that the overall expression becomes:
YE E RR-11 e E B 3w E [0, 1J" ||ll 2 n - k and Σ wi wi((Xi, 3i (x;, 2 0. (2)
Now, for fixed À and the maxinum-weigh w has very simple description: it only depends on the relative ordering of the n summands ((X, (X,, u). By classical results connected components of varieties, since the summands have only 2d 1 variables, the number of achievable orderings is only (dd(d) rather than n!, and the orderings can be enumerated efficiently (Milnor, 1964; Renegar, 1992). This allows the quantifier over [0, 1J" to be replaced by quantifier over the (2) achievable orderings, after which Theorem 4. implies that the overall expression can be decided in time n (dd)) See Appendix B for details.
Approximation via partitioning. Next, we show how to avoid the heavy algorithmic machinery used in the previous result. For Theorem 1.4, the strategy of AARTITIONANDAPPEEA the OLS solution space dR into roughly n regions, such that if we restrict to any one region, the bilinear program which defines the stability can be approximated by linear program.
is partition
Concretely, we start by writing the formulation (1) as
- --Staaillty(x.u))
sup - -
Σ (E(n)
= X (w (XA (
y))=
This has nonlinear (and nonconvex) constraint due to the pointwise product between and the residual vector XA y. Thus, i [n). This causes two issues. First, the constraint 0 S w; S becomes S 0i/((X, A) - y; S which is no longer linear. To fix this, suppose that instead of maximizing
introduce the change of variables 9i w. (X,, y)) for
all À Rd-1,
maximize over region R R where each residual () y: has constant sign og. The constraint 0 < wy < then becomes one of two linear constraints, depending on o;. Let VR denote the value of Program 3 restricted to X R. Then with the change of variables, we have that
XTg 0
R
VR
sup yERn.EER AER
0 05 S (X,,)) ys (X,, y: S. 50 [e :0
X, X)
Vk
1E
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: the region were small enough that every residual (X,, A) yi had at most (1 + e)-multiplicative variation, then the objective could be approximated within + by a linear objective.
How many regions do need? Let M |x3(0) 2/|2 be the unweighted OLSS error. If all the residuals were bounded between 6 6M/vni and M in magnitude, for À E Rbd-1 then the regions could be demarcated by O(n log1+= (n/))) hyperplanes, for total of O (1//)) regions. Of course, for some A, some residuals may be very small or very large. But (c. D--ani-conncetratioon implies that for every X, at most residuals
very small, and it can be shown that if À is 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(en).
This gives an algorithm with time complexity (ne log(1/8)) +00) To achieve the time complexity in Theorem 1.4, where the log(n/8) /8) is additive rather than multiplicative, subsampling. Every

---

Published as a conference paper at ICLR 2023
residual is still partitioned by sign, but we multiplicatively partition only random ((d/c)-size subset of the residuals. Intuitively, most residuals will still be well-approximated in any given region. This can roughly be formalized via VC dimension argument, albeit with some technical complications. See Section D for details and Appendix J for formal pseudocode of the algorithm.
Net-based approximation. The algorithm NETAPPROX for Theorem 1.5 intuitively the simplest. For any fixed À E R Program reduces to a linear program with value denoted S(A). Thus, an obvious approach is construct a net N Rd in some appropriate metric, and compute minxeN SA). This always upper bounds the stability, but to prove that it's an approximate lower bound, we need S()) to be Lipschitz under the metric.
The right metric turns out to be
XX y
XX y XX
d(,
X 2/2
Under this metric, TRd-1 essentially embeds into d-dimensional subspace of the Euclidean sphere Sn-1, and therefore has a 7-net of size 0(1/7)². Why is S(A) Lipschitz under d? First, if X equals XX - up to rescaling, then it can be seen from Program that S()) S(X). More generally, if the residuals are close up rescaling, we apply the dual formulation of ( from expression (2):
S(A) inf sup wERe (e(1)]
|well
Σ wa ( -i)( u) 2 O. 1=1
For any u, the optimal w for À and u can be rounded to some feasible w' for X' and u without decreasing the f norm too much, under strong anti-concentration. This shows that S(A) and S(X) are close. See Appendix for details and Appendix for formal pseudocode NETAPPROX.
EXPERIMENTS
In this section, we apply (modifications of) NETAPPROX and PARTITIONANDAPROO to several datasets in two and three dimensions. Due to space constraints, we defer detailed discussion of the algorithmic modifications to Appendix I.1; we simply note that the modifications are made to improve practical efficiency and usability. Most saliently, the modified algorithms do not rely on Assumptions A and B: the modified NETAPPROX provides an unconditional upper bound on stability (referred to henceforth as "net upper bound"), and the modified PARITIINNANDAPROX provides an unconditional lower bound ("LP lower bound"). As a result, we are able to experimentally verify that algorithms provide tight (and unconditional) bounds on stability for variety of datasets.
As a baseline upper bound, we implement the greedy heuristic of Kuschnig et al. (2021) which refines Broderick et al. (2020). We not aware of any prior work on lower bounding stability, so we implement simplification of our full lower bound algorithm as baseline. See Appendix for implementation details, hyperparameter choices, and discussion of error bars.
5.1 SYNTHETIC DATA
Heterogeneous data. We start with simple woo-dimeniional dataset with two disparate sub- populations, where the greedy baseline fails estimate the stability but our algorithms give tight estimates. For parameters k, and a, we generate k independent samples (X,, vi). where X; E R2 has independent coordinates Xa N(-1,0001) and X,2 N(0, 1), and yi Xn1- Then, we generate n k independent samples (X,, v:) where X, and X42 N N(0, 1). and y (0, 1). 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 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 & samples.
Indeed, this is what happens. For n = 1000 and l varying from 10 to 500, we compare our net upper bound and IP lower bound with the baselines. As seen in Figure 1, methods are always better
7

---

Published as a conference paper at ICLR 2023
ape Boadd
cuee Ead
(a) d = 2 and 1000; median (b) d = 3 and 500; one trial of 10 trials for each noise level and for each noise level and algorithm algorithm
200 pertor
Figure 1: Heterogeneous data
Figure 2: Isotropic Gaussian data
than the baselines, and certifiably approximate the stability within small constant factor. In the regime where & small, our upper bound outperforms the greedy upper bound by factor of 30.
the previous example, removing k samples caused pathological change in the sample covariance; it became singular. However, even modest, constant-factu instability in the sample covariance can cause the greedy algorithm to fail; see Appendix 1.5 for details.
Covariance shift.
Isotropic Gaussian data. Instability can arise even homogeneous data, as result of low signal- certify stability? For
to-noise ratio (Broderick et al., 2020). But when the noise level is low, broad range of noise levels, we experimentally show that this is the case. Specifically, for d {2, and noise parameter ranging from 0.1 1 to 10, generate independent samples (X,, y4 ) ) where X; ~ (0, Ia) and y/ = (X, 1) + N(0, o²). For 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 2b), the bounds looser for small noise levels but still always within a small constant factor.
and n = 1000 (Figure 2a), our LP lower
5.2 BOSTON HOUSING DATASET
The Boston Housing dataset (Harrison
Rubinfeld, 1978; Gilley et al., 1996) 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-an we designate as the response. Unfortunately the entire set of features is too large for our algorithms, for our experiments we pick various subsets of two or three features to use as covariates.
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) 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 modest positive effect: the estimated coefficient of 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 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 striking heterogeneity in the data: 73% of the samples have zn 0, and the remaining 27% of the samples (precisely those removed by the net algorithm) have crim < 0.83, i.c. 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: this case, it's between the city proper and the suburbs, and in fact the OLS regressors of these subpopulations all 13 features are markedly different
8

---

Published as a conference paper at ICLR 2023
(a) Net upper bound (y) vs greedy (b) LP lower bound (3) vs net up- upper bound (x:
per bound (x)
(c) crim (ÿ) zn (2)
Results from Boston Housing dataset. Figure (a) plots the net upper bounds on the y-axis against the greedy upper bounds on the x-axis; Figure (b) plots the LP lower bounds on the y-axis against the net upper bounds on the r-axis. In both (a) and (b), each mark corresponds to one of the 156 feature pairs. Figure (c) plots the feature zn against the feature crim (on log scale); each mark is one of the 506 datapoints.
Figure 3:
(Appendix L6).. Thus, it's important have algorithms with provable guarantees for detecting when heterogeneity causes (or doesn't 't cause) unstable conclusions.
All-feature-pairs analysis. To be thorough, we also apply our algorithms to 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 the response. We also compare to the greedy upper bound. See Figure for a perspective on the results. In each figure, each point corrssoonds 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 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 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 factor of two some are sensitive to removing less than 10 samples, and some are stable to removing even majority of the samples.
CONCLUSIONS
In this work, we studied efficient estimation of the stability OLS regressions to removing subsets of the training data. We showed that low dimensions the problem is both theoretically and experimentally tractable, whereas in high dimensions exact computation of the stability likely requires exponential time. However, this is only the beginning of the story. Most immediately, since lower bound algorithm takes time n R())) our experiments were limited no more than three dimensions. Certifying stability of OLS regressions from e.g. recent econometric studies may require additional heuristics or insights (c.g. developing fixed- perd-pareetter tractable lower bound algorithm). Beyond that, identifying reasonable assumptions under which exponential dependence on dimension can be entirely circumvented another valuable direction for future work.
Of course, machine learning extends far beyond linear regression, and for more and more complex and opaque models, stability certification is all the more crucial as tool for enhancing trustworthiness. Certainly, OLS is important in its own right, but inasmuch as it is key building block in more complex machine learning systems (from regression trees (Loh, 2011) to generative adversarial networks (Mao al., 2017) and policy iteration in linear MDPs (Lagoudakis & Parr, 2003)), our work on estimating stability of OLS is also a first step towards estimating stability for these systems.
Finally, remark that care must be taken when interpreting stability in practice. Large stability may increase trust in a model's parameters or predictions, but it does not mean that conclusions drawn from the model are "correct." Conversely, even if the stability is small, the conclusions may still be useful, with the caveat that they may be driven by a small sub- uub-population. Understanding whether this heterogeneity is problematic not is context-dependent, and is separate but important issue.

---

Published as a conference paper at ICLR 2023
REFERENCES
Joshua D Angrist, Guido W Imbens, and Donald B Rubin. Identification of causal effects using instrumental variables. Journal of the American statistical Association, 91(434):444-4555 1996.
Ainesh Bakshi and Adarsh Prasad. Robust linear regression: Optimal rates in polynomial time. In Proceedings of the 53rd Annual ACM SIGACT Symposium on Theory of Computing, pp. 102-115, 2021.
David A Belsley, Edwin Kuh, and Roy E Welsch. Regression diagnostics: Identifying inftuential data and sources of collinearity. John Wiley Sons, 1980.
Javier Bianchi and Saki Bigio. Banks, liquidity management, and monetary policy. Econometrica, 90(1):391-45 2022.
Diogo GC Britto, Paolo Pinotti, and Breno Sampaio. The effect of job loss and unemployment insurance on crime in brazil. Econometrica, 90(4): 2022.
Tamara Broderick, Ryan Giordano, and Rachael Meager. An automatic finite-sample robustness metric: Can dropping a little data change conclusions. arXiv preprint (rxXiv:2011.1499) pp. 16, 2020.
Marc-Etienne Brunet, Colleen Alkalay-Houlihar Ashton Anderson, and Richard Zemel. Understand- ing the origins of bias in word embeddings. In International conference on machine learning, pp. 803-811 1.1. PMLR, 2019.
Tony Cai and Lie Wang. Orthogonal matching pursuit for sparse signal recovery with noise. IEEE Transactions Information theory, 57(7):4680 4688, 2011.
David Card. Estimating the return schooling: Progress on some persistent econometric problems. Econometrica, 69(5): 12)-11127,1,,,, 2001.
Enrique Castillo, Ali S Hadi, Antonio Conejo, and Alfonso Fernández-Canteli. A general method for local sensitivity analysis with application to regression models and other optimization problems. Technometrics, 46(4):430-444, 2004.
Maxime Cauchois, Suyash Gupta, Alnur Ali, and John Duchi. Robust validation: Confident predictions even when distributions shift. arXiv preprint arXiv:2008.0422 2020.
Michal Černý, Jaromír Antoch, and Milan Hladík. On the possibilistic approach to linear regression models involving uncertain, indeterminate or interval data. Information Sciences, 244:26-47, 2013.
Samprit Chatterjee and Ali S Hadi. Influential observations, high leverage points, and outliers in linear regression. Statistical science, pp. 379-393, 1986.
Irene Chen, Fredrik D Johansson, and David Sontag. Why is my classifier discriminatory? Advances in neural information processing systems, 31, 2018.
Dennis Cook. Detection of influential observation in linear regression. Technometrics, 19(1):15-18, 1977.
Iliss Diakonikolas, Gautam Kamath, Daniel Kane, Jerry Li, Jacob Steinhardt, and Alistair Stewart. Sever: A robust meta-algorithm for stochastic optimization. In International Conference on Machine Learning, pp. 1596-1606... PMLR, 2019.
John Duchi and Hongseok Namkoong. Learning models with uniform performance via distributionally 2018
robust optimization. arXiv preprint
Panos Giannopoulos, Christian Knauer, and Günter Rote. The parameterized complexity of some geometric problems in unbounded dimension. In International Workshop on Parameterized and Exact Computation, pp. 198-209. Springer, 2009.
Otis W Gilley, R Kelley Pace, et al. On the harrison and rubinfeld data. Journal of Environmental Economics and Management, 31(3):403-40) 1996.
10

---

Published as a conference paper at ICLR 2023
Ryan Giordano, William Stephenson, Runjing Liu, Michael Jordan, and Tamara Broderick. swiss army infinitesimal jackknife. In The 22nd Interaational Conference on Artificial Intelligence and Statistics, pp. 1139-1147, PMLR, 2019.
Suyash Gupta and Dominik Rothenhäusler. The r-value evaluating stability with respect to distribu- tional shifts. arXiv preprint arXiv:2105.03065 2021.
Gurobi Optimization, LLC. Gurobi Optimizer Reference Manual, 2022. URL https: : /www. gurobi com.
Ali Hadi and Jeffrey S Simonoff. Procedures for the identification of multiple outliers in linear models. Journal of the American statistical association, 88(424): 1993.
David Harrison Jr and Daniel L Rubinfeld. Hedonic housing prices and the demand for clean air. Journal of environmental economics and management, 5(1):81-102, 1978.
Jamie Hayes and Olga Ohrimenko. Contamination attacks and mitigation in multi multi-party party machine learning. Advances in neural information processing systems, 31, 2018.
Hugo Hopenhayn, Julian Neira, and Rish Singhania. From population growth to firm demographics: Implications for concentration, entrepreneurship and the labor share. Econometrica, 90(4):1879 1914, 2022.
Peter Huber. Robust statistics, volume 523, John Wiley & Sons, 2004.
Andrew Ilyas, Sung Min Park, Logan Engstrom, Guillaume Leclerc, and Aleksander Madry. Datamod- els: Understanding predictions with data and data with predictions. International Conference on Machine Learning, pp. 9525-9587. PMLR, 2022.
Sookyo Jeong and Hongseok Namkoong. Robust causal inference under covariate shift via worst-case subpopulation treatment effects. In Conference on Learning Theory, pp. 2079-2084. PMLR, 2020.
Michael J Kearns and Umesh Vazirani. An introduction computational learning theory. MIT press, 1994.
Adam Klivans, Pravesh K Kothari, and Raghu Meka. Efficient algorithms for outlier-robust regression. In Conference On Learning Theory, pp. 1420-1430. PMLR, 2018
Pang Wei Koh and Percy Liang. Understanding black-box predictions via influence functions. In International conference on machine learning, pp. 1885-1894. PMLR, 2017.
Pang Wei W Koh, Kai-Siang Ang, Hubert Teo, and Percy Liang. On the accuracy of influence functions for measuring group effects. Advances in neural information processing systems, 32, 2019.
Nikolas Kuschnig, Gregor Zens, and Jes Cuaresma. Hidden in plain sight: Influential sets in linear models. 2021.
Michail G Lagoudakis and Ronald Parr. Least-squares policy iteration. The Journal of Machine Learning Research, 4:1107-1144 1149, 2003.
Edward E Leamer. Global sensitivity results for generalized least squares estimates. Journal of the American Statistical Association, 79(388):867-888 1984.
Wei-Yin Loh. Classification and regression trees. Wiley interdisciplinary reviews: data mining and knowledge discovery. 1(1):14-23, 2011.
Xudong Mao, Qing Li, Haoran Xie, Raymond YK Lau, Zhen Wang, and Stephen Paul Smolley. Least squares generative adversarial networks. In Proceedings of the IEEE international conference on computer vision, pp. 2794-2802, 2017.
John Milnor. On the betti numbers of real varieties. Proceedings of the American Mathematical Society, 15(2):255220) 1964.
11

---

Published as a conference paper at ICLR 2023
Wolfgang Polasek. Regression diagnostics for general linear regression models. Journal of the American Statistical Association, 79(386):336-33 340, 1984.
James Renegar. On the computational complexity and geometry of the first-order theory of the reals.
part Introduction. preliminaries. the geometry of semi-algebraic sets. the decision problem for the existential theory of the reals. Journal of symbolic computation, 13(3):255-299, 1992.
John D Sargan. The estimation of economic relationships using instrumental variables. Econometrica: Journal of the Econometric Sociery, pp. 393-415, 1958.
Aman Sinha, Hongseok Namkoong, Riccardo Volpi, and John Duchi. Certifying some distributional robustness with principled adversarial training. arXiv preprint vxXt::7710.10571.2000
Hideo Tanaka, Isao Hayashi, and Junzo Watada. Possibilistic linear regression analysis for fuzzy data.
European Journal of Operational Research, 40(3):389-396 1989.
Roman Vershynin. High- dighhddimessooal probability: An introduction with applications in data science, volume 47, Cambridge university press, 2018.
Jianlong Zhou, Zhidong Li, Huaiwen Hu, Kun Yu, Fang Chen, Zelin Li, and Yang Wang. Effects of influence user trust in predictive decision making. In Extended Abstracts of the 2019 CHI Conference on Human Factors in Computing Systems, pp. 1-6, 2019.
12

---

Published as a conference paper at ICLR 2023
A FURTHER RELATED WORK
Local and global sensitivity metrics. Post-hoc evaluation of the sensitivity of statistical inference to various types of model misspecification has long been recognized as an important research direction. Within this area, there is distinction between local sensitivity metrics, which measure the sensitivity of the inference to infinitesimal misspecifications of the assumed model Mo (c.g. Polasek (1984); Castillo et al. (2004); Belsley et al. (1980)), and global sensitivity metries, which measure the set of possible inferences as the model ranges in some fixed set M around Mo (e.g. Leamer (1984); Tanaka et al. (1989); Cerny et al. (2013)). For OLS in particular, there is a well-established literature on the influences of individual data points (Cook, 1977: Chatterjee Hadi, 1986), which falls under local sensitivity analysis, since deleting single data point is an infinitesimal perturbation to a dataset of size n as n + oc. 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 Belsley et al. (1980)). Existing approaches typically focus on identifying outliers in a generic sense rather than with respect to specific inference (Hadi & Simonoff, 1993), study computationally tractable variations of deletion (e.g. constant- reweighting (Leamer, 1984)).
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. Huber (2004)). While a valuable branch of research, we view robust statistics 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-trimmot datasets are not necessarily finite-sample robust (Broderick al., 2020). Rather, sample stability (along with the s-value (Gupta & Rothenhäusler, 2021)), 1)), in the regime where a constant fraction of samples removed, may be thought of measure of resilience to heterogeneity and distribution shift.
Second, 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 compu- tationally efficient estimators, for problems such as linear regression, that are robust in the strong contamination model (e.g. Klivans et al. (2018): Diakonikolas et al. (2019); Bakshi & Prasad (2021)). However, even positing that the strong contamination model is correct, estimation guarantees for these algorithms require strong, unverifiable (and unavoidable (Klivans et al., 2018) assumptions about the uncorrupted data, such as bypercontractivityy 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.
Distributionally robust optimization. A recent line of work in machine learning (Sinha et al., 2017; Duchi & Namkoong, 2018; Cauchois et al., 2020; Jeong & Namkoong, 2020) 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.
B PROOF OF THEOREM 1.2
In this section, we show how exactly compute the stability of d-dimensional dataset in time n O(dd), proving Theorem 1.2. Our main tool is Thcorem 4.1, a special case of an important result due to Renegar (1992) solving quantified polynomial systems of inequalities.
The expression Stability y) < k can indeed be written as a polynomial system of (degree 2) equations, with only 3 quantifier. Unfortunately, the number of variables in this naive formulation is n + d (n for the weights and d - 1 for the regressor), which yields an algorithm exponential in n. Thus, take advantage of the above theorem, need reformulate the expression with fewer
13

---

Published as a conference paper at ICLR 2023
variables. The following lemma rewrites the stability, via the separation theorem for convex sets, in form where the variable reduction will become apparent.
Lemma B.I. For any (X,, y and O, it holds that Stability(X..)) y) - if and only if
YE E Rd-l Vu E IR Ew e [0, ||ll 2 n - k A
Wil (X - v-),,,,,,>> u) 2 0, (4)
where X (d 1) is the matrix with columns (X*) )2, (x*)
Proof. From formulation (1) of the stability, we know that Stability (X, y) - k if and only if
YE E R Ew E [0, 1]" [|oll 2 KA X (w (XA y) O.
Fix A E pd-1 Define the set
D(n k) = { (w (XA - -)) w [0, A | 2
We are interested in the predicate D(n - k)) (anneer(X)) (²) * 0, or equivalently 0 E D(n k.) + (XT) Observe that (n - k) is convex, since w ranges over a convex set. Thus, by the separation theorem for a point and convex set, 0 E D(n k) + ker(X")) ) if and only if for every (MP=ED(n=KK)+ee((NN (v, x) 2 0. If v is not orthogonal to ker(X""), then the inner product can be made arbitrarily large. Thus, it suffices to restrict to span( X,), in which case the supremum is simply over x E D(n - k). That is, 0 E D(n - k) n ker(X") if and only if
R", we have
Vu R Bw E [0, 1]"
|wl] 2 n k
(uu (w (XA w))) 2 O.
Quantifying over X, we get the claimed expression.
The expression in Lemma B.1 still has (n) variables. However, we can now actually eliminate the variable at the cost of increasing the number of equations. This is because the optimal for fixed À and u only depends the relative order of the terms ((X,,) A) (X; u). We make the following definition:
Definition B.2. For any À E 1R and R.d let xt x(\, u) be the unique permutation on [n] such that for all < I 1,
((x,,, vn Vr u) 2 ((X ((Tm.1))) A) VT (X u).
and such that equality implies T. < M++1. Let 11 (m(A,
Rdd-1 E 12.).
Then it can be that for fixed À and u, the optimal choice of has coefficients on (A, u).,..., (A,, (n-j. and coefficient n - [n k] for ( (n= if there is any feasible w which makes the sum non-negative, then this choice of w makes the sum non-negative as well. Denoting this vector by w( u)), we have that in Equation 4 it suffices to restrict to E (w(=): π E II].
priori, the number of achievable permutations could be n!, in which case we would not have gained anything. However, because (A, u) is defined by low-degree polynomials in only 2d - 1 variables. we can actually show that [1I] is at most exponential in d, using the following result:
Theorem B.3 (Sign Partitions (Milnor, 1964; Renegar, 1992)). ge. 9m I + be arbitrary
polynomials each with total degree most d. Let SG(g) be the set of vectors E (-1, 0, 1)mh such that a is an achievable sign vector, i.e. there exists some x E R with sign ) = σ; for all i E [m]. Then (SG(g) & (md)0(n) (nn) Moreover, SG(gg can be enumerated in time (md)
Putting everything together, we have the following theorem, which proves Theorem 1.2. Theorem B.4. For any permutation on [n],, define w w(T)) E [0, 1]"h by
1 n k In
i - k ni=n-kk n +1
w( w(T).,
0 k]
otherwise
14

---

Published as a conference paper at ICLR 2023
Then for any ls [0, nl, it holds that Stability(X,y) > k if and onty if
E Bu R VT II: Σ w(T)) -w)(X,, < 0.
(5)
Moreover, II can be enumerated in time n n(d), Thus, the expression Stability (X,v) > k can be decided time n ((d))
Proof. Fix X E Fd-1 and E R If
3x e 11: Σ (((T) - 20.
(6)
then because ||{[()]]ı 2 n k; we obviously get
Ei E [0, 1)"h ||l|1 2 n k A
w wi(( (X,, 20. 0.
(7)
Conversely, if (6) is false, then in particular ((I(A, u)) produces negative Ei w((T(, ())x(((,,,, (X,, A) (X,, u). But by construction, w(((),, u)) maximizes this sum, over all w E [0" with ||||| = n k. Therefore no weight vector with & norm exactly produces nonnegative sum, and increasing the norm cannot help. Thus, (7) and (6) are equivalent. Quantifying over À and u, we have Stability v) I k if and only if
YE E gd-1 Vu E R" 3x E 11: Σ ((T)) (X u) 2 0.
Taking the negation yields expression (5). If we compute II, then this expression is -B-system of polynomial inequalities with 2d variables and II degree-2 inequalities, so by Theorem 4.1 can be decided in time [TT))())) It remains to show that II can be enumerated in time n ((d) (which bounds III..
For any i.j E [n] with define the polynomial
(A, u) = ((X,, A) (x) ((X,, (,) (X,,)) u).
For any À and u, the permutation TT (A,, u) is determined by the signs of the polynomials (fijhee at (A, u). But by Theorem B.3, the set of sign vectors can be computed in time n o(4), So II can be found in time (dd) as well.
C PROOF OF THEOREM 1.3
In this section, we prove Theorem 1.3. That is, we show that exactly computing the stability requires (d)) time under the Exponential Time Hypothesis, by simple reduction from the Maxinum Feasible Subsystem problem. This latter problem is already known to take n ((d) time in d dimensions under the Exponential Time Hypothesis:
Theorem C.1 (Theorem 13 in Giannopoulos al. (2009)). Suppose that there is an nlld) time algorithm for the following problem: given n vectors U1...., Un E Rdd real numbers Ci...., Co. E R. and integer 0 < k < n, determine whether
XERO
Then the Exponential Time Hypothesis is false.
From this, can be easily seen that (exactly) computing stability also requires n R(d) time.
15

---

Published as a conference paper at ICLR 2023
Theorem C.2. Suppose that there is no(d) algorithm for the following problem: given K1,...,Xp, E R and yi.....9 n E IR, as well as an integer 0 S k S n, determine whether
Stability(XX,2 S k.
Then the Exponential Time Hypothesis is false.
Proof. We reduce to Maximum Feasible Subsystem. Given w1,... Un E I and Ct, Cn E R, define X; (ct, Us) IR++1 and y/i G.. Then the regressor ei (1,0, ..... 0) Rd+1 perfectly fits the data set, i.c.
Σ ((X,,e)) - v1))² 0.
Thus, for any w E 1]",
OLS(X, w) B Rd+1 Σ wi ( 11)² =
Suppose that Stability (X, y) S n Then there is some w E [0, 11 and 83 = 0 with ||l| 2
and Ei=1 (X, 3) - y/ 8 = 0. Let 5 S [n] be the support of w; then [S] 2 k. Moreover for
every 2 E S, it holds that (X, 3) V/ O. Since 31 = 0, by definition of X; and ys, this implies
that (vi, (2.41) = Gi. Thus,
Σ 1k,, B2.d41) c,l k.
Conversely, suppose that there exists some A E Rd and set S C [m] of size k such that (ve, A) cs for all i E S. Define = Is E [0, , and define B = (0, A). Then it is clear that
wi (.X(( - vi)² Σ ((us. -) c2)² 0.
iE
Thus, 3 E OLS(X,XX)))) y, w), so Stability (X, X,y) S n k. This completes the reduction.
Remark C.3. Maximum Feasible Subsystem does have an en-additive approximation algorithm in time ((d/e)))) by subsampling. Thus, proving (0dd) for en-approximation of stability would require a different technique.
D PROOF OF THEOREM 1.4
In this section, we prove Theorem 1.4. The main idea of PARTTIONANDAPPROX is that under Assumption A, we can approximate the stability by partitioning d into roughly nd regions and solving linear program on each region. See Appendix J for the complete algorithm.
D.1 PARTITIONING SCHEME
Given samples (:), let M = X (() v112. where ((0) E OLS(X,, y, 1). Let [n]
uniformly random size-m subset of [n]. Let Ri...., . RG be the closed connected subsets of [gd-1 cut by the following set of equations E:
(X,, A) - yi =
Vi [>]
(X,,A) Y: oM
Vi [n], Vo (-1,1)
(K, - y: o8M/VV
Vi E [n] E (-1,1)
(X, - 9: = =
ViE E (-1,1), < < [lg (1884/(Vn/0)) vn/o)]
Formally, define region for every feasible assignment of equations to {=, <, and then replace each strict inequality by non-strict inequality, so that the region is closed. First, we observe that the number of regions is not too large, and in fact enumerate the regions efficiently.
16

---

Published as a conference paper at ICLR 2023
Lemma D.1. Each region R, is the intersection of (((E)) linear equalities inequalities. Moreover, the regions Ri,.... R.....,,Ppp... Ry be enumerated in time ((()))++0()
Proof. Order the set of equations & arbitrarily. We recursively construct the set of regions demarcated by the first equations. For each such region. we solve linear program to check if the + 1th hyperplane intersects the interior of the region. If so, we split the region according to the sign of the t + 1th hyperplane. The overall time complexity of this procedure is Opp poly(n, [E)),, where p is the final number of regions. But the number of regions which can be cut out by hyperplanes in IR most (t) by standard arguments.
Next, we argue that even though only multiplicatively partitioned small subset of the residuals, with high probability most residuals are well-approximated. More precisely, we show that for the (random) region R containing a fixed point A with high probability, for every other À in the region, most residuals at À are multiplicatively close to the corresponding residuals at X This can be proven by the generalization bound for function classes with low VC dimension:
Theorem D.2 (Theorem 3.4 in Kearns & Vazirani (1994)). Let X be set. Let C 10, be binary concept class with VC dimension d. Let 'D be distribution on X x 10, Let e, > 0 and pick m N. satisfying
1 100
!
m m Ci
for sufficiently large constant Cj. Let (is,Ei) be independent samples from D. Then with probability at least - 8, the following holds. For all h E C with h(x;) yi for all E [m],
Pr [1x)) # S.. (H
Specifically, for any X E Rd-1, we define function fx [n] {0. 1} the indicator function of samples for which the residual at À is close to the residual at X Then the region containing X is precisely the set of functions which perfectly fit S x {1), and the generalization bound implies that all such functions fit most of [n] x 1, which is what we wanted to show. The following lemma formalizes this argument, with some additional steps to deal with very small and very large residuals.
Lemma D.3. Let X E IRd-1 and let 0. Let R be the region containing X². Let BM [n] be the set of i [n] such that KX,, V/l >> M, and let Bsu S. [n] be the set of E [n]] such that 1X X²) ysl S 8M///n If > ddlog +1og log !
C
1
for an absolute constant C, then with probability at least 1 -1 over the choice of S, the following holds. For every X E R, the number o E (Bas U BBM) such that
() VEl
>C
A), 24l
is at most
Proof. Let E { 13 be defined by =iign((( X²) v.). Note that by construction of the regions, sign( (X,, X) y1) o; for every À E R.
Let S' n ([n) \ (BM BaM)). If := [in] \ (Bas B&Ms) s en, then the lemma statement is trivially true. Otherwise, by Chernoff bound,
n' m
Pr
S
2
2 1 -
2 n/3.
2n
Condition on the event | 2 n'm/(22m))) (22).. Then S' is uniform over size-m' subsets of [m] \ (BM U
B8M).. Define X' [n). and define a concept class C = (l1x X' {0. 1}/A E Rd- as the set of binary functions
(x) = [oo ((xi.)- (0::A--1++e00) ((xx)-(1++))x>) + eyj) <
17

---

Published as a conference paper at ICLR 2023
Let D be the distribution on X x {0, where (i, ~ D has i Unif([n] \, (BM U B&M))
and $ = 1. Then S' x {1} consists of at least n' (2mm/(2n) independent samples from D. For any A R. we claim that the function h, fits S' {1} perfectly. Indeed, for any i E S', we know that 0M/V1 S o,,((,,, \²) yi) S M, since i & Bas U BsM. So by construction of the regions, holds that
1 ((XiA) y:)
S +1
I+ ox((Xxx*) y/)
From the right-hand side of the equation, precisely get that hx(2) 1, as claimed. By Lemma D.4, we have vc(C)= S Thus, since
2 n'm
Cd (2n/n')e
C log 1) (2n/n')e
IS
log
2n
C' is a sufficiently large absolute constant, we can apply Theorem D.2 with failure parameter é' := nee(2n') to get that with probability at least - n/3.
sup Pr (1) + 8) Sé' (t.))-
Equivalently, for every À E R. the number of i E [m] \ (BM U BBM) such that or (, - ys) > (1+ (),,)((,,, - ys) is at most e'n' = =n/2. An identical argument with the binary functions
(x = 1 [o: - (1-c)1) -cy) 2 Of
proves that with probability least 1 1-7/3, for every X E R. the number of E [n]] (BM such - y:) most cn/2. Overall, by the union bound (taking into account the event that [S"| < n (2n)), with probability at least - 17 it holds that for every
that of ((X,, A) y1) (1 (1-c) ().cc))oXX,))))
A E R, the number of i E [m] ( such that either o; (Xx( A) yt) > (1+c)o; (x,( A) v:) (1 ((X, X) yi) is most as claimed.
vi
It remains bound the VC dimension of the function class. Lemma D.4. For any d and any (x,Bi)ke()) C Rd x IR, the VC dimension of the concept
class C (hx [m] + (0, 1}11 E Rdd, where
hx(1) 1[(X,, + 3 S O].
is at most d + 1.
Proof. Note that extending the domain of the concepts cannot decrease the VC dimension. Thus, if we define C' {h's Rd
+ {0. : R by
kx = ][ + + + 0],
then vc(C') 2 ve(C). But all of the concepts in C' are affine halfspaces in d + 1 dimensions, so ve(C') < 11
D.2 ALGORITHM
For every region R we can identify some arbitrary representative (R E R and a sign pattern E 1-1,1}" such that = implies (X,, A) ys 2 0 for all À E R, and a; = -1 implies (X,, A) y: < 0 for all X E R. Ideally, we want to compute V = maxæ VR where
x g 0
VR = sup 9ER" AER
sEim
X)
Vi
0 9i $ (X,, A)
(8)
(X,, 2) S 50 Wi [m]:a==- =
However, this is not linear program. Let BM BM [n] be the set of [n] such
that (o(R) vil > M, and let BaM Bas (R) [n] be the set of [n] such that
18

---

Published as a conference paper at ICLR 2023
AX,, - yal < 8M/vn. For each region R. we compute
((1))
(9)
X" =0
0 9 (X,, yi
Σ
9t (xx,,0(R)) yi :1)
= sup 9ER AER
min
ViE [n] :0 =
(B
(X,, )) y1 9 S 0
Vi E [n] : = -1
(10)
and define
g(R); Xx X(R)) y4
V(R) E
iE
with the convention that 0/0- 1. Note that although there is min in the objective of Program 10, it is equivalent to a linear program by a standard transformation: we can introduce variables s1,..., with the constraints s; S and s; S g;X(((; (t) y:), and change the objective maximize Eie(y) 8t- That is, the following program also computes (g( Rt), (RR)::
(0(R),((R))
(11)
X 9 0 0 9: S (X,, A) - 1/ (X,,)) 3: $9:500
Wi E [m] o; =
E Bi
arg sup gER" aER
Wi [n] :0 =-1
8, 1
Vi E [n] \ (BM U BBM)
8, 9t/ /))((x((;,,))) )o(() yi)
Vi [n] (BM U BBM)
(12)
Lemma D.5. Let (w" ,X)) be an optimal solution to Program 1. Suppose that the event of Lemma D.3
holds (with respect to X*), Then either max Bsas (R) > en,
12en 1 S max V (R) SV.
Proof. For the RHS. observe that for any region R, if we define w E R" by wi = (XX, y:), with the convention that 0/0 = 1. then the second and third constraints of Program 10 ensure that w [0, 1)a, Additionally, the first constraint ensures that
XT Σ mu(((X, vs) 0.
XE(M
Thus, (w, ((R)) is feasible for the original problem. This means that V || S V.
To prove the lower bound, suppose that max Bas R)) S en. Consider the specific region R containing the optimal parameter vector X (if there multiple choose any), and let Bag B R)) and BBM = BaM (It). Define g" E Rn by g = w/ ((X,, y1). Let Bapx C [n] (BM U BBM) be the set of E [m]] (BM U BBM) such that
I+e (LK,, - y:) (X, X² - vi v (xx,, ys > (1 + (XX, Ao - vn)
By Lemma D.3 (applied specifically to À = (R)) we know |Bapa < cn. By assumption, we know
that Baar en. And know that if i By then (o yil > M, so the same holds for
19

---

Published as a conference paper at ICLR 2023
all X E R. and in particular for X*. Thus
Σ
Σ w ((X, A²) v2
S
M
iEBM
1
Σ
W ((X; X²) -2)²
M
iE[n]
W ((x;, 3(0)) - vi)²
M
(E[m]
1
< X(3(0)
M
=1
where the third inequality is because X E OLS(X, y, w*) and the equality by definition of M.
Finally, for any i E [n]] \ (BM U Baas U Bapx), we have
(Xi, Ao())) - yi :1 2 min
Di
1 (xi,1*) - Vii (5.mm11) 11) 1+ (X,,*)
gt
g.
min
y/
Therefore
wi
iE(m)
=
Σ
(X x²)
9i
+
Σ iEBM
(X,, A²) V1
g
Σ
w +
Σ
Vi
UBa
iE Ba
iE Plys
< (1 c))
Σ
min
Si
.1 + 1 + Zen
(Xi, do (R)) Mi
$ <1 c)
Σ min (X,, Ao (R)) - yi 1
gf
+ 1 + 2en.
The sum in the last line above is precisely the objective of Program 10 at ,²). Moreover, (9", is feasible for Program 10 because XT9 XT("" ((X,, )²)
vs)) - 0 and 9; /(X,((X, X*) -vi
wi E [o. 1] for all [n]. Thus, the optimal solution ()(R), X(R)) satisfies the inequality
min
Ao(R)) - Vi 1)
gi
Σ
min (Zx, ((1))
(X,, - yi :1)
ve[n]\(ParUBan)
Finally, let Bupx S [n] be the set of i E [M] \ (BAr U BAM)) such that
c1+ ((),, - vi) > (x,,)(R))) yı v (xx,, - y; > (1 +
- wn]
By Lemma D.3, the residuals AR) multiplicatively approximate the residuals X except for en samples, and the residuals at X(R) also multiplicatively approximate the residuals at A except for en samples. Thus, [apx| I 2en, so we have
O(R)
9(R)
Σ min
- yi 1)=
(BN (o,,Ao(11)) y/i 1)
Σ
min
(X,, (I))
((12)
Σ
min
(o,,)o(R)) - yi
BEB
g(R)
Σ
min
.1 + 2en
(X,, do(RR) y/-
(BM Bap
O(R)
< <1++c)
Σ
+ 2en
(X,, ((R)) yi
(Br
( + 2cn.
20

---

Published as a conference paper at ICLR 2023
We conclude that V Σ (1 + ()) (1 c)2en + 2en. Since ((R) n, simplifying gives V S (R) 12mn 1 as claimed.
Using the above lemma in conjunction with Lemma D.3 immediately gives the desired theorem (from which Theorem 1.4 is a direct corollary).
Theorem D.6. For any c, 6, O, there an algorithm PARTITIONANDAPPRO, with time com- plexity
Cd 22
1 log
d+0(1)
n +
log
en
which, given 6, 8, T), and arbitrary samples (Xi, Vili=1 either oulputs 1 or estimate S. If the output is L., then the samples do not satisfy (e, (----a-atnconee ca-anti-concenraatinn (Assumption A), Moreover, the probability that the output is some such that
$ S-SttbilityyX, v)l 12en +
at most T).
Proof. The algorithm PARTITIONANDAPPRO does the following (see Appendix J for pseudocode).
+ log where is the constant specified in Lemma D.3, and let M = || - yl12 where ((0) E OLS(X, y, 1). Let & the set of equations described in Section D.I, with respect to a uniformly random subset S C [n]] of size m. Let R Rp be the closed connected regions cut out by E. By Lemma D.I, we enumerate , Rp in time 0(1)))++01) each is described by most |3] linear constraints. For each R, representative (o() E R by solving feasibility LP on R, and by solving vi LPs on Rt we can find a sign pattern a E {{-1, 1 such that implies (X,, A) 3--1 - 0 for all À E R. We also compute Ba = {i E [m] 1Xx, Ao (R)) yel > M} and B&M (i E [m] 1(x,, (R)) va < 6M//u). If [Barr] > en. then return 1. Otherwise, a(R).
Let m = Ce
(dloge Ce-I(dooge
find
1 implies (X,,, A} 8s 2 0 for all X E R. and a; -1
compute (R) (()), the solution to Program 12, and compute ((R) E (8,M()))
Finally, after iterating through all regions, output max ((R.).
Correctness follows from Lemmas D.3 and D.5. For each region, the time complexity is poly((E)).
As there are O( regions, the overall time complexity ((E))d++(1) Observing that | O(n + me -1 log(n/)) completes the proof.
E PROOF OF THEOREM 1.5
In the previous section, we approximated the nonlinear program (1) by partitioning RRd-1 into regions where the program could be approximated by linear program. This approach had the disadvantage of requiring that in each region, the signs of the residuals (X;, A) -:: yı were constant (so that the program could be reparametrized to have linear constraints), which necessitates making 2(n²) regions. In this section. we instead make use of the fact that for fixed X, Program therefore efficiently solvable. Our algorithm NETAPPROX is simply to (carefully) choose a finite subset. N Rd-11 solve the linear program for each À E A,, and pick the best answer.
linear program and
d-- I in an appropriate The following lemma describes how to compute N, which will be an c-net over dd metric.
Lemma E.1. For any (X,, ) and 7 O, there is set N S Id- size I S (2Vd/7)² such that for any TRd-1 with X X # y, there is some N N with X # y, and some E (11,1). such that
X y
XX -
57.
|| | 2||,
XN
Moreover, N can be computed in time O O(Vd/7)a
21

---

Observe that Published as a conference paper at ICLR 2023
Proof. Let A n x d be the matrix with columns (X"), 12 ........ (X')s, y. Let d' rank(A) and let A = UDVV be the singular value decomposition of A, where D diag($1 7-net over the unit sphere
) is the diagonal
matrix of nonzero singular values of A. Let M be a "marginally-randon mean that M is chosen from some
S under the & metric, where by "marginally-randon distribution, and every point of the has marginal distribution uniform over S (e.g. take any fixed nenet and apply uniformly random rotation). Also define B = VD-1 Suppose Ba is not identically zero. Then define
N {
(Bm))dd1 (Bm)a
m E M
Since each m E .M a generic unit vector (by the marginally-random property) and Bą is nonzero, we have that all (Bm)a are nonzero with probability 1, so the above set is well-defined. Moreover, for any m E M,
(Bm)i:d1 (Bm)a
Bm (Bm)d
Um (Bm)a
X
y = A
Now let À Rdd1 with XAA * y, and define
DV (V 11/ | 1xxt-0/ 8/1, ER..
||l|| 1|U1|12~ = ||A[}; 1|llz =
X 21/2
So there is some m E M with ||u - mll2 S 7. Let X (Bm) E N and let sign((Bm)a). Then
X - Uv
XA y 8
and
XX'
Um/(Bm)a sign((Bm)aa)
Um oUm.
XN - y
(5mm/Bmaallz
|Uml]|
Thus,
XA XA :/
XX
= =1U0 - Umml2 = =110 -mll2 S
X
//
desired. On the other hand, if Ba is identically zero, then V22 0. so y
(AA )a 0. This
boundary case can be avoided by picking any nonzero covariate (X among (X (X")),..... )2..... ... (X")a
and replacing by y + c(X"); for generic E R; this does not change Stability (X.y)
Note that 7-net M of gd' with |M & (2vd/7)4 d be constructed in time O((v )
by discretizing the cube [-1, Thus, A can be constructed in time (((V//7)d) as well, and [W] S [M] s (2Vl/7)dd
All that is left is to show that under strong anti- coffenntnnttainn, the value of the linear program is Lipschitz (under the metric described the previous lemma). This proves Theorem 1.5.
Theorem E.2. Let 8 O. There an algorithm NETAPPROX with time complexity (2Vd/(<8>)) poly(n) which, given arhitrary samples (X,, y:) produces an estimate S satisfy- ing Stability(X.g) v) < S. Moreover, if (X,, yi satisfy (e.8o- anti-concentration, then in fact
Stability(X, < & < Stability(X. v) 3em +
Proof. For any À E Rd-1, define
v(A)=
sup we 0,1
{|lw/h x (w+ *xA- -)) =0].
22

---

Published as a conference paper at ICLR 2023
so that n-Stt--lityy y) SuPXER--1 V(A). For any fixed X. can compute ((A) in polynomial time, since it is defined an LPP with n variables and d constraints. Fix e5.. The algorithm NETAPPROX does the following (see Appendix J for pseudocode): if then output O. Otherwise, let V be the net guaranteed by Lemma E.1. Then compute ((A) for every À N, and output the estimate
y) 0 has solution,
S 5=
max ( XEN
Since the algorithm involves solving (2vd/7)² linear programs, the time complexity is (2vd/7)) poly(n) (n) as claimed. Next, since the algorithm is maximizing V()) over subset of Rd-1, it's clear that $ Stability (X,y). y). It remains to prove the upper bound on S. Recall from Lemma B. that for any À E Rd-1
V(A)
inf sup wERd wE(0.1)"
Σ wo Ewi((X,,)) vx) u) 20]
Let X be a maximizer of (()) and for notational convenience let h: - V(A)). If X - y. then the algorithm correctly outputs 0. Otherwise, by the guarantee of Lemma E. there is some À E . with XX * y and (-1,1} such that
X y
XA y
57.
X 3/12
X 311
Pick any u E Rd. Since v(A*) = ki there some w = w² (ou) E [0, 1J"h such that ||w*],, 2 k and
a
Eut((E,, - (X,, ou) 2 O. Σ
Without loss of generality, there is at most one coordinate i E [n]] such that wy is strictly between 0 and 1. Also, the above inequality implies that
w
X
A)
V/
w/
(X A) y/
(X,, )²) V/1
(Xi,
X
2
8/
Xx - 3/
X y
Xx y
2
0
1|
X1 y 3/2
X
2
2 - Xuull2
where the last two inequalities by Mauchy-Schwand and the guarantee of Lemma F. respectively. Now define w E [o, 11th by the following procedure. Initially set w == w". Iterate through [n] in increasing order of ((X,,), (x,, u) and repeatedly set the current coordinate w; := 0, until
Σ ((,, y:) u) 20 0.
Obviously, this procedure will terminate with a feasible w. It terminates after making updates, then || 2 [|w"], I - Throughout the procedure, the sum Σn= wa w,((X,,)) A) - (X,,,, u) is non- decreasing. If [ A) yı 2 (8/
A y/ and (X,, u)) 2 (8/v²) ||Xu|l, and w = 1,
x y/ ||uullz So the number of such steps is at most yn/82 < But the number of steps with K(,,, A) vil (5/vm) X² / is most by (e, -)-ooo anti-concentration. Similarly, the number of steps with (X,, w)} < (8/v ||Xull is at most en. And the number of steps with ( < w < 1 is at most 1. Thus, ||l|| 2 ||0") || 33n 1.
then the sum increases by at least (52 (62/n)
We conclude that V(2)= V(X*) 3en 1, $ Stability (X, 3
23

---

Published as a conference paper at ICLR 2023
F MOTIVATION FOR ANTI-CONCENTRATION ASSUMPTIONS
F.1 SMOOTHING IMPLIES ANTICONCCENTRATION
In this section, we show that anti-concentration (Assumption A) holds under the mild assumption that the response variable is smoothed. In the following proposition, notice that we always have the crude bound (x3( r112 S ||r|l2-
Proposition F.1. Let o > and let (x,,T:);=1 be arbitrary with , E Rd and ri.... I' E R. Suppose that Z...... Zn (0, o²) are independent Gaussian random variables, and define Yi = ri + Z: for i E [n]. Then with probability at least n 2e- the dataset (X,, ye) satisfies (2d/n, (ts//n,))-nnnccuoe cattccccnnnnnnioioonnn where
1
8 min
nvd X300 - -r112 3n /2d
and g(0) OLS(X, r, 1).
Proof. Without loss of generality, (00) XX Define 3 (X"x) ys we want to
upper bound X y/ "/2 Let Q x(X"X)'xX;; then y X3 N - Q),, Q)). But (I - Q)r = - X3(10, and o(0 Q)) is stochastically dominated by o(0,,21). So y/ - X3 $ | x ll + 20 / with probability at least 2e- " by the tail bound for x²
random variables.
Next, we claim that with high probability, for every set S S [n]] of size |S| 2d, there is no 3 R such that X, - vil a/(n 22d for all i E [S]. Fix some S, We bound the probability that for all Rd.
d² no
Σ (X, y.)² IES
because this event contains the event that residuals in S are at most a (/(n/2d) Now suffices to restrict to the OLS estimator POLS have
(X3 x ( ys. Now defining P = Xs(X3 x XX. we
Us -XsBOLs (rs Zs) P(rs Zs)
= =(1 P) + Zs)
~ [[(I-P)rs. P
Note that / P an orthogonal projection onto a space of dimension d' = |S| rank(Xs) 2 d, so there is matrix M I22dxd' be such that MM = I P and M M Ig. Then we can write
US XsBOLs oMA PIrs)
= MTrr)
for a random vector A (O, Tw). This means that ys XsSOLS 2 ||l1 M rs| Mrss But for any l E IR, note that |g + ul stochastically dominates Is) where E ~ NN 1), so (E + x)² stochastically dominates g?, and thus || | + rsssl stochastically dominates A| Therefore
a2 <
Pr us - XSBOLS
Pr
AP +
A3 n
For any S i d', since the density of A, is bounded above by 1/V2TT, we have
V2 d'
Pr[Af +
A3 1/n <
Pr Pr[lA < 1/n
n 3
$
n
-dd 3d'
24

---

Published as a conference paper at ICLR 2023
Recalling that d' 2 d, this shows that for a fixed S _ [n]] of size 2d, with probability at least 1 1- n 333 there is 3 TRd with |x vil S // //(n/ 22) for all i S. A union bound over sets S of size 2d proves the claim.
Finally, we use the above bounds show that the smoothed data satisfies (2d/n, d--antiiconcentraaionn We consider two cases. First, if || |x8(0) r|l2 ovn, then with probability 1 2e- have
X3 u112 S 30 ni, so for any 3 E R.d the number of E [n] satisfying
KX,, - vel s v
x3 2/12
x 1/ S
Vn
3n 2nd
n n
is at most
F.2 DISTRIBUTIONAL ASSUMPTIONS FOR STRONG ANTI-CONCENTRATION
In this section, show that under reasonable distributional assumptions the samples (X,, ya R strong anti concenteenttrnon (Assumption B) holds with constant e, 5 > First, it holds the samples (Xe, yi) are i.i.d. and have arbitrary multivariate Gaussian joint distribution.
Proposition F2. Let n. d, e 0 and set 8 = e/4. Let E d x d be symmetric and cesitive-definite Let 7,...., Zm ~ N(0, Σ) be independent and identically distributed. If n 2 Cd log(n)/e for some absolute constant C, then with probability at least 1 exp(-S((e (()))) samples (Z,)== yd satisfy (c, -storong anti-concentration. iee. for all B E Rd, it holds that
K E [n] (ZZ B3
8 vn
[12311a]
S on
where Z : n x d the matrix with rows
Proof. Let 5 1 2im Z:ZZ' Then by concentration of Wishart matrices, it holds with probability
at least - 2 --2eppp-(1 (n d)/2) that E = 2 (see e.g. Exercise 4.7.3 in Vershynin (2018)). In this event, || /n < 2E (21, 3) for all B Rd.
Fet F= { : (R) be the class of binary functions fa(x) ||| < (2244)/( B), Observe that every function in F is the intersection of parallel half-spaces, so F has VC dimension Odd). Moreover, for any B,
(f)s(z)) (r( (
since (Z, N(0, T Σ,), and Gaussian random variables anti-concentrated (Pr(Ig) < <) S for any e >>0 0 and & ~ (0, 1)). Thus, by the Vapnik-Chervonenkis bound and assumption on n,
Pr sup fEF
(z( > Bfa(Z1) + c/2 < <exp(O(dlogn)) n) ()(2n)) S exp(-s((en))
t=1
So with probability at least - exp( exp(-S((e2)) it holds that for all B E Rdd
{ie [n] |(Z, = (c1/4)== s en.
This means that with probability at least 1 - exp((=((((()) exp( (pp(-(x) d)/2), we have that for all E Rd.
< Ki [n] ( 2311 /// /m) sen
as desired.
Second, strong anti- nnticonncentattion holds if the samples are drawn from a mixture of centered Gaussian distributions, with arbitrary weights, so long as each covariance matrix has bounded largest and smallest eigenvalues.
25

---

Published as a conference paper at ICLR 2023
Proposition F3. Let n. d, c > 0 and let Ei...., Σx be symmetric and positive-definite matrices. Suppose that there are constants X, A > such that AI = Ei...., Σx = AI, and define R. = A/X. For arbitrary weights wi...., Wk 2 0 with + ... + w2 = 1, let Zi...., Zn ~ E,=,;;N(0 (; be independent and identically distributed. If 2 Cdlog(n) for some absolute constant C, then with probability at least exp( (ep(-Q(eEn)) samples (Z,)²=1 satisfy (c, (akk/(())sstoonn anti- concentration.
Proof. The proof is similar to that of the previous proposition. First, note that each Z; can be coupled with N (0, AT) so that Z,Z," & So with probability at least 1 2 exxp( (n
2)) it holds
that 1 - " 3 2X1 and thus || /n S 2A || 1 for all 8 R4.
Next, let F = [fe: 3 E R be the class of binary functions fa(x) = (,, < < 1 /4]. Then F has VC dimension Odd) and for any B,
(fa(Z)) = Pr[] B)1 < 2
S 2
since (2,,3) ~ Σt (O, 3d Σ,3) is stochastically dominated by N [13112). By the Vapnik- Chervonenkis bound, with probability at least 1 exp( ((((n)), we get that for all / E R.
Hi E [n] : 12, B)1² < (c /4) S en.
Combining with the upper bound on | we have with probability at least 1 - exp(=S((e'n) exp(-(n d)/2) that for all 8 E R
Ki E [n] (2, 3) < (=x)(8x)) | //2///n en
as claimed.
G EXTENSION TO IV LINEAR REGRESSION
We extend the definition of stability to measure the stability of the sign of a coefficient of the IV linear regressor:
Definition G.1. For samples (Xi, yi, Z,)_ with covariates X; E Rd, response yi E Rdd and instruments Z; E IRP, the ordinary IV estimator set with weight vector w E [0. is
MVX, y, Z, w) = {B E Rd z (w (XB y)) = -).
The finite saatl- stability of (X,. Vi, Z, then defined
IV-Stability(X, 2) = inf {n - || 31 = 0A E IV(X, y. w.
w .ER"
For any k; the expression IV-Stability(X,,, Z) 2 k is still defined as bilinear system of equations in B and Our exact algorithm for Stability (X, y) never uses that the weighted residual (w (XB - y)) is multiplied by X" in the OLS solution set, rather than some arbitrary matrix ZT, all that matters is that this matrix has most d rows. Thus, with bound on the number of instruments, the algorithm generalizes to computing IV-Stability(X (X.y. Z):
Theorem G.2. There is an O(dp(d+))) -tme algorithum which, given n arbitrary samples (X,, Vi. Z:) with XI...., Xx, E IRd and Zi.... Zn E RR and , E R, and given k 2 0, decides whether IV- IV-Stability(X, Z) < k.
The algorithm NETAPPROX also generalizes to IV regression. To state the guarantee, we define an anti-concentration assumption for IV regression data.
Assumption C. Let ô > 0. We say that samples (X,, Vis Zi) satisfy (e, -)-trog anti- concentration if for every a E Ip' and 3 E Rd it holds that
Hi in] 1l ail < (8/
Za ||)l S en
and
Hi E [m] 1(x,,3) wi)l (g/vm)
83 - // y s en.
26

---

Published as a conference paper at ICLR 2023
Then it is easy to see that the proof of Theorem 1.5 immediately extends to give the following result:
Theorem G.3. For any e, > O. there is (Vd/(c62))² poly(n)-time algorithm which, given e, 6, and samples (Xi, Mi, Z:) satisfying (e, -strong anti-concentration, returns an estimate $ satisfying
IV-Stability (Xtitttyyy,, (, Z) = S IV-Stability Z) + 3en + 1.
H HEURISTIC FOR LOWER BOUNDING STABILITY
In this section we explain the "LP lower bound" which we applied in Section 5 to provide exact lower bounds on stability of various datasets.
Given a list of thresholds T and subset size m, we randomly pick set S S [n]] of size m, and enumerate the regions R1,..., Rp defined by the hyperplanes
(X, yss =t Vi V.
Fix one such region R. Similar to PARTITINAANAPPROX, we use the change of variables
(1 - (xx) X) y1). and enforce the constraint X == y). For some of the i E [n]. the residual (X;, A) y/ will have constant sign the entire region. For these samples, we enforce the constraint 0 S 9: I (X,, X) Yi (if the residual is non-negative) or (X; A) Vi 9: (if the residual is non-positive). However, because we didn't it include a hyperplane for every sample, it's likely that for some [m]. the residual attains both signs within the region, in which case the constraint E [0, 1] is not convex in and A. Thus, we relax the constraint to
inf (X, A) y: S g (X,, A) y4-
XER
XeR
Note that this is indeed relaxation, because the interval contains O. Finally, let K+ C [n]] be the set of indices for which the residual is non-negative on R, and let K_ S [n]] be the set of indices for which the residual is non-positive. Then we minimize the objective
Σ ER QUPAER
+ SEK
Σ
9i
in infxER(Xi,
Because this objective is less than or equal to Eig(n) (=,(((,,)) - -i), and because we only relaxed constraints, this program has value at most VR (the value of the exact non-linear program restricted A R).
Compared to the provable approximation algorithm described previously, this algorithm is a heuristic because each region will typically have samples for which the residual changes sign within the region. As rough intuition, if the residual (X,, A) y/s remains fairly small throughout the region then the relaxation of the constraint on w; may not be problematic. However, if the residual can attain large magnitude in the region, then relaxing the constraint on w; may significantly change the value of the program. However, we may expect that if the partition is sufficiently fine, then region where some residual allowed to blow up may have many samples whose residuals are forced be large. This motivates the use of an additional heuristic to refine the certification algorithm: if (w*,8*) ,9*) are the optimal solution to the sensitivity problem, and g(0) = OLS(XX, y, 1) is the original regressor, then must hold that
Σ ii ((X,,8)) - y:)² Σ W ((x,,
v1) <
Σ ((Xi, ((0) -21)2.
However, if R is the region containing the optimal solution and if ||w"|] 1| 2 k. then
inf
inf
Σ inf ( -):1 XER ieS
SC n
Thus, if we compute (R) in infxeR((X,,)) )1 (by computing the interval of achievable residuals) for each i E [n). then lower bound & (conditioned on R containing the opti- mal solution) by sorting (R),..... (ft) and finding the smallest subset C [m]] such that
EieR (:( S x|| =11²-
27

---

Published as a conference paper at ICLR 2023
To summarize the algorithm, for each region compute the maximum of the LP lower bound and the residual lower bound, and the output of the algorithm is the minimum over all regions. See Algorithm 4 in Appendix J for detailed pseudocode.
I FURTHER EXPERIMENTAL DETAILS
I.1 IMPLEMENTED ALGORITHMS
Net upper bound. We implement NETAPPROX with one modification: instead of deterministically picking the net M by discretization (see Lemma E. we let M be a set of random unit vectors from sd' and then compute N as in Lemma E.1. Instead of parametrizing the algorithm by the desired approximation error c. we parametrize by [.|]. Despite this change, the algorithm still provides provable, exact upper bound on Stability(A)) (X
LP lower bound. The disadvantage of NETAPPROX is that only lower bounds the stability under an assumption that seems hard to check. The MARTITIONANDAPROO algorithm is better, because unconditionally, with high probability, outputs either an accurate estimate or failure symbol However, the (2n")) time complexity (needed so that in each region, all n residuals have constant sign) may be prohibitively slow practice. For this reason, we introduce heuristic simplification of PARTITIONANDAPPER which provably lower bounds the stability with no assumptions.
At a high level, we decrease the number of regions by ignoring the requirement that within each region all residuals have constant sign. The algorithm is parametrized by a list of thresholds L and a subset
size m, and the regions are demarcated by the hyperplanes (X,, A) y/ for random choices i [m] and all t L. Now, for samples which do not have constant residual in a particular region R. the constraints w; [0, 1] are nonconvex after the change of variables. We relax these constraints to linear constraints, and relax the objective function to skip the "bad" samples. Heuristically, this relaxation should not lose much samples where the residual remains small throughout R. but may be problematic if the residual blows up. This motivates the use of a complementary lower bound heuristic based on the minimum squared-loss achievable by any X E R. See Appendix H for details.
Baseline greedy upper bound (Broderick et al., 2020; Kuschnig et al., 2021). We implement the greedy algorithm described by Kuschnig al. (2021), which refines the algorithm of Broderick et al. (2020): iteratively remove the sample with the largest local influence until the sign of the first coefficient of the OLS is reversed. After each step, recompute the influences.
Baseline lower bound. This algorithm simply computes the squared-loss-batt lower bound (used in our full lower bound algorithm) for each region. See Appendix H for details.
1.2 HYPERPARAMETER CHOICES
The net upper bound has only one hyperparameter (the number of trials), which should be chosen as large as possible subject to computational constraints. The LP lower bound has two hyperparameters (the size of the random subsample, and the set of thresholds). We choose these ad hoc subject to computational constraints. Experiments to determine the optimal tradcoff between the subsample size and threshold set could be useful. However, we note that because the LP lower bound unconditionally lower bounds the stability, no matter what hyperparameters we choose, in practice it suffices to try several sets of hyperparameters and compute the maximum of the resulting lower bounds.
Heterogeneous data experiment. For each dataset, we applied the net upper bound with 1000 trials, the LP lower bound with L (-0.01, 0, 0.01) and m 30, and the baseline lower bound with L = -0.01, 0, 0.01} and m 1000.
Isotropic Gaussian data. For each dataset with noise level a,, we applied the net upper bound with 10dd trials, the LP lower bound with L {-o, O, of and m = 30, and the baseline lower bound with L (-o, 0, a) and 1000.
28

---

Published as a conference paper at ICLR 2023
Boston Housing dataset. For the two-dimensional datasets, we applied the net upper bound with 100 trials and the LP lower bound with L = {0} and m = 100. For the three-dimensional dataset we applied the net upper bound with 1000 trials and the LP lower bound with L {0} and 30.
I.3 INDEPENDENT TRIALS & ERROR BARS
For randomized algorithms with two-sided errors, standard practice is to run the algorithm multiple times on the same dataset and report c.g. the median and eπror bars. However, our algorithms provide unconditional upper/lower bounds on the stability, so this is not necessary. If we did run the LP lower bound multiple times, we could simply report the maximum outcome, and this would be a valid lower bound on the stability; similarly for the net upper bound we could report the minimum outcome.
However, for our synthetic wwo-dimensional datasets (Figures 1 and 2a), because the data itself is random, we construct each dataset 10 independent times. For each algorithm and each dataset we compute the stability bound. Then, we report the median bound and error bars across the 10 wuw.eett-----ieeeiill-----oo
datasets. All error bars are 25th and 75th percentiles.
There are no error bars in our synthetic three-dimensional experiment (Figure 2b) because we only conducted one trial per noise level (rather than 10) due to computational constraints.
1.4 COMPUTATIONAL DETAILS
All experiments were done in Python on a Microsoft Surface Laptop, using GUROBI (Gurobi Optimization, LLC, 2022) with Academic License to solve the linear programs. Each plot took most 30 hours to generate (specifically, the three-dimensiona Gaussian isotropic data experiment took 3 hours for each of the 10 datasets, dominated by the time required for the LP lower bound algorithm; other plots were faster)
1.5 OMITTED EXPERIMENT: COVARIANCE SHIFT
Consider a dataset with n samples (Xa, yt) drawn from X; (0, and yi -X11 + X12-
1 -1]
where E -1 2 Additionally, there are ls outliers, in two types: the i type-I outliers
(X,, 9:) have X; = c(1,-3) and = -C large and negative; the one type-Il outlier (X,, B:) has
X, (n(1,1) 1) and y chosen to exactly lie on the OLS best-fit hyperplane y (x. 3). Initially, 3, > and clearly removing the last k samples suffices flip the sign. However, the initial
2 0
sample covariance is roughly Σ Σ + 11
3 The influence of sample (X,, Be) on
0
As result, the type-I outliers have negative influence on 91, so the greedy algorithm initially does not remove them. The influence only becomes positive after removing the type-II outlier, because this shifts the sample covariance to E, and therefore flips the sign of ((E (X)
coordinate of the OLS regressor 3 is ((--))) X.)
(x,, B)).
Constructing this example experimentally requires some care in the choices of k, c, and C: we need the total influence of type-A outliers (proportional to kcC) large enough that 3 > O, and need kc small enough that the type-I outliers don't affect the sample covariance much. Moreover, the number of trials needed by the net algorithm roughly scales with C. We take 1000, k 30, 0.2, and C = 300. Applying the greedy baseline and the net upper bound (with 1000 trials), we find that the former removes 97 samples while the latter removes roughly 16.7 samples. In this example, the failure of the greedy baseline can be attributed to the constant- factor shift in the sample covariance achieved by removing the type-B outlier.
29

---

Published as a conference paper at ICLR 2023
1.6
OMITTED FIGURES
15
50
trggtabe sidioed
Lenfrate Sildoodd
Figure 4: Coefficients of OLS regression for the suburb/city split (Figure (a)) and random split with the same ratio (Figure (b)). In both cases, we first rescale all covariates to mean 0 and variance 1, and add a constant variable (the last coefficient in both plots). Figure (a), we then plot the OLS regressor all samples (in blue), the OLS regressor all 134 zn 0 samples (in orange), and the OLS regressor on all 472 zn = O samples (in green). In Figure (b), we pick a subset of samples of expected size 134. We plot the OLS regressor all samples (in blue), the OLS regressor S (in orange), and the OLS regressor on S (in green).
J FORMAL PSEUDOCODE FOR ALGORITHMS
Algorithm NETAPPROX
procedure NETAPPROX ;
2:
X [(x) ..: ".....,('')) (x")
3:
if XT XXX - y) = has a solution then
4:
return 0
5:
7 6-22
6:
A [(X)) :2....; (X²) - d' rank(A)
7:
U.D.V t SVD(A) (so that A = UDVT and D is diagonal matrix of the d' nonzero
singular values)
8:
B <- VD-1
9:
Define y-net for Sdl -1 by by
M + ( im, m,...., my E /0.1///v// 2
where R d' x d' uniformly random rotation matrix Define 7-net for "residual space" (Lemma E.1) by
(Bm)) E M (Bm)a
N
11:
$ -
for N do
13:
Solve linear program
[Ill|h x (XA v) = =0]
V(A)= = sup we(o.1)
S v()))
15:
return
30

---

Published as a conference paper at ICLR 2023
Algorithm 2 GENERATEREGIONS
1: procedure (ENERATERECIOOS((SS ())
2: R R }
List of regions, each described by linear constraints
3:
for (Us,t,) do
4:
Rnext f- R
5:
for R E do
6:
Solve linear programs
Clow
inf (ve, x) zeit
and
Ohigh
sup (Ui, x) sER
7:
Clow < c < Ohigh then Add regions Rn ((e,, x) S ci} and Rn ((0,,)) 2 ci to Rnew
8:
9:
else
10:
Add region R to Rnew
R < R
return R
31

---

Published as a conference paper at ICLR 2023
Algorithm 3 PARTITIONANDAPPROX
procedure PARTTIIONANDAPPROXXXXXX
2: B(0) < arg min ||X - vll
x n matrrx with rows X1,... Xn
3: M < |x3(0) 1/12
4:
X + [X*) ) (X*) [4]
5:
Cc Cc-l log c ll
C 0 is universal constant
6:
Sample i,...., im E [n] without replacement
7:
H 1++)
8:
for E [m] do
9:
Add (X; vi) to H
Describes hyperplane (X,A) y:
10:
Add (Xi; + M), (Xi,y: - M), (X.,2 8M/Vn) and (X., vi 5M/ vT) H
for E [m]] do
12:
for 0 S k S [ (1/5) do
13:
Add (Xi, Big + M)) and (X,, + c) H
14:
R < GENERATEREGIONS((A)
for R E R do
16:
Find feasible point Ao RR) E R
17:
Вм +- {i [n] [XX, Ao (R)) yil M}
18:
BaM E [n] (X,, do (R)) yıl < 6M/ vn}
19:
BBMl > en then
20:
return
21:
for i E [m] do if inf, A) - y 2 then
22:
24:
else
0j Solve linear program
26:
(0(R),A(())
X =0
0 9i S (X,, 2) y
Vi [n] or
Σ
arg sup gER" R
8i
(N,, S 9: 50
Vi E [m] o =-1
(BA
8: 51
Vi E [n] (BM U BBM) We [N]
8; S g./ ((8, do( - vs)
(BM U BBM)
return
0(R);
max RER
(X,, ((1T)) y/
sE(m)
32

---

Published as a conference paper at ICLR 2023
Algorithm 4 LP Lower Bound
1: procedure L.P.OWERBOCND((X, ;, 9;) L,,,L,, m)
2: ((0) + arg min ||X - yll2
3: X [XTT ..: (XT) )a]
4: Sample 1 , im [n] without replacement
5: H
x-( } for j E [m] do
6:
7:
for t E do
8:
Add (X,, Di, + t) to H.
9:
R <- GENERATEREGIONS(H))
for R. E R do
11:
For each i [n]. compute
4 inf (X,, A) y/
and
sup Xi, A)- yr.
AER
xER
Set K_ fie [m] 2 0} and K_ = fiE [n] r; 0} K_
13:
Solve linear program
x" v) 0s 91 E (X,,
5(R)
inf AER
9i
Σ
9i
ViE K_
7
(X,, v1 $ 50 0
ViE K_
KEK
EK
4. 91 S:
Vi E [M] (K+ (K-)
14:
For each i E K4 U let Q. = min(?? +,); for each E [n]
(K₁ UK_ let Q; =
15:
Let Q(1) < ... < Q(1) be the sorted numbers Q .I,..., Qu and compute
X - 112 Q
N(R) - sup k
Σ
16:
return
min (a)) NER
33


Reviews:
[{"review_id": "yixWEXnIxmM", "reviewer": "Reviewer_BXKi", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The paper provides algorithms for estimating stability of OLS to dropping a small fraction of samples. The theoretical results focus on low-dimensional regime and give algorithms with run-time exponential in $d$. ", "summary_of_the_review": "Given the strengths of the paper, I recommend accept. ", "strength_and_weaknesses": "The paper is clearly strong at a technical level. The results are sound and complete. The upper bound for exact estimation, though somewhat unsatisfying in terms of run-time, is complemented by the faster approximation algorithm and the lower bound. \n\nThe author(s) provide an extensive set of experiments, which is great for a work that is theoretical by nature. \n\nThat said, I can see that the paper can be strengthened in a few ways. \n\n1. Can any of these results by extended to regularized (ridge) regression? \n2. Is it possible to identify certain distributional assumptions so that the run-time lower bound can be bypassed? It would be great to demonstrate poly-time algorithms even though it has be to under further conditions. \n3. The prior work of Broderick et al., 2020 (https://arxiv.org/abs/2011.14999) considers a broader setting, where there's a general quantity of interest $\\phi : \\mathbb R^d → \\mathbb R$ (mapping from the $d$-dimensional parameter $\\beta$ to a real) and the algorithm is asked to check whether $\\phi$ is stable wrt the dataset. The setting studied in this work is strictly more restricted, since it's only concerned with a single coordinate.  I am wondering if the results can be extended further. ", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works.", "empirical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works."}, "extra_text": {}}, {"review_id": "wrAHBgkVOr", "reviewer": "Reviewer_hNxp", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "Given a matrix $X\\in\\mathbb{R}^{n \\times d}$ and vector $y \\in \\mathbb{R}^{n}$, we are interested in understanding the _stability_ of the least squares problem $\\min_{\\beta} \\|\\|X \\beta - y\\|\\|_2$. In the prior work before this paper, _stability_ is defined as the smallest number of rows of $X$ that need to be remove to make the (wlog) first entry of the least squares solution $\\beta$ flip sign. This is computationally hard, and so the prior work relies on heuristic or brute-force algorithms.\n\nThis paper considers a fractional notion of stability, where we are able to reweigh each row of $X \\beta - y$ by some $w_i \\in [0,1]$, where $w_i=0$ means entirely removing the row and $w_i=1$ means keeping the row, and ask what is the smallest sum of $(1-w_i)$ needed to make the first entry of the least squares solution flip sign. The paper focuses on the computational hardness of this problem, and shows three key results:\n1. Exactly computing fractional stability is possible in $n^{O(d^3)}$ time.\n1. Exactly computing fractional stability requires $n^{\\Omega(d)}$ time, by reduction to the Exponential Time Hypothesis.\n1. Under a \"mild anti-concentration\" assumption that few rows of $[X \\beta]$ are nearly linearly independent from other rows, a $\\tilde{O}(n + \\frac{d}{\\varepsilon^2})^{d + O(1)}$ time algorithm approximates fractional stability to additive error\n1. Under a stronger near-incoherence assumption, a $O((\\frac{\\sqrt d}{\\varepsilon})^d \\text{poly}(n))$ time algorithm approximates fractional stability to additive error\n\nCompelling experimental results show that these algorithms can be practical on low-dimensional datasets, giving provable stability in potentially tractable runtime.", "summary_of_the_review": "The setup and results are cool, but I don't completely understand the technical proof sketches for correctness. I understand a lot of it though, so I'm still confident enough to give the paper a pass.", "strength_and_weaknesses": "I like this paper. I recommend it be published.\n\nAt a high level, it takes an intractable problem (the prior work's notion discrete of stability), keeps the intuitive underlying structure, and extends it to a fractional problem that's relatively tractable.\n\nThe paper is decently well written, with the high level idea of the algorithms being broadly clear, but some parts of the proof sketches being hard to understand. The results form a decently satisfying narrative: a terrible runtime without any assumptions, a lower bound to boot, a mild assumption to improve the terrible runtime with some approximation error, and a stronger assumption to get a potentially tractable runtime with some approximation error.\n\nThe results seem significant. The lower bound from the exponential time hypothesis makes it clear that this problem is hard, so fighting to get a runtime that's just polynomial in $n$ makes sense. Further, since verifying the stability of models is high-stakes in high-stakes applications, we really don't want to rely on heuristic algorithms. Provably correct algorithms are important. So, while the fastest runtime presented is pretty horrendous for even small $d$, it's nevertheless a good win.\n\nThe results certainly have weaknesses though, though I'm inclined to be lenient towards these weaknesses given the hardness of the problem. The key weaknesses are:\n1. The exact algorithm runs in $n^{O(d^3)}$, the lower bound for exact algorithms is $n^{\\Omega(d)}$, and the approximate algorithm with mild-anticoncentration assumption runs in $\\approx n^{O(d)}$ time. In order to just match the lower bound, the paper both relaxes the exactness of the answer and assumes structure on the input.\n\n    $\\phantom{.}$\n\n    It's not clear _where there's looseness_. Can an exact algorithm run faster, at $n^{O(d)}$? Does mild-anticoncentration suffice to beat the $n^{\\Omega(d)}$ lower bound, just using a different algorithm? Does approximation error allow us to beat the $n^{\\Omega(d)}$ lower bound? Is there a $n^{\\Omega(d)}$ lower bound for approximate algorithms given inputs that satisfy mild-anticoncentration? There's some gap here that I don't really understand, nor have any intuition for.\n\n    $\\phantom{.}$\n\n2. The proof sketches are a bit unclear. I can follow 90% of the intuition, but then some language about _demarcating regions_ or _connected component of varieties_ appears to resolve a conceptual technical block, and I don't really get how this work. I'm sure these are elaborated upon in the appendix, but that shouldn't be necessary. This hurts my confidence of correctness of the algorithms.\n\n    $\\phantom{.}$\n\n3. The final approximation algorithm, with it's stronger assumption on the input, is similar to the _incoherence_ assumption made in the matrix completion literature, which intuitively says that sampling rows from a matrix uniformly-at-random suffices to preserve $\\ell_2$ norm information (see [here for a formal paper](https://arxiv.org/pdf/1408.5099.pdf) and [here for a more approachable writeup](https://randnla.github.io/leverage-subspace-embedding/#uniform-sampling)). It's certainly weaker than a real incoherence assumption, since an $\\varepsilon$-fraction of rows are allowed to violate incoherence, but we can still contrast the assumptions.\n\n    $\\phantom{.}$\n\n    For an incoherent matrix, any uniform subsampling of $O(d \\log d)$ rows preserves the least squares solution with high probability. So, we should intuitively find that such matrices have to be stable. Admittedly, this notion of \"preserving\" a solution is not obviously related to stability, so my intuition could be really betraying me here. At the same time though, it remains unclear how much of the weight of this result is owed to the algorithm, and how much is owed to the incoherence-like assumption. Either way, it's a nice result to have out there, but I think there's a nice discussion to be had about how this differs from incoherence.\n\n4. Both assumptions on the input data, the mild-anticoncentration and the near-incoherence, read sorta like assumptions on stability. It's not clear if we should be assuming some kind of stability to \n\n\n---\n### Appendix: Formalizing the Incoherence Relationship\nFor the authors, I formalize what I see at the relationship between incoherence and the assumption used for the last algorithm of their paper. First the _leverage score_ of row $i$ of a matrix $\\bar{X}$ is\n$$\n\\tau_i[\\bar{X}] = \\max_{\\beta} \\frac{([\\bar{X}\\beta]_i)^2}{\\|\\|\\bar{X}\\beta\\|\\|_2^2}\n$$\nThen $\\tau^* := \\max_i \\tau_i$ is the incoherence of the matrix. In other words,\n$$\n\\forall \\beta,\\ \\forall i\\in[n],\\ \\frac{([\\bar{X}\\beta]_i)^2}{\\|\\|\\bar{X}\\beta\\|\\|_2^2} = \\tau^*\n$$\nAssumption B in the paper, when $\\varepsilon \\rightarrow 0$, says that $\\tau^* = \\frac{\\delta^2}{n}$, so that (intuitively) uniform $\\tilde{O}(\\delta^2 \\log(d))$ rows is needed to preserve $\\ell_2$ norms (assuming we pick $\\delta$ to be as small as possible for $\\bar{X}$).", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "Tr9iAOSZUwz", "reviewer": "Reviewer_Hpqm", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The paper studies the stability of linear regression problem by defining it as the minimum samples (in $\\ell_1$-norm) to be removed to zero out the first coordinate. The paper provides a $O(n^{d^3})$ algorithm to exactly solve the problem, and also provide a lower bound stating that the computation complexity cannot be $O(n^{o(d)})$. Further, approximate algorithms are provided and experiments are carried out on real dataset to identify stability.", "summary_of_the_review": "The paper is well written and organized. The problem is clearly formulated, and the authors tackle the exact stability problem instead of using local approximations from influence functions, which is novel and interesting to me. By exactly auditing linear regression, the authors provide theoretical guarantees for such identification and the algorithm works well on real datasets. Although the complexity of the algorithm seems to suffer from curse of dimensionality, but it is a good starting point for study stability problem in its exact form. I think this is a good paper and would like to accept.", "strength_and_weaknesses": "Unlike the previous results using influence functions to approximate the local stability, the paper study the exact problem--which means it is able to identify perturbation of the solution beyond local regime. The theory and experiments are sound. However, even with the relaxed approximate algorithm, the computation cost seems still very big.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "GeeTdxgsPM", "reviewer": "Reviewer_LbgB", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "The paper presents a fast algorithm to estimate a stability metric of ordinary least squares problems (OLS). The stability metric was introduced as\nthe minimum number of samples that need to be removed so that rerunning the analysis overturns the conclusion. Naive computation of the stability metric is computationally prohibitive, whereas the proposed algorithm can efficiently obtain an approximation in the low-dimensional regime.\n", "summary_of_the_review": "The authors propose an efficient approximation algorithm for computing the stability metric introduced in (Broderick et al., 2020). Theoretical analysis is provided to justify the robustness of the algorithm.", "strength_and_weaknesses": "[Strengths]\n- Evaluating the stability of OLS could be useful for robust regression and model interpretation.\n- The paper presents polynomial-time algorithms for estimating a particular stability metric introduced in (Broderick et al., 2020). Theoretical guarantees are also provides for the approximation quality.\n- The authors illustrated the utility of the method on Boston Housing Data.\n\n[Weaknesses]\n- The proposed method seems to be limited to only linear models/OLS problems. Generalizations to non-linear models or more advanced models are not obvious.\n- The proposed method targets one particular stability metric (Broderick et al., 2020). This metric might be sensitive to noise in the data or resampling of the data. There are widely-studied candidate stability metrics such as the Cook's distance in regression analysis that could afford cheaper computations.\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}]

Author rebuttal:

Decision: Accept: poster