Statistical Query Lower Bounds for Robust Estimation of High-dimensional Gaussians and Gaussian Mixtures

Introduction

For the unsupervised estimation problems considered here, the input is a probability distribution which is accessed via a sampling oracle, i.e., an oracle that provides i.i.d. samples from the underlying distribution. Statistical Query (SQ) algorithms are a restricted class of algorithms that are only allowed to query expectations of functions of the distribution rather than directly access samples. This class of algorithms is quite broad: a wide range of known algorithmic techniques in machine learning are known to be implementable using SQs. These include spectral techniques, moment and tensor methods, local search (e.g., Expectation Maximization), and many others (see, e.g., [CKL+06, FGR+13] for a detailed discussion). Moreover, for the unsupervised learning problems studied in this paper, all known algorithms with non-trivial performance guarantees are SQ or are easily implementable using SQs.

A number of techniques have been developed in information theory and statistics to characterize the sample complexity of inference tasks. These involve both techniques for proving sample complexity upper bounds (e.g., VC dimension, metric/bracketing entropy) and information-theoretic lower bounds (e.g., Fano and Le Cam methods). On the other hand, computational lower bounds have been much more scarce in the unsupervised setting. Perhaps surprisingly, it is possible to prove unconditional lower bounds on the computational complexity of any SQ algorithm that solves a given learning problem. Given the ubiquity and generality of SQ algorithms, an SQ lower bound provides strong evidence of the problem’s computational intractability.

In this paper, we describe a general technique that yields the first Statistical Query lower bounds for a range of fundamental high-dimensional learning problems involving Gaussian distributions. Such problems are ubiquitous in applications across the data sciences and have been intensely investigated by different communities of researchers for several decades. Our main results are for the problems of (1) learning Gaussian mixture models (GMMs), and (2) robust (agnostic) learning of a single unknown Gaussian distribution. In particular, we show a super-polynomial gap between the (information-theoretic) sample complexity and the computational complexity of any Statistical Query algorithm for these problems. In more detail, our SQ lower bound for Problem (1) is qualitatively matched by known learning algorithms for GMMs (all of which can be implemented as SQ algorithms). For Problem (2), we give a new (SQ) algorithm in this paper whose running time nearly matches our SQ lower bound.

Our SQ lower bounds are attained via a unified moment-matching technique that is useful in other contexts and may be of broader interest. Our technique yields nearly-tight lower bounds for a number of related unsupervised estimation problems. Specifically, for the problems of (3) robust covariance estimation in spectral norm, and (4) robust sparse mean estimation, we establish a quadratic statistical–computational tradeoff for SQ algorithms, matching known upper bounds.

Finally, we use our technique to obtain tight sample complexity lower bounds for high-dimensional testing problems. Specifically, for the classical problem of robustly testing an unknown mean (known covariance) Gaussian, our technique implies an information-theoretic lower bound that scales linearly in the dimension. This lower bound matches the sample complexity of the corresponding robust learning problem and separates the sample complexity of robust testing from standard (non-robust) testing. This separation is surprising because such a gap does not exist for the corresponding learning problem.

Before we discuss our contributions in detail, we provide the necessary background for the Statistical Query model and the unsupervised estimation problems that we study.

A Statistical Query (SQ) algorithm relies on an oracle that given any bounded function on a single domain element provides an estimate of the expectation of the function on a random sample from the input distribution. This computational model was introduced by Kearns [Kea98] in the context of supervised learning as a natural restriction of the PAC model [Val84]. Subsequently, the SQ model has been extensively studied in a plethora of contexts (see, e.g., [Fel16b] and references therein).

A recent line of work [FGR+13, FPV15, FGV15, Fel16a] developed a framework of SQ algorithms for search problems over distributions – encompassing the distribution estimation problems we study in this work. It turns out that one can prove unconditional lower bounds on the computational complexity of SQ algorithms via the notion of Statistical Query dimension. This complexity measure was introduced in [BFJ+94] for PAC learning of Boolean functions and was recently generalized to the unsupervised setting [FGR+13, Fel16a]. A lower bound on the SQ dimension of a learning problem provides an unconditional lower bound on the computational complexity of any SQ algorithm for the problem.

Remark. We would like to emphasize here that the SQ lower bounds shown in this paper apply to the running time of an SQ algorithm and not on its sample complexity (when we simulate the SQ algorithm by drawing samples to answer its SQ queries). Specifically, for all learning problems considered in this paper, there exist straightforward SQ algorithms (that can be simulated with sample access to the distribution) with near-optimal sample complexity, albeit with exponential running time. Specifically, lower bounds on the SQ dimension of the corresponding problems establish lower bounds on the running time of any SQ algorithm for the problem – not on its sample complexity.

In the preceding paragraphs, we were working under the assumption that the unknown distribution generating the samples is exactly a mixture of Gaussians. The more general and realistic setting of robust (or agnostic) learning – when our assumption about the model is approximately true – turns out to be significantly more challenging. Specifically, until recently, even the most basic setting of robustly learning an unknown mean Gaussian with identity covariance matrix was poorly understood. Without corruptions, this problem is straightforward: The empirical mean gives a sample-optimal efficient estimator. Unfortunately, the empirical estimate is very brittle and fails in the presence of corruptions.

Agnostically learning a single high-dimensional Gaussian is arguably the prototypical problem in robust statistics [Hub64, HRRS86, HR09]. Early work in this field [Tuk75, DG92] studied the sample complexity of robust estimation. Specifically, for the case of an unknown mean and known covariance Gaussian, the Tukey median [Tuk75] achieves $O(\epsilon)$-error with $O(n/\epsilon^{2})$ samples (see, e.g., [CGR15] for a simple proof). Since $\Omega(n/\epsilon^{2})$ samples are information-theoretically necessary – even without noise – the robustness requirement does not change the sample complexity of the problem.

The computational complexity of agnostically learning a Gaussian is less understood. Until recently, all known polynomial time estimators could only guarantee error of $\Theta(\epsilon\sqrt{n})$. Two recent works [DKK+16, LRV16] made a first step in designing robust polynomial-time estimators for this problem. The results of [DKK+16] apply in the standard agnostic model; [LRV16] works in a weaker model – known as Huber’s contamination model [Hub64] – where the noisy distribution $D$ is of the form $(1-\epsilon)G+\epsilon N$, where $N$ is an unknown “noise” distribution. For the problem of robustly estimating an unknown mean Gaussian $N(\mu,I)$, [LRV16] obtains an error guarantee of $O(\epsilon\sqrt{\log n})$, while [DKK+16] obtains error $O(\epsilon\sqrt{\log(1/\epsilon)})$, independent of the dimensionThe algorithm of [LRV16] can be extended to work in the standard agnostic model at the expense of an increased error guarantee of $O(\epsilon\sqrt{\log n\log(1/\epsilon)})$..

A natural and important open problem, put forth by these works [DKK+16, LRV16], is the following:

A statistical–computational tradeoff refers to the phenomenon that there is an inherent gap between the information-theoretic sample complexity of a learning problem and its computational sample complexity, i.e, the minimum sample complexity attainable by any polynomial time algorithm for the problem. The prototypical example is the estimation of a covariance matrix under sparsity constraints (sparse PCA) [JL09, CMW13, CMW15], where a nearly-quadratic gap between information-theoretic and computational sample complexity has been established (see [BR13b, WBS16b]) – assuming the computational hardness of the planted clique problem.

For a number of high-dimensional learning problems (including the problem of robustly learning a Gaussian under the total variation distance), it is known that the robustness requirement does not change the information-theoretic sample complexity of the problem. On the other hand, it is an intriguing possibility that injecting noise into a high-dimensional learning problem may change its computational sample complexity.

Does robustness create inherent statistical–computational tradeoffs for natural high-dimensional estimation problems?

In this work, we consider two natural instantiations of the above general question: (i) robust estimation of the covariance matrix in spectral norm, and (ii) robust sparse mean estimation. We give basic background for these problems in the following paragraphs.

Is there a computationally efficient robust covariance estimator in spectral error that uses a strongly sub-quadratic sample size, i.e., $O(n^{2-c})$ for a constant $0<c<1$?

Is there a computationally efficient robust $k$-sparse mean estimator that uses a strongly sub-quadratic sample size , i.e., $O(k^{2-c})$ for a constant $0<c<1$?

It is conjectured in [Li17] that a quadratic gap is in fact inherent for efficient algorithms.

So far, we have discussed the problem of learning an unknown distribution that is promised to belong (exactly or approximately) in a given family (Gaussians, mixtures of Gaussians). A related inference problem is that of hypothesis testing [NP33, LR05]: Given samples from a distribution in a given family, we want to distinguish between a null hypothesis and an alternative hypothesis. Starting with [GR00, BFR+00], this broad question has been extensively investigated in TCS with a focus on discrete probability distributions. A natural way to solve a distribution testing problem is to learn the distribution in question to good accuracy and then check if the corresponding hypothesis is close to one satisfying the null hypothesis. This testing-via-learning approach is typically suboptimal and the main goal in this area has been to obtain testers with sub-learning sample complexity.

In this paper, we study natural hypothesis testing analogues of the high-dimensional learning problems discussed in the previous paragraphs. Specifically, we study the sample complexity of (i) robustly testing an unknown mean Gaussian, and (ii) testing a GMM.

Robust hypothesis testing is of fundamental importance and has been extensively studied in robust statistics [HR09, HRRS86, Wil97]. Perhaps surprisingly, it is poorly understood in the most basic settings, even information-theoretically. Specifically, the sample complexity of our aforementioned robust mean testing problem has remained open. It is easy to see that the tester of Appendix C fails in the robust setting. On the other hand, the testing-via-learning approach implies a sample upper bound of $O(n/\epsilon^{2})$ for our robust testing problem – by using, e.g., the Tukey median. The following question arises:

Is there an information-theoretic gap between robust testing and non-robust testing? What is the sample complexity of robustly testing the mean of a high-dimensional Gaussian?

2 Our Results

The main contribution of this paper is a general technique to prove lower bounds for a range of high-dimensional estimation problems involving Gaussian distributions. We use analytic and probabilistic ideas to construct explicit families of hard instances for the estimation problems described in Section 1.1. Using our technique, we prove super-polynomial Statistical Query (SQ) lower bounds that answer Questions 1.1 and 1.2 in the negative for the class of SQ algorithms. We also show that the observed quadratic statistical–computational gap for robust sparse mean estimation and robust spectral covariance estimation is inherent for SQ algorithms. As an additional important application of our technique, we obtain information-theoretic lower bounds on the sample complexity of the corresponding testing problems. (We note that our testing lower bounds apply to all algorithms.) Specifically, we answer Question 1.4 in the affirmative, by showing that the robustness requirement makes the Gaussian testing problem information-theoretically harder. In the body of this section, we state our results and elaborate on their implications and the connections between them.

Our first main result is a lower bound of $n^{\Omega(k)}$ on the complexity of any SQ algorithm that learns an arbitrary $n$-dimensional $k$-GMM to constant accuracy (see Theorem 4.1 for the formal statement):

At a conceptual level, Theorem 1.1 implies that – as far as SQ algorithms are concerned – the computational complexity of learning high-dimensional GMMs is inherently exponential in the dimension of the latent space – even though there is no such information-theoretic barrier in general. Our SQ lower bound identifies a common barrier of the strongest known algorithmic approaches for this learning problem, and provides a rigorous explanation why a long line of algorithmic research on this front either relied on strong separation assumptions or resulted in runtimes of the form $n^{\Omega(k)}$.

Our second main result concerns the agnostic learning of a single $n$-dimensional Gaussian. We prove two SQ lower bounds with qualitatively similar guarantees for different versions of this problem. Our first lower bound is for the problem of agnostically learning a Gaussian with unknown mean and identity covariance. Roughly speaking, we show that any SQ algorithm that solves this learning problem to accuracy $O(\epsilon)$ requires complexity $n^{\Omega(\log^{1/4}(1/\epsilon))}$. We show (see Theorem 5.1 for a more detailed statement):

Roughly speaking, Theorem 1.2 shows that any SQ algorithm that solves the (unknown mean Gaussian) robust learning problem to accuracy $O(\epsilon)$ needs to have running time at least $n^{\Omega(\log^{1/4}(1/\epsilon))}$, i.e., quasi-polynomial in $1/\epsilon$. It is natural to ask whether this quasi-polynomial lower bound can be improved to, say, exponential, e.g., $n^{\Omega(1/\epsilon)}.$ We show that the lower bound of Theorem 1.2 is qualitatively tight. We design an (SQ) algorithm that uses $O_{\epsilon}(n^{\sqrt{\log(1/\epsilon)}})$ SQ queries of inverse quasi-polynomial precision. Moreover, we can turn this SQ algorithm into an algorithm in the sampling oracle model with similar complexity. Specifically, we show (see Theorem 8.7 and Corollary 8.8):

Our second super-polynomial SQ lower bound is for the problem of robustly learning a zero-mean unknown covariance Gaussian with respect to the spectral norm. Specifically, we show (see Theorem 5.12 for a detailed statement):

Our next SQ lower bounds establish nearly quadratic statistical–computational tradeoffs for robust spectral covariance estimation and robust sparse mean estimation. We note that both these lower bounds also hold in Huber’s contamination model. For the former problem, we show (see Theorem 6.1 for the formal statement):

We note that, in order to simulate a single query of the above precision, we need to draw $\Omega(1/\gamma^{2})=\Omega(n^{2-5c})$ samples from our distribution. Roughly speaking, Theorem 1.5 shows that if an SQ algorithm uses less than this many samples, then it needs to run in $2^{\Omega(n^{c/3})}$ time. This suggests a nearly-quadratic statistical-computational tradeoff for this problem.

For robust sparse mean estimation we show (see Theorem 6.6 for the detailed statement):

Similarly, to simulate a single query of the above precision, we need to draw $\Omega(1/\gamma^{2})=\Omega(k^{2-3c})$ samples from our distribution. Hence, any SQ algorithm that uses this many samples requires runtime at least $n^{\Omega(ck^{c})}$. This suggests a nearly-quadratic statistical-computational tradeoff for this problem.

We now turn to our information-theoretic lower bounds on the sample complexity of the corresponding high-dimensional testing problems. For the robust Gaussian mean testing problem in Huber’s contamination model, we show (see Theorem 7.5 for a more detailed statement):

As stated in the Introduction, without the robustness requirement, for any constant $\epsilon>0$, the Gaussian mean testing problem can be solved with $O_{\epsilon}(\sqrt{n})$ samples. Hence, the conceptual message of Theorem 1.7 is that robustness makes the Gaussian mean testing problem information-theoretically harder. In particular, the sample complexity of robust testing is essentially the same as that of the corresponding learning problem. Theorem 1.7 can be viewed as a surprising fact because it implies that the effect of robustness can be very different for testing versus learning of the same distribution family. Indeed, recall that the sample complexity of robustly learning an $\epsilon$-corrupted unknown mean Gaussian, up to error $O(\epsilon)$, is $O(n/\epsilon^{2})$ – i.e., the same as in the noiseless case.

As a final application of our techniques, we show a sample complexity lower bound for the problem of testing whether a spherical GMM is close to a Gaussian (see Theorem 7.6 for the detailed statement):

Similarly, the sample lower bound of Theorem 1.8 is optimal, up to constant factors, and coincides with the sample complexity of learning the underlying distribution.

3 Our Approach and Techniques

In this section, we provide a detailed outline of our approach and techniques. The structure of this section is as follows: We start by describing our Generic Lower Bound Construction, followed by our main applications to the problems of Learning GMMs and Robustly Learning an Unknown Gaussian. We continue with our applications to statistical–computational tradeoffs. We then explain how our generic technique can be used to obtain our Sample Complexity Testing Lower Bounds, which rely on essentially the same hard instances as our SQ lower bounds. We conclude with a sketch of our new (SQ) Algorithm for Robustly Learning an Unknown Mean Gaussian to optimal accuracy.

The main idea of our lower bound construction is quite simple: We construct a family of distributions $\mathcal{D}$ that are standard Gaussians in all but one direction, but are somewhat different in the remaining direction (Definition 3.1). Effectively, we are hiding the interesting information about our distributions in this unknown choice of direction. By exploiting the simple fact that it is possible to find exponentially many nearly-orthogonal directions (Lemma 3.7), we are able to show that any SQ algorithm with insufficient precision needs many queries in order to learn an unknown distribution from $\mathcal{D}$.

To prove our generic SQ lower bound, we need to bound from below the SQ-dimension of our hard family of distributions $\mathcal{D}$. Roughly speaking, the SQ-dimension of a distribution family (Definition 2.11) corresponds to the number of nearly uncorrelated distributions (with respect to some fixed distribution) in the family (see Definitions 2.9 and 2.10). It is known that a lower bound on the SQ-dimension implies a corresponding lower bound on the number and precision of queries of any SQ algorithm (see Lemma 2.12).

For the sake of the intuition, we make two observations: (1) If $A$ and $N(0,1)$ have substantially different moments of degree at most $m$, for some $m$, then $\mathbf{P}_{v}$ and $N(0,I)$ can be easily distinguished by comparing their $m^{th}$-order moment tensors. Since these tensors can be approximated in roughly $n^{m}$ queries (and time), the aforementioned lower bound construction would necessarily fail unless the low-order moments of $A$ match the corresponding low-order moments of $G$. We show that, aside from a few mild technical conditions (see Condition 3.2), this moment-matching condition is essentially sufficient for our purposes. If the degree at most $m$ tensors agree, we need to approximate tensors of degree $m+1$. Intuitively, in order to extract useful information from these higher degree tensors, one needs to approximate essentially all of the $n^{m+1}$ many such tensor entries. (2) A natural approach to distinguish between $\mathbf{P}_{v}$ and $N(0,I)$ would be via random projections. As a critical component of our proof, we show (see Lemma 3.5) that a random projection of $\mathbf{P}_{v}$ will be exponentially close to $N(0,1)$ with high probability. Therefore, a random projection-based algorithm would require exponentially many random directions until it found a good one.

We now proceed with a somewhat more technical description of our proof. To bound from below the SQ-dimension of our hard family of distributions, we proceed as follows: The definition of the pairwise correlation (Definition 2.9) implies we need to show that $\int\mathbf{P}_{v}\mathbf{P}_{v^{\prime}}/G\approx 1$, where $G\sim N(0,I)$ is the Gaussian measure, for any pair of unit vectors $v,v^{\prime}$ that are nearly orthogonal. To prove this fact, we make essential use of the Gaussian (Ornstein–Uhlenbeck) noise operator and its properties (see, e.g., [O’D14]). We explain this connection in the following paragraph.

By construction of the distributions $\mathbf{P}_{v},\mathbf{P}_{v^{\prime}}$, it follows that in the directions perpendicular to both $v$ and $v^{\prime}$, the relevant factors integrate to $1$. Letting $y=v\cdot\mathbf{x}$ and $z=v^{\prime}\cdot\mathbf{x}$ and letting $y^{\prime},z^{\prime}$ be the orthogonal directions to $y$ and $z$, we need to consider the integral

Fixing $y$ and integrating over the orthogonal direction, we get

Now, if $v$ and $v^{\prime}$ are (exactly) orthogonal, $z=y^{\prime}$ and the inner integral equals $G(y)$. When this is not the case, the $A(z)$ term is not quite vertical and the $G(z^{\prime})$ term not quite horizontal, so instead what we get is only nearly Gaussian. In general, the inner integral is equal to

where $U_{t}$ is the member of the Ornstein–Uhlenbeck semigroup, $U_{t}f(z)=\mathbf{E}[f(tz+\sqrt{1-t^{2}}G)].$ We show that this quantity is close to a Gaussian, when $v\cdot v^{\prime}$ is close to (see Lemma 3.4).

The core idea of the analysis relies on the fact that $U_{t}A$ is a smeared out version of $A$. As such, it only retains the most prominent features of $A$, namely its low-order moments. In fact, we are able to show that if $A$ and $G$ agree in their first $m$ moments, then $U_{t}A$ is $O_{m}(t^{m})$-close to a Gaussian (see Lemma 3.5), and thus the integral in question is $O_{m}((|v\cdot v^{\prime}|)^{m})$-close to $1$. This intuition is borne out in a particularly clean way by writing $A/G$ in the basis of Hermite polynomials. The moment-matching condition implies that the decomposition involves none of the Hermite polynomials of degrees $1$ through $m$. However, the Ornstein–Uhlenbeck operator, $U_{t}$, is diagonalized by the basis $H_{i}G$ with eigenvalue $t^{i}$. Thus, if $A-G$ can be written in this basis with no terms of degree less than $m$, applying $U_{t}$ decreases the size of the function by a multiple of approximately $t^{m}$.

So far, we have provided a proof sketch of the following statement (Lemma 3.4): When two unit vectors $v,v^{\prime}$ are nearly orthogonal, then the distributions $\mathbf{P}_{v},\mathbf{P}_{v^{\prime}}$ are nearly uncorrelated. Since, for $0<c<1/2$, we can pack $2^{\Omega(n^{c})}$ unit vectors $v$ onto the sphere so that their pairwise inner products are at most $n^{c-1/2}$ (Lemma 3.7), we obtain an SQ-dimension lower bound of our hard family. In particular, to learn the distribution $\mathbf{P}_{v}$, for unknown $v$, any SQ algorithm requires either $2^{\Omega(n^{c})}$ queries or queries of accuracy better than $O(n)^{(m+1)(c-1/2)}$ (Proposition 3.3). This completes the proof sketch of our generic construction.

The properties of our one-dimensional distribution $A$ are summarized in Proposition 4.2. Specifically, we construct a distribution $A$ on the real line that is a $k$-mixture of one-dimensional “skinny” Gaussians, $A_{i}$, that agrees with $N(0,1)$ on the first $m=2k-1$ moments (condition (i)). For technical reasons, we require that the chi-squared divergence of $A$ to $N(0,1)$ is bounded from above by an appropriate quantity (condition (iv)). The Gaussian components, $A_{i}$, have the same variance and appropriately bounded means (condition (ii)). We can also guarantee that the components $A_{i}$ are almost non-overlapping (condition (iii)). This implies that the corresponding high-dimensional distributions $\mathbf{P}_{v},\mathbf{P}_{v}^{\prime}$ will be at total variation distance close to $1$ from each other when the directions $v,v^{\prime}$ are nearly orthogonal, and moreover their means will be sufficiently separated.

To establish the existence of a distribution $A$ with the above properties, we proceed in two steps: First, we construct (Lemma 4.3) a discrete one-dimensional distribution $B$ supported on $k$ points, lying in an $O(\sqrt{k})$ length interval, that agrees with $N(0,1)$ on the first $k$ moments. The existence of such a distribution $B$ essentially follows from standard tools on Gauss-Hermite quadrature. The distribution $A$ is then obtained (Corollary 4.4) by adding a zero-mean skinny Gaussian to an appropriately rescaled version of $B$. Additional technical work (Lemmas 4.5 and 4.6) gives the other conditions.

Our family of hard high-dimensional instances will consist of GMMs that look like almost non-overlapping “parallel pancakes” and is reminiscent of the family of instances considered in Brubaker and Vempala [BV08]. For the case of $k=2$, consider a $2$-GMM where both components have the same covariance that is far from spherical, the vector between the means is parallel to the unit eigenvector with smallest eigenvalue, and the distance between the means is a large multiple of the standard deviation in this direction (but a small multiple of that in the orthogonal direction). This family of instances was considered in [BV08], who gave an efficient spectral algorithm to learn them.

Our lower bound construction can be thought of as $k$ “parallel pancakes” in which the means lie in a one-dimensional subspace, corresponding to the smallest eigenvalue of the identical covariance matrices of the components. All $n-1$ orthogonal directions will have an eigenvalue of $1$, which is much larger than the smallest eigenvalue. In other words, for each unit vector $v$, the $k$-GMM $\mathbf{P}_{v}$ will consist of $k$ “skinny” Gaussians whose mean vectors all lie in the direction of $v$. Moreover, each pair of components will have total variation distance very close to $1$ and their mean vectors are separated by $\Omega(1/\sqrt{k})$. We emphasize once more that our hard family of instances is learnable with $O(k\log n)$ samples – both for density estimation and parameter estimation. On the other hand, any SQ learning algorithm for the family requires $n^{\Omega(k)}$ time.

In the agnostic model, there are two types of adversarial noise to handle: subtractive noise – corresponding to the good samples removed by the adversary – and additive noise – corresponding to the bad points added by the adversary. The approach of [DKK+16] does not do anything to address subtractive noise, but shows that this type of noise can incur “small” error, e.g., at most $O(\epsilon\sqrt{\log(1/\epsilon)})$ for the case of unknown mean. For additive noise, [DKK+16] uses an iterative spectral algorithm to filter out outliers.

For concreteness, let us consider the case of robustly learning $N(\mu,I)$. Intuitively, achieving error $O(\epsilon)$ in the agnostic model is hard for the following reason: the two types of noise can collude so that the first few moments of the corrupted distribution are indistinguishable from those of a Gaussian whose mean vector has distance $\Omega(\epsilon\sqrt{\log(1/\epsilon)})$ from the true mean.

To formalize this intuition, for our robust SQ learning lower bound, we construct a distribution $A$ on the real line that agrees with $N(0,1)$ on the first $m=\Omega(\log^{1/4}(1/\epsilon))$ moments and is $\epsilon/100$-close in total variation distance to $G^{\prime}=N(\epsilon,1)$ (see Proposition 5.2). We achieve this by taking $A$ to be the Gaussian $N(\epsilon,1)$ outside its effective support, while in the effective support we add an appropriate degree-$m$ univariate polynomial $p$ satisfying the appropriate moment conditions. By expressing this polynomial as a linear combination of appropriately scaled Legendre polynomials, we can prove that its $L_{1}$ and $L_{\infty}$ norms within the effective support of $G^{\prime}$ are much smaller than $\epsilon$ (see Lemma 5.6). This result is then used to bound from above the distance of $A$ from $G^{\prime}$, which gives our SQ lower bound.

We use a similar technique to prove our SQ lower bound for robust covariance estimation in spectral norm. Specifically, we construct a distribution $A$ that agrees with $N(0,1)$ on the first $m=\Omega(\log(1/\epsilon))$ moments and is $\epsilon/100$-close in total variation distance to $G^{\prime}=N(0,(1-\delta)^{2})$, for some $\delta=O(\epsilon)$ (see Proposition 5.13). We similarly take $A$ to be the Gaussian $G^{\prime}$ outside its effective support, while in the effective support we add an appropriate degree-$m$ univariate polynomial $p$ satisfying the appropriate moment conditions. The analysis proceeds similarly as above.

For robust covariance estimation in spectral norm, our one-dimensional distribution is selected to be $A=(1-\epsilon)N(0,\sigma)+\epsilon N_{1}$, where $N_{1}$ is a mixture of $2$ unit-variance Gaussians with opposite means. By selecting $\sigma$ appropriately, we can have $A$ match the first $3$ moments of $N(0,1)$, see Theorem 6.1. For robust sparse mean estimation, it suffices to take $A=(1-\delta)N(\epsilon,1)+\delta N_{1}$, where $N_{1}$ is a unit-variance Gaussian selected so that $\mathbf{E}[A]=0$. An important aspect of both these constructions is that the chi-squared distance $\chi^{2}(A,N(0,1))$ needs to be as small as possible. Indeed, since we only match a small number of moments, our bound on $\chi^{2}(A,N(0,1))$ crucially affects the accuracy of our SQ queries (Proposition 3.3).

Note that the inner integral was bounded from above by roughly $(1+(v\cdot v^{\prime})^{m})$. A careful analysis of the distribution of the angle between two random unit vectors allows us to show that, unless $N=\Omega(n)$, the chi-squared divergence is close to $1$, and thus that this testing problem is impossible.

We give an SQ algorithm with $O(\epsilon)$-error for robustly learning an unknown mean Gaussian, showing that our corresponding SQ lower bound is qualitatively tight. Our algorithm builds on the filter technique of [DKK+16], generalizing it to the more involved setting of higher-order tensors.

As is suggested by our SQ lower bounds, the obstacle to learning the mean robustly, is that there are $\epsilon$-noisy Gaussians that are $\Omega(\epsilon)$-far in variation distance from a target Gaussian $G$, and yet match $G$ in all of their first $O(\log^{1/4}(1/\epsilon))$ moments. For our algorithm to circumvent this difficulty, it will need to approximate all of the $t^{th}$-order tensors for $t\leq k=\Omega(\log^{1/4}(1/\epsilon))$. Note that this already requires $n^{k}$ SQ queries.

The first thing we will need to show is that $k$ moments suffice, for an appropriate parameter $k$. Because of our lower bound construction, we know that $k$ needs to be at least $\Omega(\log^{1/4}(1/\epsilon))$. We show that $k=O(\log^{1/2}(1/\epsilon))$ suffices. Specifically, we prove a one-dimensional moment-matching lemma (Lemma 8.1) establishing the following: If an $\epsilon$-noisy one-dimensional Gaussian approximately matches a reference Gaussian $G$ in all of its first $k$ moments, where $k=\Theta(\log^{1/2}(1/\epsilon))$ (i.e., quadratically larger than our lower bound), then it must be $O(\epsilon)$-close to $G$ in variation distance. We note that it suffices to prove this statement in the one-dimensional case, as we can just project onto the line between the means.

We now proceed to describe our algorithm: Using the basic filter algorithm from [DKK+16], we start by learning the true mean to error $O(\epsilon\sqrt{\log(1/\epsilon)})$. By translating, we can assume that the mean is this close to . We need to robustly approximate the low-order moments of our target Gaussian $G^{\prime}$. This is complicated by the fact that even a small fraction of errors can have a huge impact on the moments of the distribution. However, any large errors are easily detectable. In particular, if any $t^{th}$ moment tensor differs substantially from that of the standard Gaussian, it will necessarily imply the presence of errors. In particular, it will allow us to construct a polynomial $p$ so that $\mathbf{E}[p(X)]-\mathbf{E}[p(G^{\prime})]$ (where $X$ is a noisy version of $G^{\prime}$) is much larger than $\epsilon\|p(G^{\prime})\|_{2}$. If this is the case, then many of our errors, $x$, must have $p(x)$ very far from the mean. By standard concentration inequalities, this will allow us to identify these points as almost certainly being errors. This in turn lets us build a filter to clean-up our distribution $X$, making it closer to $G^{\prime}$.

Repeatedly applying filters as necessary, we can reduce to the case where the higher-order moments of $X$ are close to the higher-order moments of $G$. This will tell us that, in almost all directions, the first $k$ moments of $X$ match the corresponding moments of $G$. By our moment-matching lemma, this will imply that the mean of $G^{\prime}$ is close to in these directions. We will then only need to approximate the mean of the projection of $G^{\prime}$ onto the low-dimensional subspace $V$ in which these moments fail to match. This approximation can be done in a brute-force manner (in time exponential in $\dim(V)$, which is still relatively small), completing the description of the algorithm.

4 Related Work

This work studies learning and testing high-dimensional structured distributions. Distribution learning and testing are two of the most fundamental inference tasks in statistics with a rich history (see, e.g., [NP33, BBBB72, DG85, Sil86, Sco92, DL01, LR05]) that date back to Karl Pearson. The main criteria to evaluate the performance of an estimator are its sample complexity and its computational complexity. Despite intensive investigation for several decades by different communities, the (sample and/or computational) complexity of many learning and testing problems is still not well-understood, even for some surprisingly simple high-dimensional settings. In the past few decades, a long line of work within TCS [KMR+94, Das99, FM99, AK01, VW02, CGG02, MR05, BV08, KMV10, MV10, BS10, DDS12a, DDS12b, CDSS13, DDO+13, CDSS14a, CDSS14b, ADLS17, DDS15, DDKT16, DKS16b, DKS16a] has focused on designing efficient estimators in a variety of settings. We have already mentioned the most relevant references for the specific questions we consider in Section 1.1.

With respect to computational lower bounds for unsupervised estimation problems, the most relevant references are the works [FGR+13, FPV15, KV16] that show SQ lower bounds for the planted clique and related planted-like problems. It should be noted that, beyond the fact that we also use the concept of SQ dimension, our techniques are entirely different than theirs. Prior work by Feldman, O’Donnell, and Servedio [FOS08] implicitly showed an SQ lower bound of $n^{\Omega(\log k)}$ for the problem of learning $k$-mixtures of product distributions over $\{0,1\}^{n}$. This was obtained by a straightforward reduction from the problem of learning $k$-leaf decision trees over $n$ Boolean variables. Our lower bound construction for learning GMMs is entirely different from [FOS08] that relied on the obvious combinatorial structure of the discrete setting.

A related line of work gives statistical-computational tradeoffs for sparse PCA [BR13a, BR13b, MW15, WBS16a], based on various computational hardness assumptions. These results are of similar flavor as our statistical–computational tradeoffs for SQ algorithms (Theorems 1.5 and 1.6). An important difference between these tradeoffs and the super-polynomial SQ lower bounds we prove in this paper (Theorems 1.1, 1.2, and 1.4) is that the aforementioned sparse problems are known to be tractable if we increase the sample size by a quadratic factor beyond the information-theoretic limit. In contrast, our main SQ lower bound results establish a super-polynomial gap between the information-theoretic limit and the computational complexity of any SQ algorithm.

Finally, we remark that in the supervised setting of PAC learning Boolean functions, a number of hardness results are known based on various complexity assumptions, see, e.g., [KKMS08, KS06, FGKP06, KK14, DLS14, Dan16] for the problems of learning halfspaces and learning intersections thereof.

5 Discussion and Future Directions

The main contribution of this paper is a technique that gives essentially tight SQ lower bounds for a number of fundamental high-dimensional learning problems, including learning GMMs and robustly learning a single Gaussian. To the best of our knowledge, these are the first such lower bounds for high-dimensional distribution learning problems in the continuous setting. As a corollary, we provide a rigorous explanation of the observed (super-polynomial) gap between the sample complexity of these problems and the runtime of the best known algorithms.

6 Organization

The structure of this paper is as follows: In Section 2, we introduce basic notation, definitions, and a number of useful facts that will be required throughout the paper. Our SQ lower bounds are established in Sections 3–6. Specifically, in Section 3, we give our generic high-dimensional SQ lower bound construction, assuming the existence of a one-dimensional density satisfying the necessary moment conditions. In Sections 4, 5, and 6, we construct the appropriate one-dimensional densities, thereby establishing our SQ lower bounds. Specifically, Sections 4 and 5 give our super-polynomial SQ lower bounds for the problems of learning GMMs and robustly learning an unknown Gaussian. Section 6 gives our quadratic statistical–computational tradeoffs (for SQ algorithms) for the problems of robust covariance estimation in spectral norm and robust sparse mean estimation. Section 7 gives our (information-theoretic) sample complexity lower bounds for high-dimensional testing. Finally, in Section 8 we present our (SQ) algorithm for robustly learning an unknown mean Gaussian with optimal accuracy, whose runtime qualitatively matches our SQ lower bound from Section 5.

This project evolved over a number of years. We would like to thank Vitaly Feldman for answering numerous questions about the Statistical Query model; Andy Drucker for useful discussions on Question 1.1; Anup B. Rao for asking a question that motivated Theorem 6.1; Weihao Kong and Gregory Valiant for useful discussions on Question 1.4; and Ankur Moitra and Eric Price for feedback on a previous version of this paper.

Definitions and Preliminaries

Our basic object of study is the Gaussian (or Normal) distribution and finite mixtures of Gaussians:

Throughout the paper, we will make extensive use of the pdf of the standard one-dimensional Gaussian $N(0,1)$, which we will denote by $G(x)$.

2 Formal Problem Definitions

We record here the formal definitions of the problems that we study. Our first problem of interest is learning a mixture of $k$ arbitrary high-dimensional Gaussians:

Our next question is the problem of robustly learning a Gaussian in the standard agnostic model:

Our SQ lower bounds apply to two special cases of this problem: when $\mu$ is unknown and $\Sigma=I$, and when $\mu=0$ and $\Sigma$ is unknown. For the latter case, our lower bound applies even for learning with respect to the spectral norm (which is weaker than approximation in variation distance).

We now define the problem of robustly testing a Gaussian in Huber’s model. We remind the reader that our tight sample complexity lower bound applies to this weaker model as well.

Finally, our problem of testing GMMs is the following:

3 Basics on Statistical Query Algorithms over Distributions

We begin by recording the necessary definitions of Statistical algorithms for problems over distributions. All the definitions and facts in this section are from [FGR+13]. We start by defining a general search problem over distributions.

For general search problems over a distribution, we define SQ algorithms as algorithms that do not see samples from the distribution but instead have access to an SQ oracle. We consider two types of SQ oracles from the literature.

The first oracle was defined by Kearns [Kea98] and the second was introduced in [FGR+13]. These oracles are known to be polynomially equivalent [FGR+13]. Also note that these oracles can return any value within the given tolerance, and therefore can make adversarial choices.

The main technical tool that allows us to prove unconditional lower bounds on the complexity of SQ algorithms is an appropriate notion of Statistical Query (SQ) dimension. Such a notion was defined in the context of PAC learning of Boolean functions in [BFJ+94], and subsequently generalized to search problems over distributions in [FGR+13]. We will require the simpler definition from Section 3 of that work that relies on pairwise correlations:

We will also need the following definition:

We are now ready to define our notion of dimension:

Our lower bounds proceed by bounding from below the statistical query dimension of the considered distribution learning problems. The corresponding lower bounds on the complexity of SQ algorithms for these problems are a corollary of the following result from [FGR+13]:

Statistical Query Lower Bounds: From One–Dimension to High–Dimensions

All our statistical query lower bounds are shown in two steps: We first construct a one-dimensional density $A$ satisfying certain technical conditions, and then use $A$ to construct a high-dimensional distribution which is Gaussian in all but one directions. The second step is the same for all the problems that we consider. To formally define it, we require the following construction:

That is, $\mathbf{P}_{v}$ is the product distribution whose orthogonal projection onto the direction of $v$ is $A$, and onto the subspace perpendicular to $v$ is the standard $(n-1)$-dimensional normal distribution.

Suppose that we have constructed a one-dimensional distribution $A$ satisfying the following condition:

Note that Condition 3.2-(ii) above implies that the distribution $A$ has a probability density function (pdf), which we will denote by $A(x)$. We will henceforth blur the distinction between a distribution and its pdf. The main result of this section is the following:

An intuitive interpretation of the proposition is as follows: If we do not want our SQ algorithm to use a number of queries exponential in $n^{\Omega(1)}$, then we would need $\Omega(n)^{\Omega(m+1)}$ samples to simulate a single statistical query.

The rest of this section is devoted to the proof of Proposition 3.3. In Section 3.1, we prove a correlation bound which is the main technical ingredient for the proof. In Section 3.2, we show a simple packing for unit vectors over the sphere and put the pieces together to complete the proof.

The main technical result of this section is the following:

Note that we may assume that $\chi^{2}(A,N(0,1))$ is finite, otherwise the lemma statement is trivial. Hence, we can henceforth assume that Condition 3.2 is satisfied. In particular the distributions $A$, $\mathbf{P}_{v}$ and $\mathbf{P}_{v^{\prime}}$ all have probability density functions.

To prove Lemma 3.4 we proceed as follows: We start by bounding the $\chi^{2}$-divergence between the one-dimensional projection of $\mathbf{P}_{v}$ onto $v^{\prime}$ and $N(0,1)$. As well as being a critical component towards the proof of Lemma 3.4, this fact can be used to show that random projections of $\mathbf{P}_{v}$ are close to $N(0,1)$ with high probability. Specifically, we show:

Let $\mathbf{Q}$ be the distribution of $v^{\prime}\cdot X$, for $X\sim\mathbf{P}_{v}$. Then, we have that

Let $\theta$ be the angle between $v$ and $v^{\prime}$. Let $x,y$ be orthogonal coordinates for the plane spanned by $v$ and $v^{\prime}$, with the $x$-axis in the $v^{\prime}$ direction. Note that $\mathbf{P}_{v}$ is a product of a distribution on this plane and a standard Gaussian perpendicular to it. On this plane, $\mathbf{P}_{v}$ is a product of $A$ and $N(0,1)$. Thus, we have that

so that $\mathbf{Q}=U_{\theta}(A)$. We will show that we can expand $A$ as a linear combination of eigenfunctions of $U_{\theta}$.

Using the orthogonality of $He_{i}(x)$ we can extract these coefficients, since

Since $A$ agrees with the first $m$ moments of the standard Gaussian, for $0\leq i\leq m$, we have that

This implies that $a_{0}=1$ and $a_{1},\dots,a_{m}=0$. Thus, we have

Now we consider the effect of $U_{\theta}$ on this orthogonal family. From the definition of $U_{\theta}$, we have

We will use the well-known fact that $He_{i}(x)G(x)$ is an eigenfunction of $U_{\theta}$:

We have that: $U_{\theta}(He_{i}G)(x)=\cos^{i}(\theta)He_{i}(x)G(x)\;.$

For completeness, we include a proof in Appendix D.

We can now use these eigenfunctions and eigenvalues with Equation (5) to get an expression for $U_{\theta}A$:

which can be used to express its $\chi^{2}$-divergence:

where the last line uses (6). Recalling that $\mathbf{Q}=U_{\theta}A$ and $\cos\theta=v\cdot v^{\prime}$, this completes the proof. ∎

We first show that the correlation between the high-dimensional densities $\mathbf{P}_{v}$ and $\mathbf{P}_{v^{\prime}}$ needed for Lemma 3.4 can be reduced to a one-dimensional correlation. Just as in the proof of Lemma 3.5, let $\theta=\arccos(v\cdot v^{\prime})$ and let $x,y$ be coordinates for the plane spanned by $v$ and $v^{\prime}$ with the $x$-axis in the $v^{\prime}$ direction. Each of $\mathbf{P}_{v}$ and $\mathbf{P}_{v^{\prime}}$ is a product of a distribution on this plane and a standard Gaussian perpendicular to it. On this plane, they are both products of $A$ and $N(0,1)$ with different rotations applied. Thus, we have that

Now we can bound from above this correlation in terms of the $\chi^{2}$-divergences of both distributions from $N(0,1)$, one of which we can bound using Lemma 3.5:

where the first line follows by triangle inequality, the second inequality is Cauchy-Schwarz, and the last line uses Lemma 3.5. The proof of Lemma 3.4 is now complete. ∎

2 Proof of Proposition 3.3

If $v=v^{\prime}$, then $\chi_{N(0,I)}(\mathbf{P}_{v},\mathbf{P}_{v})=\chi^{2}(\mathbf{P}_{v},N(0,I))=\chi^{2}(A,N(0,1))$. We thus have that, for

$\mathcal{D}_{D}$ is $(\gamma,\beta)$ correlated with respect to $D=N(0,I)$.

oracle to solve $\mathcal{Z}.$ From our assumption that $n\geq m^{\Omega(1/c)}$, it follows that $n\geq\Omega((m+1)\log n)^{2/c}$, and therefore $2^{\Omega(n^{c/2})}\geq n^{m+1}$. Hence, the total number of required queries is at least $2^{\Omega(n^{c/2})}\geq n^{m+1}$. This completes the proof. ∎

SQ Lower Bound for Learning Gaussian Mixtures

The main result of this section is the following:

Remark. We remark that the well-conditioned assumption in Theorem 4.1 (i.e., that the distances between the means and the largest and smallest eigenvalues of any covariance matrix are bounded) guarantees that an SQ algorithm with a bounded number of SQ queries is possible.

The proof of Theorem 4.1 follows by an application of the framework developed in Section 3 and the following proposition:

$A$ agrees with $N(0,1)$ on the first $2k-1$ moments.

Each Gaussian component $A_{i}$ has variance $\Theta\left(\frac{1}{k^{2}\log^{2}(k+1/\epsilon)}\right)$ and mean of magnitude $O(\sqrt{k})$.

We have $\chi^{2}(A,N(0,1))\leq\exp(O(k))\log(1/\epsilon)$.

Given Proposition 4.2, the proof of Theorem 4.1 follows easily.

It remains to show that $\mathbf{P}_{v}$ is a mixture of $k$ Gaussians that satisfies the necessary conditions. Note that $\mathbf{P}_{v}$, when expressed in an appropriate basis, is a product of the mixture of $k$ univariate Gaussians and the standard $(n-1)$-dimensional normal distribution. Recall that the product of two Gaussians is a Gaussian. If $A=\sum_{i=1}^{k}w_{i}N(\mu^{\prime}_{i},\delta)$, where $\delta=\Theta\left(\frac{1}{k^{2}\log(k+1/\epsilon)}\right)$ (by Proposition 4.2 (ii)), then we have that $\mathbf{P}_{v}=\sum_{i=1}^{k}w_{i}N\left(v\mu^{\prime}_{i},I-(1-\delta)vv^{T}\right)$. We can bound the variation distance between two components by:

by Proposition 4.2 (iii). Also, we have that

There is a discrete distribution $B$ on the real line, supported on $k$ points, that agrees with $N(0,1)$ on the first $2k-1$ moments. All points $x$ in the support of $B$ have $|x|=O(\sqrt{k})$.

This lemma essentially follows from standard techniques for Gaussian quadrature [AS72]. Given a (possibly infinite) interval $[a,b]$, a weighting function $\omega(x)$, and an integer $k>0$, we can find $x_{i}$ and $w_{i}$ for $1\leq i\leq k$ such that

for all polynomials $p(x)$ of degree at most $2k-1$. The Gauss-Hermite quadrature is a standard implementation of this general scheme on the interval $(-\infty,\infty)$ with $\omega(x)=e^{-x^{2}}$. Here, we take the $x_{i}$’s to be the roots of the $k$-th (physicist’s) Hermite polynomial $H_{k}(x)$. Then, we have that $w_{i}=\frac{2^{k-1}k!\sqrt{\pi}}{k^{2}H_{k-1}(x_{i})^{2}}$.

for all polynomials $p(x)$ of degree at most $2k-1$.

Note that all the weights are nonnegative by definition. Also note that $\sum_{i=1}^{k}w^{\prime}_{i}=\int_{-\infty}^{\infty}1\cdot G(x)=1$. We take $B$ to be the probability distribution with probability $w^{\prime}_{i}$ of being $x^{\prime}_{i}$, for each $1\leq i\leq k$. Then we have

It is known (see, e.g., [Sze89]) that all roots of $H_{k}(x)$ have absolute value $O(\sqrt{k})$, and so all roots of $He_{k}$. Hence, all points $x$ in the support of $B$ have $|x|=O(\sqrt{k})$. This completes the proof. ∎

On the other hand, if we want $\chi^{2}(A,N(0,1))$ to be finite, we need to have a mixture of Gaussians each with positive variance $\delta>0$.

By rescaling the distribution $B$ given by Lemma 4.3, we can find a discrete distribution $B^{\prime}$ supported on $k$ points with absolute value no bigger than $O(\sqrt{k})$ that agrees with the first $2k-1$ moments of $N(0,1-\delta)$. The rescaled distribution $B^{\prime}$ assigns probability mass $w^{\prime}_{i}$ to the points $\sqrt{1-\delta}x^{\prime}_{i}$, for $1\leq i\leq k$. Let $X\sim B^{\prime}$, $X^{\prime}\sim N(0,1-\delta)$, and $Y\sim N(0,\delta)$ that is independent of $X,X^{\prime}$. We take $A$ to be the distribution of $X+Y$. Then, we have

for all integers $1\leq j\leq{2k-1}$. By standard facts about Gaussians, $X^{\prime}+Y$ is distributed as $N(0,1)$. Finally, note that the distribution of $X+Y$ is a mixture of $k$ Gaussians $N(\sqrt{1-\delta}x^{\prime}_{i},\delta)$ with weights $w^{\prime}_{i}$. ∎

To appropriately set the parameter $\delta$, we need to consider the high-dimensional construction (Definition 3.1):

We write $A_{i}$, for $1\leq i\leq k$, for the Gaussians $N(\mu_{i},\delta)$ that $A$ is a mixture of. Fix $\epsilon>0$. By a Chernoff bound, $A_{i}$ is within the interval $[\mu_{i}-a,\mu_{i}+a]$, where $a=2\sqrt{\delta\log(1/\epsilon)}$ with probability at least $1-\epsilon$.

We again consider the plane spanned by $v$ and $v^{\prime}$. Let $x,y$ be the orthogonal coordinates with $v$ in the direction of the $x$-axis. Similarly, let $x^{\prime},y^{\prime}$ be the orthogonal coordinates with $v^{\prime}$ in the direction of the $x^{\prime}$-axis. Let $\theta$ be the angle between $v$ and $v^{\prime}$. We have that

Taking $\epsilon=\sqrt{\delta}$, we obtain that $\int_{\mathbf{x}}\min\{\mathbf{P}_{v}(\mathbf{x}),\mathbf{P}_{v^{\prime}}(\mathbf{x})\}d\mathbf{x}\leq O(k\sqrt{\delta}\log(1/\delta)\csc\theta)$. On the other hand,

This gives an upper bound on $\delta$. We don’t want $\delta$ to be too small, because of the following lemma:

We have that $\chi^{2}(A,N(0,1))\leq\exp(O(k))/\sqrt{\delta}$.

Each component $A_{i}$, for $1\leq i\leq k$, satisfies the following:

The following simple lemma helps us enforce the condition that the Gaussian components are well-separated:

We now have all the necessary tools to prove Proposition 4.2. We take $\delta=C/(k^{2}\log^{2}(k+1/\epsilon))$ for a sufficiently small constant $C>0$. Combined with Corollary 4.4, this gives condition (ii). For condition (i), note that, by Corollary 4.4, $A$ agrees with $N(0,1)$ on the first $2k-1$ moments. Since $\delta$ was selected to be smaller than $O(1/{k}\log(1/\epsilon))$, Lemma 4.7 gives condition (iii). Lemma 4.6 gives condition (iv). Finally, by our choice of $\delta$ and Lemma 4.5, we get condition (v). This completes the proof. ∎

SQ Lower Bounds for Robust Learning of a Gaussian

In this section, we prove our super-polynomial SQ lower bounds for robustly learning a high-dimensional Gaussian. In Section 5.1, we show our lower bound for robustly learning an unknown mean spherical Gaussian. In Section 5.2, we give our lower bound for robustly learning a zero mean unknown covariance Gaussian with respect to the spectral norm.

In this subsection, we use the framework of Section 3 to prove the following theorem:

The theorem will follow from the following proposition:

$A$ and $N(0,1)$ agree on the first $m$ moments.

Before we prove Proposition 5.2, we show how Theorem 5.1 easily follows from it using the machinery developed in Section 3.

Therefore, for any unit vectors $v,v^{\prime}$ with $|v\cdot v^{\prime}|\leq 1/8$ we have that:

where we used the assumption that $C$ is sufficiently large. This completes the proof. ∎

The rest of this section is devoted to the proof of Proposition 5.2. We start by describing the outline of the proof. We then provide a number of intermediate useful lemmas that we subsequently combine to complete the proof.

The proof plan proceeds as follows. For some $C=\Theta(\sqrt{\log(1/\delta)})$, we define the one-dimensional distribution $A$ to be:

For $x\notin[-C,C]$, we define $A(x)=G(x-\delta)$.

For $x\in[-C,C]$, we define $A(x)=G(x-\delta)+p(x)$, where $p(x)$ is the degree-$m$ polynomial with $\int_{-C}^{C}p(x)dx=0$ and

for $1\leq i\leq m$. (We note that $p$ is unique after fixing $m$, $C$ and $\delta$.)

We need to show that we can find appropriate values for the parameters $m$, $C$, and $\delta$ such that the $L_{1}$-norm of $p(x)$ is at most $O(\delta m^{2}/\sqrt{\log(1/\delta)})$ and that $A(x)$ is non-negative. To achieve that, we will express $p(x)$ as a linear combination of (appropriately scaled) Legendre polynomials, a family of orthogonal polynomials on $[-C,C]$. Rather than directly showing that the first $m$ moments agree, we will instead want that the expectations of the first $m$ scaled Legendre polynomials agree. Bounds on the coefficients of the Legendre polynomials in $p(x)$ allow us to obtain bounds on the $L_{1}$ and $L_{\infty}$ norms of $p(x)$ on $[-C,C]$. Choosing $m$, $\delta$, and $C$ appropriately will complete the proof of the proposition.

We start by recording the properties of Legendre polynomials that we will need:

$P_{k}(x)$ is a degree-$k$ polynomial, $P_{0}(x)=1$, and $P_{1}(x)=x$.

$\int_{-1}^{1}P_{i}(x)P_{j}(x)dx=(2/(2i+1))\delta_{i,j}$ for all $i,j\geq 0.$

$P_{k}(x)=(1/2^{k})\sum_{i=0}^{\lfloor k/2\rfloor}{k\choose i}{2k-2i\choose k}x^{k-2i}.$

As a simple corollary we obtain the following lemma:

$|P_{k}(x)|\leq(4|x|)^{k}$ for all $|x|\geq 1$.

$\int_{-1}^{1}|P_{k}(x)|dx\leq O(1/\sqrt{k}).$

We are now ready to proceed with the formal proof. The main technical result of this section is the following lemma:

We can write $p(x)=\sum_{k=0}^{m}a_{k}P_{k}(x/C)$, where $|a_{k}|=O(\delta k^{3/2}/C^{2})$, for $0\leq k\leq m$.

Before we give the proof of Lemma 5.5, we deduce two corollaries that will be useful in the proof of Proposition 5.2. First, we can obtain bounds on the $L_{1}$ and $L_{\infty}$ norms of $p(x)$ on $[-C,C]$. As an immediate corollary of Lemma 5.5 and the aforementioned properties of Legendre polynomials, we deduce:

We have that: $\int_{-C}^{C}|p(x)|dx\leq O(\delta m^{2}/C)$ and $|p(x)|\leq\delta m^{5/2}/C^{2}$, for all $x\in[-C,C].$

We now bound from above the desired $\chi^{2}$-divergence:

$\chi^{2}(A,N(0,1))=O\left(\delta^{2}+\delta m^{5/2}/C^{2}\cdot(C^{2}\delta^{{2}}+\max_{|x|\leq C}|p(x)|/G(x))\right)$.

We bound the second term from above as follows:

where the last lines uses Lemma 5.5 and Corollary 5.6. Finally, for the third term we have:

where the inequality follows from Corollary 5.6. This completes the proof of Lemma 5.7. ∎

We first note that we can express $p(x)$ as a linear combination of scaled Legendre polynomials whose coefficients are explicitly given by integrals:

We can write $p(x)=\sum_{k=0}^{m}a_{k}P_{k}(x/C)$, where $a_{k}=((2k+1)/2C)\int_{-C}^{C}P_{k}(x/C)p(x)dx$.

It follows from Fact 5.3 (ii) and a change of variables that $\int_{-C}^{C}P_{i}(x/C)P_{j}(x/C)dx=(2C/(2i+1))\delta_{i,j}$, for all $i,j\geq 0.$ We can use this to extract the $a_{k}$’s. For $1\leq k\leq m$, we have

Since the first $m$ moments of $p$ are fixed, via 10, we obtain:

Since we will apply this with the parameter $1/\delta$ exponential in $m$ and $C$, we will be able to ignore $O(\delta^{2})$ terms. We use Taylor’s theorem to expand $(G(x)-G(x-\delta))$ up to second order terms:

$G(x)-G(x-\delta)=xG(x)\delta+(\xi(x)^{2}-1)/2\cdot G(\xi(x))\delta^{2}$, for some $x\leq\xi(x)\leq x+\delta.$

By (11) and Fact 5.9, to bound the magnitude of the $a_{k}$’s, it suffices to bound the terms $\int_{-\infty}^{\infty}P_{k}(x/C)xG(x)dx$ and $\int_{-\infty}^{\infty}P_{k}(x)(\xi(x)^{2}-1)/2\cdot G(\xi(x))dx$. This is done in the following two lemmas.

For $k\leq 4C$, we have that $\int_{-\infty}^{\infty}P_{k}(x/C)xG(x)dx\leq O(\sqrt{k}/C)$.

When $k$ is even, using Fact 5.3 (iv), we have that $P_{k}(x/C)xG(x)=-(P_{k}(-x/C)\cdot(-x)G(-x))$, and so the integral is zero. When $k$ is odd, we can rewrite Fact 5.3 (v) in ascending order of terms, by using the change of variables $j=(k+1)/2-i$, as

By standard results about the moments of Gaussians, for all $j\geq 1$, we have that $\int_{-\infty}^{\infty}x^{2j}G(x)dx=(2j+1)!!:=\prod_{i=1}^{j}(2i+1).$ Thus, we can write

Note that this quantity is non-negative. We can bound it from above as follows:

The proof of Lemma 5.10 is now complete. ∎

We separate this integral into the interval $[-C,C]$ and the tails. We can use Fact 5.3 (iii) to bound the integral on $[-C,C]$, as follows:

For the tails, we need Corollary 5.4(i). For the right tail, $[C,\infty)$, we have

A similar bound holds for the left tail, which completes the proof. ∎

Putting everything together, gives Lemma 5.5. ∎

To prove Proposition 5.2, we need to set $C$ appropriately and check the bounds on $m$ needed for $A(x)$ to satisfy the necessary properties.

Note that unless $\delta$ is sufficiently small and $m^{2}\leq O(\sqrt{\log(1/\delta)})$, taking $A=N(0,1)$, instead of using our construction, satisfies the proposition. We will take $C=\Theta(\sqrt{\log(1/\delta)})$, and so we can assume that $m\leq\sqrt{C}$.

Recall that $A(x)$ is defined to be $G(x-\delta)+p(x)$ on $[-C,C]$ and $G(x-\delta)$ outside of $[-C,C]$. Firstly, $A(x)$ needs to be the pdf of a distribution. Since

for $k=0$, when $P_{k}(x)=1$, we have that $\int_{-\infty}^{\infty}A(x)dx=1$. We also need that $A(x)$ is non-negative, i.e., that $A(x)=G(x-\delta)+p(x)\geq 0$ for all $x\in[-C,C]$. Note that

using Corollary 5.6. Since $m^{2}\leq C$, we need $G(C+\delta)\geq\delta C^{3/4}$. This holds when $C=\sqrt{\ln(1/\delta)}-\delta$, since then we have $G(C+\delta)=\sqrt{\delta/2\pi}\geq\delta\sqrt{\ln(1/\delta)}^{3/4}$ for sufficiently small $\delta$. Note that this also implies that $A(x)\leq 2G(x-\delta)$ for all $x$, and $|p(x)|\leq G(x)$ for all $-C\leq x\leq C$.

The second of these and Lemma 5.7 imply (iii). For (i), by construction, we have that the first $m$ moments agree.

The proof of Proposition 5.2 is now compete. ∎

2 Robust Learning Lower Bound for Unknown Covariance Gaussian

In this subsection, we prove an SQ lower bound for robustly learning the covariance matrix of a high-dimensional Gaussian with known mean. We note that our lower bound applies even for spectral norm approximation. In particular, we show:

The theorem will follow from the following proposition:

$A$ and $N(0,1)$ agree on the first $m$ moments.

$\chi^{2}(A,N(0,1))=O(1+m^{8}\delta^{3/2}/{\log(1/\delta)^{5/2}})$.

As in the previous subsection, Theorem 5.12 follows easily from Proposition 3.3 and Proposition 5.13.

Note that we cannot directly apply Proposition 3.3, since we are not aiming to learn within small total variation distance. Instead, we are interested in a different search problem, that of finding an approximation $\widetilde{\Sigma}$ to the covariance $\Sigma$ with $\|\widetilde{\Sigma}-\Sigma\|_{2}\leq\epsilon\log(1/\epsilon){/M^{4}}$, where $\Sigma$ is the covariance of a mean Gaussian within $\epsilon$ total variation distance. Note that for $\mathbf{P}_{v}$, we have that $\Sigma=I-(1-(1-\delta)^{2})vv^{T}$. We need to argue that this search problem has at most one solution in $S$, i.e., that for any $\widetilde{\Sigma}$, the set $\mathcal{Z}^{-1}(\widetilde{\Sigma})=\{\mathbf{P}_{v}:v\in S\text{ and }\|\widetilde{\Sigma}-I-(1-(1-\delta)^{2})vv^{T}\|_{2}\leq\epsilon\log(1/\epsilon){/M^{4}}\}$ has $|\mathcal{Z}^{-1}(\widetilde{\Sigma})|\leq 1$, where $S$ is as in Lemma 3.7.

For $S$ as in Lemma 3.7 with $c=1/6$ and with $n$ larger than a sufficiently large constant, $|\{\mathbf{P}_{v}:v\in S\textrm{ and }\|\widetilde{\Sigma}-I-(1-(1-\delta)^{2})vv^{T}\|_{2}\leq\epsilon\log(1/\epsilon){/M^{4}}\}|\leq 1$, for all $\widetilde{\Sigma}$.

Suppose for a contradiction that this set has size at least $2$ for some $\widetilde{\Sigma}$ and let $v,v^{\prime}$ be distinct elements. Let $\Sigma_{v}=I-(1-(1-\delta)^{2})vv^{T}$ and define $\Sigma_{v^{\prime}}$ similarly. Then we have $\|\widetilde{\Sigma}-\Sigma_{v}\|_{2}\leq\epsilon\log(1/\epsilon){/M^{4}}$ and $\|\widetilde{\Sigma}-\Sigma_{v^{\prime}}\|_{2}\leq\epsilon\log(1/\epsilon){/M^{4}}$. By the triangle inequality, we have $\|\Sigma_{v}-\Sigma_{v^{\prime}}\|_{2}\leq{2}\epsilon\log(1/\epsilon){/M^{4}}$. However, we also have that $|v\cdot v^{\prime}|\leq O(n^{-1/3})\leq 1/2$. Now we get that $v^{T}\Sigma v=(1-\delta)^{2}$, but $v^{T}\Sigma_{v^{\prime}}v=(1-|v\cdot v^{\prime}|^{2})\cdot 1+|v\cdot v^{\prime}|^{2}\cdot(1-\delta)^{2}\geq 3/4+(1/4)(1-\delta)^{2}$. We thus obtain

where the last inequality assumes that $\epsilon$ is at most an appropriately small universal constant. Since $\delta=2\epsilon\ln(1/\epsilon){/M^{4}}$, this leads to a contradiction. ∎

Similarly to the previous subsection, we choose to define the univariate distribution $A$ to have probability density function given by

where $C$ is a sufficiently small multiple of $\sqrt{\log(1/\delta)}$ and $p(x)$ is the unique degree-$m$ polynomial that causes $A$ and $G(x)$ (the pdf of $N(0,1)$) to have the same first $m$ moments. Once again, we may write $p(x)=\sum_{k=0}^{m}a_{k}P_{k}(x/C)$, where $a_{k}=((2k+1)/2C)\int_{-\infty}^{\infty}(G(x)-G(x/(1-\delta))/(1-\delta))P_{k}(x/C)dx$. The bulk of our proof will now be in bounding the $a_{k}$’s.

The first thing to note is that since $G(x)-G(x/(1-\delta))$ is even, $a_{k}$ is for $k$ odd. For $k$ even, we will need to compute this expression using Fact 5.3 (v). In particular, we have that

where in the last step we assume that $k$ is less than a sufficiently small multiple of $C$.

It is now clear that $A$ is a pseudo-distribution that matches its first $m$ moments with $N(0,1)$. Firstly, in order to check that $A$ is a distribution, it is clear that $A(x)>0$ for $|x|>C$. For $|x|\leq C$ we have that $|A(x)-G(x/(1-\delta))/(1-\delta)|\leq\sum_{k=0}^{m}|a_{k}|\leq\delta 10m^{4}C^{-3}$. Since this is smaller than $\delta^{1/2}<G(x/(1-\delta))/(1-\delta)$, we have that $A(x)\geq 0$ everywhere.