Insufficient Statistics Perturbation: Stable Estimators for Private Least Squares
Gavin Brown Jonathan Hayase Samuel Hopkins Weihao Kong Xiyang Liu Sewoong Oh Juan C. Perdomo Adam Smith


Stellar atmospheric parameters of ∼ 11,000 RR Lyrae stars from LAMOST Spectra
    Chunqian Li
    April 29, 2024
=============================================================================




    We present a sample- and time-efficient differentially private algorithm for ordinary least squares, with error that depends linearly on the dimension and is independent of the condition number of X^⊤ X, where X is the design matrix.
    All prior private algorithms for this task require either d^3/2 examples, error growing polynomially with the condition number, or exponential time. 
    Our near-optimal accuracy guarantee holds for any dataset with bounded statistical leverage and bounded residuals.
    Technically, we build on the approach of <cit.> for private mean estimation, adding scaled noise to a carefully designed stable nonprivate estimator of the empirical regression vector.





§ INTRODUCTION


We present a sample- and time-efficient differentially private algorithm for ordinary least squares (OLS) regression.
Central throughout the theory and practice of data science, OLS is used in numerous domains, ranging from causal inference, to control theory, to (of course) supervised learning. 

Given covariates X∈ℝ^n× d and responses y∈ℝ^n, the OLS estimator is defined as

    β_ols = *X^ X^-1 X^ y .

Among the many reasons for the popularity of OLS is the fact that it is a statistically and computationally efficient way of solving linear regression. Speaking informally, OLS has low excess error whenever the number of samples n is as large as the problem dimension d. Crucially, its statistical performance does not depend on the condition number κ(X^ X), the ratio between the maximum and minimum eigenvalues. 
Furthermore, it can be computed in closed-form using only basic linear-algebraic operations, with no need for the subtle hyperparameter tuning often inherent in first-order methods.

Given its widespread use in the analysis of personal data, there is a long line of work giving differentially private algorithms to approximate OLS.
However, designing practical and efficient algorithms for this problem has been a particularly challenging endeavor; so far there are no clear answers even in the d=1 case when x and y are both scalars <cit.>. Existing algorithms for DP regression suffer from one of three limitations: they either have poor dimension dependence in their sample complexity, place unnaturally restrictive assumptions on the geometry of the data, or run in exponential time. 

In terms of private algorithm design, one natural and well-established approach is sufficient statistics perturbation, which privately produces separate estimates of X^ X and X^ y and then combines them to produce a single parameter estimate.
Such approaches are often efficient and some versions come with formal accuracy guarantees. 
An exemplar in this line is the AdaSSP algorithm of <cit.>.
The central drawback in all these algorithms, however, is the sample complexity as d grows: privately producing an accurate estimate of X^ X requires roughly d^3/2 samples <cit.>.
Furthermore, many approaches within this class add noise proportional to the worst-case sensitivity of X^ X and X^ y <cit.>. 
To deal with the fact that this sensitivity is unbounded in the case of real-valued data, these results assume uniform norm bounds on the covariates x and responses y (e.g., x≤ B_y, |y| ≤ B_y). While conceptually simple, they fail to capture the intrinsic complexity of the problem and do not satisfy natural properties like scale invariance.






Another approach comes via private optimization, searching for a parameter estimate that approximately minimizes the sum of squared errors.
Despite a wealth  private convex optimization methods that can be applied directly to linear regression, 
off-the-shelf approaches again require d^3/2 samples for accurate estimates.
A notable exception is the recent work from <cit.>, whose algorithm based on private gradient descent succeeds with only roughly d samples.
However, its error grows with the square of the condition number, a high price to pay for many problems.
A polynomial dependence on κ(X^ X) is inherent in private first-order optimization for linear regression, as the smoothness of the optimization task is directly linked to the condition number.

The only known approach that avoids these two issues is the exponential-time algorithm of <cit.>, which comes from the framework they call high-dimensional propose-test-release (HPTR), after the propose-test-release (PTR) framework of <cit.>.

We see the mirror of this story in private mean estimation, where <cit.> and <cit.> recently gave the first sample- and time-efficient private algorithms with error guarantees that adapt to the covariance of the data.
All prior private algorithms achieving this guarantee require d^3/2 examples, error growing polynomially with the condition number of the covariance, or exponential time.

In this work, we build on the work of <cit.> and present the first computationally efficient (in fact, practically implementable) differentially private estimator for linear regression with sample complexity independent of κ(X^ X) and the optimal linear dependence on the dimension d.
Furthermore, we make no use of norm bounds. We establish its  utility  under the “textbook” conditions one would typically require to run OLS in the non-private setting. More specifically, the algorithm is accurate as long no observation has high statistical leverage or a large residual, formalized in <ref>.





 §.§ Our Results


In this work, we introduce a new algorithm, , for differentially private linear regression. At a high level, works in two main phases. In the first, we search for a reweighting of the dataset such that running OLS on this reweighted version is roughly stable. Having successfully found this set of weights, we simply compute the OLS solution on this weighted version of the data and add appropriately  shaped Gaussian noise to the solution. While the approach is conceptually very simple, establishing its correctness requires several significant technical advances.








Our estimator satisfies differential privacy (DP), the gold standard for privacy protection in statistical data analysis. DP requires that an algorithm provides approximately the same output on any datasets that differ in only one entry.

Let 𝒳 and 𝒴 be sets. An algorithm 𝒜: 𝒳^n →𝒴
is (, δ)-differentially private if for every x = (x_1, …, x_n) ∈𝒳^n and x' = (x_1',…, x_n') ∈𝒳^n such that x, x' agree on all but one coordinate and for all Y⊆𝒴, 

    [](x) ∈ Y     ≤   e^ [](x') ∈ Y + δ .




 One of the core advances we make, in light of most previous results on DP regression, is that we do not require any norm bounds on the data. We only assume the types of conditions a circumspect statistician would always verify to ensure that OLS is a sensible procedure.
In particular, we establish the utility of our estimator whenever the dataset is free of outliers, or “good.” 

definition]def:goodness_set
    Fix parameters ,  > 0.
    A dataset X, y∈^n × d×^n is called ,-good if X^ X is invertible and the following conditions hold for all i ∈n.
    
        
  * Bounded leverage: x_i^ (X^  X)^-1 x_i ≤.
        
  * Bounded residuals: []x_i - y_i≤.
    




Note that both of these conditions hold in various natural, well-studied settings. For instance, when x,y are both subgaussian and drawn from a well-specified linear model with true parameter β^*, these conditions hold with high probability when L ≈ d/n and R ≈σ where σ^2 = (y - ⟨β^*, x⟩)^2 (see <ref>). 
This idea, of outliers being observations with high leverage or large residuals,
is quite classical and found across standard texts.
For instance, a standard rule-of-thumb identifies high-leverage points as those with leverage greater than 2d/n or 3d/n <cit.>.
The precise forms of stability we need to ensure privacy, however, are far from classical.
They require a carefully designed algorithm, which we elaborate in <ref>. 
Such a stable estimator, in turn, implies that we can achieve differential privacy with small amounts of noise. 



In these well-behaved instances, 
our mechanism takes a most straightforward form:
it returns the OLS solution plus a small amount of Gaussian noise. 

claim]obs:output_on_good_data 
    If (X,y) is (,)-good for parameters   R > 0 and ≤  
    c'ε^2 log^-2 (1/δ), for some constant c', then 
    (X,y;,δ,,) 
    releases a sample drawn from [], c^2 (X^ X)^-1, where c^2=Θ( ^2 log (1/δ)/^2).

Attentive readers will detect a modest sleight-of-hand: (,)-goodness is a property of the data and a priori unknown, yet the algorithm gets  and  as inputs!
Nevertheless, an analyst with beliefs about the data generation process can set these parameters appropriately.
The maximum leverage score does not depend on the scale of the data, only its concentration properties.
Since it lies within [0,1], one might pick the  hyperparameter adaptively by calling a handful of times.
Similarly, if the analyst believes the labels are generated by a process such as y_i x_iβ^* + (0,σ^2), they can privately produce an accurate estimate σ using standard tools (see <ref>).
We believe alternative standardized or studentized definitions could remove this need for prior knowledge about σ.
These alternatives would likely increase the complexity of our proofs.

The difficulty in our work lies in proving that is differentially private (<ref>); reasoning about utility is simple once we have <ref>.
More specifically, seeing how the output distribution on good data matches standard statistical practice (and classical CLT-like analyses of OLS), we can quickly derive error bounds. 
For instance, in the simplest case of fixed design, where we consider only the randomness of the labels generated from a well-specified linear model y_i = x_iβ^* + (0,σ^2), we have ∼(β^*, σ^2· (X^ X)^-1). Hence, from <ref>, we see that, relative to the empirical OLS solution, the private estimator is just a slightly noisier version of the true parameter (and has the same kind of error covariance). 

More formally, we can analyze the mean squared error (MSE) of our algorithm on any good dataset.

corollary]cor:fixed_design_closeness
Fix (X,y), >0, and δ∈[0,1].
If (X,y) is (,)-good for ≤  c'ε^2 log^-2 (1/δ) for some constant c' and R>0, then X,y;,δ,,, releases β̃ such that, for some absolute constant c',
    
    *1/n[]y - Xβ̃^2   =  1/n[]y - X^2 + c' ^2 d/n log 1/δ/^2.



By <ref>, we have - β̃= c· (X^ X)^-1/2 u for u∼(0,). 
We then expand:

    _u []y - Xβ̃^2    = _u []y - X + X(-β̃)^2 
       = _u []y - X + X · c (X^ X)^-1/2 u^2 
       =y - X^2 + c^2 ·_u [u^ (X^ X)^-1/2 X^ X (X^ X)^-1/2 u],

where the cross terms drop out as u is independent and mean-zero. 
The matrices cancel and we are left with [u^ u], which is exactly d.


We emphasize that this result holds without any assumption that the data arises from a specific family of distributions. 
It assumes (X,y) is fixed and (L,R)-good to bound the difference from the empirical OLS solution on (X,Y). 
However, if we do add such distributional assumptions, it is easy to show that our algorithm produces a private parameter estimate that closely approximates the true regression parameter. We state this fact as part of the following theorem, our main result.







    Fix ε, η∈0, 1, δ∈0, ε/10, and n, d ∈.
    takes a dataset X, y∈^n × d×^n, privacy parameters ε, δ, and outlier thresholds _0, _0.
    
        
  * is ε, δ-differentially private.
        
  * Let X ∈^n × d be drawn i.i.d. from a d-dimensional subgaussian distribution 𝒟 with mean 0 
        and covariance Σ≻ 0.
        Let y_i = β^ x_i + z_i where the z_i are drawn i.i.d. from a subgaussian distribution with mean 0 and variance σ^2
        (see <ref> in <ref>).
        [We state the utility guarantees of our estimator for the case where data is drawn from a well-specified linear model to simplify the presentation and enable direct comparisons to previous work. However, as per <ref>, on good data our algorithm is always close to the OLS solution. Hence, we can prove closeness to the population quantity whenever the OLS solution concentrates.]
        If 
        _0 = Θ*d/n, _0=Θσ,
        and
        
    n   =  Ω*d/α^2 + d√(log1/δ)/αε + dlog 1/δ^2/ε^2 ,

        with a large enough constant for some α>0, then  returns β̃ such that, with high probability, 
         
    []β̃- β_Σ  ≤  σα.

        Here, Θ and Ω hide logarithmic factors in 1/α, log1/ε, and log1/δ as well as polynomial factors of the subgaussian parameters.
        
  * <ref> can be implemented to require
             one product of the form A^⊤ A for A∈^n× d,
             one product of the form A B for A∈ R^n× d and B∈^d× d,
             one inversion of a positive definite matrix in ^d× d; and
             further computational overhead of Õ(nd/).
    

Informally, this running time corresponds to O(nd^ω-1 + nd/), where ω<2.38 is the matrix multiplication exponent.
For modest privacy parameters, the running time of our algorithm is dominated by the time needed to compute the nonprivate OLS solution itself.

This is the first computationally efficient algorithm whose sample complexity is linear in d and has no dependence on the condition number κ(X^⊤ X). 
This almost matches the best known sample complexity of an exponential-time algorithm from <cit.>; we have an additional d(log(1/δ))^2/^2 term, but this term does not depend on the final accuracy α. 


[h]
InputInput
Requirerequire
dataset X, y; privacy parameters ε, δ, outlier thresholds _0, _0
k*12 log 3/δ/ε + 8;  c^2  56448exp[]432 k^2 _0_0_0^2·log12/δ/ε^2 
_0 > 1/96k or _0 > 3ε/56log 12/δ
 

_1, w X, _0, k*<ref>
_2, v X, y, w, _0, _0, k*<ref>
ε/3δ/3max_1, _2 =
     

_v  X^v X 
β̂ (_v)^-1 X^v y*OLS weighted by v 
β̃∼[]β̂, c^2_v^-1 

InSufficient Statistics Perturbation ()



We now briefly sketch the steps of the proof and discuss the paper's organization.
We establish <ref>'s subclaims in <ref>.
As outlined above, the utility analysis is straightforward once we have <ref> in hand.
The full analysis is presented in <ref>.
It is easy to see that <ref> runs in polynomial time. 
In <ref>, we analyze a careful implementation.

The bulk of the work comes in the privacy analysis. 
In <ref>, we analyze the greedy residual thresholding algorithm, with the main result about that algorithm being <ref>, the “intertwining” property.
Then, in <ref>, we establish our guarantees for .
The main results about are <ref>, which says that the score is low-sensitivity, and <ref>, which says that the weights are stable.
<ref> pulls these together to establish the privacy of .

<ref> covers additional related work.
<ref> provides necessary preliminaries.
<ref> contains proofs deferred from the main text.
<ref>, via standard tools, shows how to privately estimatie .
<ref> contains details on the lower bound of <cit.>.



 §.§ Optimality


For modest values of the privacy parameters, the error of our algorithm is dominated by the empirical error of OLS. 
Informally speaking, we obtain privacy “for free.”

Formally, our error guarantees are close to tight for random-design regression with subgaussian covariates and subgaussian label noise.
Suppressing constants and logarithmic factors other than log 1/δ, <ref> says that we can achieve β̃- β_Σ≤σα with high probability with

    n ≈d/α^2 + d√(log 1/δ)/α + d(log 1/δ)^2/^2.

Known lower bounds imply this task requires

    n ≳d/α^2 + d/α + log 1/δ/.

The first term corresponds to the classical analysis of OLS.
The second term was established by <cit.> and holds even for parameter estimation in ℓ_2 norm; see <ref> for a more detailed discussion.
The third term, the minimal number of samples required to produce any estimate, is from <cit.> and holds even for one-dimensional mean estimation with known variance.
The exponential-time algorithm of <cit.> nearly matches all three terms.
For constant  and δ = 1/poly(n), our algorithm's error guarantee in this setting is tight up to logarithmic factors.
An exciting topic for future work is determining the existence or impossibility of efficient algorithms with error matching <ref> up to constant factors.



 §.§ Techniques



At a high level, our algorithm follows the blueprint for private mean estimation laid out by <cit.> and made computationally efficient by <cit.> and <cit.>.
Our approach closely follows that of <cit.>, henceforth BHS.
We now sketch our algorithm, discuss how our analysis differs from that of BHS, and investigate how the notions we use are, in a sense, “correct” for the task of private least squares.



  
Overview of  Perhaps the most natural approach for private estimation of regression coefficients is to perturb the ordinary least square estimator, . However, without restrictions on the data, the sensitivity of  is unbounded. 
Our key observation is that, on datasets with bounded leverage and bounded residuals,
the OLS solution is actually quite stable. 
If we could restrict our inputs to only such outlier-free data sets, we might hope to release  plus noise with shape (X^⊤ X)^-1.

While this would provide accuracy, it fails on privacy: we must accommodate worst-case data. 
We use the PTR framework of <cit.> to test if our input contains a large good subset. 
We propose a greedy pruning algorithm which, in each iteration, removes the data point with the largest residual and recomputes OLS  on the remaining data. 
Similar approaches abound in the literature on robust statistics, but we prove key new properties about how this algorithm behaves across adjacent data sets and different outlier thresholds.



  
Adaptively selecting outlier thresholds
Our algorithm takes as input target bounds  and  for the leverage and residuals, respectively.
This simplifies our analysis but is not strictly necessary.
The maximum leverage can only lie within the interval [0,1], so one could imagine calling repeatedly within this space (via a well-chosen grid or binary search) to find an appropriate setting, perhaps via a small validation set.
Independently, one could privately learn an appropriate value for  directly through standard techniques; we give a complete description in <ref>.



  
Proof techniques 
While our work builds on a long line of research connecting robust statistics and differential privacy, it especially relies on the 
recent algorithmic approach of BHS, who gave improved algorithms for private mean and covariance estimation.
At a bird's eye view, our recipe for private linear regression follows the main ideas behind the mean estimator of BHS. However, key parts of the implementation and analysis differ significantly in the more complicated linear regression setting. 

We start by discussing the ways in which our main proof strategy is similar to BHS.
As mentioned previously, we introduce a notion of “good” outlier-free datasets for linear regression.
We repeatedly call a greedy algorithm to find a series of good weight vectors 
across a range of carefully chosen outlier thresholds. 
We use these vectors to privately test that our input data is sufficiently close to the good set and to finally produce a vector of weights  over the input.
Crucially, this weight-finding procedure is stable: if run on any adjacent dataset, it would produce a vector that is close in ℓ_1 distance.



When adapting their analysis, however, we run into immediate issues. 
For both their definitions of good set, BHS prove a number of strong properties that are false in the context of regression.
For mean and covariance estimation, the good sets are unique (i.e., for any dataset and outlier threshold, there exists a unique largest good set), are directly found by the natural greedy algorithm, and enjoy a form of monotonicity (e.g., introducing a new point to the dataset cannot alter the good set very much). 
For the definition we use, there is no unique “largest good set” which introduces significantly more complexity in the analysis.
What's more, adding a single point may significantly affect the downstream choices made by the greedy algorithm which further complicate the relevant stability calculations. 

In more detail, a key step in the BHS analysis establishes the following “intertwining” property in the context of mean and covariance estimation. 
Suppose we call the greedy algorithm on a dataset D with outlier threshold B and find a largest good subset S⊆ [n].
If we then call the same algorithm on an adjacent dataset D' with a slightly larger outlier threshold B', the largest good subset T will satisfy the property that S⊆ T (ignoring the index that differs between D and D').

We establish an analogous statement (<ref>) about the output of our greedy residual thresholding, <ref>, even though we cannot prove the same uniqueness and monotonicity statements.
More specifically, we develop a novel regression-specific argument  that uses the closed-form expressions describing how the least squares solution changes when an observation is added or removed.
The exact arguments are formalized in <ref> and <ref>, but we  sketch the ideas here.

For a dataset (X,y) and index i∈ [n], let ŷ_i = x_i^ be the fitted value.
We denote by e_i the i-th residual: e_i ŷ_i - y_i.
Recall the hat matrix:

    H  X (X^ X)^-1 X^ ,

so called because it maps the true labels to their “hat” values: ŷ = H y.
The leverage scores (also known as sensitivities or self-influences) form its diagonal entries, while its off-diagonal entries will be called (by us) the cross-leverage scores:

    h_i = H_i,i   = x_i^ (X^ X)^-1 x_i    and    
        H_i,j = x_i^ (X^ X)^-1 x_j.

Note that, by Cauchy–Schwarz, the cross-leverages are no larger in magnitude than the leverages.

What happens if we remove an observation, say, (x_j,y_j), from the dataset?
This takes the form of a rank-one update. Applying the Sherman-Morrison formula we can derive closed-form expressions for the changes in the OLS solution as well as (for any i∈[n]) the leverage score and residual of point i after removing j.
Using the subscript “(-j)” to denote the quantity after removal, we have

    h_i - h_i(-j)   =  - H_i,j^2/1 - h_j
      - β_ols(-j)   = (X^ X)^-1 x_j/1-h_j·[]x_j - y_j
    
        e_i - e_i(-j)   = H_i,j·e_j/1-h_j.

These well-known formulas have elementary derivations; the second and third correspond to the DFBETA (“difference in β”) and DFFITS (“difference in fits”) regression diagnostics <cit.>.
All three seamlessly generalize to the case where the points are weighted. 
We can reuse them to reason about what happens when we add points to the dataset.

Beyond their use in our formal arguments, these formulas show how our goodness definition in <ref> is essentially the “right” one to analyze stability. The leverage score and the magnitude of the residual exactly determine the sensitivity of the least-squares solution to adding or removing that data point. 
To see this in more detail, consider the effect of dropping a point from a typical dataset:

    [](X^ X)^1/2 ( - β_ols(-j))^2 
               = *(X^ X)^1/2·(X^ X)^-1 x_j/1-h_j· (y_j - x_j)^2 
       = h_j· e_j^2/(1- h_j)^2 = Δ.

Arguing heuristically for now, if removing a point changes the OLS solution by Δ, to ensure privacy one must ensure noise of magnitude at least Δ. It is impossible to do any better. Note that by working with (,)-good sets we can guarantee that the noise we add for privacy,  (0,c^2 (X^ X)^-1) where c^2 ≥^2 ·log 1/δ/^2,  has magnitude roughly Δ. 
This insight shows our accuracy guarantees are sharp.




 §.§ Notation



We use n to denote the set 1, …, n and = 1, 2, ….
For a vector v∈^n its support is v = i ∈n| v_i ≠ 0.
If we have a set S⊆n, then vS∈^n has vS_i = v_i for i ∈ S and vS_i = 0 otherwise.
Also we define S = n∖ S.
We use vv_2 and v_^1/2v.
If M ∈^n × n is a matrix, then M_2 denotes its spectral norm.







§ ANALYSIS OF GREEDY RESIDUAL THRESHOLDING




Inputinput
Requirerequire
dataset X, y; outlier threshold R; starting weights w
𝚃𝚁𝚄𝙴
    β_w (X,y,w) 
    i^* _i∈w[]y_i - x_i^β_w  
    y_i^*- x_i^*^β_w > R
        w_i^* 0 
    
        w 
    




In this section we establish the key properties of our greedy residual-thresholding algorithm.
This analysis contains the bulk of the technical novelty in our work.
The main result is <ref>, the “intertwining” property that relates the outputs of on adjacent datasets.

Since we will be dealing extensively with weighted sets from now on, we expand the definition of good sets in <ref> to vectors of weights.
[(,)-goodness, weighted]defnngoodnessdef
definition]def:goodness
    Fix a dataset (X,y)∈^n× d×^n and parameters ,>0.
    A vector w∈ [0,1]^n is (,)-good for (X,y) if, denoting W=diag(w), X^ W X is invertible and the following two conditions hold for all i∈supp(w).
    
        
  * Bounded leverage: x_i^ (X^ W X)^-1 x_i ≤.
        
  * Bounded residuals: []x_iβ_w - y_i≤, where β_w = (X^ WX)^-1X^ W y.
    
    Furthermore, we will say that w is L,∞-good for X, y if (<ref>) holds, but not (<ref>).




 §.§ Stability and Goodness for Ordinary Least Squares


Our analyses rely on how goodness is affected when adding and or removing mass.
As discussed in <ref>, closed-form expressions characterize the effects of removing a single point <cit.>.
The following claim generalizes these results to removing multiple weighted points, or adding weight to points already included in the regression.
In addition, it shows how these results interact with goodness.
We defer the proof to <ref>.
[Changing Weight Within Support]
    Let w, w'∈[0,1]^n satisfy
        w'⊆w and w-w'_1 ≤1/2. If w is (L, ∞)-good for (X,y), then for all i ∈w,
    
    x_i^ (X^(w') X)^-1 x_i    ≤ (1 + 2w-w'_1) .

If, in addition to the previous conditions, it also holds that w is (L,R)-good for (X,y), then
    
    x_i^β_w - x_i^β_w'   ≤ 2 w - w'_1 .

    In particular, since w'⊆w,  Equations (<ref>) and (<ref>) apply to all i∈w'.
    Consequently, w' is (η,η)-good for (X,y), where η = 1 + 2w -w'_1.




We next present a claim about adding a point to existing good weights: either the 
expanded weights are good or the new point has a large residual (in which case our greedy algorithm, presented later, will identify it).
We illustrate these cases in <ref>.
Such a claim also holds when we add sets of points.

Mathematically, this proof contains little innovation beyond <ref>.
However, it provides a key conceptual bridge.
We see that it connects directly to our greedy algorithm, which removes large residuals.
[Adding Weight Outside Support]
    Let w'∈[0,1]^n be an (,)-good vector for a dataset (X,y) and 
    let v∈[0,1]^n be a vector such that w'∩v=∅.
    Define w=w'+v and η = 1 + 8v_1. Assume the following two conditions hold:

    
        
  * The matrix X^(w) X is invertible and for all j ∈w, 
    x_j^ (X^(w) X)^-1x_j≤ 2.

        
  * The weights v satisfy v_1 ·≤1/8  
    
    If max_i ∈w[]y_i - x_i^β_w > η, then _i ∈w[]y_i - x_i^β_w⊆v.


    We prove the contrapositive: if there exists j^* ∈_i∈wy_i - x_i^β_w with j^*∉v, then for all i  ∈w, y_i - x_i^β_w≤η.

    
    Note that since w = w'∪v and w'∩v=∅, j^*∉v implies j^*∈w'. 
    We first produce a lower bound on the j^* residual under w'.
    By the triangle inequality, 
    
    []e_j^*' = []y_j^* - x_j^*^β_w'   = []y_j^* - x_j^*^β_w + x_j^*^β - x_j^*^β_w'
       ≥[]y_j^* - x_j^*^β_w - []x_j^*^β - x_j^*^β_w'
       = []e_j^* - []x_j^*^β_w - x_j^*^β_w'.

    Note that by assumption, w is (2L, |e_j^*|)-good for (X,y). Since w'⊆w, and w-w'_1=v_1 ≤1/8, we can apply 
    <ref> to get that x_j^*^β_w - x_j^*^β_w'≤ 2 w-w'_1 (2L) |e_j^*|. Using this upper bound, we get that:
    
    []e_j^*'≥[]e_j^* -  4 []w-w'_1 []e_j^*

        To complete the proof, we use the upper bound e_j'≤, which holds by the assumption that j∈w' (and that w' is (L,R)-good). Rearranging our previous inequality, we get that for all i ∈w, 
        
    y_i - x_i^β_w≤e_j^*≤e_j^*'/1-4w-w'_1≤1 + 8 w-w'_1·,

    where we have used the inequality (1-z)^-1≤ 1 + 2z, which holds for all z ∈ (0,1/2]. 




 §.§ Guarantees for Leverage Filtering


receives as input a vector w, the “starting weights,” and iteratively zeros out any weights corresponding to residual outliers, 
recomputing the weighted OLS solution as it goes.
These starting weights w will come from the subroutine of BHS, which filters out high-leverage outliers.
The exact algorithm we use differs superficially from the version in BHS, who use it for covariance estimation and call it “Stable Covariance.”
Our application only needs its properties as a leverage-filtering procedure.
We give a complete description of our variant in  <ref>.

The filtering algorithm has “goodness” guarantees (when the score is modest, many points receive weight and no point with weight has high leverage), utility guarantees (on outlier-free data, all points receive full weight), and stability guarantees on adjacent datasets (the score is low-sensitivity and the weights are stable).
We now give the formal statement.
    
stabcovcount

    There is a deterministic algorithm receiving as input a list of vectors X∈^n× d, a leverage threshold , and a discretization parameter k∈ and returning as output an integer  and a vector w∈ [0, 1]^n.
    Let W = (w).
    Assume kL≤ 1.
    If <k the following hold.
    
        
  * w_1≥ n-k.
            As a consequence, w≥ n - k.
        
  * For all i∈w, we have x_i^ (X^ W X)^-1 x_i ≤.
        
    
    On “outlier-free” data as defined below, the algorithm's output is as follows.
    
        
        
  * If x_i^ (X^ X)^-1 x_i ≤/2e^2 for all i∈[n] then =0 and w=1⃗. 
        
    
    To present the stability guarantees, let X and X' be datasets that differ in one entry.
    For any values of k and , consider
    
    , w    (X,,k) 
    ', w'    (X',,k).

    We have the following sensitivity bounds. 
    
        
        
  * -'≤ 2.
        
  * If ,'<k then w-w'_1≤ 2.
    




 §.§ Properties of 


The first claim we prove says that, when we run followed immediately by , the returned weights are good.

    Let (X,y) be a dataset, k∈ be a discretization parameter, and 
    ,>0 be outlier thresholds. 
    Assume k≤ 1/2.
    Consider the outputs of the following calls, where the latter uses the output of the former:
    
    , w    (X,, k)  , 
    
            u    (X,y, , w).

    If < k and u_1≥ n - k 
    then u is (2, )-good for X,y.


    By the guarantees of , <ref>, when <k the weights w give us a bound of  on leverage. That is, for all i ∈w, 
    
    x_i^(X^ W X)^-1 x_i ≤ L.

Furthermore, since only alters w by setting some entries to 0, we have that u_1 = w_1 - w-u_1. Using the assumption that u_1 ≥ n-k and the trivial bound w≤ n, we get that w-u_1≤ k.
    Thus, setting U  =(u), by the first part of †<ref>, we have for all i∈u that
    
    x_i^ (X^ U X)^-1 x_i    ≤[]1+ 2w-u_1·≤ 2,

    where we used the assumption that k ≤ 1/2.
    Since only returns a vector when the largest absolute residual is no greater than , we are done.


Our next claim relates the runs of residual thresholding on adjacent datasets at nearby residual thresholds.
This is the main result about our thresholding procedure.

[Intertwining]
    Let (X,y) and (X',y') be adjacent datasets that differ on index i^*.
    Let k∈ be a discretization parameter.
    Let ,, and ' >0 be any  outlier thresholds such that 
    k≤1/96 and '≥exp108kL.
    Consider the outputs of the following calls:
    
    w,    (X,, k) 
    
            w', '    (X',, k),

    which we feed into:
    
    u    (X,y, , w) 
    
            u'    (X',y', ', w').

    Define I = u∩w'∖i^*.
    If ,' < k and u_1 ≥ n-k then 
    I⊆u'.


As we will see in <ref>, where uses to obtain stable weights, many indices of the weights are easily accounted for.
This includes i^*, which can be handled as a special case, as well as w∖w' and w'∖w whose stability is established by <ref>.
Ignoring those cases for now, we wish to show that any point that is not filtered under (X,y) will also not be filtered under (X',y') provided that the threshold used to filter (X',y') is sufficiently large.
We illustrate these cases in <ref>.
Now we are ready to state the proof of <ref>.









    Our first goal will be to show that w'I is sufficiently good.
    (Recall our notation: w'I∈ [0,1]^n takes the value w_i' for i∈ I and 0 elsewhere.)
    First, we see that u is 2,-good for X, y by noting that < k, u_1 ≥ n-k, and k≤ 1/2 and applying <ref>.
    Next, we show that w'I is close to u in ℓ_1 distance.
    In particular, by definition of , if i ∈u then u_i = w_i. Hence,
    
    *u - w'I_1 
           = ∑_i=1^n *u_i - *w'I_i
       = ∑_i∈u*u_i - *w'I_i + ∑_i∉u*u_i - *w'I_i

If i ∉u, then u_i = 0 and hence *w'I_i = 0 since w'I is by definition only nonzero outside the support of u. Hence, the second term in the last equation is 0. Moving on, by definition of the set I,

    *u - w'I_1     = ∑_i∈u[]w_i - *w'I_i
       = ∑_i∈u[]w_i - *w'w'∖i^*_i
       ≤[]w_i^* - *w'w'∖i^*_i^* + ∑_i∈u[]w_i - *w'w'_i
        = w_i^* - 0 + ∑_i∈uw_i - w'_i
       ≤w_i^* + w - w'_1
       ≤ 3 .

The the last inequality follows from the last part of <ref>. 

    Since I ⊆u, and L ≤ kL≤ 1/12 by assumption, it holds that *u - w'I_1 ≤ 3 ≤ 1/(2L), and we can apply <ref> to show that w'I is 2η_1 L, η_1-good for X, y where η_1 = 1 + 12L.
    Furthermore, w'I is 2η_1 L, η_1-good for X', y' because w'I_i^* = 0.

    Now, we will show that during the execution of (X',y', ', w') we will never discard any i ∈ I.
    Let w'^j denote the weights obtained in the jth iteration of the while-loop out of m total iterations, such that w'^0 = w' and w'^m = u'.
    We proceed to show the loop invariant w'^jI = w'I for j ∈0, …, m.
    Since w' = w'^0, the invariant holds initially.
    
    In each iteration, by the loop invariant we can decompose w'^j-1 as
    
    w'^j-1 = w'^j-1I + w'^j-1I = w'I + w'^j-1I.

    Now, we note that we have u≥u_1 ≥ n-k by assumption and since ' < k, we also have w'≥ n - k by <ref>.
    Therefore by inclusion-exclusion, *I≤ 2k + 1 and so *w'^j-1I≤ 2k + 1.
    Now since 96kL ≤ 1 by assumption, we note that
    
    2[]w'^j-1Iη_1 L≤ 22k+1η_1 L ≤ 22k+11 + 12L L ≤ 12kL ≤1/8

    and for η_2 = 1 + 162k + 1η_1 L,
    
    η_2η_1R 
               = 1 + 162k + 11 + 12L L1 + 12LR
       ≤1 + 96k L1 + 12LR
       ≤exp96k Lexp12LR
       ≤exp108k LR
       ≤ R'.

    Thus by the 2η_1 L, η_1-goodness of w'I and <ref>, if maxy_i - x_i^β_w'^j-1≥' then y_i - x_i^β_w'^j-1∉I and so w'^jI = w'I.
    Finally, it follows that u'I = w'I.
    Therefore, since I ⊆w', we have I ⊆u'.


We now observe that our greedy residual thresholding subroutine only removes more points when run with smaller thresholds.
We now state this simple fact for future reference.

    Fix a dataset (X,y) and starting weights w.
    For outlier thresholds ≤', consider running
    
    u    (X,y,,w) 
    
            u'    (X,y,',w).

    For all i∈ [n], u_i ≤ u_i'.

    


§ ANALYSIS OF 




Inputinput
Requirerequire
dataset X, y; base outlier thresholds _0,_0; weights w; discretization parameter k
∀ j∈ [2k], _j  ()^j ·_0 
j∈{0,…,2k}
    u^(j)(X,y,_j, w) 
    ^(j)min{k, n - u^(j)_1 + j} 

min_j∈{0,…,k}^(j) 
v1/k∑_j=k+1^2k u^(j) 
, v



In this section, we prove the stability guarantees for <ref>, our new regression estimator.
<ref> repeatedly calls <ref>, , over a range of slowly increasing outlier thresholds.
These thresholds are indexed by a number j∈{0,1,…,2k}, where k is a discretization parameter.
(Later, we connect this discretization to the privacy parameters, setting k≈log(1/δ)/.)
The key lemma used in these proofs is <ref>, which relates the weights found on a dataset (X,y) at level j to the weights found on an adjacent dataset (X',y') at level j+1.

We start by showing that the  value and  the weight vector, v, returned by <ref> are low-sensitivity.
[Score is Low-Sensitivity]
    Let (X,y) and (X',y') be adjacent datasets. 
    Fix outlier thresholds , and discretization parameter k.
    Assume kL ≤1/96.
    Let 
    
    _1, w    (X,,k) 
    _1', w'    (X',,k)

    and
    
    _2, v    (X,y,w,,,k) 
    _2', v'    (X',y',w',,,k).

    If _1,_1' < k, then _2-_2'≤ 4.

 
 We observe that all  variables are at most k by construction.
    Without loss of generality, assume _2 ≤_2'.
    First, we consider the case when _2 = k. In this setting, since _2' ≤ k we must have _2 = _2', hence _2 - _2'=0   and we are done.
    
    Now, consider the case when _2 < k.
    Then, by definition of the  algorithm, there must exist a j^* ∈0, …, k such that n - u^j^*_1 + j^* = _2, where u^j^* are the weights returned by the  subroutine (run within ) at outlier threshold R_j^*.
    
    Let u = u^j^* and u' = u'^j^*+1 denote the weights returned by on dataset and outlier thresholds X, y, _j^* and X', y', _j^*+1 respectively.
    Defining I as in <ref>, we note that 
    
    u'_1 ≥*u'I_1 = *u - u - u'I_1 ≥*u_1 - u - u'I_1.

    Now, seeking to bound the last term using <ref>, we note that
    
        
  * _j^*+1/_j^*≥exp108kL by the definition in <ref>,
        
  * kL ≤1/96 by assumption
        
  * Since, u^(j^*)_1 = u_1 and _2 = n - u^(j^*)_1 + j^* < k, it holds that 
    u_1 > n- k + j^* ≥ n - k.
 
    
    Therefore, we can apply <ref>, which implies that u'I = w'I.
    This gives
    
    u - u'I   = *wu - w'I
       = *wu - w'u∖i^*
       ≤*w - w'u + 1
       ≤ 3.

    Finally, combining the previous results gives
    
    _2' 
               ≤ n - *u'^j^*+1_1 + j^*+1
       ≤ n - u^j^*_1 - 3 + j^*+1
       = n - u^j^*_1 + j^* + 4
       = _2 + 4.

The first inequality in the calculation about holds by definition of . The second one uses our previous two calculations.

    


[Weights are Stable]
    Let (X,y) and (X',y') be adjacent datasets. 
    Fix outlier thresholds , and discretization parameter k.
    Assume kL ≤1/96.
    Let 
    
    _1, w    (X,,k) 
    _1', w'    (X',,k)

    and
    
    _2, v    (X,y,w,,,k) 
    _2', v'    (X',y',w',,,k).

    If _1,_1',_2,_2' < k, then v - v'_1≤ 5.


    Consider the execution of resulting in a weight vector u.
    We observe that receives weight vector w as input and modifies it by setting a subset of the weights to zero. 
    Thus we can write the weight u_i = w_i · 1{u_i ≠ 0}.
    This (rather trivial) modification allows us to write the output of in terms of counts: letting c_i = ∑_j=k+1^2k 1{u_i^(j)≠ 0}, we have
    
    v_i = 1/k∑_j=k+1^2k u_i^(j) = w_i/k∑_j=k+1^2k 1{u_i^(j)≠ 0} = w_i c_i/k.

    Now note that by <ref>, for j ∈k+1, …, 2k we can have u_i^j-1≠ 0 only if u_i^j≠ 0.
    This implies that u_i^2k-c_i≠ 0 and  u_i^2k-c_i-1 = 0.
    
    Now, since _2 < k, we know that there exists some j^* ∈0, …, k such that n - *u^j^*_1 + j^* < k. 
    Applying, <ref> again, we see that u^j_1 ≥u^j^*_1 > n-k for all j ≥ j^*.
    From this, we can conclude that i | c_i ≠ k < k.
    
    Now, define c' analogously as the counts under X', y' and note that all of the previous observations apply under X', y' as well.
    Consider some i ∈w∩w'∖i^* and suppose without loss of generality that c_i ≥ c'_i.
    Our goal will be to show that c'_i ≤ c_i ≤ c'_i + 1.
    If we have c_i = k, then c_i = c'_i, so we turn our attention to the case where c_i < k.
    We know that u_i^2k-c_i≠ 0.
    Now, seeking to show that u_i'^2k-c_i+1≠ 0 using <ref>, we note that
    
        
  * u_i^2k-c_i and u_i'^2k-c_i+1 were computed using outlier thresholds _2k-c_i and _2k-c_i+1 which satisfy _2k-c_i+1/_2k-c_i≥exp108kL by the definition in <ref>,
        
  * we have kL ≤1/96 by assumption, and
        
  * we have []u_i^2k-c_i_1 > n - k as we observed previously.
    
    Therefore we can apply <ref>, which implies that u_i'^2k-c_i+1≠ 0.
    Recalling our previous observation, we obtain c_i ≤ c'_i + 1 as desired.
    In summary, if i ∈w∩w'∖i^* then we can write c'_i = c_i + Δ_i where Δ_i≤ 1.

    
    Now, define D = i ∈w∩w'∖i^*|Δ_i ≠ 0 and note that D≤ 2k since c_i and c'_i both contain at most k elements not equal to k.

    Now, we are ready to complete the proof by noting that we can decompose the quantity we wish to bound into four terms
    
    kv - v'_1 
               = *c_i^*w_i^* - c'_i^*w'_i^* + ∑_i ∈w
     i∉w'
     i ≠ i^*c_iw_i + ∑_i ∉w
     i∈w'
     i ≠ i^*c'_iw'_i + ∑_i ∈w
     i∈w'
     i ≠ i^*c_iw_i - c'_iw'_i


    This is valid because c_iw_i - c'_iw'_i appears exactly once on the right for each i.
    Now we will consider each term separately.
    The first term is at most k because c_i, c'_i are bounded by k and w_i, w'_i are bounded by 1.
    The summands of the second and third terms can be rewritten as c_iw_i - w'_i and c'_iw_i - w'_i and are thus bounded by kw_i - w'_i and kw_i - w'_i respectively. 
    Now focusing on the last term, we have
    
    ∑_i ∈w
     i∈w'
     i ≠ i^*c_iw_i - c'_iw'_i   ≤∑_i ∈w
     i∈w'
     i ≠ i^**c_iw_i - c_iw'_i + Δ_iw'_i
       ≤∑_i ∈w
     i∈w'
     i ≠ i^*c_iw_i - c_iw'_i + ∑_i ∈w
     i∈w'
     i ≠ i^*
    i ∈ DΔ_iw'_i
       ≤∑_i ∈w
     i∈w'
     i ≠ i^*kw_i - w'_i + 2k.

    Combining the bounds for each term, we have
    
    kv - v'_1    ≤ 3k + ∑_i ∈w
     i∉w'
     i ≠ i^*kw_i - w'_i + ∑_i ∉w
     i∈w'
     i ≠ i^* kw_i - w'_i + ∑_i ∈w
     i∈w'
     i ≠ i^*kw_i - w'_i
       ≤ 3k + k∑_iw_i - w'_i
       = 3k + kw - w'_1
       ≤ 5k

    where the last line is an application of <ref>.



The previous claim shows that the weights produced by on adjacent datasets are close in ℓ_1 (when  is less than k).
In the next claim, we prove that the weights are good.
Note that this does not follow immediately from <ref>, which says that the weights returned by are good.
Since returns an average of  the vectors returned by , we have to argue that the average of good sets is good.

[Weights are Good]
    Fix a dataset (X,y), outlier thresholds ,, and discretization parameter k.
    Assume k ≤1/4.
    Consider the following calls:
    
    _1, w    (X, , k)
    _2, v    (X,y, w, ,,k).

    Then the vector v is (4,2_2k)-good for (X,y), where _2k= exp(216 k^2 ) ·.


    calls repeatedly, producing a vector u^(j) for each residual threshold _j.
    Recall from <ref> that for all i≤ j, since R_i ≤ R_j (within ) we have that u^(i)≤ u^(j) elementwise.
    This implies that the support of u^(2k) contains all other supports, including that of the average v.
    Additionally, since _2<k (by construction within the algorithm), there exists some j^*∈{0,…,k} with u^(j^*)_1 ≥ n-k, so the same lower bound holds for all j≥ k.
    Together, these facts imply that u^(2k) - u^(j)_1 ≤ k for all j≥ k.
    Since the ℓ_1 norm is convex and v=_j[u^(j)], by Jensen's inequality we have u^(2k)-v_1 ≤ k as well.

    To finish the proof, we apply <ref>: since v⊆u^(2k), u^(2k) - v_1 ≤ k, and u^(2k) is (2, _2k)-good, we conclude that v is (',')-good for
    
    '    ≤ (1 + 4 k)· 2≤ 4
    '    ≤ (1 + 4 k)·_2k≤ 2 R_2k.

    Recalling that _2k= ()^2k·_0, we finish the proof.





§ PRIVACY ANALYSIS OF <REF>



Our privacy analysis follows the blueprint established by <cit.>.
We use the well-known propose-test-release (PTR) framework of <cit.> and first privately check (via our low-sensitivity ) if it is safe to proceed.
If this check passes, we compute a vector of weights v∈ [0,1]^n.
We use this vector to compute a weighted covariance S_v and weighted least squares solution β̂_v.
The output is then drawn from 𝒩(β̂_v, c^2 S_v^-1) for some appropriate constant c.

On adjacent datasets, we may compute different weights v,v'.
We know that, when the PTR checks pass, these vectors are close in ℓ_1. 
The main work in this section, then, lies in connecting this stability of weights to stability of parameters, which in turn implies 𝒩(β̂_v, c^2 S_v^-1)≈_(,δ)𝒩(β̂_v', c^2 S_v'^-1).
Note that this is more complicated than the standard Gaussian mechanism, since both the shape and location of the noise change.

Before proving <ref>, our main privacy claim, we collect the necessary statements.
First, we recall the privacy check of BHS, which (in place of the standard Laplace-noise-and-threshold) simplifies our analysis.
[PTR Mechanism]
    Fix 0<≤ 1, 0< δ≤/10, and 0<Δ.
    There is an algorithm : →, that satisfies the following conditions:
    
        
  * Let  be a set and g:^n →_≥ 0 a function.
            If, for all x,x'∈^n that differ in one entry,  g(x)-g(x')≤Δ, then (g(·)) is (,δ)-DP.
        
  * (0)=.
        
  * For all z≥Δlog 1/δ/+2Δ, (z)=.
    


The next claim relates bounded leverage, ℓ_1 closeness, and covariance closeness.
This statement comes directly from BHS, Lemma 23; similar claims were used in <cit.>.
We use the notation _1, _2 to denote the maximum of  *_1^-1/2_2_1^-1/2- _ and *_2^-1/2_1_2^-1/2- _.
Recall that _1,_2=_1^-1,_2^-1 (<ref>).


    Let ∈ (0,1) and let X,X'∈^n× d be adjacent (i.e., they differ in one out of n rows).
    For vectors v,w∈[0,1]^n, let _v = X^(v) X and _w = (X')^(w) X'.
    Suppose v and w both have bounded leverage: for all i∈v, x_i^_v^-1x_i ≤ and for all j∈w, x_j^_w^-1 x_j≤.
    Then _v and _w are positive definite and, 
    if (1+v-w_1)≤1/2, satisfy
    
    _v, _w≤ 2[]2+v-w_1.





An analogous claim says that, if we have two vectors v and v' that are (,)-good on adjacent datasets and are close in ℓ_1, then the regression parameters they induce are close.
We defer the proof to <ref>, as similar claims appear in the robust statistics literature <cit.>.
claim]claim:parameter_stability_from_weight_stability_different_datasets
    Let (X,y) and (X',y') be datasets differing in one entry.
    Let vector v be (,)-good for (X,y) and let vector w be (,)-good for (X',y').
    Set V=(v) and likewise W.
    Let S_v=X^ V X, β_v = (X^ VX)^-1X^ Vy, and β_w = ((X')^ WX')^-1(X')^ Wy'.
    Assume (v-w_1+2)≤1/4.
    We have S_v^1/2(β_v-β_w)^2 ≤ 4(v-w_1 + 2)^2 ^2.



We use the following relationship between the closeness of covariance matrices and the indistinguishability of their induced Gaussians <cit.>.

    Fix ∈ (0,1) and δ∈ (0,1/10] and let _1, _2 ∈^d× d be positive definite matrices.
    If
        _1,_2≤/3log 2/δ
    then (0, _1)≈_(,δ)(0, _2).






We also need two standard privacy facts: privacy of the Gaussian mechanism and the DP almost-triangle inequality.
[Gaussian Mechanism]
    Fix ,δ∈(0,1) and let u,v be vectors.
    If u-v_2≤Δ, then for any c^2 ≥Δ^2 ·2 log 2/δ/^2 we have  (u,c^2)≈_(,δ)(v,c^2).


[See <cit.>]

    Suppose for some ε and δ that distributions p_1, p_2, and p_3 satisfy p_1 ≈_ε, δ p_2 and p_2 ≈_ε, δ p_3.
    Then p_1 ≈_2ε, 1 + e^εδ p_3.


We are now ready to prove our main privacy claim.

lemma]lemma:main_privacy
    For ε∈0, 1, δ∈0, ε/10, and _0, _0 > 0,
     <ref> is ε, δ-differentially private.


    Consider the execution of <ref> on two adjacent datasets X, y and X', y', yielding _1, _2, v, β̂ and _1', _2', v', β̂' respectively.
    Note that in order to not immediately fail, we must have 
    
    _0 ≤min*1/96k, 3ε/56log 12/δ

    where k = *12 log 3/δ/ε + 8.

    

  
Privacy of the test First, we will show that 
    
    *max_1, _2 - max_1', _2'≤ 4.

    By <ref>, we have _1 - _1'≤ 2.
    Without loss of generality, assume that _1 ≥_1'.
    
    Considering the case where _1 = k, we have max_1, _2 = k and
    
    max_1', _2'≥_1' ≥_1 - 2 ≥ k - 2

    so in this case,
    *max_1, _2 - max_1', _2'≤ 2.
    
    Now if _1 < k then _1' < k as well, so we can apply <ref> to get _2 - _2'≤ 4.
    Then by noting that max is 1-Lipschitz in the ∞-norm, we have
    
    *max_1, _2 - max_1', _2'
       ≤max*_1 - _1', _2 - _2'
       ≤max2, 4
       ≤ 4.

    Finally we see that ε/3δ/3max_1, _2 is ε/3, δ/3-DP by <ref>.
    
    

  
Privacy of the parameter estimate Now we will proceed under the assumption that 
    
    ε/3δ/3max_1, _2 = ε/3δ/3max_1', _2' = ,

    with the goal of showing that β̂, c^2S_v^-1≈_2ε/3, 2δ/3β̂', c^2S_v'^-1.
    Since the PTR checks passed, <ref> says that
    
    _1, _2, _1', _2' < k

    where k = *12 log 3/δ/ε + 8, which matches the assignment in <ref>.
    Now we can apply <ref> to obtain v - v'_1 ≤ 5 and observe that v, v' are both 4_0, 2 exp216 k^2 _0_0-good by <ref>. 
    We will use the stability and goodness of the weights to establish the stability of both β̂ and S_v.
    
    <ref> requires 28L_0 ≤7/24k≤1/4, which is true by assumption.
    The claim implies that  S_v^1/2β̂ - β̂^2 ≤Δ^2 where Δ^2 = 3136exp432 k^2 _0_0_0^2.
    Next, we see that transforming β, β' by S_v^-1/2 allows us to apply <ref>, giving
    
    β̂, c^2S_v^-1≈_ε/3, δ/6β̂', c^2S_v^-1.

     as long as c^2 ≥Δ^2 ·18log 12/δ/ε^2, which is satisfied by construction in <ref>.
    
    Then, since 24_0 ≤ 1/4k≤ 1/2 by assumption, 
    <ref> tells us that
        S_v,S_v'≤ 56_0.
    We apply <ref> and <ref> to obtain
    
    β̂', c^2S_v^-1≈_ε/3, δ/6β̂', c^2S_v'^-1,

    since 56_0 ≤ 3ε/log 12/δ, which we assumed to be true.
    Finally, we apply <ref> to combine the two results, observing that e^ε < e < 3, to complete the proof.





§ UTILITY ANALYSIS OF <REF>

 

Given the privacy guarantee of <ref>, we analyze the utility of <ref> under the standard subgaussian linear model. The definition of subgaussian variables and necessary concentration inequalities are provided in 
<ref>. 
We first note that data from the standard subgaussian linear model is good with high probability.
lemma]lemma:subgaussian_good
    Let X ∈^n × d be drawn i.i.d. from a d-dimensional subgaussian distribution 𝒟 with mean 0, (full-rank) covariance Σ, and subgaussian parameter K_𝒟.
    Let y_i = β^ x_i + z_i where the z_i are drawn i.i.d. from a subgaussian distribution with mean 0, variance σ^2, and subgaussian parameter K_σ.
    There exists constants K_, K_, K_n > 0 such that for any η∈0, 1, if n ≥ K_nK_𝒟^4d + log3/η then X, y is ,-good, where
    
    = K_ K_𝒟^2 ·d + log3n/η/n    and    =K_ K_σσ√(log3n/η),

     with probability at least 1 - η.


Identical calculations about leverage appeared in <cit.>.

Recall that we can write the vector of residuals e=ŷ - y = (H-)z, where H is the hat matrix and z the noise vector. 
Denote by r_i the i-th row of H-.
Then e_i = r_i^ z, which (for any fixed H) implies e_i has subgaussian norm r_i K_σ.
We know that r_i≤ 1, since -H is idempotent and symmetric: 1≥ (-H)_i,i = r_i^ r_i.
Thus all e_i will be bounded with high probability.





As we noted in <ref>, when the input data is good returns the OLS estimate plus noise.
We use a bound on the error of the OLS estimate from prior work.

lemma]lemma:ols_error
Under the distributional assumption of Lemma <ref>, there exists an absolute constant K_OLS such that, for any δ∈ (0, 1), if n>K_OLSK_𝒟(d+log(1/δ)), then with probability 1-δ, we have

    []β̂_OLS-β^2_Σ≤K_OLSK_σ^2σ^2(d+log(1/δ))/n.



Now, we are ready to prove the main accuracy lemma by bounding the norm of the added noise.

lemma]lemma:accuracy_main
    Let X ∈^n × d be drawn i.i.d. from a d-dimensional subgaussian distribution 𝒟 with mean 0, (full-rank) covariance Σ, and subgaussian parameter K_𝒟.
    Let y_i = β^ x_i + z_i where the z_i are drawn i.i.d. from a subgaussian distribution with mean 0, variance σ^2, and subgaussian parameter K_σ.
    There exists constants K_, K_ > 0 such that for any η∈0, 1, if
    
    _0 = K_ K_𝒟^2 ·d + log3n/η/n,        _0=K_ K_σσ√(log3n/η),

    and
    
    n = Ω*K_𝒟^4*d+log*1/εηlog 1/δ^2/ε^2,

    then with probability at least 1 - η <ref> successfully returns β̃ such that
    
    []β̃ - β_Σ≤*K_σσ√(d+log1/η/n) +  K_𝒟K_σσ·d + logn/η√(logn/ηlog1/δ)/ε n,

    where Ω hides log factors in K_D and log1/δ.


    We begin by determining how many samples are needed to ensure that (i) the algorithm does not fail immediately and (ii) the data is (_0,_0)-good (for the specified values) with high probability.
    
    In order to not fail, we require 
    _0 = *ε/log1/δ.
    Meanwhile, in order to apply <ref>, we require _0 = ε/log1/δ^2.
    It is clear that the second requirement implies the first.
    Thus, we can expand our choice of L_0 to get
    
    _0 = K_ K_𝒟^2 ·d + log3n/η/n = *ε^2/log 1/δ^2.

    Using the fact that a/log a = Ωb implies a = Ωb log b, this translates to
    
    n = Ω* K_𝒟^2 log 1/δ^2/ε^2*d+log*K_𝒟^2 log 1/δ^2/ε^2η.

    We note that this implies n = Ω[]K_𝒟d + log1/δ as required by <ref>.
    Now, in order to apply <ref>, we additionally require n = Ω*K_𝒟^4d + log1/η. 
    As in <ref>, this requirement gives us []Σ^-1/2ΣΣ^-1/2 - _2 ≤ 1/2 by <ref>, which we will use later.
    Combining the two outstanding requirements and dropping lower order terms gives
    
    n = Ω*K_𝒟^4*d+log*K_𝒟^2 log 1/δ^2/ε^2ηlog 1/δ^2/ε^2.


    
    This ensures that X, y is _0, _0-good with probability 1 - η.
    When this happens the PTR check passes deterministically.
    Thus, we now turn to evaluating the accuracy of our regression estimate.
    We apply the triangle inequality about :
    
    β - β̃_Σ   ≤β - _Σ +  - β̃_Σ.

    We analyze these terms separately.

    The first term in <ref> is solely about the empirical quantity.
    By <ref>, with probability at least 1 - η we have
    
    - β_Σ = *K_σσ√(d+log1/η/n).


    To bound the second term in <ref>, we apply <ref>, which states that on good data β̃ is drawn from 
        (, c^2 X^ X^-1) where c^2=Θ_0  _0^2 log1/δ/^2.
    Equivalently, we draw z ∼(0,) and set β̃ + c (X^ X)^-1/2 z.
    Plugging this in, we have 
    
    - β̃_Σ   =   c (X^ X)^-1/2 z_Σ
       =  c ·Σ^1/2 (X^ X)^-1/2 z_2.

    We plug in Σ̂= 1/n X^ X the empirical covariance and apply Cauchy–Schwarz:
    
    - β̃_Σ   ≤c/√(n)·Σ^1/2Σ̂^-1/2_2 ·z_2.

    By <ref> the matrix norm is at most a constant, and by <ref> we can bound z_2^2= Od +log 1/η with probability at least 1-O(η).
    Plugging these in, along with our expressions for c, _0, and _0, we arrive at the expression in the lemma.
    Applying a union bound over the three failure cases finishes the proof.





§ RUNNING-TIME ANALYSIS OF <REF>


In this section we prove the following guarantee about the running time of , whose computational requirements are quite lightweight.
The core ideas in this proof appeared in the analogous claim of BHS.

lemma]lemma:main_time
    <ref> can be implemented to require
    
        
  * one product of the form A^⊤ A for A∈^n× d,
        
  * one product of the form A B for A∈ R^n× d and B∈^d× d,
        
  * one inversion of a positive definite matrix in ^d× d;
        and
        
  * further computational overhead of Õ(nd/).
    

Ignoring bit complexity, this corresponds to time Õ(nd^ω-1 + nd/), where ω<2.38 is the matrix multiplication exponent.
For modest privacy parameters, the running time of our algorithm is dominated by the time needed to compute the nonprivate OLS solution itself.

To establish this claim, we provide a second version of , <ref>, which is more computationally efficient.
We show that this alternative algorithm is functionally equivalent.


Inputinput
Requirerequire
dataset X, y; base outlier thresholds _0,_0; weights w; discretization parameter k
∀ j∈ [2k], _j  ()^j ·_0 
COUNT 0 
j∈{2k,2k-1,…,0}
    𝚃𝚁𝚄𝙴
        check for large residuals
        β_w (X,y,w)*via rank-one update
        i^* _i∈w[]y_i - x_i^β_w  
        y_i^*- x_i^*^β_w≤_j or COUNT≥ k
            break*too many outliers or no large residuals
        
        w_i^* 0*otherwise, remove weight
        COUNTCOUNT + 1 
    
    
    COUNT≥ k            
            too many outliers
            ∀ i ≤ j, ^(i) k 
            ∀ i ≤ j, u^(i) 0^n 
            break 
    
    store result and move to next threshold
    u^(j) w 
    ^(j)min{k, n - ‖ u^(j)‖_1 + j} 

min_j∈{0,…,k}^(j) 
v1/k∑_j=k+1^2k u^(j) 
, v
, More Efficient Implementation



    From BHS, Lemma 20 in Section 2.3, we see that we can implement using one product A^⊤ A, one product A B, one matrix inversion, and at most O(log(1/δ)/) additional operations, each of which requires Õ(nd) time.
    We need two additional conclusions from their analysis: can be implemented to return the inverse weighted covariance (X^ W X)^-1 in the same asymptotic running time and we can update all leverage scores in time Õ(nd) when removing a single observation.

    With the weights w and inverse covariance in hand, we call . 
    The initial regression parameter can be computed in Õ(nd) time, as we compute the vector X^⊤ W y with a matrix-vector product (since W is diagonal) and multiply it with the inverse covariance.
    Computing all residuals is linear-time.
    
    Each outlier removal and associated set of updates can also be implemented in Õ(nd) time. 
    This is because the removal of a single point corresponds to a rank-one update, which can be done efficiently.
    Recall from <ref> the equation for updating the least squares solution after removing a data point:
    
    β_ols(-j)   =  + (X^ X)^-1 x_j/1-h_j·[]y_j - x_j.

    (A nearly identical formula applies when the data are weighted.)
    Since we have the previous leverage scores and inverse covariance, this update can be performed in time O(nd).
    As before, with the new regression parameter all the residuals can be recalculated in linear time.

    Setting these details aside, we turn to the crux of the analysis: is functionally equivalent to <ref>, our efficient version.
    
    <ref> iterates through the residual thresholds in decreasing order.
    This is identical to independently calling the greedy algorithm repeatedly from scratch, since the removal process is deterministic (we can break ties in a consistent manner, e.g., using the index of the points). 
    Formally, for any > ' and any fixed X,y,w, the result of (X,y,',w) is identical to calling u(X,y,,w) and then (X,y,',u).

    <ref> also tracks a count of observations it removes and halts if that number reaches k.
    If it halts at level ℓ, for all j≤ℓ it sets u^(j)=0^n and ^(j)=k.
    To show that this has no effect on the outcome of the algorithm, we suppose <ref>'s count reaches k and analyze two cases.
    Let ℓ∈{0,…,2k} be the residual threshold index at which the count k was reached.
    We know that u^(ℓ) returned by (X,y,_ℓ,w) in <ref>, the main version of , satisfies ‖ u^(ℓ)‖_1 ≤ n-k, since after any removal the weight is zero.
    This also holds for all u^(j) with j≤ℓ.

    Case 1: 
    Suppose ℓ>k, i.e., ℓ falls among the indices used to compute the weights.
    Then for all j∈{0,…,k}, the indices used to compute the scores, <ref> computes u^(j) with ‖ u^(j)‖_1 ≤ n-k.
    This means <ref> computes =k, as does <ref> (since it sets _j =k for all j≤ k).
    (Recall that this causes to fail deterministically, so the weights do not impact the output.)

    Case 2: 
    If ℓ≤ k, then ℓ falls among the indices used to compute the score.
    (Thus, <ref> return the same weights.)
    <ref> sets ^(j)=k for all j≤ℓ.
    We claim that <ref> also computes ^(j)=k for all j≤ℓ.
    To see this, recall that on these indices <ref> computes u^(j) with u^(j)_1≤ n-k.

    
    



    
    


    To finish the proof, we note that the final β_v and S_v^-1 computed by can be computed with at most k rank-one updates from their initial values.
    Since k = O(log(1/δ)/), we are done.






§ ACKNOWLEDGEMENTS

This work is supported in part by the National
Science Foundation under grant no. 2019844, 2112471, 2238080, and 2229876, NSF Graduate Research Fellowships Program, the Machine Learning Alliance at MIT CSAIL, and Microsoft Grant for Customer Experience Innovation.
Part of this work was done while GB was at Boston University. AS and GB (while at Boston University) were supported in part by NSF awards CCF-1763786 and CNS-2120667 as well as Faculty Awards from Google and Apple. JCP was supported in part by the Harvard Center for Research on Computation and Society.










§ ADDITIONAL RELATED WORK

 



  
Private Mean Estimation
Many of the developments in private linear regression have analogs in private mean estimation, albeit rearranged chronologically. 
Consider the canonical mean estimation problem with “covariance-adaptive” error guarantees, which respect the shape of the dataset: the error is measured in Mahalanobis norm with respect to the covariance matrix Σ of the data X, Σ^-1/2(μ̂-μ). This scales each direction according to the directional variance, providing a more relevant measure of utility. This is closely related to how linear regression error corresponds to the Σ-norm, β̂-β_Σ. 

For non-private mean estimation, the distinction between Euclidean and Mahalanobis norm error are minor. For example, the empirical mean achieves  such a geometry-aware guarantee and, like the OLS estimator, is accurate with roughly d samples with no dependence on the condition number of the covariance of the data.

For private estimation, geometry-aware estimation is significantly more challenging, since the privately learning the geometry, i.e., the covariance matrix, is more sample expensive than the primary task of mean estimation. At the same time, it seemed like it was necessary to design the privacy noise that matches the shape of the data covariance. The standard Gaussian mechanism privately estimates the mean with only d samples but its error depends polynomially on the condition number, a high price to pay when the estimator does not respect the geometry. 
The work of <cit.> allows us to privately learn the covariance and apply the Gaussian mechanism on whitened data, but, as with the SSP approaches for linear regression, this requires d^3/2 samples. A long line of work that follows either makes more strict assumptions on the geometry or pays a price in the sample complexity <cit.>. 
Of particular relevance for our work are frameworks introduced by <cit.> and <cit.>, which remove outliers in a way that depends on the rest of the dataset (e.g., asking that inliers be close to a large number of other examples). 
These frameworks bear some similarity to our techniques and to those of <cit.>, especially their “Stable Mean” estimator.
Informally, our approach improves over theirs in the ability to adapt the definition of outlier and the resulting geometry as points are removed.

<cit.> was the first to address this geometry-aware challenge in private mean estimation. 
They termed it the covariance estimation bottleneck and gave two exponential-time approaches for avoiding it, achieving accurate estimation with Õ(d) samples and no dependence on the condition number. 
The first, which combined the exponential mechanism with PTR (Propose-Test-Release), served as a direct inspiration to the HPTR (High-dimensional PTR) framework of <cit.>.
The concurrent works of <cit.> and <cit.>  
built on the second algorithm of <cit.>, giving time-efficient algorithms matching the guarantees of the exponential-time approaches. The sample complexity has linear dependence on the dimension d and no dependence on the condition number κ(X^⊤ X). 
As in <cit.> and <cit.>, our goal is to achieve the same for linear regression.   





  
Private Linear Regression.
Commensurate with its centrality in statistical theory and practice,  significant effort has gone into producing differentially private algorithms for least squares <cit.>. 

One standard theme in many of these works is the class of assumptions that directly enable global sensitivity analysis.
Prime examples include assuming that the covariates satisfy an ℓ_2 norm bound or that the true parameter lies in some ball about  the origin.
Such guarantees are incomparable with our definition of goodness (for example, our definition allows arbitrarily large covariates, but covariates with bounded norms may still have high leverage).
Under some collections of these assumptions, state-of-the-art guarantees are achieved in <cit.>, which in our setting translates into a sample complexity of n= Ω(d^1.5/(αε)) to achieve (1/σ) β̂-β_Σ≤α. 
Both these prior algorithms and ours analyze accuracy under the assumption that the input data is “outlier-free." 
The prior work uses conditions on the norm of the covariates or the magnitude of the labels. These assumptions lend themselves handily to global sensitivity calculations. In contrast, our work uses a notion of outlier-freeness which is more in line with standard statistical practice: we ask that the dataset have no high-leverage or high-residual points. 

When applied to data from the standard sub-Gaussian linear models, is the first computationally efficient algorithm to achieve linear dependence in the dimension d and no dependence on the condition number κ(X^⊤ X) (see <ref>). This nearly matches the best known sample complexity of <cit.> that relies on an exponential time approach of HPTR: n=Õ(d/α^2 + (d+log(1/δ))/(αε)) samples suffice to achieve an error of (1/σ)β̂- β_Σ≤α. 
Existing computationally efficient approaches based on gradient descent either assume the covariance matrix is close to identity <cit.> or have polynomial dependence on the condition number κ( X^⊤ X) <cit.>. The best known sample complexity of an efficient algorithm is by <cit.>: n=Õ(d/α^2 + (κ^1/2 dlog(1/δ))/(αε)).  










Iterative Thresholding.
Our algorithm is a special case of the family of iterative thresholding algorithms, a longstanding heuristic for robust linear regression that dates back to Legendre. 
Its theoretical properties in the non-asymptotic regime have been extensively studied recently in <cit.>. <cit.> and <cit.> studied the iterative trimmed estimator under generalized linear models and <cit.> studied the mixed linear regression setting. It is worth noting that most iterative thresholding algorithm in the robust linear regression setting will alternate between finding the OLS solution of the current set and finding the set with the smallest residual under the current regression coefficient, and no data point is permanently removed in each iteration. In contrary, our algorithm will permanently remove one data point in each iteration before recomputing the OLS solution.



§ PRELIMINARIES


We collect here known preliminary results that we use in our analyses. 




    Let _1, _2 be positive-definite matrices and define 
    
    _1, _2 =  max**_1^-1/2_2_1^-1/2- _, *_2^-1/2_1_2^-1/2- _.

    Then _1, _2 = _1^-1, _2^-1.


    Note that _1^-1/2_2_1^-1/2 and _2^1/2_1^-1_2^1/2 are similar and likewise, _2^-1/2_1_2^-1/2 and _1^1/2_2^-1_1^1/2 are similar.
    Thus,
    
    (_1, _2) 
               = max**_1^-1/2_2_1^-1/2- _, *_2^-1/2_1_2^-1/2- _
       = max**_2^1/2_1^-1_2^1/2- _, *_1^1/2_2^-1_1^1/2- _
       = _2^-1, _1^-1
       = _1^-1, _2^-1

    as desired.




 §.§ Subgaussian Random Variables and Concentration Inequalities


For formal proofs of claims and further discussion we refer to <cit.>.


    Let ∈ be a random variable. 
    The subgaussian norm of , denoted _ψ_2, is 
        _ψ_2 = inft>0 : exp^2/t^2≤ 2.


definition]def:subgaussian_vector
    Let ∈^d be a random variable with mean μ and covariance Σ. 
    Call  subgaussian with parameter K if there exists K≥ 1 such that for all v∈^d we have
    
    -μv_ψ_2  ≤   K √( v^Σ v ).


For example, the Gaussian distribution (μ,Σ) is subgaussian with parameter K=O(1).

[Concentration of Norm]
    Let _1,…,_n be drawn i.i.d. from a d-dimensional subgaussian distribution with parameter K_>0, mean μ, and (full-rank) covariance Σ.
    There exists a constant K_1>0 such that, with probability at least 1-β, we have both
    
    *Σ^-1/2_1 - μ^2 ≤ K_1 K_^2 d + log 1/β  and  *Σ^-1/2*1/n∑_i=1^n _i - μ^2 ≤ K_1 K_^2·d + log 1/β/n.



[Concentration of Covariance]
    Let _1,…, _n be drawn i.i.d. from a d-dimensional subgaussian distribution with parameter K_>0, mean μ=0, and (full-rank) covariance Σ.
    Let Σ̂= 1/n∑_i=1^n _i_i^ be the empirical covariance.
    There exist positive absolute constants K_1 and K_2 such that, for any β∈ (0,1), if n≥ K_2d + log 1/β, then
    with probability at least 1-β we have
    
    Σ^-1/2Σ̂Σ^-1/2 - _2   ≤   K_1 K_^2√(d + log 1/β/n).






 §.§ Details on 



As a preprocessing step, performs a leverage-score filtering routine introduced by BHS.
The algorithm we use differs only superficially from their version.
(For instance, they compute a set of weights w ∈ [0,1]^n and a weighted covariance estimate, while we only care about the weights themselves.)
For completeness, we now state the version we use here.
Recall that <ref> contains the relevant guarantees proved by BHS.
[h]
InputInput
dataset X∈^n× d; outlier threshold _0; discretization parameter k∈
A  [n] 
j = 2k,2k-1,…,0
     _j = exp{j /k}·_0 
     𝙾𝚄𝚃 = ∅
        S_A ∑_i∈ A x_i x_i^ 
        𝙾𝚄𝚃i ∈ A : x_i^S_A^-1 x_i > _j 
        A  A ∖𝙾𝚄𝚃
    
    A_j  A 

min{k, min_0≤ j ≤ k{n - A_j + j}} 
i=1,…,n
    w_i 1/k∑_j=k+1^2ki∈ A_j 

, w 

Stable Leverage Filtering (), BHS









§ DEFERRED PROOFS









We now give the proof for <ref>, which characterizes the effect of removing weighted points from a least squares model.
This is a natural generalization of standard results <cit.>.

We start by setting up notation and bounding a term useful in proving both (<ref>) and (<ref>).

    

  
Setup
    Assume without loss of generality that w= [n], as any points outside the support of w are irrelevant. 
    Let v = w' - w, 
     v_1 = ρ, 
    W=(w) (likewise W' and V), and C = X^ W X.
    Decompose V = P N (for “positive” and “negative”) where P,N are diagonal matrices with P_i,i = √(v_i) and N_i,i = sign(v_i)·√(v_i).


    Let Δ = I + Y where Y = NX C^-1 X^ P. If Y_2≤<1 then I+ Y ≽(1- )I. Consequently, Δ is invertible and Δ^-1_2 ≤ 1 / (1-). To prove that Y_2 ≤, we use the fact that Y_2 ≤Y_F and compute:
    
    NX C^-1 X^ P_F^2    = ∑_i,j N_i,i P_j,j· x_i^ C^-1 x_j  ^2 
       ≤∑_i,j N_i,i^2 P_j,j^2 ^2 
       = ^2 ∑_i,jv_iv_j
       = ^2 ∑_i v_i∑_j v_j
       = ^2 ·w - w'_1^2.

    Ultimately, we arrive at Δ^-1_2 ≤ (1- w-w'_1 )^-1. In the first line of the calculation above, we used the fact that,  by our initial assumption on the goodness of (X,w),

    x_i^ C^-1 x_j_2 ≤C^-1/2 x_i_2 C^-1/2 x_j_2  ≤ L^2.

Since, w-w'_1 ≤ 1/2, we conclude that Δ^-1_2≤ 2.

    

  
Bounding the Leverage Scores
    Consider an index i∈w 
    Write out its leverage under w' and plug in our notation:
    
    h_i'    = x_i^X^ W' X^-1 x_i 
       = x_i^X^ W X + X^ V X^-1 x_i 
       = x_i^C + (PX)^ NX^-1 x_i.

    Recalling the fact that, x_i^ C^-1 x_i = h_i and Δ = I + NXC^-1X^ P, we can apply the Woodbury matrix identity to arrive at the following relationship between h_i' and h_i: 
    
    h_i'    = x_i^[C^-1 - C^-1 (PX)^I + NX C^-1 (PX)^^-1 NX C^-1]x_i 
       = h_i - x_i^ C^-1 X^ P Δ^-1 NX C^-1x_i.

    We want to upper bound this leverage, h'_i, so we take the absolute value of the right-hand term and apply Cauchy–Schwarz:
    
    x_i^ C^-1 X^ P Δ^-1 NX C^-1x_i   ≤Δ^-1_2 ·P X C^-1 x_i_2 ·NX C^-1 x_i_2.

    We already argued that Δ^-1_2≤ 2, so we turn to the second term in the product:
    
    PX C^-1 x_i_2^2    = ∑_j=1^n P_j,j^2 (x_j^ C^-1 x_i)^2 
       ≤∑_j=1^n v_i·^2 = w-w'_1 ·^2.

    Note that an identical bound also holds for NXC^-1 x_i_2, hence we have
    
    h_i'    ≤ h_i + 2 w-w'_1 ^2 
       ≤+ 2 w-w'_1 ^2 = 1 + 2 w-w'_1·.

    
    
    

  
Bounding the Residuals
    As above, we use the Woodbury matrix identity to derive an expression for the regression line β_w'.
    
    β_w'   = X^ W' X^-1 X^ W' y 
       = X^ W X + X^ V X^-1X^ W y + X^ V y
       = C + (PX)^ NX^-1X^ W y + X^ V y
       = C^-1 - C^-1 (PX)^Δ^-1 NX C^-1X^ W y + X^ V y.

    This expression expands into four terms, which simplify nicely:
    
    β_w'   = β_w 
          + C^-1 X^ V y 
          - C^-1 X^ P Δ^-1 N X C^-1 X^ W y 
          - C^-1 X^ P Δ^-1 N X C^-1 X^ V y 
       = β_w 
          + C^-1 X^ P Δ^-1 N (P Δ^-1N)^-1 V y 
          - C^-1 X^ P Δ^-1 N X β_w 
          - C^-1 X^ P Δ^-1 N X C^-1 X^ V y 
       =β_w   + C^-1 X^ P Δ^-1 N N^-1Δ P^-1 V y - X β_w - X C^-1 X^ V y.

    We can further simplify the expression in the parentheses. In particular, by plugging in the definition of Δ, we have that N^-1Δ P^-1 V y - X β_w - X C^-1 X^ V y can be rewritten as:
    
    N^-1I + NX C^-1 X^ P P^-1 V y - X β_w - X C^-1 X^ V y 
       = N^-1 P^-1Vy + N^-1 N X C^-1 X^ P P^-1 Vy 
               - X β_w -  X C^-1 X^ Vy 
       = V^-1 V y + XC^-1 X^ Vy  - X β_w - X C^-1 X^ Vy 
       = y - Xβ_w.

    Plugging back in, we arrive at 
    
    β_w'   = β_w + C^-1 X^ P Δ^-1 N y - Xβ_w.

To finish the proof, we consider the residual on a point x_i∈w:
    
    []y_i - x_i^β_w'   = []y_i - x_i^β_w + x_i^β_w - x_i^β_w'
       ≤[]y_i - x_i^β_w + []x_i^β_w - x_i^β_w'
       ≤ + []x_i^ C^-1 X^ P Δ^-1 N y - Xβ_w,

    The bound []y_i - x_i^β_w follows from our initial assumption that w is (L,R)-good for (X,y).
    We can bound the second term in an analogous way to how we bounded the leverage scores. In particular, using the fact that N_2 ≤√(w-w'_1),
    
    []x_i^ C^-1 X^ P Δ^-1 N y - Xβ_w   ≤*P X C^-1 x_i_2 ·*Δ^-1_2 ·*N (y - X β_w)_2 
       ≤ 2 w - w'_1.

    This completes the proof. 

    












    


    


    


    

    









We now prove <ref>, which says that, if two weight vectors on adjacent datasets are both good and close in total variation distance, then their least-squares solutions are close as well.
We start by considering the setting where the vectors correspond to the same dataset.

    Assume v,w are both (,)-good for dataset (X,y).
    Define _v = X^diag(v) X, β_v = S_v^-1X^ V y, and likewise _w, β_w.
    If v-w_1≤1/4, then
        _v^1/2(β_v- β_w)^2 ≤ 4 v -w _1^2 ^2.


From this claim, <ref> is an easy corollary.

    Vectors v,w∈[0,1]^n on adjacent datasets (X,y) and (X',y'), respectively, correspond to vectors v',w'∈ [0,1]^n+1 over the datasets' union. 
    We have v'-w'_1≤v-w_1 + 2.


We start by expanding out the squared norm and substituting in the definition of _v:
    
    []_v^1/2(β_v- β_w)^2
                   = []β_v -β_w, _v (β_v- β_w)
       = *β_v -β_w, ∑_i∈ [n] v_i · x_i x_i^ (β_v- β_w).

    Next we expand the sum across β_v -β_w and add and subtract ∑_i v_i · x_i y_i, which makes the right-hand side of the inner product look like a pair of gradients.
    
    []_v^1/2(β_v- β_w)^2
                   = *β_v -β_w, ∑_iv_i · x_i x_i, β_v -∑_iv_i · x_i x_i, β_w
       = *β_v -β_w, ∑_i v_i · x_i []x_i, β_v -y_i - ∑_iv_i · x_i []x_i, β_w - y_i.

    By definition, β_v is the vector that sets the first gradient to zero, so we have
    
    []_v^1/2(β_v- β_w)^2
                   = *β_v -β_w, 0 - ∑_iv_i · x_i []x_i, β_w - y_i.

    We now add and subtract the gradient at β_w weighted by w, which leaves a gradient term (also zero by definition) and the differences v_i - w_i:
    
    []_v^1/2(β_v- β_w)^2
                   = *β_v -β_w, - ∑_iw_i · x_i []x_i, β_w - y_i + ∑_i(w_i-v_i) · x_i []x_i, β_w - y_i
       = *β_v -β_w, 0 - ∑_i(w_i-v_i) · x_i []x_i, β_w - y_i
       = ∑_i[]β_v -β_w, (w_i-v_i) · x_i []x_i, β_w - y_i.

    We now insert _v^1/2_v^-1/2 in the middle of each inner product.
    We apply Cauchy–Schwarz to each term and pull out the scalars (recall that w_i and v_i are scalars, x_i is a vector):
    
    []_v^1/2(β_v- β_w)^2
                   = ∑_i*_v^1/2(β_v -β_w), (w_i-v_i) ·_v^-1/2 x_i []x_i, β_w - y_i
       ≤∑_i[]_v^1/2(β_v -β_w)·[](w_i-v_i) ·_v^-1/2x_i []x_i, β_w - y_i
       = ∑_iv_i - w_i·[]_v^1/2(β_v -β_w)·[]_v^-1/2x_i ·x_i,β_w-y_i .

    Both sides of the equation have a []_v^1/2(β_v -β_w) term; these cancel.
    We apply ℓ_1/ℓ_∞ on the weight differences, leverage scores, and residuals.
    There is some subtlety here: define set U = v∪w.
    Then we have
    
    []_v^1/2(β_v- β_w)   ≤∑_iv_i - w_i·_v^-1/2 x_i·[]x_i,β_w-y_i
       ≤v - w_1 ·max_i∈ U_v^-1/2x_i·[]x_i,β_w - y_i.

    By the definition of goodness, if i∈v then we have _v^-1/2 x_i≤√().
    Similarly, if i∈v, we have y_i - x_i^β_w≤.
    
    However, these bounds may not hold for points outside the relevant support.
    <ref> allows us to bound the leverage: for all i∈w∪v we have _v^-1/2x_i_2^2 ≤ 2, since we assumed (1+v-w_1)≤1/2.
    
    Similarly, bounding the residual involves a simple trick alongside <ref>.
    Let w̌ be the entry-wise minimum of {w,v}, so w̌_i = min{w_i, v_i}.
    We have w̌ - w_1, w̌-v_1≤w-v_1 and, furthermore, the support of w̌ is contained in both the support of w and that of v.
    Thus, we can apply <ref>: assuming i∈w (since otherwise i∈v and we have a bound on the residual)
    
    y_i - x_i^β_v   = y_i - x_i^β_w + x_i^β_w - x_i^β_w̌ + x_i^β_w̌ - x_i^β_v
       ≤y_i - x_i^β_w + x_i^β_w - x_i^β_w̌ + x_i^β_w̌ - x_i^β_v
       ≤ + 2 w̌ - w_1  + 2 w̌ - v_1 
       ≤[]1 + 4 v - w_1 .

    Since v-w_1≤1/4, this is at most 2.





§ ESTIMATION OF Σ^2


[H]    
   Private σ^2 Estimator 
    
   	S={(x_i, y_i)}_i=1^n, target privacy (ε_0,δ_0), target failure probability ζ.
	Fn:
	 
	Partition S into k=⌊ C_1log(1/(δ_0ζ))/ε⌋ subsets of equal size and let G_j be the j-th partition, where each dataset is of size b=|G_j|=⌊ n/k⌋.

	
	
	
	For each j∈ [k], denote ψ_j= min_β(1/|G_j|) ∑_i∈ G_j(y_i-β^⊤ x_i)^2.
	
	Partition [0, ∞) into bins of geometrically increasing intervals Ω:= {…,[ 2^-2/4, 2^-1 / 4),[ 2^-1 / 4, 1),[1,2^1 / 4),[2^1 / 4, 2^2/4), …}∪{[0,0]}

	Run (ε_0, δ_0)-DP histogram learner of Lemma <ref> on {ψ_j}_j=1^k over Ω 

	if all the bins are empty then 
	Return ⊥

	Let [ℓ, r] be a non-empty bin that contains the maximum number of points in the DP histogram

	Return  ℓ
	 





	Algorithm <ref> is (ε_0, δ_0)-DP. Let S={(x_i, y_i)}_i=1^n be a dataset of i.i.d. samples with  x_i ∼(0, Σ), y_i=x_i^⊤β^*+z_i and z_i∼(0, σ^2)  
for some unknown true parameter β^*=Σ^-1𝔼[y_i x_i] ∈ℝ^d and unknown Σ and σ^2.  Suppose 
	
    n=Õ(dlog(1/(δ_0ζ))/ε_0) ,

	with a large enough constant then 
Algorithm <ref> returns ℓ such that, with probability 1-ζ, 

    1/√(2)σ^2≤ℓ≤√(2)σ^2 ,

where Õ(·) hides logarithmic factors in log(1/ε_0),log(1/δ_0).

We provide a proof in Appendix. <ref>.




 §.§ Proof of Lemma <ref> on the private σ^2 estimation


The privacy proof follows from the DP-histogram from Lemma <ref>. We provide proof for utility.

For each of the partition G_j, we show a_i:=1/b∑_i∈ G_j(y_i-β̂_j^⊤ x_i)^2 concentrates around the true parameter β^* where β̂_j:=_β(1/|G_j|) ∑_i∈ G_j(y_i-β^⊤ x_i)^2. Let f(β)=1/b∑_i∈ G_j(y_i-β^⊤ x_i)^2. We know f(β̂_j) = min_β1/b∑_i∈ G_j(y_i-β^⊤ x_i)^2 ≤1/b∑_i∈ G_j(y_i-β^*⊤ x_i)^2=1/b∑_i∈ G_j z_i^2

Since z_i^2 are sub-exponential, from Bernstein bound, we know  there exists constant c_1>0 such that with proability 1-ζ,
1/b∑_i=1^b z_i^2≤σ^2(1+c√(log(1/ζ)/b)+clog(1/ζ)/b).

Now we show lower bound of f(β̂_j). For any β, we also have

    f(β) = 1/b∑_i∈ G_j(y_i-w^⊤ x_i)^2= 1/b∑_i∈ G_j(z_i+x_i^⊤(w^*-w))^2 .


Let β̃:=(Σ^1 / 2(β^*-β), σ) ∈ℝ^d+1 and x̃_i:=(Σ^-1 / 2 x_i, z_i / σ) ∈ℝ^d+1 for i ∈ [n]. By definition, we can see that x̃_i is zero-mean sub-Gaussian with covariance 𝐈_d+1.


    f(β) = 1/b∑_i∈ G_j(β̃^⊤x̃_i)^2.


Following Lemma 9 from <cit.>, we know for any vector β̃, there exists c_2>0 such that with probability 1-ζ, 


    |1/b∑_i∈ G_j(β̃^⊤x̃_i)^2-β̃^2| =  |β̃^⊤(1/b∑_i∈ G_jx̃_ix̃_i^⊤-𝐈_d) β̃|
    	 ≤ c_2√(d+1+log(1/ζ)/b)+c_2d+1+log(1/ζ)/b


This means for any w, we have


    f(β) = 1/b∑_i∈ G_j^b(β̃^⊤x̃_i)^2   ≥ (1-c_2√(d+1+log(1/ζ)/b)-c_2d+1+log(1/ζ)/b) (Σ^1/2(β-β^*)^2+σ^2) 
       ≥ (1-c_1√(d+1+log(1/ζ)/b)-c_2d+1+log(1/ζ)/b) σ^2 .


Together with the upper bound, this implies that there exists constant c_3>0, such that with probability 1-ζ,


    |1/b∑_i∈ G_j a_i-σ^2|=|f(β̂_j)-σ^2|≤ c_3(√(d+1+log(1/ζ)/b)+d+1+log(1/ζ)/b)σ^2 .


By union bound, there exists a constant c_4>0 such that if b≥ c_4(d+log(k/ζ)), then for all j∈ [k],


    |ψ_j-σ^2| = |1/b∑_i∈ G_j a_i-σ^2|≤ 2^1/8σ^2 .
 

With probability 1-ζ, {ψ_j}_j=1^k  lie in interval of size 2^1/4σ^2. Thus, at most two consecutive bins are filled with {ψ_j}_j=1^k. Denote them as I=I_1∪ I_2. 


 Our analysis indicates that ℙ(ψ_i∈ I)≥ 0.99. By private histogram in Lemma <ref>, if k≥ c_5log(1/(δ_0ζ))/ε_0, |p̂_I-p̃_I|≤ 0.01 where p̂_I is the empirical count on I and p̃_I is the noisy count on I. Under this condition, one of these two intervals are released. This results in multiplicative error of √(2). 	


 For every K∈ℕ∪{∞}, domain Ω, for every collection of disjoint bins B_1,…, B_K defined on Ω, n∈ℕ, ≥ 0,δ∈(0,1/n), β>0 and α∈ (0,1) there exists an (,δ)-differentially private algorithm M:Ω^n→ℝ^K such that for any set of data X_1,…,X_n∈Ω^n


  * p̂_k = 1/n∑_X_i∈ B_k1

  * (p̃_1,…,p̃_K) M(X_1,…,X_n), and

  * 
    n≥min{8/βlog(2K/α),8/βlog(4/(αδ))}


then,

    ℙ(|p̃_k-p̂_k|≤β)≥ 1-α






§ LOWER BOUND



Series of advances have been made in designing tools for lower bounds in statistical estimation.

Fingerprinting <cit.>. 

Our lower bound is a direct corollary of a similar lower bound on linear regression from 
<cit.>. 


 
 	Let _Σ,  σ^2 be a class of distributions over (x_i, z_i)∈ℝ^d ×ℝ, where x_i are i.i.d. samples from a d-dimensional subgaussian distribution with mean 0 
        and covariance Σ≻ 0, and z_i are i.i.d. samples from a subgaussian distribution with mean 0 and variance σ^2
        (see <ref> in <ref>). We observe labelled examples from linear model:  y_i = β^ x_i + z_i with 𝔼[x_i z_i]=0.
        Let _ε,δ be a class of (ε, δ)-DP estimators that are functions over the datasets S={(x_i, y_i)}_i=1^n. Then if 0<ε<1, d≲ nε, δ≲ n^-(1+γ) for some γ>0, there exists constant C>0 such that 
        
    inf _M ∈ℳ_ε, δsup__Σ,  σ^2, β𝔼M(y, x)-β_Σ^2 ≥ C σ^2(d/n+d^2/n^2 ε^2).

 
 
 
 
 	We will apply the lower bound below from <cit.>.
 	
 	
Consider i.i.d. observations {(y_1, x_1), ⋯,(y_n, x_n)} drawn from the Gaussian linear model:

    f_β(y | x)=1/√(2 π)σexp(-(y-x^⊤β)^2/2 σ^2) ; x ∼ f_x .

Suppose 𝔼[xx^⊤] is diagonal, and λ_max(𝔼[xx^⊤])< C'<∞,  X_2 ≲√(d) almost surely. If d  ≲ n ε, 0<ε<1 and δ≲ n^-(1+γ) for some γ>0, then

    inf _M ∈ℳ_ε, δsup _β∈ℝ^d𝔼M(y, x)-β_2^2 ≳σ^2(d/n+d^2/n^2 ε^2) .

 
Note that this lower bound holds for every construction of x_i that satisfies the assumption. We construct one instance of joint distribution P∈_Σ, σ^2 such that it also satisfies the assumptions in Theorem <ref>. Let {x_i}_i=1^n be i.i.d. samples from (0, 𝐈_d). And let x̃_i=x_i·𝐈[x_i≤√(d)].  Clearly, {x̃_i}_i=1^n satisfies that 𝔼[x̃x̃^⊤] is diagonal, λ_max(𝔼[x̃x̃^⊤])< 1 and x̃_i are bounded by √(d). Let z_i be independent Gaussian distribution with variance σ^2. By Theorem <ref>, we know there exists constant C such that

    inf _M ∈ℳ_ε, δsup__Σ,  σ^2, β𝔼M(y, x)-β_Σ^2 ≥ C σ^2(d/n+d^2/n^2 ε^2).

 	
 


