⁶

1 Introduction↩︎

We study regression tasks in settings that deviate from the ideal i.i.d. assumption, motivated by the prevalence of imperfect or incomplete data in modern machine learning. Missing values can arise from sensor failures, survey nonresponse, or data aggregation processes such as crowdsourcing [1] and peer grading [2], [3]. Moreover, datasets may include corrupted entries, which may stem from data poisoning attacks [4]–[8], out-of-distribution examples [9], or biological anomalies [10]–[12]. These two types of impurities (incompleteness and corruption) have given rise to two largely distinct lines of research: learning with missing data and robust statistics.

1.0.0.1 Incomplete data

Rubin’s seminal work [13] classifies missing data mechanisms based on their dependence on observed values: if missingness is independent of the data, it is Missing Completely At Random (MCAR); if it depends only on observed values, it is Missing At Random (MAR); otherwise, it is Missing Not At Random (MNAR). A common approach to handling missing data is imputation of missing entries [14]–[16], though not all methods rely on it. Early work in the area focused on parameter estimation [17]–[20], whereas more recent research has shifted toward prediction [14], [21], [22], a more challenging task since missing features during test time can prevent model evaluation.

1.0.0.2 Corrupted data

The robust statistics, initiated in the 1960s [23], [24], focuses on parameter estimation from fully observed data with a small fraction of arbitrary corruptions. A recent resurgence in this area [25], [26] has extended classical one-dimensional results to high dimensions while preserving sample and computational efficiency (see [27]). While this literature can also handle missing data by naïvely imputing erased locations with arbitrary values, this simple strategy ignores the fact that the erasure locations are known and can potentially be leveraged by a more sophisticated algorithm. Additional related work on missing data and robust statistics is given in 5.5.

In this work, we study the relationship between these two types of data impurities in the context of parameter estimation in Gaussian linear regression.

We model data impurities via a coordinate-wise adversary: for each coordinate⁷, at most an \(\eta\)-fraction of samples may be modified by erasure or replacement. Unlike much of the robust statistics literature—which treats each sample as fully clean or corrupted—partial corruptions are common, especially when coordinates correspond to different sensors or measurement methods [28]. [def:model] is a strong MNAR variant, as missing-entry locations are chosen adversarially based on the dataset.

The most relevant prior work on [def:model,def:model2] is [29], which focuses on mean estimation. While linear regression has been studied under various missing data settings [19], [30], [31], existing work typically assumes much milder missingness than the fully adversarial pattern in [def:model]. To the best of our knowledge, even for the most “textbook” and standard setting of Gaussian covariates used in [def:model], linear regression under the aforementioned corruption models remains unexplored, and our understanding is limited even from an information-theoretic standpoint. Specifically, we ask:

Due to the adversarial nature of [def:model], the error need not vanish as the number of samples grows. Instead, the best achievable error is a positive function of \(\eta, \sigma, d, \beta\) that we seek to characterize. We provide matching upper and lower bounds, with the upper bounds attained by algorithms that are both sample and computationally efficient.

Another question arises from comparing [def:model,def:model2]. When coordinates are missing rather than replaced, their locations are visible and may be leveraged for improved error, thus:

Prior algorithms for related problems with missing data [29], [32], like mean estimation under [def:model] and low-rank structure [29], rely heavily on knowing corruption locations and only handle special cases under replacement corruption. Surprisingly, for linear regression, the answer to the above question is negative.

1.1 Our Results↩︎

We now present the upper and lower bounds that address the two stated questions. The upper bounds follow from existing computationally-efficient estimators, while the lower bounds are the main contributions of this work.

We start with three estimators, all applicable to the stronger corruption model of [def:model2] (i.e., replaced data) and efficient in runtime and sample complexity. The first, \(\mathcal{A}_1\), is simply using an estimator from existing robust statistics literature, designed for the setting where \(\varepsilon\)-samples are (entirely) corrupted with \(\varepsilon=\eta d\). The second, \(\mathcal{A}_2\), is based on the fact that under 1, \(\operatorname*{\mathbb{E}}[yX] = \beta\), thus one can estimate each coordinate of \(\beta\) using existing 1D robust Gaussian mean estimators tolerant to \(2\eta\) corruption (like trimmed mean or median). The third, \(\mathcal{A}_3\), is the trivial algorithm returning always the zero vector. See 5 for details on these algorithms. As a note, \(\mathcal{A}_1,\mathcal{A}_2,\mathcal{A}_3\) are known to generalize naturally to other covariate distributions beyond Gaussians, but for concreteness and for clear comparison with our lower bounds, we focus on the most basic setting stated below.

Depending on \(\eta\), \(\sigma\), and \(\|\beta\|\), a different algorithm among \(\mathcal{A}_1,\mathcal{A}_2,\mathcal{A}_3\) achieves the lowest error. [tab:eta-main,table:medium-eta,table:small-eta] in 5 present this by partitioning the parameter space into regimes based on which algorithm’s error is best (5.2 has the explicit construction of the partition). In 5.3, we show that \(\mathcal{A}_1\), \(\mathcal{A}_2\), and \(\mathcal{A}_3\) can be unified into a single algorithm that adaptively switches between the three strategies to consistently get the best error (up to constant factors).

We now present the lower bounds (the main results of this paper) in a combined statement in 1. Since some parameter regimes in 1 overlap, we show only the strongest lower bound in 1. A more expanded version is given in [tab:eta-main,table:medium-eta,table:small-eta] which might be more convenient to look at, but due to space constraints are moved to 5. 1 provides information-theoretic lower bounds, applying to any estimator, independent of runtime and for any sufficiently large sample size.⁸

A few comments are in order regarding the theorem statement. First, the failure probability in the theorem statement is close to a half (instead of one) because the lower bound is obtained from a reduction to a hypothesis testing problem that asks for distinguishing between two parameter vectors \(\beta, \beta'\) (see 2 for more details), thus a trivial algorithm that outputs each of \(\beta, \beta',\) with probability \(1/2\) succeeds with the same probability. Second, parts [it:thm-big-eta] and [it:thm-interm-eta] hold under arbitrary additive noise, not necessarily Gaussian. This is because the conclusion as stated in the theorem holds in the noiseless case (\(\sigma=0\)), thus by a straightforward reduction this implies hardness under any noise distribution: given an estimator for the noisy model, one could add the same noise to noiseless hard instance and then apply the estimator. Finally, in terms of the covariate distribution, 1 shows that the problem is hard even for the most restricted and basic setting of Gaussian covariates; the hardness trivially extends to any family of distributions that contains Gaussians.

We now discuss the conceptual takeaways and how 1 addresses our main research questions.

1.1.0.1 Missing vs. Replaced Data

The upper bounds hold under the stronger model of data replacement ([def:model2]), while the lower bounds assume the weaker model of missing data ([def:model]). Since the bounds match (up to constants) for essentially the entire range of \(\eta\), \(\sigma\), and \(\|\beta\|\), both models are equally hard in terms of optimal estimation error. While this may be unsurprising for the robust statistics framework where samples are either entirely corrupted or uncorrupted, our more delicate lower bound constructions show that the fact remains true even under coordinate-wise erasure/corruption.

1.1.0.2 Corruption Threshold for Non Trivial Estimation

In the classical robust statistics model with \(\varepsilon\)-fraction of (entirely) corrupted samples, estimation becomes impossible at \(\varepsilon= 1/2\), since half the samples could come from two different models and there is no way to tell which are the inliers. 1[it:thm-big-eta] shows that the analogous transition point for [def:model,def:model2] is \(\eta = \Theta(1/\sqrt{d})\). This is interesting, as there is no equally obvious explanation for this threshold in this model.

1.1.0.3 \(\|\beta\|\)-Dependence on Error is Necessary

Ignoring the regime \(\eta \geq 1/\sqrt{d}\), where meaningful estimation is impossible, our results show that \(\|\beta\|\)-dependence in the error is unavoidable in general —in stark contrast to the classical robust model (with each sample either fully clean or corrupted), where optimal error is independent of \(\|\beta\|\). Here is the intuition: in that setting, a variant of \(\mathcal{A}_2\) gives error \(O(\varepsilon(\|\beta\| + \sigma))\) (with \(\varepsilon\) corruption rate). But by shifting labels via \(y \gets y - \widehat{\beta}^\top X\), one reduces the regression target to \(\beta - \widehat{\beta}\), which has smaller norm, and can recursively apply the algorithm to eliminate the \(\|\beta\|\)-dependence. This technique fails under [def:model2]: even one corrupted coordinate in \(X\) corrupts \(\widehat{\beta}^\top X\), and the transformation \(y \gets y - \widehat{\beta}^\top X\) can amplify the label corruption rate from \(\eta\) to \(\eta d\) in a single step. Our results show that no workaround exists.

1.1.0.4 Role of Label Corruptions

As can be inferred from the proofs, all of our lower bounds hold even if \(y\) is never masked. Thus, the difficulty in estimating \(\beta\) comes only from the missing \(X\)-coordinates.

1.1.0.5 Additional Discussion

1[it:thm-big-eta] may extend to other activation functions. Since 1[it:thm-big-eta] and [it:thm-interm-eta] hold for \(\sigma=0\), they imply that for any estimator, the mean squared prediction error can be bounded below for a broad class of GLMs, such as a class of GLMs with non-decreasing Lipschitz activations, including ReLUs in particular. See 5.4 for the relevant discussion.

1.2 Overview of Techniques↩︎

The general approach for our lower bounds is the following. We first reduce estimation to a hypothesis testing problem. If no algorithm can distinguish between the regressor being \(\beta\) vs. \(\beta'\), then no estimator can achieve error less than \(\|\beta - \beta'\|/2\). To show testing hardness, it suffices to construct a coupling between the distributions of the labeled examples in 1 with regressors \(\beta\) and \(\beta'\) and with small coordinate-wise disagreements, that is, a pair of jointly distributed labeled examples \((X,y), (X',y')\) such that \(\operatorname*{\mathbb{P}}[X_i \neq X_i'] \leq \eta\) for all \(i\in [d]\) and \(\operatorname*{\mathbb{P}}[y \neq y'] \leq \eta\). This is because an adversary can then delete differing coordinates, making the datasets indistinguishable (cf. [fact:coupling]). This leaves us with the (highly) non-trivial task of constructing that coupling, outlined below.

1.2.0.1 Basic construction

The high-level idea of the coupling construction boils down to two main steps: (1) induce the same \(y\)-marginal distribution for some sufficiently distinct \(\beta\) and \(\beta'\), and (2) for each value of \(y\), construct a masking scheme to make \(X \, | \, y\) and \(X' \, | \, y\) indistinguishable. Recall that, since \(X\) (and \(X'\)) and the linear regression noise are all Gaussian, the conditional distributions \(X \, | \, y\) and \(X' \, | \, y'\) are also Gaussian; denote them by \(\mathcal{N}(\mu, \Sigma)\) and \(\mathcal{N}(\mu',\Sigma')\). In all considered regimes of \(\eta\), our choices of \(\beta\) and \(\beta'\) are such that the covariances \(\Sigma\) and \(\Sigma'\) are the same. Then, to obtain a masking scheme that makes the two conditional distributions the same while fully-leveraging the coordinate-wise masking capabilities of the adversary, in [lem:hybrid], we construct a sequence of Gaussians \(N^{(0)}, N^{(1)}, N^{(2)}, \ldots, N^{(d)}\) with the same covariance matrix, which interpolates between \(X \, | \, y\) and \(X' \, | \, y\). Namely, \(N^{(0)} = X \, | \, y\) and \(N^{(d)} = X' \, | \, y\), and each consecutive pair \(N^{(i)}, N^{(i+1)},\) \(i \in \{0, \dots, d-1\},\) has means that differ over a single coordinate. We then construct a masking scheme for the pair \(N^{(i)}\) and \(N^{(i+1)}\) using a maximal coupling⁹ that masks only the \(i^\text{th}\) coordinate, with probability \(D_\mathrm{TV}(N^{(i)}, N^{(i+1)})\) (cf. [cl:one95step]), and that makes the pair of Gaussians indistinguishable.

1.2.0.2 Case \(\|\beta\| \leq \sigma\)

The above basic construction is sufficient for \(\|\beta\| \leq \sigma\) (1[it:thm-small-beta]). In this case, we choose \(\beta = (r, (b/\sqrt{d}) \mathbb{1}_{d-1})\) and \(\beta' = -\beta\), where \(b \leq \eta \sqrt{d}\sigma\), \(\mathbb{1}_\ell\) denotes the dimension-\(\ell\) vector of all ones, and \(r\) is chosen to set the norms of \(\beta, \beta'\) to the target value.

1.2.0.3 Handling the case \(\sigma/\|\beta\| \to 0\) and \(\eta \gg 1/\sqrt{d}\)

¹⁰ The basic construction described above is only meaningful (both in terms of the error and the proof itself) when \(0 < \|\beta \| \leq \sigma\). As \(\sigma/\|\beta\| \to 0\), the conditional distributions \(X\,|\,y\) and \(X'\,|\,y\) approach Gaussians with singular covariances, making the TV-distances between the conditionals equal to one. Thus, the distributions cannot be made indistinguishable by deleting only an \(\eta < 1\) fraction of coordinates.

Our strategy in this case is to effectively reduce the data dimension by one and treat one of the coordinates as noise. For this discussion, fix \(\|\beta\| = \sqrt{d}\), let \(\beta=\mathbb{1}_d\) and \(\beta' = -\beta\); the same ideas transfer to other values of \(\|\beta\|\). The core task again is to couple \(X \, |\, y\) with \(X' \,|\, y'\). Now, given the particular form of \(\beta\) and \(\beta'\) and that \(\beta = -\beta'\), the core task can be re-phrased as coupling \(X \sim \mathcal{N}(0,I)\) conditioned on \(\sum_i X_i = t\) with \(X'\) conditioned on \(\sum_i X_i' = t'\) where \(t = -t'\). Since there is no noise (\(\sigma = 0\)), the key idea is to treat the last coordinate as “noise”: we couple the first \(d-1\) coordinates using the basic coupling method, which gives \(\operatorname*{\mathbb{P}}[X_i \neq X_i'] \leq |t - t'|/d\) for \(i \in [d-1]\). We then set the final coordinate to ensure the total sums are \(t\) and \(t'\). This makes the last coordinate differ with probability 1 between the two scenarios, but randomizing the coordinate order leads to coupled \(X, X'\) where each coordinate differs with probability \((|t - t'| + 1)/d\). This idea and construction is formally stated in [lem:impr-coupling] (which gives a construction for general \(t\) and \(t'\) instead of only \(t = -t'\)). Finally, since in expectation \(|t - t'| = 2|t| = 2|y|\) is \(O(\sqrt{d})\), this yields a bound of \(\Omega(\|\beta\|)\) as long as \(\eta \gg 1/\sqrt{d}\), leading to 1[it:thm-big-eta].

1.2.0.4 Handling \(1/d \ll \eta \ll 1/\sqrt{d}\)

When \(1/d \ll \eta \ll 1/\sqrt{d}\), the approach described so far fails because \(\operatorname*{\mathbb{E}}[|t - t'|]\) is too large, making the probability of coordinate-disagreement larger than the assumed \(\eta\) upper bound. Instead, we construct a new hard instance with \(\beta\) and \(\beta'\) proportional to \((\varepsilon\mathbb{1}_{d/2}, \mathbb{1}_{d/2})\) and \((-\varepsilon\mathbb{1}_{d/2}, \mathbb{1}_{d/2})\) respectively. The coupling proceeds as follows: sample the labels \(y = y'\) from the distribution \(\mathcal{N}(0,(1+\varepsilon^2)/d))\), which is the distribution of the labels. Now consider coupling \(t := \sum_{i=1}^{d/2} X_i\) and \(t' := \sum_{i=1}^{d/2} X_i'\), which are the sums of the first half of the coordinates under each hypothesis. From calculations, the marginals \(t\) and \(t'\) are both univariate Gaussians with the same variance, with difference in mean equal to \(\Theta(\varepsilon y)\). We thus sample \(t\) and \(t'\) in a coupled way where \(t - t' = \Theta(\varepsilon y)\) always (this is possible as the distribution of \(t\) is the same as the distribution of \(t'\) shifted by \(\Theta(\varepsilon y)\)). Now it remains to couple \(X \,|\, y\) and \(X'\,|\,y'\), which we do by first coupling the first half of the coordinates using the coupling between \(t\) and \(t'\), then coupling the second half of the coordinates based on the values of \(y,t\) and \(y',t'\). More concretely, we use [lem:impr-coupling] to couple \((X_1,\ldots,X_{d/2}) \,|\, t\) and \((X'_1,\ldots,X'_{d/2}) \,|\, t'\). Then, we use [lem:impr-coupling] again to couple \((X_{d/2+1},\ldots,X_d)\) and \((X'_{d/2+1},\ldots,X'_d)\) so their sums are \(y - \varepsilon t\) and \(y + \varepsilon t'\), respectively. Going through the calculations, this strategy yields an expected number of disagreeing coordinates bounded by \(\Theta(1 + \varepsilon\sqrt{d})\), which is at most \(\eta d\) as long as \(\varepsilon\ll \eta\sqrt{d}\) (recall that we assume \(\eta \gg 1/d\) in this paragraph). This translates to an estimation error lower bound of \(\Omega(\eta \sqrt{d} \|\beta\|)\), proving 1[it:thm-interm-eta].

1.2.0.5 The case of small \(\eta\)

The remaining regime is \(\eta \ll 1/d\) (1[it:thm-small-eta]). The (already complicated) construction in the last paragraph does not handle that case, given that the expected number of disagreeing coordinates there is \(\Omega(1) \ge 1/d \gg \eta\). For this final regime, we use a similar hard hypothesis test instance, showing the impossibility of distinguishing between \(\beta = (B\mathbb{1}_{d/2},E \mathbb{1}_{d/2})/\sqrt{d/2}\) and \(\beta' = (B\mathbb{1}_{d/2},-E\mathbb{1}_{d/2})/\sqrt{d/2}\) for appropriately chosen parameters \(B\) and \(E\). The main technical innovation here is a deconstruction of label values into the contribution from the first and second half of the coordinates and a more refined way to sample them in the coupling construction. This is sketched in 4 and detailed in 8.

1.3 Preliminaries↩︎

We include only the essential preliminaries here; the full version appears in 6. For \((X,Y) \sim P\), we write \(X\) for the marginal of \(X\) and \(X|(Y=y)\) for the conditional. We denote expectations as \(\operatorname*{\mathbb{E}}_{X \sim D}[X]\). For \(X,Y\) with pdfs \(P_X\) and \(P_Y\), we write \(D_\mathrm{TV}(X,Y)\) or \(D_\mathrm{TV}(P_X,P_Y)\) for the total variation distance, defined as \(\tfrac{1}{2}\int |P_X(u) - P_Y(u) |\mathrm{d}u\). We write \(D_\mathrm{KL}(X \;||\; Y)\) or \(D_\mathrm{KL}(P_X \;||\; P_Y)\) for the Kullback–Leibler divergence, defined as \(D_\mathrm{KL}(P_X \;||\; P_Y) = \int_x P_X(z) \log (P_X(z)/P_Y(z))\mathrm{d}z\). We use \(\mathbb{R}_+\) for non-negative reals, \(\mathbb{Z}_+\) for non-negative integers, and \([n] = {1,\dots,n}\). For a vector \(x\), \(\|x\|\) denotes its Euclidean norm. Let \(I_d\) be the \(d \times d\) identity matrix, and \(\mathbb{1}_d\) the all-ones vector in \(\mathbb{R}^d\). We use \(\top\) to denote matrix transposes. We write \(a \lesssim b\) (or \(O(b)\)) to mean \(a \leq Cb\) for an absolute constant \(C > 0\), independent of the parameters involved (similarly for \(\Omega(\cdot)\)). We write \(\widetilde{O}\) and \(\widetilde{\Omega}\) to hide polylog factors.

If \(P\) and \(Q\) are distributions over \(\mathcal{X}\), a coupling \(\Pi\) of \(P\) and \(Q\) is any distribution over \(\mathcal{X}\times \mathcal{X}\) with marginals \(P\) and \(Q\) such that under \((X,Y) \sim \Pi\), the marginals are \(X \sim P\) and \(Y \sim Q\).

2 Warm-up: Coupling via Hybrid Argument and Proof of 1[it:thm-small-beta]↩︎

We show that our goal reduces to finding a coupling with bounded coordinate-wise disagreements. We also present a method (the “hybrid argument”) for building such couplings. This alone suffices for 1[it:thm-small-beta], but stronger results require a modification that we give in 3.

2.0.0.1 Estimation Hardness via Hypothesis Testing

To show linear regression hardness with error \(<\varepsilon\), it suffices to show an indistinguishability result for vectors \(\beta^{(0)},\beta^{(1)}\) with \(\|\beta^{(0)} - \beta^{(1)}\| > 2\varepsilon\):

• (Null Hypothesis) The data follows the model in [def:linear_model,def:model] with \(\beta = \beta^{(0)}\).

• (Alternative Hypothesis) The data follows the model in [def:linear_model,def:model] with \(\beta = \beta^{(1)}\).

This is sufficient because if an estimator \(\widehat{\beta}\) exists with error less than \(\varepsilon\), the test that rejects the null hypothesis when \(\|\widehat{\beta} - \beta^{(0)}\| > \|\widehat{\beta} - \beta^{(1)}\|\) will succeed, as shown by a simple triangle inequality.

2.0.0.2 Couplings with Small Coordinate-Wise Disagreements

We now argue that constructing a suitable coupling suffices for proving hardness of the testing problem. One can think of data generation as follows: a pair of two samples \(((X,y), (X',y'))\) is drawn from a coupling \(\Pi\) between the distributions of 1 with \(\beta = \beta^{(0)}\) and \(\beta = \beta^{(1)}\). Depending on the hypothesis, either \((X,y)\) or \((X',y')\) is added to the dataset, after which the adversary may modify it ([def:model]). If \(\operatorname*{\mathbb{E}}[\mathbb{1}(X_i \neq X_i')] \leq \eta\) for all \(i \in [d]\), then by deleting disagreeing coordinates, the adversary can make the dataset look identical under both hypotheses, making testing impossible. [fact:coupling] formalizes this. This lemma actually uses a slightly relaxed coordinate-wise disagreement condition: \(\operatorname*{\mathbb{E}}[ \sum_i \mathbb{1}(X_i \neq X_i)] \leq \eta d\). At first glance, this may appear to correspond to a slightly different definition of the contamination model, in which the adversary is allowed to corrupt \(\eta d\) coordinates per sample in expectation. However, [def:model] instead imposes a budget of \(\eta n\) corruptions per coordinate. We can overcome this issue by assuming that the distribution of \((X,y)\) is invariant to permutations of the coordinates of \(X\), because this would imply that after a random permutation each coordinate is deleted with equal probability (at most \(\eta\)). Moreover, because the adversary’s budget in [def:model] is deterministic (rather than in expectation), there is a small probability—dictated by the Chernoff–Hoeffding bound—that our randomized adversary does not satisfy [def:model]. In this case, we may simply assume that the testing problem is solvable. This completes the proof sketch of [fact:coupling]. Finally, we note that we will actually require the permutation invariance condition in the statement of [fact:coupling] to hold separately over the first and second halves of the coordinates, rather than over the entire vector. This modification ensures that the lemma remains applicable in the more involved construction of 1[it:thm-small-eta], which treats the two halves of the coordinates independently. The formal proof is given in 7.1.

2.0.0.3 Hybrid Argument for Coupling Construction

We discuss an argument for constructing couplings satisfying 1 . Instead of coupling the distributions of 1, we will first show, in [lem:hybrid] below, how to construct such couplings between two \(d\)-dimensional Gaussians \(D\) and \(D'\) with same covariance and different means. As we will see in the next subsection, this directly leads to a coupling between the distributions of the linear regression model.

The idea is the following: if the two means differ only in the \(i\)-th coordinate, we can adapt the standard maximal coupling ([fact:maximal]) so that any disagreement comes exclusively from the \(i\)-th coordinate.

[cl:one95step] holds because the marginals over \([d] \setminus \{i\}\) are identical under \(Q\) and \(Q'\), allowing us to match these coordinates, and apply [fact:maximal] to the distribution of the \(i\)-th coordinate conditioned on the rest.

To prove [lem:hybrid] for means differing in more coordinates, we can consider the distributions \(Q_i\) from the statement and couple each consecutive pair using [cl:one95step]. Chaining together the couplings completes the proof of [lem:hybrid]. The formal proofs are in 7.1.

2.0.0.4 Proof Sketch of 1[it:thm-small-beta]

We present a proof sketch with routine calculations omitted here; the full proof is in 7.2. We first focus on the regime \(\|\beta\| \leq \eta \sqrt{d} \sigma/(2\sqrt{2})\). The remaining regime can be handled by an argument that essentially reduces to the \(\|\beta\| \leq \eta \sqrt{d} \sigma/(2\sqrt{2})\) regime and will be explained at the end. For any \(b \in [0,\eta \sqrt{d}\sigma/(2\sqrt{2})]\), we consider the hypothesis testing problem of distinguishing between the regression vectors \(\beta^{(0)} = (b/\sqrt{d},\ldots,b/\sqrt{d})\) and \(\beta^{(1)} = - \beta^{(0)}\) from \(n\) samples from the models in [def:linear_model,def:model]. Note that \(\|\beta^{(0)}\| = \| \beta^{(1)} \| = b\) and \(\|\beta^{(0)} - \beta^{(1)}\| = \sqrt{2} b\). By the reduction from estimation to hypothesis testing, showing that no algorithm can solve this testing problem implies that no estimator has Euclidean error smaller than \(b/\sqrt{2}\).

To use [fact:coupling], we need a coupling \((X,y),(X',y')\) where \((X,y)\) follows 1 with \(\beta = \beta^{(0)}\) and \((X',y')\) follows the same model with \(\beta = \beta^{(1)}\), satisfying \(\operatorname*{\mathbb{E}}[\sum_{i=1}^d \mathbb{1}(X_i = X_i')] \leq \eta d\) and \(\operatorname*{\mathbb{P}}[y \neq y'] \leq \eta\). We construct this coupling by: (i) drawing \(t\) from \(\mathcal{N}(0,\|\beta\|^2 + \sigma^2)\), and setting \(y = y' = t\) (the label distribution is the same for both hypotheses), and (ii) drawing \(X\) and \(X'\) from an appropriate coupling of the conditional distributions given \(y = t\). We thus focus on how to do step (ii) with small coordinate-wise disagreements.

By turning this into a coupling for the full labeled examples as described in the beginning of this proof sketch, the expected sum of disagreements is bounded by simply taking expectation with respect to \(t \sim \mathcal{N}(0,\sigma^2 + b^2)\) above. By using \(\operatorname*{\mathbb{E}}[|t|] \leq \sqrt{\sigma^2 + b^2}\) and plugging in \(b \leq \eta \sqrt{d} \sigma/(2\sqrt{2})\), the expected sum coordinate-wise disagreement becomes at most \(\eta d/2\). Applying [fact:coupling] completes the proof of the lower bound \(\Omega(\|\beta\|)\) in the regime \(\|\beta\| \leq \eta \sqrt{d} \sigma/(2\sqrt{2})\).

For the remaining regime \(\|\beta\| > \eta \sqrt{d} \sigma /(2\sqrt{2})\), consider instead the testing problem in \(\mathbb{R}^{d}\) between \(\beta = (r, b/\sqrt{d}, \ldots, b/\sqrt{d})\) and \(\beta' = (r, -b/\sqrt{d}, \ldots, -b/\sqrt{d})\), where \(b \leq \eta \sqrt{d}\sigma/(2\sqrt{2})\) is as before and \(r\) is a tunable parameter allowing \(\|\beta\|\) to get any desired value larger than \(\eta \sqrt{d} \sigma/(2\sqrt{2})\). The labels can be written as \(y = r X_1 + y_0\) and \(y' = r X_1' + y_0'\), where \(y_0 = (b/\sqrt{d}, \ldots, b/\sqrt{d})^\top X_{2:d} + \xi\) and \(y_0' = (-b/\sqrt{d}, \ldots, -b/\sqrt{d})^\top X_{2:d}' + \xi'\), with \(X_{2:d}\) denoting the last \(d\) coordinates of \(X\). From earlier, there exists a coupling between \((X_{2:d}, y_0)\) and \((X_{2:d}', y_0')\) with expected disagreements at most \(\eta d/2\) per sample which extends trivially to the full \((X, y)\) and \((X', y')\) by adding shared Gaussian noise \(\mathcal{N}(0, r^2)\) to the labels. Applying [fact:coupling] as before completes the proof. 0◻

3 Improved Core Coupling Construction↩︎

The previous strategy has a key limitation that the result holds only for \(\|\beta\| \leq \sigma\). When \(\sigma = 0\), the matrix \(\Sigma\) in 2 becomes singular, causing the proof to break down and yield no meaningful bound. For some regimes in 1, we want to overcome this, for example, to establish the lower bound \(\Omega(\|\beta\|)\) whenever \(\eta \geq 7/\sqrt{d}\), regardless of the value of \(\sigma\). The key observation is that, because \(\beta\) in 2 points in the all-ones direction, the zero-eigenvalue issue disappears when we restrict our attention to the first \(d-1\) coordinates. We can thus apply the technique to those coordinates and adjust the final one to ensure the total sum is correct. This idea is formalized in [lem:impr-coupling], which provides an improved version of 2 for the case \(b = \sqrt{d}\) and \(\sigma = 0\) (the result extends trivially to any scaling of the vector and any \(\sigma > 0\)). We will later use [lem:impr-coupling] to prove our remaining main results.

4 Proofs of parts [it:thm-big-eta]-[it:thm-small-eta] of 1↩︎

The proof of 1[it:thm-big-eta] is almost the same as the proof of 1[it:thm-small-beta]: we test between linear regression setups with \(\sigma {=} 0\) and \(\beta {=} (1,\ldots,1)\) vs. \(\beta' {=} -\beta\). It suffices to use \(\sigma {=} 0\) and \(\|\beta\| {=} \sqrt{d}\), since coupling these cases implies a coupling for any scaling \(\|\beta\|\) and \(\sigma {>} 0\) (see 8 for the full proof). Using [lem:impr-coupling,fact:coupling] completes the proof.

The remaining lower bound proofs also use [lem:impr-coupling]. However, since we no longer aim to prove a lower bound of order \(\|\beta\|\) (as the upper bound is smaller), we will not use the hypothesis testing setup with \(\beta = (1,\ldots,1)\) vs.\(\beta' = (-1,\ldots,-1)\). Instead, we adjust it as follows. The full proofs for this section are provided in 8.

4.0.0.1 Proof Sketch of 1[it:thm-interm-eta]

For this lower bound we need a slightly different hypothesis testing problem: Our two hypotheses will use regressors \(\beta = s(\varepsilon,\ldots,\varepsilon,1,\ldots,1)\) and \(\beta' = s(-\varepsilon,\ldots,-\varepsilon,1,\ldots,1)\) respectively, where \(s > 0\) can be any scaling factor (so that our lower bound holds for every possible norm of \(\|\beta\|\)) and \(\varepsilon\leq 1\) will be chosen shortly. We claim that [lem:impr-coupling] implies a coupling \(\Pi\) between the distributions of 1 for \(\beta\) and \(\beta'\) such that \(\operatorname*{\mathbb{P}}_{((X,y)(X',y')) \sim \Pi}[y \neq y']=0\) and \[\label{eq:coupling-combined} \operatorname*{\mathbb{E}}\limits_{((X,y),(X',y')) \sim \Pi}\left[\sum_{i=1}^d \mathbb{1}(X_i \neq X_i') \right] \leq 2 + O(\varepsilon\sqrt{d} ) .\tag{4}\] We first complete the theorem proof sketch given 4 . Setting \(\varepsilon:= (\eta - 4/d)\tfrac{\sqrt{d}}{2C }\) (for \(C\) the constant in the big-O), the RHS of 4 becomes \(\eta d/2\), making [fact:coupling] applicable. By the reduction from estimation to testing, this implies every estimator has error \(\Omega(\|\beta-\beta'\|)\), which is \(\|\beta - \beta'\| = \varepsilon s\sqrt{d/2} = \tfrac{\varepsilon}{\sqrt{1+\varepsilon^2}}\|\beta\| \gtrsim \varepsilon\|\beta\| \gtrsim (\eta - 4/d)\sqrt{d}\|\beta\| \gtrsim \eta \sqrt{d}\|\beta\|\), where the last step used \(\eta \geq 5/d\). In Appendix we use better constants so only require \(\eta \geq (2{+}o(1))/d\).

We now prove 4 . By scaling, it suffices to do this only for \(s=1\) and \(\sigma=0\) (i.e., no additive noise). The coupling is the following: (i) Sample \(z\) from the distribution of labels (which is \(\mathcal{N}(0,\sigma^2 + \|\beta\|^2)= \mathcal{N}(0,(1+\varepsilon^2)d/2)\) for both hypotheses) and set \(y=y'=z\). (ii) Sample \(t\) and \(t'\) from the distribution of the sum of the first half of the coordinates for the two hypotheses. By using standard facts for Gaussians, it turns out that these two distributions are \(\mathcal{N}(\varepsilon z/(1+ \varepsilon^2), d/(2(1+\varepsilon^2)),\mathcal{N}(-\varepsilon z/(1+ \varepsilon^2), d/(2(1+\varepsilon^2))\). We can also choose to have that \(t-t' = 2\varepsilon z/(1+ \varepsilon^2)\) always. (iii) Sample the first half of the coordinates using the coupling of [lem:impr-coupling] such that they sum to \(z-\varepsilon t\). (iv) Sample the second half using [lem:impr-coupling] such that they sum to \(z+\varepsilon t'\). The expected disagreements from the two applications of [lem:impr-coupling] are \[\begin{align} \operatorname*{\mathbb{E}}\left[ \textstyle \sum_{i=1}^d \mathbb{1}\left( X_i {\neq} X_i'\right) \right] \leq 2 + O\big(\operatorname*{\mathbb{E}}_{t,t'}\left[|t {-} t'| \right] {+} \varepsilon\operatorname*{\mathbb{E}}_{t,t'}\left[|t {+} t'| \right] \big). \end{align}\] Using \(t-t' = 2\varepsilon z/(1 + \varepsilon^2)\) with \(z \sim \mathcal{N}(0,(1{+}\varepsilon^2)d/2)\), \(\varepsilon\leq 1\) and that \(t,t'\) have variance \(O(d)\) proves 4 .0◻

The constant additive term in 4 breaks the proof of 1[it:thm-interm-eta] when \(\eta < 2/d\). To handle this, we improve the coupling below. Roughly speaking, the idea is to break the variables we want to couple into smaller parts and couple the parts separately. When the first parts match, we force the rest to match as well, which (in expectation) reduces the coordinate-wise disagreements.

4.0.0.2 Proof Sketch of 1[it:thm-small-eta]

Let \(B,E \in (0,1)\) be parameters, let \(\mathbb{1}_{d}\in \mathbb{R}^d\) denote the all-ones vector, and define the problem of distinguishing between regressor \(\beta := (\tfrac{B}{\sqrt{d/2}}\mathbb{1}_{d/2}, \tfrac{E}{\sqrt{d/2}}\mathbb{1}_{d/2})\) and \(\beta' := (\tfrac{B}{\sqrt{d/2}}\mathbb{1}_{d/2}, \tfrac{-E}{\sqrt{d/2}}\mathbb{1}_{d/2})\). The claim now is existence of a coupling \(\Pi\) between labeled examples of the linear regression models of 1 with \(\beta\) and \(\beta'\) such that, \(\operatorname*{\mathbb{P}}_{((X,y)(X',y')) \sim \Pi}[y \neq y']=0\) and \[\begin{align} \label{eq:coupling-combined-2} \operatorname*{\mathbb{E}}_{((X,y),(X',y')) \sim \Pi}\left[\sum_{i=1}^d \mathbb{1}(X_i \neq X_i') \right] \lesssim \frac{E}{\sigma} + \sqrt{d} \frac{E}{B} . \end{align}\tag{5}\] As before, we aim to apply [fact:coupling] to conclude that the hypothesis testing problem is hard. To do that, we need to ensure that the RHS above is at most \(\eta d/2\). Equivalently \(E/\sigma \ll \eta d\) and \(\sqrt{d}E /B \ll \eta d\). We thus choose \(E:= \frac{1}{2C } \min(\eta d \sigma, \eta \sqrt{d}B)\). Noting that \(\|\beta - \beta'\| \gtrsim \min(\eta d \sigma, \eta \sqrt{d} \|\beta\|)\) concludes the proof of 1[it:thm-small-eta].

It remains to argue 5 . The coupling for that is the following. We will use the notation \(X=(X_1,X_2)\) to denote the first and second half of coordinates of \(X\) and \(S_1 = \mathbb{1}_{d/2}^\top X_1 /\sqrt{d/2}\), \(S_2 = \mathbb{1}_{d/2}^\top X_2 /\sqrt{d/2}\) to denote the scaled sums of these coordinates (and we will use similar primed letters notation for the corresponding parts of \(X'\)). The coupling is the following: (i) Draw \(y\) from the distribution of labels and let \(y'=y\). (ii) Draw \((S_2,S_2')\) from a maximal coupling between the distributions under the null and alternative hypotheses conditioned on \(y\). (iii) If \(S_2 = S_2’\): draw \(X_2\) from the null distribution conditioned on the value of \(S_2\) from the previous step, and set \(X_2’=X_2\). Draw \((S_1,S_1’)\) from a maximal coupling between null and alternative distributions conditioned on \((S_2, y)\). (iii a) If \(S_1=S_1’\) draw \(X_1\) from the null distribution conditioned on \(S_1\) and set \(X_1’=X_1\). (iii b) Otherwise use [lem:impr-coupling] to jointly sample \((X_1,X_1’)\) conditioned on the values of \(S_1\) and \(S_1’\). (iv) If \(S_2 \neq S_2’\): couple the conditional distributions \(S_1\) under the null and \(S_1’\) under the alternative given \((S_2, y)\) and \((S_2', y)\), respectively. Then use [lem:impr-coupling] to jointly sample \((X_2, X_2')\) conditioned on \(S_2, S_2'\), and similarly for \((X_1, X_1')\) conditioned on \(S_1, S_1'\).

It is true (and straightforward) that this is indeed a coupling, since the labeled examples are constructed in steps and each step uses the correct conditional distributions under the null/alternative hypotheses, conditioned on all previous steps.

For the disagreement analysis, one must derive the precise form of all the conditional distributions mentioned above and apply [fact:maximal,lem:impr-coupling] appropriately. This requires some technical work, which we defer to 8. At a high level however, the key point is that the case-splitting in the coupling construction implies that in some cases there are no disagreements, and [lem:impr-coupling] is applied only to the remaining cases, which occur with probability strictly less than \(1\). As a result, the additive \(1\) term in the disagreement bound from [lem:impr-coupling] is scaled by the probability of the corresponding case. Ultimately, this ensures that the expected total number of disagreements is strictly less than \(1\), which in turn allows our theorem to apply in the regime \(\eta < 1/d\).

Acknowledgments↩︎

The authors would like to thank the anonymous COLT reviewer for identifying a minor bug in an earlier version of this paper.

Appendix↩︎

The structure of the appendix is as follows: In 5, we provide material omitted from 1 of the main body, including the detailed description of upper bounds for our learning task, implications for other corruption models, a detailed summary of related work, and open problems. 6 records the notation and mathematical background required for our technical results. Finally, 7 and 8 give the technical proofs omitted from 2 and 4 respectively.

5 Omitted Details from 1↩︎

5.1 Discussion on [fact:known-algorithms]↩︎

We provide a more detailed presentation of the three algorithms for linear-regression under coordinate-wise corruptions with the specific references to prior robust statistics literature below.

5.2 Derivation of Regimes in 1↩︎

With [fact:known-algorithms] in hand, we discuss how the parameter space can be partitioned based on which algorithm from [fact:known-algorithms] achieves the best error (up to absolute constant factors). Due to space constraints in 1.1, the regimes were summarized there in a single, somewhat condensed table (1). In this section, we present a more detailed version of this partitioning in [tab:eta-main,table:medium-eta,table:small-eta] and explicitly explain how the different regimes are derived. 1 in 1 is a combined version of [tab:eta-main,table:medium-eta,table:small-eta] into a single table.

5.2.0.1 Regime \(7/\sqrt{d} \leq \eta \leq 1\)

Algorithm \(\mathcal{A}_1\) from [fact:known-algorithms] is not applicable in this regime. Among the other two algorithms, \(\mathcal{A}_3\) always has the best error (because the error of \(\mathcal{A}_2\) is in the order of \(\eta \sqrt{d}\sqrt{\|\beta\|^2 + \sigma^2} \geq \eta \sqrt{d} \|\beta\| \geq 7\|\beta\|\)). This explains the last column in ¿tbl:tab:eta-main?.

5.2.0.2 Regime \(0.49/d \leq \eta < 7/\sqrt{d}\)

Again, algorithm \(\mathcal{A}_1\) is not applicable in this regime. Between \(\mathcal{A}_2\) and \(\mathcal{A}_3\), we note that whenever \(\|\beta\| \geq \eta \sqrt{d} \sigma\), the error of \(\mathcal{A}_2\) is \(O(\eta \sqrt{d}\sqrt{\|\beta\|^2 + \sigma^2})=O(\eta \sqrt{d}\|\beta\| + \eta \sqrt{d}\sigma) = O(\|\beta\|)\), i.e., \(\mathcal{A}_2\) has the best error (up to constant factors). Moreover, when \(\|\beta\| < \sigma\), the term that dominates in the error of \(\mathcal{A}_2\) is \(\eta\sqrt{d}\sigma\) and when \(\|\beta\| \geq\sigma\) the dominating \(\mathcal{A}_2\) error term is \(\eta \sqrt{d} \|\beta\|\). This explains the regimes shown in [table:medium-eta].

5.2.0.3 Regime \(0 \leq \eta < 0.49/d\)

In this regime, all three algorithms are applicable. For the same reason as before, the error of \(\mathcal{A}_3\) is better than that of \(\mathcal{A}_2\) when \(\|\beta\| < \eta \sqrt{d} \sigma\) and \(\mathcal{A}_2\) is better than \(\mathcal{A}_3\) otherwise. It also outperforms \(\mathcal{A}_1\), whose error is on the order of \(\eta d \sigma\). In the sub-regime \(\eta \sqrt{d} \sigma \leq \|\beta\| < \sqrt{d} \sigma\), the error of \(\mathcal{A}_2\) is at most \(O(\eta \sqrt{d} \sigma + \eta \sqrt{d} \|\beta\|)\), with both terms being smaller than \(O(\eta d \sigma)\). Thus, \(\mathcal{A}_2\) achieves the best error in this sub-regime. In the final sub-regime \(\|\beta\| \geq \sqrt{d}\sigma\), the error of \(\mathcal{A}_1\) is of the order \(\eta d \sigma\), which is smaller than that of \(\mathcal{A}_3\) (which is of the order \(\|\beta\| \geq \sqrt{d} \sigma \geq \sigma \geq \eta d \sigma\)). Comparing with \(\mathcal{A}_2\), the dominating term in the error of \(\mathcal{A}_2\) is the \(\eta \sqrt{d} \|\beta\|\) term which is larger than \(\eta d \sigma\) whenever \(\|\beta\| > \sqrt{d}\sigma\). Thus \(\mathcal{A}_1\) has the best error (up to a constant factor) in the sub-regime \(\|\beta\| > \sqrt{d}\sigma\). This explains [table:small-eta].

5.3 Unified Algorithm with Error Guarantees as in [tab:eta-main,table:medium-eta,table:small-eta]↩︎

We further claim that there exists a single algorithm that can automatically adapt to the best choice among \(\mathcal{A}_1, \mathcal{A}_2, \mathcal{A}_3\) from [fact:known-algorithms], and achieve the error upper bounds in [tab:eta-main,table:medium-eta,table:small-eta], without requiring a priori knowledge of the regime in which \(\|\beta\|\) lies. The algorithm only needs to know \(\eta\) (which could also be avoided by applying standard techniques like Lepski’s method, but is beyond the scope of this work).

5.4 Other Implications of Our Results↩︎

In this brief subsection, we discuss some immediate implications of our results.

5.4.0.1 Implications on Models with a Restricted Set of Corrupted Coordinates

In this paper, we primarily discuss the case in which every coordinate is subject to an \(\eta\) fraction of corruptions. However, another reasonable model would be to consider the case where only a subset \(\mathcal{S} \subset \{1, 2, \dots, d\}\) of \(|\mathcal{S}| = r \leq d\) coordinates is subject to adversarial erasures as in [def:model] or more general corruptions as in [def:model2]. It is straightforward to verify that all our results generalize to this model, with estimation error bounds now depending on \(r\) in place of \(d\) and on \(\|\beta_{\mathcal{S}}\|_2\) in place of \(\|\beta\|_2,\) where \(\beta_{\mathcal{S}} = \sum_{i\in \mathcal{S}}\beta_i e_i\) and \(e_i\) denotes the \(i^{\rm th}\) standard basis vector.

5.4.0.2 Implications for Generalized Linear Models

While our primary focus in this work is on lower and upper bounds for the attainable estimation error \(\|\beta - \widehat{\beta}\|_2\), some of our lower bounds have broader implications on the attainable squared error of generalized linear models (GLMs)—where the “ground truth” labels follow \(f(\beta^\top X)\) for some possibly more general function \(f\) called the activation, instead of \(\beta^\top X\) corresponding to the case of linear regression (where \(f(t) =t\) is the identity function). In this case, the focus is usually on bounding the mean squared error for a predictor \(\widehat{\beta}\), defined by \(\mathcal{L}(\widehat{\beta}) = \operatorname*{\mathbb{E}}_{(X, y)\sim \mathcal{D}}[(y - f(\widehat{\beta}^\top X))^2]\) for \((X, y)\) jointly distributed according to a distribution \(\mathcal{D}\). It is immediate that in the case of linear regression model from 1, we have \[\begin{align} \mathcal{L}(\widehat{\beta}) &= \sigma^2 + \operatorname*{\mathbb{E}}_{X\sim \mathcal{N}}[(\beta^\top X - \widehat{\beta}^\top X)^2]\\ &= \sigma^2 + \operatorname*{\mathbb{E}}_{X\sim \mathcal{N}}[(\beta - \widehat{\beta})^\top XX^\top (\beta - \widehat{\beta})]\\ &= \sigma^2 + \|\beta - \widehat{\beta}\|_2^2. \end{align}\] Thus, our lower bounds on \(\|\beta - \widehat{\beta}\|_2\) directly imply lower bounds on the mean squared error for linear regression, using the above derivation.

Concerning GLMs in broader generality, let us consider the realizable case, where all labeled examples \((X, y)\) satisfy \(y = f(\beta^\top X)\) for some fixed vector \(\beta\) and activation function \(f.\) Notice first that the lower bounds in 1[it:thm-big-eta] and 1[it:thm-interm-eta] apply even when \(\sigma = 0\).¹¹ Our lower bounds for linear regression are proved by constructing a coupling between the distribution of \(X\) conditioned on the value of \(y=t\) under two distinct hypotheses corresponding to sufficiently different prediction/weight vectors \(\beta\) and \(\beta'\). Since for \(\sigma = 0\) in the case of linear regression this is equivalent to conditioning on the value of \(\beta^\top X = (\beta')^\top X = t\), the same lower bounds (implying impossibility of distinguishing between \(\beta\) and \(\beta'\)) generalize to the case where \(y = f(t)\). This is because considering more general functions \(f\) that are not identity (as in the case of linear regression) can only further obscure information if \(f\) is not invertible (as is the case, for example, for ReLU activations where \(f(t) = \max\{0, t\}\)).

The above discussion implies that the same lower bounds as in 1[it:thm-big-eta] and 1[it:thm-interm-eta] apply for the value of \(\|\widehat{\beta} - \beta\|_2\), where \(\beta\) is the “ground truth” predictor and \(\widehat{\beta}\) the predictor that can be constructed by an algorithm. A consequence is that the mean squared error can be bounded below as a function of \(\|\widehat{\beta} - \beta\|_2^2\) for a broad class of GLMs. In particular, for the class of \((a, b)\)-unbounded activations—namely, monotonically non-decreasing functions \(f\) that are \(b\)-Lipschitz, satisfy \(f(0) = 0,\) and are such that \(f'(t) \geq a > 0\) for \(t > 0\)—known results (cf.[35], [36]) imply that for \(X\) distributed according to the standard normal distribution, we have \(\mathcal{L}(\widehat{\beta}) = \Omega(\mathrm{poly}(b/a))\|\widehat{\beta} - \beta\|_2^2.\) In particular, since ReLU activation is \((a, b)\)-unbounded with \(a = b = 1,\) we have that for the ReLU activation, \(\mathcal{L}(\widehat{\beta}) = C\|\beta - \widehat{\beta}\|_2^2\) for a universal constant \(C > 0.\) As a consequence, the mean squared error is of the order \(\Omega(\|\beta\|_2^2)\) for \(\eta \in [7/\sqrt{d}, 1]\) and of the order \(\Omega(c^2\eta^2 d \|\beta\|_2^2)\) for \(\eta \in [(2+c)/d, 7/\sqrt{d}].\) In other words, unless both \(\eta\) and \(\|\beta\|_2\) are small, no meaningful prediction—as measured by the mean squared error—is possible.

These observations reinforce the argument that linear regression is the most basic model to study when considering limitations imposed by adversarial coordinate-wise corruptions studied in our work.

5.5 Related Work↩︎

There is a vast literature on learning with missing data in various forms. To the best of our knowledge, this work is the first to address the fully adversarial missingness model of [def:model] in the context of high-dimensional regression. The most relevant prior work is [29], which considers fully adversarial missingness, but for the simpler task of mean estimation. [37] also considers mean estimation and uses a milder setting for the missingness model. We provide a discussion of these works below, followed by a broader (though not exhaustive) review of related literature, with an emphasis on settings most closely aligned with ours.

5.5.0.1 Comparison with [29]

This work studies the optimal error of Gaussian mean estimation in the same coordinate-wise adversary of [def:model,def:model2] (among others). Concretely, it considers a slightly more general problem than the vanilla Gaussian mean estimation where data has some additional low-dimensional structure: each sample \(x_i\) follows \(\mathcal{N}(\mu,\Sigma)\) but can be written as \(x_i = A z_i\), where \(A \in \mathbb{R}^{d \times r}\) and \(z \in \mathbb{R}^{r}\) is some lower dimensional Gaussian. The goal is to estimate \(\mu\). Their lower bounds on the estimation error are based on the approach of constructing couplings with small coordinate-wise disagreements (as in the present work). In terms of algorithms, the approach of [29] for the missing data of [def:model] is a two-step algorithm that first estimates the missing values and then runs existing estimators on the imputed dataset. They show that the first step becomes NP-hard when data is replaced ([def:model2]) instead of missing. However, they are able to get some preliminary results by a randomized algorithm for special cases (like when \(A\) is known). In the most general case for the mean estimation problem of [29], the answer to the question of whether missing data are easier to handle than replaced data remains unclear. In the case of the linear regression problem considered in this paper, we find that the answer is negative. On the technical side, our work builds on the coupling technique of [29], but significantly extends it to handle new challenges in the linear regression setting. The main difficulty is that the estimation error behaves differently across many distinct regimes, each requiring a separate analysis. As discussed in 1.2, a direct application of the [29] coupling fails when the additive noise \(\sigma\) is small, leaving several regimes unaddressed. We resolve this by coupling only the first \(d-1\) coordinates and treating the last one as additive noise. This regime also requires a redesigned hypothesis testing setup. The original test inspired from [29] that compares \(\beta = (1,\ldots,1)\) and \(\beta' = (-1,\ldots,-1)\), is tailored to show an \(\Omega(\|\beta\|)\) lower bound, i.e., that no non-trivial estimation is possible. In our case, however, some non-trivial estimation is possible, so we design a finer test: we split the coordinates into two halves with differing values in \(\beta\) and \(\beta'\) and apply separate couplings to each half. Finally, one particularly challenging regime remains, which we address by further splitting both the coordinates and the labels into two parts, and using a more delicate coupling between the parts of covariates and labels.

5.5.0.2 Comparison with [37]

This work also studies mean estimation under a missing data model. Their setting involves a combination of two types of erasures: up to an \(\eta\)-fraction of coordinates may be missing, and up to an \(\varepsilon\)-fraction of entire samples may be deleted. The second type of deletion (entire samples) can be fully adversarial, while the first type (coordinate-wise) is less adversarial than in our model, because in [37] the adversary must commit to the missingness pattern before seeing the data. The algorithms proposed in [37] are based on imputing missing values, and thus do not naturally extend to the contamination setting, where adversaries may replace values rather than erase them.

5.5.0.3 Literature on Missing Data

The study of missing data has a long history, motivated by concerns like those discussed in 1. In a seminal work, [13] distinguishes three types of missingness: when missingness is independent of the data, it is called missing completely at random (MCAR); when it depends only on the observed data, it is missing at random (MAR); otherwise, it is missing not at random (MNAR), where missingness depends on both observed and unobserved values. See also [38]–[41] for a more in-depth review.

Much of the literature on regression with missing data has focused on parameter estimation. Early works include [17], [18], as well as [20], which proposes a sparse variant of the Lasso estimator. More recent approaches revisit the problem from the perspectives of imputation and collaborative learning [30], [31].

Other recent efforts have shifted focus to the label prediction problem, which is a fundamentally different goal, as the test set is also expected to contain missing entries. Notably, the Bayes optimal predictor in this setting decomposes into a sum over predictors for each missingness pattern, and the number of such patterns can be exponential. This pattern-specific structure has been studied in [21], [22], [42], which characterize minimax-optimal rates under MAR and MCAR (and in [43] for classification). A body of work also focuses on the idea of using a two-step procedure: first impute the missing entries, then apply algorithms designed for complete data. This method has been shown to be Bayes optimal in various settings [14], [15], [44].

In light of the above, our work fits into the MNAR setting, but to the best of our knowledge, is not directly related to any of the aforementioned works. Unlike the prediction-focused literature (which is more challenging and often incurs exponential complexity) we study parameter estimation. Moreover, prior work on parameter estimation for linear regression either assumes MAR [19], [31] or considers much more benign forms of MNAR; for instance, [30] assumes that the missingness pattern is drawn i.i.d. from an arbitrary distribution.

5.5.0.4 Truncated Statistics

Although technically this literature falls into the Missing Not At Random category from the previous paragraph, it is sufficiently well developed to have its own place in the literature. Truncated statistics refer to the situation where samples falling outside of a fixed set are censored. Parameter estimation under this setting traces back to the early work of Galton, Pearson and Lee [45]–[47]. Despite early interest in the problem, computationally efficient algorithms for multivariate Gaussians were developed relatively recently, with some notable works being [48], [49] for mean estimation and [50]–[52] for linear regression. While the previous works operate under a setting where each sample is either entirely censored or entirely visible, [53] extends the setting to partially truncated data, which each coordinate individually visible or censored (for the problem of Gaussian mean estimation). That said, the setting considered in [53] is still milder than ours. It considers either (i) a setting where each coordinate is erased based on membership in a fixed, predetermined set (to which the algorithm has oracle access), or (ii) erasures based on a simple projection rule. Under this model, consistent estimation is possible, as shown in [53], whereas in our fully adversarial erasure setting, it is provably not.

5.5.0.5 Robust Statistics

Robust Statistics was developed in the 1960s [23], [24] to handle scenarios where a small fraction of the data points are arbitrarily corrupted. As in the setting of this work, such corruptions do not allow for consistent estimation. Early work focused on characterizing the information-theoretic error for univariate Gaussian mean estimation. Since then, a large body of work developed estimators for various tasks [32]. However, these estimators were computationally inefficient for high-dimensional tasks. The field saw a resurgence with the work of [25], [26] that tackled high-dimensional robust mean estimation in polynomial time. In the last decade, a plethora of robust estimators have been introduced for different tasks, including linear regression [54]–[57] considered in this work. For a comprehensive treatment, see the book [27]. The key difference is that the de facto contamination models in robust statistics treat each sample as either fully clean or fully corrupted, making non-trivial estimation impossible as the corruption rate approaches \(1/2\). This fails to capture cases where each coordinate is only mildly corrupted, i.e., few samples are corrupted per coordinate, even if most samples are affected overall. Such settings may still allow for meaningful estimation, but algorithms from the robust statistics literature are inapplicable.

5.6 Open Problems↩︎

While our work addresses the estimation error for the standard model of linear regression with coordinate-wise corruptions by providing matching upper and lower bounds for essentially all possible parameter regimes, there are several avenues for future research that merit further investigation.

First, it would be interesting to consider \(X \sim \mathcal{N}(0,\Sigma)\) in the linear regression problem, with \(\Sigma \neq I\) and unknown to the algorithm in the linear regression model. Does the conclusion that dealing with missing data is no easier than handling adversarially replaced data still hold in this setting?

Second, although the algorithms are efficient in terms of runtime (i.e., polynomial time) and sample complexity (using \(\widetilde{O}(d/\eta^2)\) samples), the optimal sample complexity for the problem remains unclear. In the clean setting (where all samples follow the distribution of 1), the optimal sample complexity is \(d/u^2\), where \(u\) is the target estimation error. In some regimes of [tab:eta-main,table:medium-eta,table:small-eta], the estimation error scales with \(\|\beta\|\), which can be very large. This implies that (at least in the noiseless setting) fewer samples are needed to achieve that level of error.

Third, as discussed in 5.4, some of our lower bounds directly imply lower bounds for the mean squared error of a broad class of generalized linear models that includes ReLU as a special case, even with perfect, noise-free labels (i.e., in the realizable case). It is however unclear if those lower bounds are tight or if they can be strengthened further, since non-invertible activations like ReLU can further obscure information available to the algorithm (because the label is zero whenever the argument \(\beta^\top X\) is negative). Relatedly, it would be interesting to develop polynomial-time algorithms that can match the information-theoretic lower bounds (or prove no such algorithms exist, even when restricted to classes such as Statistical Query algorithms).

Finally, our results demonstrate that adversarial coordinate-wise deletions make both estimation and prediction challenging, even in the most basic model of Gaussian linear regression and with infinite samples. On the other hand, a fully random model of deletions (where deletions are independent of the observed data) allows for the diminishing estimation error as the sample size is increased. These two extremes beg the question of what lies in between, when the deletions are neither fully adversarial nor fully random.

6 Preliminaries↩︎

We provide the full version of the preliminaries here.

6.1 Notation↩︎

6.1.0.1 Basic Probability Notation

We write \(X \sim D\) to denote a random variable \(X\) that is distributed according to the distribution \(D\). We use \(P_X(x)\) to denote the pdf of \(X\). For multiple variables, we use \((X,Y) \sim D\) to denote that \(X\) and \(Y\) are jointly distributed according to \(D\), and write \(P_{X,Y}(x,y)\) for the pdf of that joint distribution. If \((X,Y) \sim D\) we will also use the notation \(X|Y=y\) to denote the random variable distributed according to the conditional distribution of \(X\) given the occurrence of the value \(y\) for \(Y\), i.e., the distribution whose pdf is the conditional density function \(P_{X|Y}(x|y) = P_{X,Y}(x,y)/P_Y(y)\), and we denote by \(X\) the random variable distributed according to the marginal distribution of \(D\), i.e., the distribution with pdf \(P_X(x) = \int P_{X,Y}(x,y) \mathrm{d}y\). We use \(\operatorname*{\mathbb{E}}_{X \sim D}[X]\) for the expectation of \(X\). If \(X \sim P_X\) and \(Y \sim P_Y\) are continuous random variables over some domain \(\mathcal{X}\), with pdfs \(P_X(x)\) and \(P_Y(y)\), we denote by \(D_\mathrm{TV}(X,Y) = \frac{1}{2} \int_{x \in \mathcal{X}} |P_X(x) - P_Y(x)| \mathrm{d}x\) the total-variation distance between the distributions \(P_X\) and \(P_Y\) (by slightly abusing notation, \(D_\mathrm{TV}(X,Y)\) may some times be also denoted as \(D_\mathrm{TV}(P_X,P_Y)\)). We denote the Kullback–Leibler (KL) divergence between distributions \(X \sim P_X\) and \(Y \sim P_Y\) by \(D_\mathrm{KL}(X \,\|\, Y) = \int_{x \in \mathcal{X}} P_X(x) \log\left(\frac{P_X(x)}{P_Y(x)}\right) \mathrm{d}x\), assuming \(P_X\) is absolutely continuous with respect to \(P_Y\).

6.1.0.2 Basic Notation

We use \(\mathbb{Z}_+\) for the set of positive integers. We denote \([n]=\{1,\ldots,n\}\). For a vector \(x\) we denote by \(\| x \|\) its Euclidean norm. Let \(I_d\) denote the \(d\times d\) identity matrix (omitting the subscript when it is clear from the context). We use \(\mathbb{1}_d\) for the all-ones vector in \(\mathbb{R}^d\). We use \(\top\) for the transpose of matrices and vectors. We use \(|A|\) to denote the determinant of matrix \(A\). We use \(a\lesssim b\) to denote that there exists an absolute universal constant \(C>0\) (independent of the variables or parameters on which \(a\) and \(b\) depend) such that \(a\le Cb\). In our notation \(a = O(b)\) has the same meaning as \(a \lesssim b\) (similarly for \(\Omega(\cdot)\) notation) We use \(\widetilde{O}\) and \(\widetilde{\Omega}\) to hide polylogarithmic factors.

6.1.0.3 Couplings

If \(P\) and \(Q\) are probability distributions over \(\mathcal{X}\), then a coupling \(\Pi\) of \(P\) and \(Q\) is any distribution over \(\mathcal{X}\times \mathcal{X}\) such that the marginals of \(\Pi\) coincide with \(P\) and \(Q\). We write \((X,Y) \sim \Pi\) to denote random variables that are distributed according to the coupling; \(X\) is marginally distributed according to \(P\) and \(Y\) according to \(Q\).

6.2 Useful Probability and Linear Algebraic Facts↩︎

We start with some facts about Gaussian distributions, namely the inner product of a Gaussian and a fixed vector is another Gaussian, and the distribution conditioned on one of the variables is also Gaussian.

The following fact provides a closed-form formula for the KL-divergence between two multivariate Gaussians. The formula can be derived by direct computation and using properties from Section 8.2 of [58].

We will bound total variation distances using the previous fact combined with Pinsker’s inequality, stated below.

Combining the above two facts we obtain the following bound on the TV distance between univariate Gaussians.

The following fact states that there exists a coupling between two distributions such that the probability of disagreement is at most their total variation distance.

We also require the multiplicative version of the Chernoff-Hoeffding bound for binary random variables:

Finally, the following linear algebraic fact provides a useful formula for inverting matrices of the form identity minus a rank-one matrix.

7 Omitted Details from 2↩︎

7.1 Couplings with Small Coordinate-Wise Disagreements↩︎

We restate and prove the following statements.

Before presenting the proof, we note that although the statement is written in terms of Gaussian distributions (which is the setting we will apply it to later), the argument applies to any distributions \(D\) and \(D'\), where \(D'\) is a shifted version of \(D\) (in the sense that a random variable from \(D'\) can be written as a random variable from \(D\) plus a deterministic vector).

7.2 Proof of 1[it:thm-small-beta]↩︎

We restate and prove the following lower bound.

8 Omitted Details from 4↩︎

In this section we restate and prove the lower bounds corresponding to parts [it:thm-big-eta],[it:thm-interm-eta] and [it:thm-small-eta] of 1. Since these bounds depend on [lem:impr-coupling], we restate that lemma below for convenience:

We start with 1[it:thm-big-eta]

We now move to 1[it:thm-interm-eta] which is restated and proved below.

We restate and prove 1[it:thm-small-eta] below.

References↩︎