Faster Algorithms for High-Dimensional Robust Covariance Estimation

Introduction

In this paper, we study the outlier robust setting when a small constant fraction of our samples can be arbitrarily corrupted. We work in the following model of corruptions (see, e.g., [DKK+16]) that generalizes Huber’s contamination model ([Hub64]):

Recent work in TCS ([DKK+16, LRV16]) gave the first polynomial time robust estimators for a range of high-dimensional statistical tasks, including mean and covariance estimation. Since these initial papers ([DKK+16, LRV16]), a growing body of subsequent works have obtained polynomial-time robust learning algorithms for a variety of unsupervised and supervised high-dimensional models. (See Section 1.3 for more related work.)

It should be noted that the aforementioned robust estimators have already been useful in exploratory data analysis. Specifically, [DKK+17] evaluated the robust covariance estimators of [DKK+16] and [LRV16] to detect patterns in a well-known genetic dataset ([NJB+08]) in the presence of corruptions. Perhaps surprisingly, it was found that the robust algorithms developed for $\mathcal{N}(\mathbf{0},\Sigma)$ outperformed all previous approaches on this real dataset, essentially matching the setting where there are no corruptions at all.

Once a polynomial-time algorithm for a computational problem has been discovered, the next step is to focus on designing asymptotically faster algorithms for the problem – with linear time as the ultimate goal. We note that the aforementioned robust estimators ([DKK+16, LRV16]) are significantly slower than their non-robust counterparts (e.g., computing the empirical mean/covariance), hence may not be scalable when the dimension is very high. This raises the following natural question:

Can we design robust estimators that are as efficient as their non-robust analogues?

This direction was initiated in [CDG19] who gave a robust mean estimation algorithm with runtime $\widetilde{O}(Nd)/\poly(\epsilon)$, The $\widetilde{O}(\cdot)$ notation hides logarithmic factors in its argument. nearly matching the runtime of computing the empirical mean (when $\epsilon$ is constant).

For the sake of direct comparison, we note that the filtering-based robust covariance estimator of [DKK+16] has runtime $\Omega(N^{2}d)=\Omega(d^{5})/\poly(\epsilon)$. On the other hand, the recursive dimension-halving estimator of [LRV16] requires $\Omega(\log d)$ SVD computations of a $d^{2}\times d^{2}$ “covariance” matrix, hence has runtime $\Omega(d^{2\omega})$, where $\omega$ is the exponent of matrix multiplication. (Plugging in the best-known value for $\omega$ ([Gal14]) gives a runtime of $\Omega(d^{4.74})$.)

We note that the runtime of our algorithm, while being super-linear, essentially matches the best-known runtime to compute the empirical covariance matrix (see Section 6.1). Moreover, we provide evidence (Section 6.2) that this runtime may be a bottleneck even for the weaker task of obtaining an implicit representation to the output (by reweighing the input samples). It should be noted that all known computationally efficient robust estimators fit in this framework.

Our first algorithmic result states that we can robustly estimate the covariance matrix of a high-dimensional Gaussian within multiplicative, dimension-independent error, with running time that almost matches that of computing the empirical covariance matrix.

We also develop a related robust covariance estimation algorithm (Algorithm 3) with additive error guarantee, whose running time does not depend on the condition number of $\Sigma$.

We will prove Theorem 1.2 in Section 3.1 and Theorem 1.3 in Appendix B.1.

In Section 6, we provide evidence that the runtime of our algorithm may be difficult to improve. Specifically, we show that the best-known runtime for computing (or even approximating) the empirical covariance is $\Omega(d^{3.25})$. Moreover, even outputting a set of weights for the samples (such that the weighted empirical covariance works) seems to require $\Omega(d^{3.25})$ time with current methods.

2 Our Approach and Techniques

Iterative Refinement.

The second difficulty is that existing robust mean estimation algorithms rely on the assumption that the distribution of the good data either has known covariance or unknown bounded covariance. By reducing covariance estimation to the mean estimation of $X\otimes X$, we run into the difficulty that the covariance of such vectors would correspond to the fourth order moments of the original variables $X$. So, directly applying the mean estimation algorithms does not give our desired strong guarantees.

We establish two guarantees for the robust estimation of the mean of $Z$. In the first phase, we only use the fact that the covariance of $Z$ is bounded. Repeating this will give a rough estimation of the covariance of $X$. In the second phase, we use the fact that the current estimate $\Sigma_{t}$ is already very close to $\Sigma$, therefore $Y_{i}=\Sigma_{t}^{-1/2}X_{i}$ is very close to a standard Gaussian $N(0,I)$ for the good samples. In this case, we need to open up the algorithm in [CDG19] and prove stronger robust estimation guarantees tailored for the specific distribution of $Z=Y\otimes Y$. Our main algorithm (Algorithm 1) combines these two refinement steps to match the best-known robustness guarantees for covariance estimation.

Evidence of Hardness.

A natural way of improving the running time of non-robust covariance estimation is to try to approximate the product $X^{T}X$ using oblivious sketching (see, e.g., [Woo14] for a survey) which works roughly as follows. One samples a random $S$ from a certain family of random matrices, and computes $S\cdot X$ where $S$ has much fewer than $N$ rows. For structured random families of matrices, like fast JL matrices, $S\cdot X$ can be computed very quickly. Then one instead computes $X^{T}S^{T}SX$ with the guarantee that $\|X^{T}S^{T}SX-X^{T}X\|_{F}^{2}$ is small. Note that the matrix product $(X^{T}S^{T})\cdot(SX)$ can be performed more quickly if $S$ has a small number of rows. Unfortunately, for the guarantees we want, all known constructions of $S$ require $\Omega(N)$ rows. In fact, we prove an information-theoretic result that any oblivious sketching matrix $S$ must have $\Omega(N/\log N)$ rows in Lemma 6.1, thus ruling out this approach for achieving faster running time. Our proof uses arguments from communication complexity, arguing that such a family of sketching matrices would imply a better protocol for solving multiple copies of the Gap-Hamming communication problem.

Another way of trying to improve the runtime is to give an alternative definition of the problem: Instead of outputting a $d\times d$ matrix that is close to $\Sigma$, the algorithm outputs a set of nonnegative weights $w$ such that $\|w\|_{1}=1$, $\|w\|_{\infty}\leq\frac{1}{(1-\epsilon)N}$, and $\|\sum_{i=1}^{N}w_{i}X_{i}X_{i}^{T}-\Sigma\|_{F}=O(1)$. This bypasses the arguments above, since in the case of no corruption, we do not have to actually output $X^{T}X$ and can just set $w_{i}=1/N$ for all $i$. However, even for this relaxed version of the problem, we show that unless one can solve a certain “column norm” distinguishing problem faster than rectangular matrix multiplication, one cannot solve this problem faster than $(d,d^{2},d)$ matrix multiplication time. This problem can be intuitively stated as follows: the good samples are drawn from $\mathcal{N}(0,I)$, and the corrupted samples are drawn from a mean-zero Gaussian distribution with a very slight and known perturbation to the identity covariance matrix. One needs to identify a large fraction of these corrupted samples. Even though the covariance of the perturbed Gaussians is known, it is so slight that the norms of the corrupted samples are very similar to the uncorrupted ones. Therefore, one needs to measure these norms along certain directions, which requires computing a matrix product of the samples with a worst-case covariance matrix. We show that outputting the weights described above requires solving this problem, which we conjecture to be hard.

3 Related and Prior Work

Learning in the presence of outliers is an important goal in statistics and has been studied in the robust statistics community since the 1960s ([Hub64]). After several decades of work, a number of sample-efficient and robust estimators have been discovered. The reader is referred to ([DKK+16, LRV16]) for a detailed summary of this line of work. Until recently, all known computationally efficient high-dimensional estimators could only tolerate a negligible fraction of outliers. Recent work ([DKK+16, LRV16]) gave the first efficient robust estimators for basic high-dimensional unsupervised tasks. Since these works, there has been a flurry of research activity on robust learning algorithms in both supervised and unsupervised settings ([BDLS17, CSV17, DKK+17, DKS17, DKK+18, SCV18, DKS18b, DKS18a, HL18, KSS18, PSBR18, DKK+19, KKM18, DKS19, LSLC18, CDKS18]).

The most relevant prior work is that of [CDG19], initiating the direction of obtaining fast algorithms for robust high-dimensional estimation. For the problem of robust mean estimation, [CDG19] proposed a primal-dual approach – building on the convex programming approach of [DKK+16] – yielding an algorithm with runtime $\widetilde{O}(Nd)/\poly(\epsilon)$. This improved on the $\widetilde{O}(Nd^{2})$ runtime of the iterative filtering method in [DKK+16]. Our algorithm uses the same primal-dual framework; however, we emphasize that a standard application of their framework would only lead to a runtime of $\widetilde{O}(Nd^{2})/\poly(\epsilon)$. To obtain our improved runtime, we need to overcome a number of technical obstacles, as we explained in Section 1.2.

Preliminaries

Connections Between the Second and Fourth Moments of a Gaussian.

Fast Rectangular Matrix Multiplication.

We will frequently use fast rectangular matrix multiplication in our algorithms. Let $(a,b,c)$-matrix multiplication time denote the time it takes to multiply an $a\times b$ matrix with a $b\times c$ matrix. It is forklore (see, e.g., [LR83]) that $(a,a,b)$, $(a,b,a)$, and $(b,a,a)$ matrix multiplications require the same number of arithmetic operations. More specifically, we use the algorithm proposed in [Gal12]. They obtained new upper bounds on $(n,n^{\alpha},n)$ matrix multiplication time. We use a special case of their result ($\alpha=2$).

We can compute $(n,n^{2},n)$-matrix multiplication in time $O(n^{3.26})$. This implies that for $d>0$ and $N=\widetilde{O}(d^{2}/\epsilon^{2})$, we can multiply a $d\times N$ matrix with an $N\times d$ matrix in time $\widetilde{O}(d^{3.26}/\epsilon^{2})$.

To multiply a $d\times N$ matrix with an $N\times d$ matrix when $N=\widetilde{O}(d^{2}/\epsilon^{2})$, we split the matrices into blocks of size $d\times d^{2}$ or $d^{2}\times d$, multiply each pair of matrices and then add the results together. The total running time is $\frac{N}{d^{2}}\cdot O(d^{3.26}+d^{2})=\widetilde{O}(d^{3.26}/\epsilon^{2})$.

Transposition Principle for Matrix-Vector Multiplication.

The transposition principle (see, e.g., [Bor57, Fid73]) plays a central role in our faster implementation of positive SDP solvers. It states that matrix-vector multiplication by $A^{\top}$ has almost exactly the same computational complexity as matrix-vector multiplication by $A$.

Estimating the Covariance of a Gaussian Distribution

In this section, we present our key structural and computational lemmas, and use them to prove our main algorithmic results (Theorems 1.2 and 1.3).

We first present our algorithm (Algorithm 1) for robustly estimating the covariance of Gaussian distributions with multiplicative error guarantees. Algorithm 1 starts with an upper bound $\Sigma_{0}$ on the true covariance matrix $\Sigma$, and iteratively compute more and more accurate upper bounds $\Sigma_{t}\succeq\Sigma$.

First we need a reasonable starting point before we can run any iterative refinement steps.

Consider the same setting as in Theorem 1.2. We can compute a matrix $\Sigma_{0}$ in $\widetilde{O}(d^{3.26}/\epsilon^{2})$ time such that, with high probability, $\Sigma\preceq\Sigma_{0}\preceq(\kappa\poly(d))\Sigma$ and $\left\|\Sigma_{0}\right\|_{2}\leq\poly(d)\left\|\Sigma\right\|_{2}$.

We use two different iterative refinement steps (Lemmas 3.2 and 3.3), which correspond to the two loops in Algorithm 1. In the first phase (Lemma 3.2), we only have a crude upper bound on $\Sigma$.

The first phase can only converge to a matrix $\Sigma_{T_{1}}$ with $\Sigma\preceq\Sigma_{T_{1}}\preceq(1+O(\sqrt{\epsilon})\Sigma$. In the second phase (Lemma 3.3), we already have a fairly accurate estimate of $\Sigma$, so the refinement steps converge faster and eventually we can get to a matrix $\Sigma_{T_{2}}$ with $\Sigma\preceq\Sigma_{T_{2}}\preceq(1+O(\epsilon\log(1/\epsilon))\Sigma$.

Consider the same setting as in Theorem 1.2. Let $0<\tau_{t}<\tau_{0}$ for some universal constant $\tau_{0}$. Given $\tau$ and $\Sigma_{t}$ with $\Sigma\preceq\Sigma_{t}\preceq(1+\tau_{t})\Sigma$, we can compute in time $\widetilde{O}(d^{3.26}/\epsilon^{8})$ an upper bound matrix $\Sigma_{t+1}$ and a hypothesis matrix $\widehat{\Sigma}$ such that, with high probability, for $\tau_{t+1}=O(\sqrt{\epsilon\tau}+\epsilon\log(1/\epsilon))$,

We defer the proof of Lemma 3.1 to Appendix B, and the proofs of Lemmas 3.2 and 3.3 to Section 3.2. We first use these three lemmas to prove Theorem 1.2 (correctness and runtime of Algorithm 1).

For any integer $t\geq 0$, given the upper bound matrix $\Sigma_{t}$, we can use Lemma 3.2 to obtain a better upper bound $\Sigma_{t+1}$ such that $\Sigma_{t+1}\preceq\Sigma+O(\sqrt{\epsilon})\Sigma_{t}$. Since $\Sigma_{0}\preceq\kappa\poly(d)\Sigma$, after $T_{1}=O(\log\kappa+\log d)$ iterations, we have a matrix $\Sigma_{T_{1}}$ with $\Sigma\preceq\Sigma_{T_{1}}\preceq\left(1+O(\sqrt{\epsilon})\right)\Sigma$.

At this point we have a pretty accurate upper bound $\Sigma_{T_{1}}$ with $\Sigma\preceq\Sigma_{T_{1}}\preceq(1+\tau_{T_{1}})\Sigma$, where $\tau_{T_{1}}=O(\sqrt{\epsilon})$. For any integer $t\geq T_{1}$, given $\Sigma_{t}$ and $\tau_{t}$, we can use Lemma 3.3 to obtain a better upper bound matrix $\Sigma_{t+1}$ such that $\Sigma\preceq\Sigma_{t+1}\preceq\Sigma+\tau_{t+1}\Sigma_{t}$, where $\tau_{t+1}=O(\sqrt{\epsilon\tau_{t}}+\epsilon\log(1/\epsilon))$. Similar to the previous step, after $O(\log\log(1/\epsilon))$ iterations, we have a matrix $\Sigma_{T_{2}}$ such that $\Sigma\preceq\Sigma_{T_{2}}\preceq(1+\tau_{T_{2}})\Sigma$, where $\tau_{T_{2}}=O(\epsilon\log(1/\epsilon))$.

Finally, using Lemma 3.3 one more time with $\Sigma_{T_{2}}$ and $\tau_{T_{2}}$, we can get a matrix $\widehat{\Sigma}$ with

We note that both Lemmas 3.2 and 3.3 hold with high probability, so we can take a union bound over the failure probabilities and conclude that with high probability all iterative refinement steps are successful, and therefore, we can return $\widehat{\Sigma}$ as our final answer.

Now we analyze the running time of Algorithm 1. We call Lemma 3.1 once to compute $\Sigma_{0}$, which takes time $\widetilde{O}(d^{3.26}/\epsilon^{2})$. After that, we use two iterative refinement steps. The total number of iterations is $O(\log\kappa+\log d+\log\log(1/\epsilon))$. In each iteration, we invoke either Lemma 3.2 or 3.3. Since both lemmas run in time $\widetilde{O}(d^{3.26}/\epsilon^{8})$, the overall running time is $\widetilde{O}(d^{3.26}/\epsilon^{2})+O(\log\kappa+\log d+\log\log(1/\epsilon))\cdot\widetilde{O}(d^{3.26}/\epsilon^{8})=\widetilde{O}(d^{3.26}/\epsilon^{8})$. ∎

2 Implementing the Iterative Refinement Steps

In the first phase (Lemma 3.2), we only have a crude upper bound on $\Sigma$. For any $\Sigma_{t}\succeq\Sigma$, we have $\Sigma_{Y}=\Sigma_{t}^{-1/2}\Sigma\Sigma_{t}^{-1/2}\preceq I$, which implies $\Sigma_{Z}\preceq 2I$ (Lemma 2.1). Because $Z$ has bounded covariance, we can use the following robust mean estimation algorithm from [CDG19].

In the second phase (Lemma 3.3), we use the fact that the current estimate $\Sigma_{t}$ is already very close to $\Sigma$, therefore $Y=\Sigma_{t}^{-1/2}X$ is very close to $\mathcal{N}(0,I)$. In this case we need an algorithm with stronger robust estimation guarantees tailored for the specific distribution of $Z=Y\otimes Y$.

It is worth noting that we cannot use these mean estimation algorithms in a black-box manner. This is because writing down the input explicitly takes $\Omega(Nd^{2})=\widetilde{\Omega}(d^{4}/\epsilon^{2})$ time, and these algorithms run in time $\widetilde{\Omega}(Nd^{2})/\poly(\epsilon)$ in $d^{2}$ dimensions. One of our main contributions is to show that it is possible to open up these algorithms and take advantage of the additional structures of our inputs (they all have the form $Y_{i}\otimes Y_{i}$) to implement both algorithms to run in time $\widetilde{O}(d^{3.26})/\poly(\epsilon)$.

We give a description of the algorithm for Lemma 3.4 in Appendix B.2 (Algorithm 5). We prove Lemma 3.5 and present the corresponding algorithm (Algorithm 2) in Section 4. In Section 5, we show that both algorithms can be implemented to run in time $\widetilde{O}(d^{3.26})/\poly(\epsilon)$ (Proposition 3.6). We first use Lemmas 3.4, 3.5, and Proposition 3.6 to prove the iterative refinement lemmas.

Given an upper bound $\Sigma_{t}$ on the true covariance matrix, we can rotate the input samples to compute $Y_{i}=\Sigma_{t}^{-1/2}X_{i}$. Let $Z_{i}=Y_{i}\otimes Y_{i}$. Note that when $X\sim\mathcal{N}(0,\Sigma)$, the random variable $Y=\Sigma_{t}^{-1/2}X$ is drawn from a Gaussian distribution with covariance $\Sigma_{Y}=\Sigma_{t}^{-1/2}\Sigma\Sigma_{t}^{-1/2}\preceq I$. Lemma 2.1 implies that, if $\Sigma_{Y}\preceq I$, then $\Sigma_{Z}\preceq 2I$. Therefore, $(Z_{i})_{i=1}^{N}$ is an $\epsilon$-corrupted set of samples drawn from a distribution with bounded covariance, so we can apply Algorithm 5 to robustly estimate its mean.

Let $\widehat{\Sigma}=\Sigma_{t}^{1/2}M\Sigma_{t}^{1/2}$. Using ${\left\|AB\right\|}_{F}\leq\left\|A\right\|_{2}{\left\|B\right\|}_{F}$, we can prove the first part of the lemma,

As a result, we have that $-O(\sqrt{\epsilon})I\preceq\Sigma_{t}^{-1/2}(\widehat{\Sigma}-\Sigma)\Sigma_{t}^{-1/2}\preceq O(\sqrt{\epsilon})I$, or equivalently

Now $\Sigma_{t+1}=\widehat{\Sigma}+O(\sqrt{\epsilon})\Sigma_{t}$ is a better upper bound, which satisfies $\Sigma\preceq\Sigma_{t+1}\preceq\Sigma+O(\sqrt{\epsilon})\Sigma_{t}$.

Since all the $Z_{i}$’s have the form $Y_{i}\otimes Y_{i}$, Proposition 3.6 shows that Algorithms 5 has running time $\widetilde{O}(d^{3.26}/\epsilon^{8})$. Given the output of Algorithms 5, we can compute the new upper bound $\Sigma_{t+1}$ in time $O(d^{\omega})$ using a constant number of $d\times d$ matrix additions and multiplications. ∎

Let $Y_{i}=\Sigma_{t}^{-1/2}X_{i}$ and $Z_{i}=Y_{i}\otimes Y_{i}$. We know that

By Lemma 2.1, we have $\left\|\Sigma_{Z}-2I\right\|_{2}\leq O(\tau)$. By standard concentration results, $Z=Y\otimes Y$ has exponential concentration about its mean in any direction. Therefore, $(Z_{i})_{i=1}^{N}$ is an $\epsilon$-corrupted set of samples drawn from a distribution that satisfies the conditions in Lemma 3.5, and we can apply Algorithm 2 to robustly learn the mean of $Z$. By Lemma 3.5, we can compute a matrix $M$ such that

Let $\widehat{\Sigma}=\Sigma_{t}^{1/2}M\Sigma_{t}^{1/2}$. Using ${\left\|AB\right\|}_{F}\leq\left\|A\right\|_{2}{\left\|B\right\|}_{F}$, we can prove the first part of the lemma,

This gives a better upper bound $\Sigma_{t+1}=\widehat{\Sigma}+O(\sqrt{\epsilon\tau}+\epsilon\log(1/\epsilon))\Sigma_{t}$ such that $\Sigma\preceq\Sigma_{t+1}\preceq\Sigma+O(\sqrt{\epsilon\tau}+\epsilon\log(1/\epsilon))\Sigma_{t}$.

We omit the running time analysis because it is identical to the one in the previous proof. ∎

Robust Mean Estimation Subroutines

We first present the robust mean estimation algorithm (Algorithm 2) that achieves Lemma 3.5.

There are two obstacles for applying the algorithmic framework of [CDG19] to our setting. First, the input samples $Z_{i}=X_{i}\otimes X_{i}$’s are $d^{2}$-dimensional vectors. Writing down these vectors explicitly takes time $\Omega(Nd^{2})=\Omega(d^{4})$. Therefore, we want to solve the SDPs (1) and (2) on input $Z_{i}$ without computing them explicitly. We resolve this issue in Section 5 (Proposition 3.6).

Second, their algorithms have error $O(\sqrt{\epsilon})$ for bounded-covariance distributions, and error $O(\epsilon\sqrt{\log(1/\epsilon)})$ for sub-gaussian distributions with identity covariance matrix. While we can directly use their result for bounded-covariance distributions for Lemma 3.4, we need to develop a new algorithm for Lemma 3.5. In Lemma 3.5, we have a distribution with exponential decaying tails, and we know its covariance is $\tau$-close to the identity matrix. We want to robustly estimate its mean, with optimal error guarantees that depend on both $\epsilon$ and $\tau$. We generalize the analysis of [CDG19] to handle this case. Lemma 3.5 is proved in Appendix B.3.

Faster Implementation of Robust Mean Estimation with Tensor Inputs

The bottleneck of both Algorithms 5 and 2 are solving SDPs (1) and (2). In this section, we prove Proposition 3.6, which states that when all input samples have the tensor-product form $Y\otimes Y$, we can solve these SDPs in time $\widetilde{O}(d^{3.26})/\poly(\epsilon)$.

We first convert the SDPs (1) and (2) into packing/covering SDPs as follows.

Here $\rho$ is a binary search parameter that is between $\frac{1}{d}$ and $1$. At the core of nearly-linear time width-independent SDP solvers (e.g., [ALO16, PTZ16]) is an application of matrix multiplicative weight update, where the algorithm maintains a weighted sum $\Psi$ of the matrices. In iteration $t$, we have $\Psi^{t}=\sum_{i=1}^{n}w_{i}A_{i}$, and we will update the weights based on the values of $A_{i}\bullet\frac{\exp(\Psi^{t})}{\tr(\exp(\Psi^{t}))}$.

Let $A_{1},\ldots,A_{n}$ be $m\times m$ PSD matrices given in factorized form $A_{i}=C_{i}C_{i}^{\top}$. Consider the following pair of packing and covering SDPs:

We will approximate each $\exp(\Psi)\bullet A_{i}$ and $\tr(\exp(\Psi))$ separately. Observe that $\Psi=\sum_{i=1}^{n}w_{i}A_{i}$ and $\exp(\Psi)$ have the same block structure as the $A_{i}$’s. Due to the special structure of the bottom-right block, we can compute its contribution to $\tr(\exp(\Psi))$ and $\exp(\Psi)\bullet A_{i}$ exactly. Therefore, we can focus on the top-left block. Moreover, because the goal is to compute a multiplicative approximation of the top-left block’s contribution to $\tr(\exp(\Psi))$ and $\exp(\Psi)\bullet A_{i}$, we can ignore the scalar $\rho$. We prove the following lemma.

$(Z-\nu{\bf 1}^{\top})s=Zs-({\bf 1}^{\top}s)\nu$ and we can compute $Zs=\left(YD_{s}Y^{\top}\right)^{\flat}$ via fast rectangular matrix multiplication (Lemma 2.2),

matrix-vector multiplication with $(Z-\nu{\bf 1}^{\top})$ or $(Z-\nu{\bf 1}^{\top})^{\top}$ has the same running time by the Transposition Principle (Lemma 2.3 in Section 2)), and

multiplication with a diagonal matrix $D_{w}$ can be done in time $O(N)$.

We approximate $\exp(\Psi)\bullet(Z_{i}-\nu)(Z_{i}-\nu)^{\top}$ using a similar approach: $\exp(\Psi)\bullet(Z_{i}-\nu)(Z_{i}-\nu)^{\top}=\left\|\exp(\Psi/2)(Z_{i}-\nu)\right\|_{2}^{2}\approx_{\epsilon/4}\left\|M(Z_{i}-\nu)\right\|_{2}^{2}\approx_{\epsilon/4}\left\|QM(Z_{i}-\nu)\right\|_{2}^{2}$. Notice that the last line is precisely the squared norm of the $i$-th column of $QM(Z-\nu{\bf 1}^{\top})$. For the same reasons as in the previous case, we can compute this matrix in time $\widetilde{O}(d^{3.26}/\epsilon^{5})$. ∎

We can approximate $\exp(A)$ with a matrix polynomial of $A$, whose degree depends on the spectral norm of $A$ and the desired precision (see, e.g., [AK16]).

By Lemma 5.1, we only need to show that the required oracle algorithm can be implemented in time $\widetilde{O}(d^{3.26}/\epsilon^{5})$. We approximate $\tr(\exp(\Psi))$ and each $\exp(\Psi)\bullet A_{i}$ separately. Given $\Psi=\sum_{i=1}^{N}w_{i}A_{i}$, we will compute the contribution from bottom-right block explicitly, and use Lemma 5.2 for the top-left block. The bottom-right block adds $\sum_{i=1}^{N}\exp(w_{i}(1-\epsilon)N)$ to $\tr(\exp(\Psi))$, and for every $i$, it adds $w_{i}\exp(w_{i}(1-\epsilon)N)$ to $\exp(\Psi)\bullet A_{i}$. ∎

Evidence of Hardness

In this section, we provide some evidence which suggests that the running time of our algorithm has near-optimal dependence on $d$. We start by noting that our sample complexity $N=\widetilde{\Omega}(d^{2}/\epsilon^{2})$ is tight up to polylogarithmic factors, and this holds even when there is no corruption. For the rest of this section, we will assume both $\epsilon$ and $\kappa$ are constants, and focus on the dependence on $d$ in the running time. Since the running time of our algorithm is dominated by $(d,d^{2},d)$-matrix multiplication time, faster matrix multiplication algorithms time will improve our running time.

In Section 6.1, we show that even when there are no corrupted samples, it is not known how to compute the empirical covariance matrix faster than $(d,d^{2},d)$-matrix multiplication time. We give a communication complexity lower bound that rules out all oblivious matrix sketching approaches.

In Section 6.2, to circumvent the difficulty raised in Section 6.1, we consider a weaker problem where the algorithm only need to find a set of good weights (instead of a $d\times d$ matrix). We give a reduction to show that this problem is still at least as hard as some basic matrix computation question, which we do not know how to solve faster than $(d,d^{2},d)$-matrix multiplication time.

Our algorithm matches the running time of the best non-robust covariance estimation algorithm. When there are no corrupted samples and $N=\widetilde{\Omega}(d^{2}/\epsilon^{2})$, with high probability, the empirical second-moment matrix $\frac{1}{N}\sum_{i=1}^{N}X_{i}X_{i}^{\top}$ is $\epsilon$-close to the true covariance matrix in Frobenius norm. However, it is not known how to (approximately) compute this empirical second-moment matrix faster than $(d,d^{2},d)$ matrix multiplication time.

For approximate matrix product of an $N\times d$ matrix $A$ with $\left\|A\right\|_{2}=O(1)$, we want to choose a sketching matrix $S$ so that ${\left\|A^{\top}S^{\top}SA-A^{\top}A\right\|}_{F}^{2}=O(1)$. Known results for approximate matrix product state that if $S$ has $s$ rows, then ${\left\|A^{\top}S^{\top}SA-A^{\top}A\right\|}_{F}^{2}=O\left(\frac{{\left\|A\right\|}_{F}^{4}}{s}\right)$ with probability at least $9/10$, see, e.g., Section 2.2 of [Woo14] for a survey. In the context of Problem 1, letting $A=\frac{1}{\sqrt{N}}X$, we have $\left\|A\right\|_{2}=O(1)$ and ${\left\|A\right\|}_{F}^{4}=O(d^{2})$. The error is then $O(d^{2}/s)$, and consequently $S$ must have $s=\Omega(d^{2})$ rows for the error to be at most $O(1)$.

We can show that the argument above is almost tight for all oblivious sketches.

Let $N=d^{2}$. There is no distribution over $t\times N$ matrices $S$, oblivious to the underlying input $N\times d$ matrix $A$, where $t=o(d^{2}/\log d)$, such that with probability at least $2/3$, it holds that $\|A^{\top}S^{\top}SA-A^{\top}A\|_{F}^{2}\leq C_{1}\frac{\|A\|_{F}^{4}}{d^{2}}$, where $C_{1}=4\cdot 25\cdot 2000^{2}$ is a positive constant.

Suppose, to the contrary, there were such a distribution on matrices $S$ satisfying $t=o(d^{2}/\log d)$.

For $N=d^{2}$, consider a uniformly random $N\times d$ matrix $A\in\{-\frac{1}{d},\frac{1}{d}\}^{N\times d}$. Then $\|A\|_{F}^{2}=d$, and so for a random matrix $S$ from our family and a random input $A$ from this family of inputs, it holds that with probability at least $\frac{2}{3}$, $\|A^{\top}S^{\top}SA-A^{\top}A\|_{F}^{2}\leq C_{1}\frac{\|A\|_{F}^{4}}{d^{2}}=C_{1}$. By anti-concentration of the binomial distribution, with probability at least $\frac{99}{100}$, at least a $\frac{99}{100}$-fraction of the off-diagonal entries of $A^{\top}A$ have absolute value at least $\frac{1}{1000d}$.

Consequently, for at least a $\frac{24}{25}$-fraction of the entries in the bottom left $\frac{d}{2}\times\frac{d}{2}$ submatrix of $A^{\top}A$, we have the property that the entry has the same sign as in $A^{\top}S^{\top}SA$, and also the entries in $A^{\top}S^{\top}SA$ are at least $\frac{1}{2000d}$. Indeed, otherwise we would have $\|A^{\top}S^{\top}SA-A^{\top}A\|_{F}^{2}>(\frac{d}{2})^{2}\cdot\frac{1}{25}\cdot(\frac{1}{1000d}-\frac{1}{2000d})^{2}>C_{1}$ with probability at least $\frac{99}{100}$ over the choice of $A$ and $S$, and in particular there exists a fixed $A$ for which this holds with probability at least $\frac{99}{100}$ over the choice of $S$, contradicting our assumption on the family of matrices $S$.

Now consider the following two-player communication game with public shared randomness. Alice has the first $\frac{d}{2}$ columns of $A$, denoted $A_{L}\in\{-1,1\}^{N\times d/2}$ while Bob has the remaining $\frac{d}{2}$ columns of $A$, denoted $A_{R}\in\{-1,1\}^{N\times d/2}$. The entries in the in the bottom left $\frac{d}{2}\times\frac{d}{2}$ submatrix of $A^{\top}A$ are exactly the inner products between all columns of Alice and all columns of Bob. Suppose there were such a family of matrices $S$ as described above. Alice and Bob use the public coin to agree upon $S$ with no communication. Alice then computes $S\cdot A_{L}$, and rounds each entry to the nearest power of $(1+\frac{1}{\poly(d)})$. Note all entries of $S\cdot A_{L}$ need to be at most $\poly(d)$ and rounding preserves $\|A^{\top}S^{\top}SA-A^{\top}A\|_{F}^{2}$ up to additive $\frac{1}{\poly(d)}$. Therefore, we maintain the property that, at least a $\frac{24}{25}$-fraction of the entries in the bottom left $\frac{d}{2}\times\frac{d}{2}$ submatrix of $A^{\top}A$ have the same signs in $A^{\top}A$ and $A^{\top}S^{\top}SA$. Alice sends each of the rounded entries of $S\cdot A_{L}$, which is $\Theta(td\log d)=o(d^{3})$ bits. Bob then computes $S\cdot A_{R}$ and thus forms $S\cdot A$, from which he can compute $A^{\top}S^{\top}SA$. At this point, Bob can recover the sign of a uniformly random entry in the bottom left $d/2\times d/2$ submatrix of $A^{\top}A$ with probability at least $2/3-1/25-1/100>3/5$.

Notice that the sign of such an entry is the same as solving the Gap-Hamming communication problem under the uniform distribution: in this communication problem there are two players, Alice and Bob, who hold uniformly random vectors $x,y\in\{-1,1\}^{N}$, respectively, and wish to decide if $\langle x,y\rangle>0$ or $\langle x,y\rangle<0$. This problem requires $\Omega(N)$ randomized communication complexity [CR12]. Moreover, as shown by Braverman et al. [BGPW16], the information complexity of this problem is $\mathcal{I}=\Omega(N)$ bits. In our setting, we can think of Alice as having $\frac{d}{2}$ independent instances $x^{1},\ldots,x^{d/2}$, and Bob having an index $i\in\{1,2,\ldots,\frac{d}{2}\}$ as well as a vector $y$ and Bob wants to solve the Gap-Hamming problem on the pair $(x^{i},y)$. However, only Alice is allowed to speak, and she sends a single message to Bob, without knowing $i$. By standard direct sum arguments in communication complexity [BR11] (see also [PSW14] where Gap-Hamming composed with the Index problem was used), the randomized one-way communication complexity of this problem is $\Omega(d\cdot\mathcal{I})=\Omega(d^{3})$ bits. However, the communication cost of our protocol is $\Theta(td\log d)=o(d^{3})$ bits, which is a contradiction. Consequently, we must have $t=\Omega(d^{2}/\log d)$, as desired. ∎

It is worth noting that this lower bound holds for any possible algorithm one can run on $SA$ (i.e., the algorithm can do more than just computing $A^{\top}S^{\top}SA$), so it is a stronger information-theoretic statement.

2 Finding Good Weights

To circumvent the difficulty of Problem 1, we could redefine our problem so that the algorithm does not need to output a $d\times d$ matrix, instead it outputs a set of good weights $w$ such that $\|\sum_{i=1}^{N}w_{i}X_{i}X_{i}^{\top}-\Sigma\|_{F}=O(\sqrt{\epsilon})$. We will show that, even for this weaker problem of finding good weights, one still need to come up with faster algorithms for a basic matrix problem.

Consider the following instance. Let $U$ be an arbitrary $d\times\frac{d}{2}$ matrix with orthonormal columns. Let the good distribution be $D=\mathcal{N}(0,\Sigma=I)$, and the noise distribution $D^{\prime}$ is defined as

We draw $(1-\epsilon)N$ samples from $D$ and $\epsilon N$ samples from $D^{\prime}$. The empirical covariance matrix of the mixed distribution is $\widehat{\Sigma}=(1-\epsilon)\Sigma+\epsilon\Sigma^{\prime}=\left(1-\frac{c}{\sqrt{d}+1/\epsilon}\right)I+\left(\frac{2c}{\sqrt{d}+1/\epsilon}\right)UU^{\top}$. Observe that ${\left\|\widehat{\Sigma}-\Sigma\right\|}_{F}^{2}=d\left(\frac{c}{\sqrt{d}+1/\epsilon}\right)^{2}=\Omega(1)$, so the bad samples are distorting the empirical covariance matrix by more than we could tolerate.

Therefore, a natural way of distinguishing them is to compute $U^{\top}X$, which requires $(d,d^{2},d)$-matrix multiplication time. We could compute the column norms of $SU^{\top}A$, where $S$ is a Johnson-Lindenstrauss matrix. However, $S$ must have $\frac{1}{\epsilon^{2}}$ rows to obtain $(1+\epsilon)$-approximation, and therefore $S$ must have $\widetilde{\Omega}(d^{2})$ rows. Even if one uses a sparse matrix $S$, one has that $SU^{\top}$ is a dense matrix, and it is unclear how to compute $SU^{\top}A$ quickly.

Consider the same setting as in Problem 2. Given a set of weights $w$ such that $\left\|w\right\|_{1}=1$, $\left\|w\right\|_{\infty}\leq\frac{1}{(1-\epsilon)N}$, and $\|\sum_{i=1}^{N}w_{i}X_{i}X_{i}^{\top}-I\|_{F}=O(1)$, we can solve Problem 2 in $O(N)$ time.

For the rest of proof we assume the samples meet these conditions.

Let $\Sigma_{w}=\sum_{i=1}^{N}w_{i}X_{i}X_{i}^{\top}$. Let $w_{G}$ and $w_{B}$ denote the total weights on $G$ and $B$ respectively. Since ${\left\|\Sigma_{w}-I\right\|}_{F}=O(1)$, by Cauchy-Schwarz,

Putting these two inequalities together, we get that $w_{B}\leq\frac{\epsilon}{4}$. In other words, the average weight of a bad sample is $\frac{1}{4N}$.

Let $S=\{i\in[N]:w_{i}\leq\frac{1}{2N}\}$. By Markov’s inequality, we have $|S\cap B|\geq\frac{|B|}{2}$. Since $\left\|w\right\|_{\infty}\leq\frac{1}{(1-\epsilon)N}$ and $w_{G}=\left\|w\right\|_{1}-w_{B}\geq 1-\frac{\epsilon}{4}$, again by Markov’s inequality, we get that $|S\cap G|\leq\epsilon N$ and hence $|S|\leq|B|+\epsilon N=2\epsilon N$. ∎

Acknowledgments

This work was done in part while some of the authors were visiting the Simons Institute for the Theory of Computing. Ilias Diakonikolas was supported by NSF Award CCF-1652862 (CAREER) and a Sloan Research Fellowship. Rong Ge is supported by NSF Award CCF-1704656, CCF-1845171 (CAREER), a Sloan Research Fellowship, and a Google Faculty Research Award. David Woodruff was supported in part by Office of Naval Research (ONR) grant N00014-18-1-2562.

References

Appendix A Omitted Proofs from Section 2

If $\Sigma\preceq I$, then $\Sigma_{Z}\preceq 2I$.

If $0\leq\tau<1$ and $\left\|\Sigma-I\right\|_{2}\leq\tau$, then $\left\|\Sigma_{Z}-2I\right\|_{2}\leq 6\tau$.

Let $A$ be the unique matrix such that $A^{\flat}=v$. We have

Note that $\Sigma$ is a covariance matrix, so it is always symmetric and PSD. We can write $\Sigma$ as $\Sigma=\sum_{i=1}^{d}\lambda_{i}v_{i}v_{i}^{\top}$. Let $\lambda_{\max}$ and $\lambda_{\min}$ denote the maximum and minimum eigenvalues of $\Sigma$. Let $\widehat{A}=\frac{A+A^{\top}}{2}$. The right hand side is equal to

For (i), by assumption $\lambda_{\max}=\left\|\Sigma\right\|_{2}\leq 1$, so we have $\left\|\Sigma_{Z}\right\|_{2}\leq 2$.

For (ii), we know $1-\tau\leq\lambda_{\min}\leq\lambda_{\max}\leq 1+\tau$ and $0<\tau<1$. It follows that $(1-2\tau)2I\preceq\Sigma\preceq(1+3\tau)2I$, and thus $\left\|\Sigma_{Z}-2I\right\|_{2}\leq 6\tau$. ∎

Appendix B Omitted Proofs from Section 3

Let $(G_{i})_{i=1}^{N}$ be the original set of good samples drawn from $\mathcal{N}(0,\Sigma)$, and let $(X_{i})_{i=1}^{N}$ be the corrupted samples. Let $S$ denote the set of $(1-\epsilon)N$ samples with the smallest norm $\left\|X_{i}\right\|_{2}$. We define $\Sigma_{0}=2\left(\frac{1}{N}\sum_{i\in S}X_{i}X_{i}^{\top}\right)$.

We first show that $\Sigma_{0}\succeq\Sigma$ with high probability. Since the adversary corrupts at most $\epsilon N$ samples and we throw away $\epsilon N$ samples, we are left with at least $(1-2\epsilon)N$ good samples in $S$. We will use the fact that removing any $(2\epsilon)$-fraction of the good samples will not change the empirical covariance too much. Let $Y_{i}=\Sigma^{-1/2}G_{i}$ so that if $G_{i}\sim\mathcal{N}(0,\Sigma)$ then $Y_{i}\sim\mathcal{N}(0,I)$. When $N=\widetilde{\Omega}(d/\epsilon^{2})$, for any $T\subset[N]$ with $|T|=(1-2\epsilon)N$, we have that with high probability,

We set $T\subseteq S$ to be a set of $(1-2\epsilon)N$ good samples in $S$, i.e., $G_{i}=X_{i}$ for all $i\in T$. Let $M=\frac{1}{N}\sum_{i\in T}Y_{i}Y_{i}^{\top}$. We know that $M=\frac{1}{N}\sum_{i\in T}\Sigma^{-1/2}G_{i}G_{i}^{\top}\Sigma^{-1/2}$ by definition, and $M\succeq(1-O(\epsilon\log(1/\epsilon)))I\succeq\frac{I}{2}$ by the above concentration inequality. Therefore,

Next we show that $\left\|\Sigma_{0}\right\|_{2}\leq\poly(d)\left\|\Sigma\right\|_{2}$ and $\Sigma_{0}\preceq(\kappa\poly(d))\Sigma$. Let $\sigma^{2}$ denote the largest eigenvalue of $\Sigma$. Again let $Y_{i}=\Sigma^{-1/2}G_{i}$, we know that when $N=\widetilde{\Omega}(d/\epsilon^{2})$, with high probability,

We assume this condition holds for the rest of the proof. As a result, $\left\|G_{i}\right\|_{2}=\left\|\Sigma^{1/2}Y_{i}\right\|_{2}\leq O(\sigma\sqrt{d\log d})$ for all $i$. Since only corrupted samples can have larger norm, and we remove the $\epsilon N$ samples with the largest norm, all samples in $S$ have norm at most $O(\sigma\sqrt{d\log d})$. This gives an upper bound on the spectral norm of $\Sigma_{0}$,

This proves $\left\|\Sigma_{0}\right\|_{2}\leq\poly(d)\left\|\Sigma\right\|_{2}$. Moreover, by the definition of condition number we know that $\Sigma\succeq\frac{\sigma^{2}}{\kappa}I$, which implies $\Sigma_{0}\preceq(\sigma^{2}\poly(d))I\preceq(\kappa\poly(d))\Sigma$.

We conclude the proof by noting that $\Sigma_{0}$ can be computed by multiplying a $d\times|S|$ matrix with an $|S|\times d$ matrix. This can be done in time $\widetilde{O}(d^{3.26}/\epsilon^{2})$ by fast rectangular matrix multiplication (Lemma 2.2). ∎

Let us first prove the guarantee for the crude $O(\sqrt{\epsilon})$ additive estimation.

Under the same setting as Theorem 1.3, there exists universal constant $C_{0}$ such that Algorithm 4 outputs an estimate $\widehat{\Sigma}$ that satisfies ${\left\|\widehat{\Sigma}-\Sigma\right\|}_{F}=O(\sqrt{\epsilon})\left\|\Sigma\right\|_{2}$ with high probability in time $\widetilde{O}(d^{3.26})/\poly(\epsilon)$.

By Lemma 3.1, in Algorithm 4, we can compute $\Sigma_{0}\succeq\Sigma$ such that $\left\|\Sigma_{0}\right\|_{2}\leq\poly(d)\left\|\Sigma\right\|_{2}$. Lemma 3.2 allows us to iteratively compute $\Sigma_{t+1}$ such that $\Sigma\preceq\Sigma_{t+1}\preceq\Sigma+O(\sqrt{\epsilon})\Sigma_{t}$. It follows that $\left\|\Sigma_{t+1}\right\|_{2}\leq\left\|\Sigma\right\|_{2}+O(\sqrt{\epsilon})\left\|\Sigma_{t}\right\|_{2}$, and thus after $O(\log d)$ iterations we have a matrix $\Sigma_{T}$ with $\left\|\Sigma_{T}\right\|_{2}\leq 2\left\|\Sigma\right\|_{2}$. Using Lemma 3.2 with $\Sigma_{T}$, we can get a matrix $\widehat{\Sigma}$ with $\|\widehat{\Sigma}-\Sigma\|_{F}=O(\sqrt{\epsilon})\left\|\Sigma_{T}\right\|_{2}=O(\sqrt{\epsilon})\Sigma$. The running time follows from the running time of Lemmas 3.1 and 3.2. ∎

Suppose the constant hiding in the $O(\cdot)$ notation in Lemma B.1 is $C_{0}$, we have the following immediate corollary of Lemma B.1.

In Algorithm 3, the matrix $M_{0}$ satisfies ${\left\|M_{0}-\Sigma\right\|}_{F}\leq C_{0}\sqrt{\epsilon}\left\|\Sigma\right\|_{2}$.

In Algorithm 3 we will choose $C_{1}=20C_{r}$ and define $S_{1}$ to be the subspace where the eigenvalues of $M_{0}$ are at least $C_{1}\sqrt{\epsilon}$. We can then show the following lemma:

In Algorithm 3, with high probability, the matrix $M_{1}$ satisfies

We assume the calls to compute $M_{0},M_{1}$ are successful, which happens with high probability.

We continue to bound $\left\|\Sigma[S_{1}^{\perp}]\right\|_{2}$. Notice that by Corollary B.2,

Here the last step uses the fact when $\epsilon_{0}$ is small enough $\left\|M_{0}\right\|_{2}\leq 2\left\|\Sigma\right\|$. ∎

We will now show that $M_{2}$ computed by Algorithm 3 has low additive error in the subspace $S_{12}$.

In Algorithm 3, with high probability, the matrix $M_{2}$ satisfies

Moreover, if we consider $\Pi_{S_{12}}X_{i}$ as vectors of dimension equal to the dimension of $S_{12}$, the algorithm runs in time $\widetilde{O}(d^{3.26})/\poly(\epsilon)$.

We assume the calls for computing matrices $M_{0},M_{1},M_{2}$ are all successful, which happens with high probability.

Let $Y_{i}=\Pi_{S_{12}}X_{i}$, we know the covariance of these samples are exactly equal to $\Sigma[S_{12}]$.

In this case, by the guarantee of Theorem 1.2 we know

We can left and right multiply by $M_{2}^{1/2}$ and get

When $\epsilon$ is small enough this implies $\left\|M_{2}\right\|_{2}\leq 2\left\|\Sigma[S_{12}]\right\|_{2}\leq 2\left\|\Sigma\right\|_{2}$, therefore as desired we have

The only thing left to establish is the running time. To bound the running time we will show $\kappa(\Sigma[S_{12}])=O(1/\epsilon)$. Here we restrict the attention to the subspace $S_{12}$, so if $S_{12}$ as dimension $k$, $\kappa$ is the ratio of the largest eigenvalue and the $k$-th eigenvalue.

Let $S^{\star}$ be the subspace of eigenvectors of $\Sigma$ with eigenvalue at most $\epsilon\left\|\Sigma\right\|_{2}$. We first show the following claim:

For any unit vector $v\in S^{\star}$, $\left\|\Pi_{S_{12}}v\right\|_{2}^{2}\leq 1/5$.

Since $\left\|\Pi_{S_{12}}v\right\|_{2}^{2}=\left\|\Pi_{S_{1}}v\right\|_{2}^{2}+\left\|\Pi_{S_{2}}v\right\|_{2}^{2}$, we will bound the contributions separately. Notice that for any $v\in S^{\star}$, we have $v^{\top}\Sigma v\leq\epsilon$, therefore by Corollary B.2,

On the other hand, $v^{\top}M_{0}v\geq C_{1}\sqrt{\epsilon}\left\|\Sigma\right\|_{2}\left\|\Pi_{S_{1}}v\right\|_{2}^{2}$. Combining the two equations we get $\left\|\Pi_{S_{1}}v\right\|_{2}^{2}\leq 2C_{r}/C_{1}=1/10$. The proof for $\left\|\Pi_{S_{2}}v\right\|_{2}^{2}$ is exactly the same except we use Lemma B.3. ∎

For any two subspaces $U$ and $V$, one can check that

where $\angle(u,v)$ is the angle between $u,v$. Therefore we know for any vector $v\in S_{12}$, $\left\|\Pi_{S^{\star}}v\right\|_{2}^{2}\leq 1/5$. This implies for every $v\in S_{12}$,

This shows $\lambda_{k}(\Sigma_{S_{12}})\geq\frac{4\epsilon}{5}\left\|\Sigma\right\|_{2}$, so $\kappa\leq\frac{5}{4\epsilon}=O(1/\epsilon)$. ∎

Finally we give the guarantee for $M_{3}$.

In Algorithm 3, with high probability, the matrix $M_{3}$ satisfies

We assume the calls to compute $M_{0},M_{1},M_{3}$ are successful, which happens with high probability.

Set $Y_{i}=(\epsilon^{1/2}\Pi_{S_{1}}+\epsilon^{1/4}\Pi_{S_{2}}+\Pi_{S_{3}})X_{i}$. Let $\Sigma_{Y}$ denote the covariance of $Y$. It is easy to check that $\Sigma_{Y}\preceq 3\left(\epsilon\Sigma[S_{1}]+\sqrt{\epsilon}\Sigma[S_{2}]+\Sigma[S_{3}]\right)$. By Lemma B.2,

Combining these we have $\left\|\Sigma_{Y}\right\|_{2}\leq O(\epsilon)\left\|\Sigma\right\|_{2}$. Therefore by Lemma B.1 we know the estimation $M_{3}$ satisfies ${\left\|M_{3}-\Sigma_{Y}\right\|}_{F}\leq O(\epsilon^{1.5})\left\|\Sigma\right\|_{2}$.

On the other hand, it is easy to check that $\Sigma_{Y}[S_{1},S_{3}]=\frac{1}{\sqrt{\epsilon}}\Sigma[S_{1},S_{3}]$, therefore we have

Finally we are ready to combine all the steps.

We will assume the four calls to Algorithm 1 and Algorithm 4 are all successful, which happens with high probability. The running time follows from Lemma B.1 and Lemma B.4. Now the resulting matrix looks like

By Lemmas B.3, B.4, B.5 we know for each one of these nine blocks the error is bounded by $O(\epsilon\log(1/\epsilon))\left\|\Sigma\right\|_{2}$. Therefore, the entire matrix also has error at most $O(\epsilon\log(1/\epsilon))\left\|\Sigma\right\|_{2}$. ∎