{'fraction_non_alphanumeric': 0.04729391295178341, 'fraction_numerical': 0.01732001159939894, 'mean_word_length': 4.045085782683867, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 19, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Conditional Linear Regression for Heterogeneous Covariances

Often machine learning and statistical models will attempt to describe the majority of the data. However, there may be situations where only a fraction of the data can be fit well by a linear regression model. Here, we are interested in a case where such inliers can be identified by a Disjunctive Normal Form (DNF) formula. We give a polynomial time algorithm for the conditional linear regression task, which identifies a DNF condition together with the linear predictor on the corresponding portion of the data. In this work, we improve on previous algorithms by removing a requirement that the covariances of the data satisfying each of the terms of the condition have to all be very similar in spectral norm to the covariance of the overall condition.


INTRODUCTION
Linear regression is a technique frequently used in statistical and data analysis. The task for standard linear regression is to fit a linear relationship among variables in a data set. Often, the goal is to find the most parsimonious model that can describe the majority of the data. In this work, we consider the situation where only a small portion of the data can be accurately modeled using linear regression. More generally, in many kinds of real-world data, portions of the data of significant size can be predicted significantly more accurately than by the best linear model for the overall data distribution: Rosenfeld et al. (2015) showed that there are attributes that are significant risk factors for gastrointestinal cancer in certain subpopulations, but not in the overall population. Hainline et al. (2019) demonstrated that a variety of standard (real-world) regression benchmarks have portions that are fit significantly better by a different linear model than the best model for the overall data set; Calderon et al. (2020) presented further, similar findings. We will consider cases where linear regression fits well when the data set is conditioned on a simple condition, which is unknown to us. We study the task of finding such a linear model, together with a formula on the data attributes describing the condition, i.e., the portion of the data for which the linear model is accurate. This problem was introduced by Juba (2017), who gave an algorithm for conditional sparse linear regression, using the maximum residual as the objective. This was extended to the usual squared-error loss (as well as other ℓ p losses) by Hainline et al. (2019). Juba (2017) also gave an algorithm for the general (non sparse) case that could only find a small fraction of the largest such condition. All of these algorithms find conditions describing subpopulations that are a union of some basic subsets of data, selected by "terms." For example, simple families of terms may be obtained by considering the data for which a small set of categorical attributes take some specific values, or based on whether the value of some real attributes lie in specific quantiles of the distribution. Calderon et al. (2020) gave an algorithm for non sparse linear regression that matches the size of the largest condition, but only under a new assumption, that the covariances of the data satisfying each of the terms of the condition have to all be very similar in spectral norm to the covariance of the overall condition.
Uniform covariances across terms is an extremely restrictive assumption: it means that essentially the only difference between populations selected by the various terms may be in their means. Note that the problem presupposes that there is significant heterogeneity in the conditional covariances across the distribution overall, or else the same linear model would be equally accurate across the various subsets; concretely, the risk factors found by Rosenfeld et al. (2015) are a kind of correlation between a factor and the target variable that exists in the identified subpopulation, but not in the larger population. In a real-world data set, there is no reason to expect that the only relationships that exist would involve the target attribute; relationships between other pairs of attributes may appear when we consider one term or another. For example, intuitively, if the data lies on a curved manifold, then conditioning on some attribute taking values that select one portion of the curve or another would alter such relationships since the tangent space changes, and the covariance matrix in any local region of the manifold only has eigenvectors lying in the tangent space. (Note that we only require a common component in the orthogonal subspace for a linear model to exist.) So, Calderon et al's algorithm can only be guaranteed to find highly homogeneous subsets of a distribution that features significant heterogeneity overall. In their work, Calderon et al. concluded with the question of whether or not this new property was necessary to obtain a computationally efficient algorithm.

Our Contribution
In this work, we answer the questions posed by Calderon et al. and Juba (2017), solving the form of the task originally sought: we give a polynomial-time algorithm that identifies a condition that covers as much of the distribution as the optimal condition and a linear model which provides a good fit when conditioned on said condition, even if the terms of the condition feature heterogeneous covariances. The only assumptions on the data we use are bounds on the moments of the data itself (including hypercontractivity) and generalizations of the standard Gaussian noise assumption on the subset of the data described by the unknown condition. Note that in regression, the error can be arbitrarily large with arbitrary probability, so bounds on the moments of the data are necessary to empirically estimate the error.
Our algorithm is inspired by the list-decodable subspace recovery algorithm presented by Bakshi and Kothari (2021). Their work uses the sum of squares method to construct an algorithm which addresses robust subspace recovery. We present an analogous algorithm for conditional linear regression. As in Calderon et al. (2020), this is done by using a collection of subsets, which we will call "terms," in place of individual points. We thus make use of the fact that by drawing many examples per term, the noise in the data selected by a term can be better controlled, leading to more accurate estimates. We stress that in contrast to the guarantee that Bakshi and Kothari obtain for their problem, we can obtain arbitrary accuracy with an algorithm that runs in fixed polynomial time, with an exponent that does not depend on the desired accuracy; we only require a sufficient (polynomial) number of examples from the target distribution, and our running time has only a low-order polynomial dependence on the size of the data set. (The dependence on the dimension, by contrast, while fixed, is a higher degree polynomial due to our use of the sum-of-squares method; reducing this dependence is a key challenge for future work, see Section 4.)

Related Work
Our work is both technically and conceptually related to "list-decodable" linear regression. Classical work in robust statistics (Huber, 1981;Rousseeuw and Leroy, 1987) considers situations where a minority subset of the data consists of "outliers" that should be ignored. In this classical setting, it did not make sense to consider the possibility that a majority of the data could be "outliers," in part because there would no longer be a unique, dominant solution to consider. But, a recently proposed model of "list-decodable" robust statistics (Charikar et al., 2017) (similarly to classical algorithms such as RANSAC (Fischler and Bolles, 1981) for subspace discovery) overcomes this obstacle by permitting a list of possible estimates or models to be produced, provided that the list is not too long (generally, O(1/µ) estimators for a µ-fraction of the data) and that an accurate estimate appears somewhere in the list. In particular, algorithms for list-decodable linear regression have been proposed (Raghavendra and Yau, 2020;Karmalkar et al., 2019) (see also Bakshi and Kothari (2021)). Although we have formulated our problem in such a way that we only produce a single arbitrary model as output, we could have returned a list of models as well (or vice-versa, select a suitable model from such a list). The main distinction is that in this line of work, on the one hand, one does not need to produce a DNF that identifies the inliers, in contrast to our setting. Note that without this formula, we cannot tell when we should use one of the models versus another to make predictions for new data. Of course, on the other hand, in these works one is not promised that such a DNF exists, either, and so the approach used in our analysis cannot be used in these problems.
Another similar line of work to conditional linear regression is selective regression (El-Yaniv and Wiener, 2012): here as well, the objective is to identify a fraction of the data that can be fit well. But in contrast to our setting, here a linear predictor is first identified, and then a data set is given, and finally a ranking of that data is produced. The interpretation is that the top µ-fraction of the ranking comprise the data for which the predictor is expected to be most accurate. In contrast, we jointly produce a linear model and a DNF that identifies which further examples, drawn from the same distribution, will be accurately predicted by the model. In learning with rejection or abstention (Cortes et al., 2016), on the other hand, a formula that selects a subset of the data is identified, but the problem formulation assigns a penalty to each example that is "rejected" -thus, we have a default loss value that our classifier can take in place of the loss that would be incurred by this prediction, and this overall loss is minimized over the entire data set. The cost of rejection here takes the place of the probability of the subset µ in our problem.
All of these works have some similarities to classical topics such as fitting mixture models (McCulloch and Searle, 2001;Jiang, 2007). The primary difference is that in such work, first, every data point should have been drawn from some linear model in the mixture; if some large subset of the data cannot be fit well by linear models, there is no guarantee that the model will identify a small subset that can be fit well. A second difference is that such models do not provide a (DNF) rule to decide whether or not subsequent data is drawn from one of the components versus another. There are a number of topics such as regression trees (Quinlan, 1992), cluster-wise regression (Park et al., 2017), etc. that do provide such rules, but again, if the data overall cannot be fit well, they do not guarantee that small subsets of the data that can be fit well will be found.

PRELIMINARIES
We will assume that we have a data set containing N samples, from a distribution D over {0, 1} n × R d × R. Each sample consists of an n dimensional vector of Boolean attributes x, a d dimensional real valued vector of predictor variables y, and a real valued response, z, which we would like to predict. We will denote the ith sample as (x, y, z) (i) and abbreviate it as x (i) when there is no ambiguity.
For linear regression, we want to find a a vector of coefficients v such that z can be predicted by v, y . Typically, v is found using ordinary least squares, which minimizes the sum of ( v, y − z) 2 over all data points. However, since we are interested in cases where the majority of data cannot be fit, we want to find a subset of points described by condition c where there exists a good linear model. Similar to previous work, we will consider conditions represented by Disjunctive Normal Form (DNF) formulas; other natural families of formulas are either weaker or yield intractable problems (Juba, 2017). A k-DNF is defined to be a disjunction (OR) of terms where each term is a conjunction (AND) of no more than k attributes.
Throughout this paper, we will use · F to denote the Frobenius norm and · 2 to denote the L2 norm of a vector. We will also define X I as the characteristic function where X I (x) = 1 if x ∈ I and 0 otherwise. For brevity, we will use [N ] = {n ∈ N|1 ≤ n ≤ N }. For a matrix M , M 0 denotes M is positive semidefinite. Finally, we will use Π to denote projection matrices.
Definition 2.1 (Conditional Linear Regression). Given a sample of N points, (x, y, z) (i) , from a distribution D over {0, 1} n × R d × R, the task of conditional linear regression is to find a k-DNF, c, and linear predictor, v = (v 1 . . . v d ) T , such that, with high probability, | v, y − z| is bounded by ǫ when conditioned on c(x) = 1 and c(x) = 1 is satisfied by at least a µ fraction of the data.
We will present an algorithm that finds c and v that is close to the optimal values of c * and v * given that the distribution of samples conditioned on c * follows certain regularity conditions. Our algorithm is obtained by solving a sum-of-squares relaxation (Parrilo, 2000;Lasserre, 2001;Nesterov, 2000;Shor, 1987) of a polynomial optimization problem: Definition 2.2 (Sum of squares relaxation). Given a system of polynomial inequalities for polynomials in R[x 1 , . . . , x n ], g 1 (x) ≥ 0, . . . , g r (x) ≥ 0, h 1 (x) = 0, . . . , h s (x) = 0, the degree-ℓ sum-ofsquares relaxation is the following semidefinite optimization problem.
The set of program variables u is indexed by monomials over x 1 , . . . , x n of total degree at most ℓ, with u α denoting the variable for the monomial with degree vector α ∈ N n .
For a "shift" polynomial p ∈ R[x 1 , . . . , x n ] of degree t, letting p γ denote the coefficient of the monomial x γ in p, the degree-ℓ localizing matrix M ℓ (pu) is defined by M ℓ (pu) (α,β) = γ p γ u α+β+γ for α and β of total degree at most ℓ/2 − t. Now, the program has the constraints that Given a bound on the magnitudes of the values involved, the solutions to semidefinite programs can be approximated to arbitrary precision in polynomial time by various algorithms; the current state-of-the-art is due to Jiang et al. (2020).
A helpful interpretation of the sum-of-squares relaxation is that it defines a set of "pseudodistributions" that relax the moments of probability distributions supported on solutions to the system of inequalities, with an associated "pseudo-expectation" operator defined on polynomials of degree up to ℓ (we borrow the presentation from Raghavendra and Yau (2020)): A low-degree pseudo-distribution is generally not a probability distribution, and we generally cannot sample from it. The quality of a sum-of-squares relaxation is characterized by "sum-ofsquares proofs"-under mild conditions, the bounds obtained by the relaxation of a given degree match the optimal bound that can be proved via a sum-of-squares proof (Parrilo, 2000;Lasserre, 2001;Nesterov, 2000;Shor, 1987): Definition 2.5 (Sum of Squares proofs (Grigoriev and Vorobjov, 2002)). Fix a set of polynomial Moreover, a recent technique for extracting solutions from the sum-of-squares relaxation has been to use identifiability (Raghavendra et al., 2019): roughly, if there is a sum-of-squares proof that a portion of the solution is determined up to small ℓ 2 error, then we can read that portion of the solution off directly from the degree-1 variables. Conceptually, our analysis will follow this approach. Towards such an analysis, we require that the data/noise possess some niceness properties in a form that is amenable to sum-of-squares: Definition 2.6 (Certifiable Hypercontractivity, Section 2.1 of Bakshi and Kothari (2021) We say that the hypercontractivity is ℓ-certifiable if there is a degree-ℓ sum-of-squares proof of this inequality with Q as an indeterminate.
In linear regression, it is standard and convenient to assume Gaussian residuals with mean 0 and constant variance σ 2 . We stress that in order to empirically estimate the squared error of candidate linear models, we need some control on the moments of the loss function, which for various candidates corresponds to various quadratic polynomials; hypercontractivity of such quadratic polynomials gives such bounds. Since many distributions, including uniform distributions, Gaussians, products of subgaussian distributions etc., are certifiably hypercontractive (Bakshi and Kothari, 2021), it is reasonable to assume that we have samples from a certifiably hypercontractive distribution. Lemma 2.7 shows that empirical distributions from a certifiably hypercontractive distribution are also certifiably hypercontractive (Bakshi and Kothari, 2021). We also assume that degree-2 polynomials similarly have certifiably bounded variances, which also holds for such distributions (Bakshi and Kothari, 2021).

Preprocessing
For our algorithm, we will consider the data as m disjoint subsets which we will call terms, {I j } m j=1 . To find k-DNF conditions, we will create a term for each setting of each set of k distinct attributes, but the algorithm can be used with any family of sets we choose. We will weight each term by the number of points, |I j |, in term I j . We will define I good as the collection terms of the optimal k-DNF c * . From the perspective of the subspace recovery task (Bakshi and Kothari, 2021), I good represents the collection of inliers, and the remaining points represent the outliers. Since the condition c * is satisfied by a µ fraction of the data, |I good | = µN . Additionally, our algorithm will assume that each term I j is pairwise disjoint. This can be achieved by duplicating points that satisfy more than one term (Calderon et al., 2020). We have N data points and m terms, so there will be at most mN data points following the duplication procedure. Previously, |I good | = | Ij ∈I good I j | which increases to |I good | = Ij ∈I good |I j | with duplicate points. Therefore size of I good may blow up at most by a factor of the number of terms in I good . If we use N ′ to denote the number of points after duplication, notice that N ′ good /N ′ ≥ N good /mN . Thus the proportion of good points after duplication is at least a µ/m (Calderon et al., 2020). Hereafter, N , I good , {I j } m j=1 , and µ will be used to describe the data after duplication.
Our algorithm only obtains good estimates when each term in I good is large. Since the contribution of each term is weighted by its size, we can remove sufficiently small terms without compromising the quality of the estimates. Finally, we will extend y (i) by a constant, 1, to allow for an intercept in the linear model.

Main algorithm: identifying the data subspace
Recall that our goal is to find the linear predictor v such that v, y (i) = z (i) for all points satisfying condition c * . If we extend v with -1 and y (i) with z (i) , the previous equation is equivalent to v, y (i) = 0. This equation describes a hyperplane, which allows us to look at our conditional linear regression problem from the perspective of a subspace recovery problem. Note that both v and y (i) are now (d + 2)-dimensional vectors. Let us use Π to denote a projection matrix that projects onto this subspace and Π * to denote the optimal projection.
and predefined disjoint subsets, {I j } m j=1 . Let I good be an unknown collection of terms containing points that satisfy c * such that P (x ∈ I good ≥ µ). If there exists a linear predictor, v * , such that | v * , y | ≤ ǫ and a projection matrix, Π * , that projects onto the hyperplane described by the linear predictor, then we want to findΠ that approximates Π * such that Π − Π * F ≤ error.
Note that the Frobenius norm bounds the spectral norm from above, where the squared error of v is precisely v ⊤ Π * v, so a solution to this reformulation gives an empirical estimate that is off by at most v 2 2 · error. The advantage of this reformulation is that the subspace is uniquely determined by the set of terms, whereas if the subspace has lower dimension, there may be a large space of linear models, which are thus not identifiable. We remark that Karmalkar et al. (2019); Raghavendra and Yau (2020) made a stronger assumption to obtain identifiability for regression, that the data was anti-concentrated. Anticoncentration does not hold if, e.g., the data lies on a lowerdimensional manifold, so it is a restrictive assumption. By contrast, we only need hypercontractivity and bounded moments to estimate the subspace. We now state our main theorem, which is adapted from Theorem 1.4 of Bakshi and Kothari (2021).
Theorem 3.2. Let Π * be a projection matrix for a subspace of dimension r. Let D|c * be a distribution with mean 0, covariance Π * , and with 2-certifiably C-hypercontracitve degree-2 polynomials with certifiably C-bounded variance. Then, there exists an algorithm that takes n ≥ Ω (d log d/µ) 16 samples from the distribution conditioned on condition c * and outputs a list L of O(1/µ) projection matrices such that with probability at least 0.99 over the draw of the sample and randomness of the algorithm, there is aΠ ∈ L satisfying Π − Π * F ≤ O(1/µ) in polynomial time.
Let A w,v,ǫ,Π be the SoS program in Figure 1 where w, v, ǫ, and Π 1 , . . . , Π m are indeterminates. We interpret the program constraints as follows: 1. Defines y ′(i) such that the inliers have mean 0 while preserving the constant 1 for the intercept of the model. 2. The linear predictor v fits well when conditioned on c.
3. The number of samples when conditioned on c comprises a µ fraction of the data. 4. The Boolean constraint: w j ∈ {0, 1} for all j. 5. Defines Π j as a projection matrix corresponding to the distribution of points in term I j . 6. The residuals of the linear model follow a Gaussian distribution with mean 0 and standard deviation σ, which bounds the average noise on the inliers. 7. The samples are certifiably hypercontractive. 8. Similarly, the variance is certifiably bounded. 9. The second moment of each predictor variable of samples satisfying c is bounded by α. 10. Finally, the fourth moment of each predictor variable of samples satisfying c is bounded by β. We note that constraints 7 and 8 are infinite families of constraints. But since there are sum-ofsquares proofs of these constraints, similar to Bakshi and Kothari (2021), we can use the quantifier elimination technique (Fleming et al., 2019, Section 4.3.4) to rewrite these as standard constraints. We note that we substitute the cubic polynomial y ′⊤ Qy ′ for Q in the proof of the hypercontractive inequality, thus obtaining a degree-6 sum-of-squares proof from the original degree-2 proof.
For the following analysis, y (i) will denote the samples after centering the inliers, which is equivalent to y ′(i) from the program.  Proof. The squared Frobenius norm is equivalent to summing the square of each element. Thus we will bound the square of each element using Chebyshev's inequality and then compute the sum. The matrix 1 µN N i=1 w j X Ij (x (i) )(x (i) )Π j y (i) y (i)⊤ Π j can treated as a random quantity that represents an empirical estimate of w(I j )Π j . Recall The empirical covariance matrix can be rewritten as

For each
.
3. Take J to be a random multi-set formed by Let us use the notation A r,s to denote the element in row r and column s of a matrix A. Using Chebyshev's Inequality, the square of each element in the Frobenius norm can be bounded by with probability 1 − δ/(d + 2) 2 . Since y (i) has been extended with one, let us treat the first element of y (i) as 1. Therefore, when r = 1 or s = 1, For all other pairs of r and s, the variance is bounded by β + 1α 4 . Due to the same coordinate of y (i) being 1 for all i ∈ [N ], the second moment, α, must be at least 1. Therefore, Tying it all together, the square of each element can be bounded by w j |I j |w j (d + 2) 2 µ 2 N 2 δ β + α + 7α 4 σ 2 with probability 1 − δ/(d + 2) 2 . This is a (d + 2) dimensional square matrix so there are (d + 2) 2 elements. Therefore, the squared Frobenius norm is bounded by w j |I j |w j (d + 2) 4 µ 2 N 2 δ β + α + 7α 4 σ 2 and by a union bound, this holds with probability at least 1 − δ.
Lemma 3.3, which uses Chebyshev's Inequality, gives us the nice property of being inversely proportional to N 2 . Therefore as we draw more data, the right hand side goes to 0.
The proofs of Lemmas 3.4 and 3.5 are similar to Lemmas 4.5 and 4.6 in Bakshi and Kothari (2021), but by using Lemma 3.3 we are able to form a tighter bound.
with probability at least 1 − kδ where w(I j ) = 1 µN N i w j X Ij x (i) . Lemma 3.6 admits a sum-of-squares proof and is proved by using Lemmas 3.4 and 3.5. The full proofs for all new lemmas are included in the supplementary material.
Lemma 3.7 (Large weight on inliers from high-entropy constraints. Fact 4.4 in Bakshi and Kothari (2021) and Lemma 3.1 in Raghavendra and Yau (2020)). LetẼ ξ be a pseudo-distribution of degree ≥ 2 that satisfies A w,v,ǫ,Π and minimizes Ẽ ξ m j=1 The main theorem, Theorem 3.2, can be obtained by using Lemmas 3.6 and 3.7 and following a similar argument to Theorem 1.4 of Bakshi and Kothari (2021).
Proof. Given a distribution D that is certifiably hypercontractive, Lemma 2.7 implies that a large enough sample of inliers will also be certifiably hypercontractive with high probability. Algorithm 1 finds a pseudo-distribution that satisfies A w,v,ǫ,Π and minimizes k j=1 w j x (i) ∈Ij X Ij (x (i) ). Then µN δ by using Lemma 3.6. By applying Jensen's Inequality and taking the square root, we have Due to the definition ofΠ i from the algorithm, we can rewrite the inequality as . Then, from Lemma 3.7, Z ≥ µ ⇒ 1 Z ≤ 1 µ . Dividing by Z on both sides thus yields: it follows that i ∈ I good with probability at least 1 µn k j w j i∈Ij E[X Ij (x i )] = Z ≥ µ. By Markov's inequality applied to the last equation, with probability 1 2 over the choice of i conditioned on i ∈ I good , Thus, in total, with probability at least µ/2, Therefore, with probability of at least 0.99 over the draw of the random set J, the list constructed by the algorithm containsΠ such that For the running time, the SDP can be solved in polynomial time and dominates the running time. Therefore, the algorithm runs in polynomial time overall.
For a projection Π and a linear predictor, v, which satisfies v, y (i) = 0 for all i ∈ I good when disregarding noise, v ⊤ Π = 0 ⊤ . Therefore, for each candidate Π, we can recover a candidate linear predictor u by treating u as a solution to the linear system u ⊤ Π = 0 ⊤ .

Obtaining a k-DNF Condition
Once we obtain an approximation forv, we will use the method described by Calderon et al. (2020) to obtain a k-DNF condition for c. Since we have used the same reductions, specifically the reduction to disjoint terms by duplicating points, we will be able to invoke their analysis of their method. Let us define a loss function, f (i) : H → R, for each point as Let f Ij (v) be a the average loss over points in I j . Finally, let us definef (v) = E[f I good (v)] as the expectation of average loss for points in I good . We will search through the Boolean data space {x} for conditions c for each u. When we find a pair (u, c) such that f c (u) is small and there are are enough points satisfying c, then we return that pair as a solution. For some constant accuracy parameter γ and Lipschitz constant Recall that during preprocessing, we duplicated points that satisfied more than 1 term so that each term is disjoint. For m terms, each point can be copied at most m times. Let t be the number of terms in c.  2020)). Let u be such that u − v * < γ. Then |f (u)| ≤ t(γ + ǫ).
After obtaining u such that f i (u) is close to f i (v * ), Calderon et al. (2020) use a greedy setcover algorithm to find the corresponding conditions c. The algorithm greedily chooses terms I j satisfying i∈Ij f (i) (u) ≤ (1 + γ)µǫN to maximize the number of additional points in I j that did not satisfy the previously chosen terms. It iterates until the number of points satisfying the chosen terms is at least (1 − γ/2)µN (Calderon et al., 2020).
Lemma 3.9 (Lemma 3.5 in Calderon et al. (2020)). If there exists an optimal k-DNF c * that is satisfied by a µ-fraction of the points with total loss ǫ, then, the weighted greedy set cover algorithm finds a k-DNFĉ that is satisfied by a (1 − γ)µ-fraction of the points with total loss O(t log(µN )ǫ).
Thus, we can obtain a pair (u, c) that gives empirical error that is only greater than the optimal by a O(t log(µN )) factor. Given that our assumed bounds on the moments of the data distribution implies that the square of the loss (being a quadratic polynomial) is bounded, we can use the bounds of Cortes et al. (2013) to bound the generalization error of linear regression on each possible k-DNF, to thus obtain that the true generalization error is similarly bounded. For y ∈ B ⊆ R d , where B is the l 2 radius of B, this only incurs a polynomial increase in the sample complexity in B, t, and 1/γ overall.

DISCUSSION AND FUTURE DIRECTIONS
We have thus shown that the assumption of homogeneous covariances used by Calderon et al. (2020) is not needed to obtain a polynomial-time algorithm for conditional linear regression. On the other hand, although we obtain a polynomial running time and sample complexity, the exponents are quite large. In particular, the sample complexity we obtain in our analysis is far from optimal for this problem. Thus, our algorithm is impractical in its current form and our contribution is strictly theoretical. The main direction for future work is to develop a practical algorithm that does not require the homogeneous covariance assumption.
We now elaborate on the obstacles. Much of the overhead in the sample complexity arises from the use of the certifiable hypercontractivity assumption, which requires relatively high-degree polynomial expressions. Another source of sub-optimality seems to arise from the way we invoke bounds obtained via Chebyshev's inequality on the Frobenius error of the projector onto the data subspace. We conjecture that this can be improved by a more refined analysis, but it is still not clear whether or not certifiable hypercontractivity -the dominant source of overhead for most purposes -is really necessary.
The overhead in the computational complexity arises from the use of the sum-of-squares semidefinite program relaxation. Although the degree of the relaxation is moderate, the technique unfortunately yields algorithms that solve a large semidefinite program, and thus inherently tends to simply scale poorly. In many cases, however, the development of a sum-of-squares algorithm has led to the subsequent development of a spectral algorithm that can be practical. Examples of this sequence include a number of robust estimation tasks Hopkins et al., 2019;Diakonikolas et al., 2020;Depersin, 2020) and tasks to identify hidden structures in a data set (Hopkins et al., 2016(Hopkins et al., , 2017, that are variously of similar flavor to our problem.

A Sum of Squares
Lemma A.1 (SoS Hölder's Inequality. Fact 3.11 in Bakshi and Kothari (2021)). Let w 1 , . . . , w n be indeterminates and let f 1 , . . . , f n be polynomials of degree m in vector valued variable x. Let k be a power of 2. Then, Lemma A.2 (SoS Almost Triangle Inequality. Fact 3.8 in Bakshi and Kothari (2021)). Let a, b be indetermnates. Then, for any t ∈ N, Lemma A.3 (Cancellation within SoS. Fact 3.12 in Bakshi and Kothari (2021)). Let a be an indeterminate. Then, Lemma A.4 (Soundness. Fact 3.4 in Bakshi and Kothari (2021)). If D r A for a level-ℓ. pseudodistribution D and there exists a sum-of-squares proof A r ′ B, then D r·r ′ +r ′ B.
Lemma A.5 (SoS Cauchy Schwarz. Fact 2.4 in Raghavendra and Yau (2020)). Let x 1 , . . . , x n , y 1 , . . . , y n be indeterminates, then  (2021) and Lemma 3.1 in Raghavendra and Yau (2020)). LetẼ ξ be a pseudodistribution of degree ≥ 2 that satisfies A w,v,ǫ,Π and minimizes Ẽ ξ m j=1 Proof. (This proof is the same as Lemma 3.1 in Raghavendra and Yau (2020).) For the sake of of simplicity, let w i = i∈Ij w j X Ij x (i) and note that w i ∈ {0, 1}. LetẼ P denote a pseudodistribution corresponding to the actual assignment {w ′ i } i∈ [N ] and letẼ D be the pseudo-expectation that minimizes E D [w] . For a constant κ ∈ [0, 1], define the pseudo-expectationẼ R as a mixture ofẼ P andẼ D .Ẽ We can use the definition ofẼ R to expand the left hand side.
By rearranging the terms, we get As κ → 0, the right hand side tends to µ 2 N .

Lemma B.2 (Lemma 3.3 -Frobenius Closeness of Empirical and True Covariances).
Let ǫ i = 0 d+1 ǫ i and w(I j ) = 1 µN w j X Ij (x (i) ). Then with probability 1 − δ, Proof. The squared Frobenius norm is equivalent to summing the square of each element. Thus we will bound the square of each element using Chebyshev's inequality and then compute the sum. The matrix 1 µN N i=1 w j X Ij (x (i) )(x (i) )Π j y (i) y (i)⊤ Π j can treated as a random quantity that represents an empirical estimate of w(I j )Π j . Recall from the program. Thus The empircal covariance matrix can be rewritten as Let us use the notation A r,s to denote the element in row r and column s of a matrix A. Using Chebyshev's Inequality, the square of each element in the Frobenius norm can be bounded by due to nonnegativity. By using the SoS Cauchy Schwarz and SoS AM-GM inequalities, Lemmas A.5 and A.6, Since y (i) has been extended with one, let us treat the first element of y (i) as 1. Therefore, when r = 1 or s = 1, For all other pairs of r and s, the variance is bounded by β + 1α 4 . Due to the same coordinate of y (i) being 1 for all i ∈ [N ], the second moment, α, must be at least 1. Therefore, Tying it all together, the square of each element can be bounded by wj |Ij |wj (d+2) 2 µ 2 N 2 δ β + α + 7α 4 σ 2 with probability 1 − δ/(d + 2) 2 . This is a (d + 2) dimensional square matrix so there are (d + 2) 2 elements. Therefore, By taking the union bound, this holds with probability at least 1 − δ.  (2021)).
Proof. For a (d + 2) × (d + 2) matrix-valued indeterminate Q and using the SoS Hölder's Inequality, Lemma A.1, we have Using certifiable hypercontractivity combined with the bounded variance constraints, we have By combining Equations 1 and 2 and substituting Q = 1 By using Lemma B.2, we have with probability at least 1 − δ.
Proof. This follows the same proof as Lemma B.3.
with probability at least 1 − kδ where w( Proof. Define w * (I j ) = By dividing both sides of the inequality by w(I j ) 2 , we have By using Cancellation within SoS, Lemma A.3 and multiplying both sides of the inequality by w(I j ) 1/2 , we have Since w(I j ) = 1 µN N i=1 w j X Ij (x (i) ) and is bounded above by w * (I j ), then with probability 1 − 2δ: Define Π and Π * as a weighted average of Π j and Π j * respectively, where the weights are proportional to |I j |. Thus Π = m j=1 w(I j )Π j and Π * = m j=1 w(I j )Π j * . By summing both sides of the inequality over all j ∈ [m], we have By using the Cauchy-Schwarz Inequality and the triangle inequality to rewrite the left hand side, we have Lemma B.6 (Lemma 2.7 -Certifiable Hypercontractivity Under Sampling. Lemma 6.11 in Bakshi and Kothari (2021)). Let D be a 1-subgaussian, 2h-certifiably C-hypercontractive distribution over R d . Let S be a set of n = Ω((hd) 8h ) i.i.d. samples from D. Then, with probability at least 1 − 1/poly(n), the uniform distribution on S is h-certifiably (2C)-hypercontractive.

B.1 Proof of Main Theorem
Theorem B.7 (Main Theorem -Theorem 3.2). Let Π * be a projection matrix for a subspace of dimension r. Let D be a distribution with mean 0, covariance Π * , and 2-certifiably C-hypercontracitve degree-2 polynomials. Then, there exists an algorithm that takes n ≥ Ω (d log d/µ) 16 samples from the distribution conditioned on condition c * and outputs a list L of O(1/µ) projection matrices such that with probability at least 0.99 over the draw of the sample and randomness of the algorithm, there is aΠ ∈ L satisfying Π − Π * F ≤ O(1/µ) in polynomial time.
Proof. This follows the same proof as Theorem 1.4 in Bakshi and Kothari (2021). Since D is certifiably C-hypercontractive, Fact B.6 implies that ≥ n = Ω(d log d/µ) 16 samples suffice for the uniform distribtution on the inliers, I good , to have 2-certifiably C-hypercontractive degree 2 polynomials with probabiliy at least 1 − 1/d. Let ξ 1 be the event that this succeds, and condition on it.
Letμ be a pseudo-distribution of degree-24 satisfying A w,v,ǫ,Π , minimizing m j=1 w j N i=1 X Ij (x (i) ) as described in Algorithm 1. Observe that such a pseudo-distribution is guaranteed to exist: take the pesudo-distribution supported on a single point, (w, Π) such that w i = 1 iff i ∈ I good and Π = Π * . It is straight forward to check that Π * is indeed a rank r projection matrix and rank d + 2 projection matrix and m j=1 N i=1 X Ij (x (i) ) = µN . Conditioned on ξ 1 , the hypercontractivity constraint is also satisfied by the inliers.
Since Lemma B.5 admits a sum-of-squares proof, it follows from Fact A.4 that the polynomial inequality is preserved under pseudo-expectations.
]. Then, from Lemma B.1, Z ≥ µ ⇒ 1 Z ≤ 1 µ . Dividing by Z on both sides thus yields: Since each index i ∈ [N ] is chosen with probability it follows that i ∈ I good with probability at least 1 µn k j w j i∈Ij E[X Ij (x i )] = Z ≥ µ. By Markov's inequality applied to the last equation, with probability 1 2 over the choice of i conditioned on i ∈ I good , Π i − Π * F ≤ 2 µ mC 2 5 w j (d + 2) 4 (β + α + 7α 4 σ 2 ) µN δ .
Therefore, with probability of at least 0.99 over the draw of the random set J, the list constructed by the algorithm containsΠ such that Π i − Π * F ≤ 2 µ mC 2 5 w j (d + 2) 4 (β + α + 7α 4 σ 2 ) µN δ . Now to account for the running time of the algorithm, the SDP for the program can be solved in polynomial time, so the algorithm runs in polynomial time overall.