Robust Sub-Gaussian Principal Component Analysis and Width-Independent Schatten Packing

Introduction

We study two natural, but seemingly unrelated, problems in high dimensional robust statistics and continuous optimization respectively. As we will see, these problems have an intimate connection.

Problem 1: Robust sub-Gaussian principal component analysis. We consider the following statistical task, which we call robust sub-Gaussian principal component analysis (PCA). Given samples $X_{1},\ldots,X_{n}$ from sub-Gaussian See Section 2 for a formal definition. distribution $\mathcal{D}$ with covariance $\bm{\Sigma}$, an $\epsilon$ fraction of which are arbitrarily corrupted, the task asks to output unit vector $u$ with $u^{\top}\bm{\Sigma}u\geq(1-\gamma)\left\lVert\bm{\Sigma}\right\rVert_{\infty}$ Throughout we use $\left\lVert\mathbf{M}\right\rVert_{p}$ to denote the Schatten $p$-norm (cf. Section 2 for more details). for tolerance $\gamma$. Ergo, the goal is to robustly return a $(1-\gamma)$-approximate top eigenvector of the covariance of sub-Gaussian $\mathcal{D}$. This is the natural extension of PCA to the robust statistics setting.

There has been a flurry of recent work on efficient algorithms for robust statistical tasks, e.g. covariance estimation and PCA. From an information-theoretic perspective, sub-Gaussian concentration suffices for robust covariance estimation. Nonetheless, to date all polynomial-time algorithms achieving nontrivial guarantees on covariance estimation of a sub-Gaussian distribution (including PCA specifically) in the presence of adversarial noise require additional algebraic structure. For instance, sum-of-squares certifiably bounded moments have been leveraged in polynomial time covariance estimation [HL18, KSS18]; however, this is a stronger assumption than sub-Gaussianity.

In many applications (see discussion in [DKK+17]), the end goal of covariance estimation is PCA. Thus, a natural question which relaxes robust covariance estimation is: can we robustly estimate the top eigenvector of the covariance $\bm{\Sigma}$, assuming only sub-Gaussian concentration? Our work answers this question affirmatively via two incomparable algorithms. The first achieves $\gamma=O(\epsilon\log\epsilon^{-1})$ in polynomial time; the second achieves $\gamma=O(\sqrt{\epsilon\log\epsilon^{-1}\log d})$ in nearly-linear time under a mild gap assumption on $\bm{\Sigma}$. Moreover, both methods have nearly-optimal sample complexity.

Problem 2: Width-independent Schatten packing. We consider a natural generalization of packing semidefinite programs (SDPs) which we call Schatten packing. Given symmetric positive semidefinite $\mathbf{A}_{1},\ldots,\mathbf{A}_{n}$ and parameter $p\geq 1$, a Schatten packing SDP asks to solve the optimization problem

Here, $\left\lVert\mathbf{M}\right\rVert_{p}$ is the Schatten-$p$ norm of matrix $\mathbf{M}$ and $\Delta^{n}$ is the probability simplex (see Section 2). When $p=\infty$, (1) is the well-studied (standard) packing SDP objective [JY11, ALO16, PTZ16], which asks to find the most spectrally bounded convex combination of packing matrices. For smaller $p$, the objective encourages combinations more (spectrally) uniformly distributed over directions.

Learning with adversarial outliers. The study of estimators robust to a small fraction of adversarial outliers dates back to foundational work, e.g. [Hub64, Tuk75]. Following more recent work [LRV16, DKK+19], there has been significant interest in efficient, robust algorithms for statistical tasks in high-dimensional settings. We focus on methods robustly estimating covariance properties here, and defer a thorough discussion of the (extensive) robust statistics literature to [Ste18, Li18, DK19].

There has been quite a bit of work in understanding and giving guarantees for robust covariance estimation where the uncorrupted distribution is exactly Gaussian [DKK+17, DKK+18, DKK+19, CDGW19]. These algorithms strongly use relationships between higher-order moments of Gaussian distributions via Isserlis’ theorem. Departing from the Gaussian setting, work of [LRV16] showed that if the distribution is an affine transformation of a 4-wise independent distribution, robust covariance estimation is possible. This was extended by [KSS18], which also assumed nontrivial structure in the moments of the distribution, namely that sub-Gaussianity was certifiable via the sum-of-squares proof system. To the best of our knowledge it has remained open to give nontrivial guarantees for robust estimation of any covariance properties under minimal assumptions, i.e. sub-Gaussian concentration.

All aforementioned algorithms also yield guarantees for robust PCA, by applying a top eigenvector method to the learned covariance. However, performing robust PCA via the intermediate covariance estimation step is lossy, both statistically and computationally. From a statistical perspective, $\Omega(d^{2})$ samples are necessary to learn the covariance of a $d$-dimensional Gaussian in Frobenius norm (and for known efficient algorithms for spectral norm error [DKS17]); in contrast, $O(d)$ samples suffice for (non-robust) PCA. Computationally, even when the underlying distrubution is exactly Gaussian, the best-known covariance estimation algorithms run in time $\Omega(d^{3.25})$; algorithms working in more general settings based on the sum-of-squares approach require much more time. In contrast, the power method for PCA in a $d\times d$ matrix takes time $\widetilde{O}(d^{2})$ We say $g=\widetilde{O}(f)$ if $g=O(f\log^{c}f)$ for some constant $c>0$.. Motivated by this, our work initiates the direct study of robust PCA, which is often independently interesting in applications.

We remark there is another problem termed “robust PCA” in the literature, e.g. [CLMW11], under a different generative model. We defer a detailed discussion to [DKK+17], which experimentally shows that algorithms from that line of work do not transfer well to our corruption model.

Width-independent iterative methods. Semidefinite programming (SDP) and its linear programming specialization are fundamental computational tasks, with myriad applications in learning, operations research, and computer science. Though general-purpose polynomial time algorithms exist for SDPs ([NN94]), in practical settings in high dimensions, approximations depending linearly on input size and polynomially on error $\epsilon$ are sometimes desirable. To this end, approximation algorithms based on entropic mirror descent have been intensely studied [WK06, AK16, GHM15, AL17, CDST19], obtaining $\epsilon$ additive approximations to the objective with runtimes depending polynomially on $\rho/\epsilon$, where $\rho$ is the “width”, the largest spectral norm of a constraint.

2 Our results

Robust sub-Gaussian principal component analysis. We give two algorithms for robust sub-Gaussian PCA We follow the distribution and corruption model described in Assumption 1.. Both are near-sample optimal, polynomial-time, and assume only sub-Gaussianity. The first is via a simple filtering approach, as summarized here (and developed in Section 3).

Under Assumption 1, let $\delta\in$, and $n=\Omega\left(\tfrac{d+\log\delta^{-1}}{(\epsilon\log\epsilon^{-1})^{2}}\right)$. Algorithm 6 runs in time $O(\tfrac{nd^{2}}{\epsilon}\log\tfrac{n}{\delta\epsilon}\log\tfrac{n}{\delta})$, and outputs $u$ with $u^{\top}\bm{\Sigma}u>(1-C^{\star}\epsilon\log\epsilon^{-1})\|\bm{\Sigma}\|_{\infty}$, for $C^{\star}$ a fixed multiple of parameter $c$ in Assumption 1, with probability at least $1-\delta$.

Our second algorithm is more efficient under mild conditions, but yields a worse approximation $1-\gamma$ for $\gamma=O(\sqrt{\epsilon\log\epsilon^{-1}\log d})$. Specifically, if there are few eigenvalues of $\bm{\Sigma}$ larger than $1-\gamma$, our algorithm runs in nearly-linear time. Note that if there are many eigenvalues above this threshold, then the PCA problem itself is not very well-posed; our algorithm is very efficient in the interesting setting where the approximate top eigenvector is identifiable. We state our main algorithmic guarantee here, and defer details to Section 5.

We remark that $\Omega(d\epsilon^{-2})$ samples are necessary for a $(1-\epsilon)$-approximation to the top eigenvector of $\bm{\Sigma}$ via uncorrupted samples from $\mathcal{N}(0,\bm{\Sigma})$, so our first method is sample-optimal, as is our second up to a $\widetilde{O}(\epsilon^{-1})$ factor.

Width-independent Schatten packing. Our second method crucially requires an efficient solver for Schatten packing SDPs. We demonstrate that Schatten packing, i.e. (1) for arbitrary $p$, admits width-independent solvers. We state an informal guarantee, and defer details to Section 4.

Preliminaries

For $\mathbf{M}=\sum_{j\in[d]}\lambda_{i}v_{i}v_{i}^{\top}$, the satisfying $\mathbf{N}$ is $\tfrac{\sum_{j\in[d]}\pm\lambda_{i}^{p-1}v_{i}v_{i}^{\top}}{\left\lVert\mathbf{M}\right\rVert_{p}^{p-1}}$, so $\mathbf{N}\mathbf{M}$ has spectrum $\tfrac{|\lambda|^{p}}{\left\lVert\mathbf{M}\right\rVert_{p}^{p-1}}$.

Multivariate $\mathcal{D}$ has sub-Gaussian proxy $\bm{\Gamma}$ if its restriction to any unit $v$ is $\left\lVert v\right\rVert_{\bm{\Gamma}}^{2}$-sub-Gaussian, i.e.

We consider the following standard model for gross corruption with respect to distribution a $\mathcal{D}$.

As we only estimate covariance properties, the assumption that $\mathcal{D}$ is mean-zero only loses constants in problem parameters, by pairing samples and subtracting them (cf. [DKK+19], Section 4.5.1).

Robust sub-Gaussian PCA via filtering

In this section, we sketch the proof of Theorem 1, which gives guarantees on our filtering algorithm for robust sub-Gaussian PCA. This algorithm obtains stronger statistical guarantees than Theorem 2, at the cost of super-linear runtime; the algorithm is given as Algorithm 6. Our analysis stems largely from concentration facts about sub-Gaussian distributions, as well as the following (folklore) fact regarding estimation of variance along any particular direction.

In other words, we show that using corrupted samples, we can efficiently estimate a $1+O(\epsilon\log\epsilon^{-1})$-multiplicative approximation of the variance of $\mathcal{D}$ in any unit direction Corollary 5 gives a slightly stronger guarantee that reusing samples does not break dependencies of $u$.. This proof is deferred to Appendix B for completeness. Algorithm 6 combines this key insight with a soft filtering approach, suggested by the following known structural fact found in previous work (e.g. Lemma A.1 of [DHL19], see also [SCV17, Ste18]).

Let $\{a_{i}\}_{i\in[m]}$, $\{w_{i}\}_{i\in[m]}$ be sets of nonnegative reals, and $a_{\max}=\max_{i\in[m]}a_{i}$. Define $w^{\prime}_{i}=\left(1-\frac{a_{i}}{a_{\max}}\right)w_{i}$, for all $i\in[m]$. Consider any disjoint partition $I_{B}$, $I_{G}$ of $[m]$ with $\sum_{i\in I_{B}}w_{i}a_{i}>\sum_{i\in I_{G}}w_{i}a_{i}.$ Then, $\sum_{i\in I_{B}}w_{i}-w_{i}^{\prime}>\tfrac{1}{2a_{\max}}\sum_{i\in[m]}w_{i}a_{i}>\sum_{i\in I_{G}}w_{i}-w_{i}^{\prime}$.

Our Algorithm 6, $\mathsf{PCAFilter}$, takes as input a set of corrupted samples $\{X_{i}\}_{i\in[n]}$ following Assumption 1 and the corruption parameter $\epsilon$. At a high level, it initializes a uniform weight vector $w^{(0)}$, and iteratively operates as follows (we denote by $\mathbf{M}(w)$ the empirical covariance $\sum_{i\in[n]}w_{i}X_{i}X_{i}^{\top}$).

$u_{t}\leftarrow$ approximate top eigenvector of $\mathbf{M}(w^{(t-1)})$ via power iteration.

Compute $\sigma_{t}^{2}\leftarrow\mathsf{1DRobustVariance}(\{X_{i}\}_{i\in[n]},u_{t},\epsilon)$.

If $\sigma_{t}^{2}>(1-O(\epsilon\log\epsilon^{-1}))\cdot u_{t}^{\top}\mathbf{M}(w^{(t-1)})u_{t}$, then terminate and return $u_{t}$.

Sort indices $i\in[n]$ by $a_{i}\leftarrow\left\langle u_{t},X_{i}\right\rangle^{2}$, with $a_{1}$ smallest.

The analysis of Algorithm 6 then proceeds in two stages.

As Lemma 2 then applies, the procedure always removes more mass from bad points than good, and thus can only remove at most $2\epsilon$ mass total by the corruption model. Thus, the weights $w^{(t)}$ are always roughly uniform (in $\mathfrak{S}_{O(\epsilon)}^{n}$), which by standard concentration facts (see Appendix A) imply the quality of the approximate top eigenvector is good. Moreover, the iteration count is bounded by roughly $d$ because whenever the algorithm does not terminate, enough mass is removed from large spectral directions. Combining with the termination criteria imply that when a vector is returned, it is a close approximation to the top direction of $\bm{\Sigma}$. Details can be found as Lemma 15 and in the proof of Theorem 1.

Schatten packing

The following result is shown in [MRWZ16].

$\mathsf{PackingLP}$ (Algorithm 1) solves Problem 1 in $O(\textup{nnz}(\mathbf{A})\cdot\tfrac{\log(d)\log(nd/\epsilon)}{\epsilon^{2}})$ time.

Oour interpretation of the analysis of [MRWZ16], combines two ingredients: a potential argument and mirror descent, which yields a dual feasible point if $\left\lVert w_{t}\right\rVert_{1}$ did not grow sufficiently.

Potential argument. The potential used by [MRWZ16] is $\log(\sum_{j\in[d]}\exp([\mathbf{A}w_{t}]_{j}))-\left\lVert w_{t}\right\rVert_{1}$, well-known to be a $O(\log d)$-additive approximation of $\left\lVert\mathbf{A}w_{t}\right\rVert_{\infty}-\left\lVert w_{t}\right\rVert_{1}$. As soon as $\left\lVert\mathbf{A}w_{t}\right\rVert_{\infty}$ or $\left\lVert w_{t}\right\rVert_{1}$ reaches the scale $O(\tfrac{\log d}{\epsilon})$, by nonnegativity this becomes a multiplicative guarantee, motivating the setting of threshold $K$. To prove the potential is monotone, [MRWZ16] uses step size $K^{-1}$ and a Taylor approximation; combining with the termination condition yields the desired claim.

Mirror descent. To certify that $w_{t}$ grows sufficiently (e.g. the method terminates in few iterations, else dual feasibility holds), we interpret the step $w_{t+1}\leftarrow w_{t}\circ(1+\eta g_{t})$ as approximate entropic mirror descent. Specifically, we track the quantity $\sum_{0\leq t<T}\left\langle\eta g_{t},u\right\rangle$, and show that if $\left\lVert w_{t}\right\rVert_{1}$ has not grown sufficiently, then it must be bounded for every $u\in\Delta^{n}$, certifying dual feasibility. Formally, for any $g_{t}$ sequence and $u\in\Delta^{n}$, we show

The last inequality followed by $g_{t}$ being an upwards truncation. If $\left\lVert w_{T}\right\rVert_{1}$ is bounded (else, we have primal feasibility), we show the entire above expression is bounded $O(\log\tfrac{nd}{\epsilon})$ for any $u$. Thus, by setting $T=O(\tfrac{\log(nd/\epsilon)}{\eta\epsilon})$ and choosing $u$ to be each coordinate indicator, it follows that the average of all $v_{t}$ is coordinatewise at least $1-\epsilon$, and solves Problem 1 as a dual solution.

In all iterations $t$ of Algorithm 2, defining $\Phi_{t}:=\left\lVert\mathbf{A}w_{t}\right\rVert_{p}-\left\lVert w_{t}\right\rVert_{1}$, $\Phi_{t+1}\leq\Phi_{t}$.

We now prove our main result, which leverages the potential bound following the framework of Section 4.1. In the proof, we assume that entries of $\mathbf{A}$ are bounded by $n\epsilon^{-1}$; this does not incur more loss than a constant multiple of $\epsilon$ in the guarantees, and a proof can be found as Lemma 16.

Algorithm 2 runs in time $O(\textup{nnz}(\mathbf{A})\cdot\tfrac{p\log(nd/\epsilon)}{\epsilon})$. Further, its output solves Problem 2.

The runtime follows from Line 7 (each iteration cost is dominated by multiplication through $\mathbf{A}$), so we prove correctness. Define potential $\Phi_{t}$ as in Lemma 3, and note that as $w_{0}=\tfrac{\epsilon}{n^{2}d}\mathbf{1}$,

The second inequality followed from our assumption on $\mathbf{A}$ entry sizes (Lemma 16). If Algorithm 2 breaks out of the while loop of Line 4, we have by Lemma 3 that for $x$ returned on Line 11,

Thus, primal feasibility is always correct. We now prove correctness of dual feasibility. First, let $V_{x}(u)=\sum_{i\in[n]}u_{i}\log(\tfrac{u_{i}}{x_{i}})$ be the Kullback-Leibler divergence from $x$ to $u$, for $x$, $u\in\Delta^{d}$. Define the normalized points $x_{t}=\tfrac{w_{t}}{\left\lVert w_{t}\right\rVert_{1}}$ in each iteration. Expanding definitions,

The only inequality used the bounds, for $g,\eta\in$,

Telescoping (4) over all $T$ iterations, and using $V_{x_{0}}(u)\leq\log n$ for all $u\in\Delta^{n}$ since $x_{0}$ is uniform, we have that whenever Line 4 is not satisfied before the check on Line 7 (i.e. $t\geq T$),

The last inequality used $\left\lVert w_{T}\right\rVert_{1}\leq\epsilon^{-1}$ by assumption, and $\left\lVert w_{0}\right\rVert_{1}=\tfrac{\epsilon}{nd}$. Next, since each $g_{t}\geq\mathbf{1}-\mathbf{A}^{\top}(v_{t})^{p-1}$ entrywise, defining $\bar{z}=\tfrac{z}{T}$,

Combining (5) and (6), and rearranging, yields by definition of $T$,

3 Schatten-norm packing semidefinite programs

We generalize Algorithm 2 to solve Schatten packing semidefinite programs, which we now define.

We assume that $p$ is an odd integer for simplicity (sufficient for our applications), and leave for interesting future work the cases when $p$ is even or noninteger. The potential used in the analysis and an overall guarantee are stated here, and deferred to Appendix C. The proofs are simple modifications of Lemma 3 and Theorem 4 using trace inequalities (similar to those in [JLL+20]) in place of scalar inequalities, as well as efficient approximation of quantities in Line 5 via the standard technique of Johnson-Lindestrauss projections.

In all iterations $t$ of Algorithm 3, defining $\Phi_{t}:=\left\lVert\sum_{i\in[n]}[w_{t}]_{i}\mathbf{A}_{i}\right\rVert_{p}-\left\lVert w_{t}\right\rVert_{1}$, $\Phi_{t+1}\leq\Phi_{t}$.

Let $p$ be odd. Algorithm 3 runs in $O(\tfrac{p\log(nd/\epsilon)}{\epsilon})$ iterations, and its output solves Problem 3. Each iteration is implementable in $O(\textup{nnz}\cdot\tfrac{p\log(nd/\epsilon)}{\epsilon^{2}})$, where nnz is the number of nonzero entries amongst all $\{\mathbf{A}_{i}\}_{i\in[n]}$, losing $O(\epsilon)$ in the quality of Problem 3 with probability $1-\textup{poly}((nd/\epsilon)^{-1})$.

We remark that the framework outlined in Section 4.1 is flexible enough to handle mixed-norm packing problems. Specifically, developments in Section 5 require the following guarantee.

for $\mathcal{A}(x):=\sum_{i\in[n]}x_{i}\mathbf{A}_{i}$. Given estimate of OPT exponentially bounded in $\tfrac{nd}{\epsilon}$, there is a procedure calling Algorithm 7 $O(\log\frac{nd}{\epsilon})$ times giving $x\in\Delta^{n}$ with $\left\lVert x\right\rVert_{\infty}\leq\tfrac{(1+\alpha)(1+\epsilon)}{n}$, $\left\lVert\mathcal{A}(x)\right\rVert_{p}\leq(1+\epsilon)\textup{OPT}$. Algorithm 7 runs in $O(\tfrac{\log(nd/\epsilon)\log n}{\epsilon^{2}})$ iterations, each implementable in time $O(\textup{nnz}\cdot\frac{p\log(nd/\epsilon)}{\epsilon^{2}})$.

Our method, found in Appendix C, approximately solves (7) by first applying a standard binary search to place $\mathcal{A}(x)$ on the right scale, for which it suffices to solve an approximate decision problem. Then, we apply a truncated mirror descent procedure on the potential $\Phi(w)=\log(\exp(\left\lVert\mathcal{A}(w)\right\rVert_{p})+\exp(\tfrac{n}{1+\alpha}\left\lVert w\right\rVert_{\infty}))-\left\lVert w\right\rVert_{1}$, and prove correctness for solving the decision problem following the framework we outlined in Section 4.1.

Robust sub-Gaussian PCA in nearly-linear time

We give our nearly-linear time robust PCA method, leveraging developments of Section 4. Throughout, we will be operating under Assumption 1, for some corruption parameter $\epsilon$ with $\epsilon\log\epsilon^{-1}\log d=O(1)$; $\epsilon=O(\frac{1}{\log d\log\log d})$ suffices. We now develop tools to prove Theorem 2.

Algorithm 4 uses three subroutines: our earlier $\mathsf{1DRobustVariance}$ method (Lemma 1), an application of our earlier Proposition 2 to approximate the solution to

and a method for computing approximate eigenvectors by [MM15] (discussed in Appendix D).

Here, $\lambda_{j}(\mathbf{A})$ is the $j^{th}$ largest eigenvalue of $\mathbf{A}$. The total time required by the method is $O(\textup{nnz}(\mathbf{A})\tfrac{tp\log d}{\varepsilon})$.

Algorithm 4 is computationally bottlenecked by the application of Proposition 2 on Line 2 and the call to $\mathsf{Power}$ on Line 4, from which the runtime guarantee of Theorem 2 follows straightforwardly. To demonstrate correctness, we first certify the quality of the solution to (8).

The proof of this is similar to results in e.g. [DKK+19, Li18], and combines concentration guarantees with a union bound over all possible corruption sets $B$. This implies the following immediately, upon applying the guarantees of Proposition 2.

Let $w$ be the output of the solver. Recall that $\mathbf{M}=\sum_{i=1}^{n}w_{i}X_{i}X_{i}^{\top}$. Additionally, define

Notice in particular that $\mathbf{M}=\mathbf{M}_{G}+\mathbf{M}_{B}$, and that all these matrices are PSD. We next prove the second, crucial fact, which says that $\mathbf{M}_{G}$ is a good approximator to $\bm{\Sigma}$ in Loewner ordering:

The proof combines the strategy in Lemma 5 with the guarantee of the SDP solver. Perhaps surprisingly, Corollary 1 and Lemma 6 are the only two properties about $\mathbf{M}$ that our final analysis of Theorem 2 will need. In particular, we have the following key geometric proposition, which carefully combines trace inequalities to argue that the corrupted points $\mathbf{M}_{B}$ cannot create too many new large eigendirections.

be sorted eigendecompositions of $\mathbf{M}$ and $\bm{\Sigma}$, so $\lambda_{1}\geq\ldots\geq\lambda_{d}$, and $\sigma_{1}\geq\ldots\geq\sigma_{d}$. Let $\gamma$ be as in Theorem 2, and assume $\sigma_{t+1}<(1-\gamma)\sigma_{1}$. Then,

For concreteness, we will define the parameters

Suppose for contradiction that all $v_{j}^{\top}\bm{\Sigma}v_{j}<(1-\gamma)\sigma_{1}$ for $j\in[t]$. By applying the guarantee of Corollary 1 and Fact 2, it follows that

Let $s\in[d]$ be the largest index such that $\sigma_{s}>\left(1-\tfrac{\gamma}{4}\right)\sigma_{1}$, and note that $s\leq t$. We define

That is, $\mathbf{N}$ is the restriction of $\mathbf{M}^{p-1}$ to its top $s$ eigendirections. Then,

In the proof of Proposition 4, we used the following facts.

Let $\mathbf{A}\succeq\mathbf{B}\succeq 0$ be symmetric matrices and $p$ a positive integer. Then we have

By the Courant-Fischer minimax characterization of eigenvalues,

The guarantees of Proposition 4 were geared towards exact eigenvectors of the matrix $\mathbf{M}$. We now modify the analysis to tolerate inexactness in the eigenvector computation, in line with the processing of Line 5 of our Algorithm 4. This yields our final claim in Theorem 2.

In the setting of Proposition 4, and letting $\{z_{j}\}_{j\in[t]}$ satisfy (9), set for all $j\in[t]$

Then with probability at least $1-\delta$,

Assume all $y_{j}$ have $y_{j}^{\top}\bm{\Sigma}y_{j}\leq(1-\gamma)\sigma_{1}$ for contradiction. We outline modifications to the proof of Proposition 4. Specifically, we redefine the matrix $\mathbf{N}$ by

Because $\sum_{j\in[s]}z_{j}z_{j}^{\top}$ is a projection matrix, it is clear $\mathbf{N}\preceq\mathbf{M}^{p-1}$. Therefore, by combining the derivations (13) and (14), it remains true that

We now bound these two terms in an analogous way from Proposition 4, with negligible loss; combining these bounds will again yield a contradiction. First, we have the lower bound

Here, the last inequality applied the assumption (9) with respect to $\mathbf{M}^{p}$. Next, we upper bound

Finally, we prove Theorem 2 by combining the tools developed thus far.

Correctness of the algorithm is immediate from Corollary 2 and the guarantees of $\mathsf{1DRobustVariance}$. Concretely, Corollary 2 guarantees that one of the vectors we produce will be a $(1-\gamma)$-approximate top eigenvector (say some index $j\in[t]$), and $\mathsf{1DRobustVariance}$ will only lose a negligible fraction $O(\epsilon\log\epsilon^{-1})$ of this quality (see Lemma 1); the best returned eigenvector as measured by $\mathsf{1DRobustVariance}$ can only improve the guarantee. Finally, the failure probability follows by combining the guarantees of Lemmas 1, 5, and 6.

We now discuss runtime. The complexity of lines 2, 4, and 5, as guaranteed by Proposition 2, Proposition 3, and Lemma 1 are respectively (recalling $p=\widetilde{O}(\epsilon^{-0.5})$)

Throughout we use that we can compute matrix-vector products in an arbitrary linear combination of the $X_{i}X_{i}^{\top}$ in time $O(nd)$; it is easy to check that in all runtime guarantees, nnz can be replaced by this computational cost. Combining these bounds yields the final conclusion. ∎

We thank Swati Padmanabhan and Aaron Sidford for helpful discussions.

References

Appendix A Concentration

We use the following concentration facts on sub-Gaussian distributions following from standard techniques, and give an application bounding Schatten-norm deviations.

Under Assumption 1, there are universal constants $C_{1}$, $C_{2}$ such that

By observing (3), it is clear that the random vector $\widetilde{X}=\bm{\Sigma}^{-\frac{1}{2}}X$ for $X\sim\mathcal{D}$ has covariance $\mathbf{I}$ and sub-Gaussian proxy $c\mathbf{I}$. For any fixed unit vector $u$, by Lemma 1.12 of [RH17], the random variable $(u^{\top}\widetilde{X})^{2}-1$ is sub-exponential with parameter $16c$, so by Bernstein’s inequality (Theorem 1.13, [RH17]), defining $\widetilde{X}_{i}=\bm{\Sigma}^{-\frac{1}{2}}X_{i}$ for each $X_{i}\sim\mathcal{D}$,

Next, by a standard application of the triangle inequality (see e.g. Exercise 4.3.3, [Ver16])

with probability at least $1-\exp\left(C_{1}d-C_{2}n\min(t,t^{2})\right)$ for appropriate $C_{1}$, $C_{2}$. The conclusion follows since its statement is scale invariant, so it suffices to show as we have that

Let $p\geq 2$. Under Assumption 1, there are universal constants $C_{1}$, $C_{2}$ with

Suppose the event in Lemma 9 does not hold, which happens with probability at least $1-\exp(C_{1}d-C_{2}n\min(t,t^{2}))$. Define for shorthand $\mathbf{M}:=\tfrac{1}{n}\sum_{i\in G^{\prime}}X_{i}X_{i}^{\top}-\bm{\Sigma}$ and let its spectral decomposition be $\sum_{j\in[d]}\lambda_{j}v_{j}v_{j}^{\top}$. By the triangle inequality and Fact 2,

A.2 Concentration under weightings in 𝔖ϵn\mathfrak{S}_{\epsilon}^{n}

We consider concentration of the empirical covariance under weightings which are not far from uniform, in spectral and Schatten senses.

Under Assumption 1, let $\delta\in$, $p\geq 2$, and $n=\Omega\left(\tfrac{d+\log\delta^{-1}}{(\epsilon\log\epsilon^{-1})^{2}}\right)$ for a sufficiently large constant. Then for a universal constant $C_{3}$,

Because the vertices of $\mathfrak{S}_{\epsilon}^{n}$ are uniform over sets $S\subseteq G^{\prime}$ with $|S|=(1-\epsilon)n$ (see e.g. Section 4.1, [DKK+19]), by convexity of the Schatten-$p$ norm it suffices to prove

For any fixed $S$, and recalling $|S^{c}|=\epsilon n$, we can decompose this sum as

By applying Corollary 3, it follows that by setting $t=\tfrac{1-\epsilon}{2}\cdot\epsilon\log\epsilon^{-1}$ and our choice of $n$ that

Moreover, for any fixed $S^{c}$, setting $t=\tfrac{1-\epsilon}{2}\cdot C_{3}\log\epsilon^{-1}$ where $C_{3}$ is a sufficiently large constant, so that for sufficiently small $\epsilon$, $t=\min(t,t^{2})$,

Here, we used that $\log\binom{n}{\epsilon n}=O\left(n\epsilon\log\epsilon^{-1}\right)$. Finally, union bounding over all possible sets $S^{c}$ imply that with probability at least $1-\tfrac{\delta}{2}$, the following events hold:

Combining these bounds in the context of (18) after applying the triangle inequality, we have with probability at least $1-\tfrac{\delta}{2}$ for all $S$ the desired conclusion,

Under Assumption 1, let $n=\Omega\left(\tfrac{d+\log\delta^{-1}}{(\epsilon\log\epsilon^{-1})^{2}}\right)$ for a sufficiently large constant. For universal $C_{3}$ and all $w\in\mathfrak{S}_{\epsilon}^{n}$, with probability at least $1-\tfrac{\delta}{2}$,

Therefore, again using the formula (18) and the triangle inequality yields the desired conclusion for all directions $v$, which is equivalent to the spectral bound of the lemma statement. ∎

Appendix B Deferred proofs from Section 3

In this section, we prove Lemma 1, which allows us to robustly estimate the quadratic form of a vector in the covariance of a sub-Gaussian distribution from corrupted samples. Algorithm 5 is folklore, and intuitively very simple; it projects all samples onto $u$, throws away the $2\epsilon$ fraction of points with largest magnitude in this direction, and takes the mean of the remaining set.

We have by Hölder’s inequality that for any $p,q\geq 1$ with $p^{-1}+q^{-1}=1$,

The second inequality is Lemma 1.10 [RH17]. Setting $p=\log\epsilon^{-1}$ yields the result. ∎

The runtime claim is immediate; we now turn our attention to correctness. We follow notation of Assumption 1, and in a slight abuse of notation, also define $a_{i}=\left\langle X_{i},u\right\rangle^{2}$ for $i\in G^{\prime}$. First, for $X\sim\mathcal{D}$, then $\left\langle u,X\right\rangle^{2}-u^{\top}\bm{\Sigma}u$ is sub-exponential with parameter at most $16cu^{\top}\bm{\Sigma}u$ (Lemma 1.12, [RH17]). By Bernstein’s inequality, we have that if $X\sim\mathcal{D}$, then for all $t\geq 1$,

Using this in a standard Chernoff bound, we have that with probability $1-\tfrac{\delta}{2}$,

Define the interval $I=[0,T]$ and let $S$ be the set of points in $[n]$ that survive the truncation procedure, so that $\sigma^{2}_{u}=\tfrac{1}{|S|}\sum_{i\in S}a_{i}$. Given event (24), $a_{i}\in I$ for all $i\in S$, since there are at most $\epsilon n$ points in $G$ outside $I$, and $|B|\leq\epsilon n$. We decompose the deviation as follows:

Here we overloaded $i\in I$ to mean that $a_{i}$ lies in the interval $I$, and conditioned on $S$ lying entirely in $I$. We bound each of these terms individually. First, for all $i\in G^{\prime}\cap I$, conditioning on (24) (i.e. all $a_{i}\in I$), $a_{i}-u^{\top}\bm{\Sigma}u$ is an independent sample from $Y$. Thus, by (25) and Bernstein’s inequality,

with (conditional) probability at least $1-\tfrac{\delta}{2}$. By a union bound, both events occur with probability at least $1-\delta$; condition on this for the remainder of the proof. Under this assumption, we control the other three terms of (26). Observe that $|B\cap S|\leq\epsilon n$, $|(G^{\prime}\setminus G)\cap I|\leq\epsilon n$, and $|(G\cap I)\setminus S|\leq\epsilon n$. Further, by definition of $I$, every summand is at most $64c\log\epsilon^{-1}\cdot u^{\top}\bm{\Sigma}u$. Thus,

Combining (27), (28), (29), and (30) in derivation (26) and dividing by $|S|$ yields the claim. ∎

Finally, we also give an alternative set of conditions under which we can certify correctness of $\mathsf{1DRobustVariance}$. Specifically, this assumption will be useful in lifting indpendence assumptions between $u$ and our samples $\{X_{i}\}_{i\in[n]}$ in repeated calls within Algorithm 6.

Under Assumption 1, let the following conditions hold for universal constant $C_{4}$:

Note that (32) is a factor $\epsilon$ weaker in its guarantee than Corollary 4, and is over weights in a different set $\mathfrak{S}_{1-\epsilon}^{n}$. Standard sub-Gaussian concentration (i.e. an unweighted variant of Corollary 4) and modifying the proof of Corollary 4 to take the constraint set $\mathfrak{S}_{1-\epsilon}^{n}$ and normalizing over vertex sets of size $\epsilon n$ yield the following conclusion.

Let $n=\Omega\left(\tfrac{d+\log\delta^{-1}}{(\epsilon\log\epsilon^{-1})^{2}}\right)$ for a sufficiently large constant. Assumption 2 holds with probability at least $1-\tfrac{\delta}{2}$.

We give a variant of Lemma 1 with slightly stronger guarantees for $\mathsf{1DRobustVariance}$; specifically, it holds for all $u$ simultaneously for a fixed set of samples satisfying Assumption 2.

Under Assumption 2, Algorithm 5 outputs $\sigma_{u}^{2}$ with $|u^{\top}\bm{\Sigma}u-\sigma_{u}^{2}|<Cu^{\top}\bm{\Sigma}u\cdot\epsilon\log\epsilon^{-1}$, for $C$ a fixed multiple of the parameter $c$ in Assumption 1, and runs in time $O(nd+n\log n)$.

We discuss how to modify the derivations from Lemma 1 appropriately in the absence of applications of Bernstein’s inequality. First, note that appropriately combining (31) and (32) in a derivation such as (18) yields the following bound (deterministically under Assumption 2):

Now, consider the decomposition (26). We claim first that similarly to (28), (29), (30) we can bound each summand in the latter three terms by $O(u^{\top}\bm{\Sigma}u\log\epsilon^{-1})$; to prove this, it suffices to show that at least one filtered $a_{i}$ attains this bound, as then by definition of the algorithm, each non-filtered $a_{i}$ will as well. Note that a fraction between $\epsilon$ and $2\epsilon$ of points in $G\subset G^{\prime}$ is filtered (since there are only $\epsilon n$ points from $B$). The assumption (32) then implies precisely the desired bound on some filtered $a_{i}$ by placing uniform mass on filtered points from $G$, and applying pigeonhole. So, all non-filtered $a_{i}$ are bounded by $O(u^{\top}\bm{\Sigma}u\log\epsilon^{-1})$, yielding analogous statements to (28), (29), (30).

Finally, an analogous derivation to (27) follows via an application of the bound (33), where we place uniform mass on the set $G^{\prime}\cap I$ and adjust constants appropriately, since the above argument shows that under the assumption (32), we have that at most $2\epsilon n$ indices $i\in G^{\prime}$ have $a_{i}\not\in I$. ∎

B.2 Preliminaries

For convenience, we give the following preliminaries before embarking on our proof of Theorem 1 and giving guarantees on Algorithm 6. First, we state a set of assumptions which augments Assumption 2 with one additional condition, used in bounding the iteration count of our algorithm.

Under Assumption 1, let Assumption 2 hold, as well as the following additional condition for the same universal constant $C_{4}$:

Standard sub-Gaussian concentration inequalities and a union bound, combined with our earlier claim Lemma 11, then yield the following guarantee.

Let $n=\Omega\left(\tfrac{d+\log\delta^{-1}}{(\epsilon\log\epsilon^{-1})^{2}}\right)$ for a sufficiently large constant. Assumption 3 holds with probability at least $1-\delta$.

B.3 Analysis of 𝖯𝖢𝖠𝖥𝗂𝗅𝗍𝖾𝗋\mathsf{PCAFilter}

For this section, for any nonnegative weights $w$, define $\mathbf{M}(w):=\sum_{i\in[n]}w_{i}X_{i}X_{i}^{\top}$. We now state our algorithm, $\mathsf{PCAFilter}$. At all iterations $t$, it maintains a current nonnegative weight vector $w^{(t)}$ (initialized to be the uniform distribution on $[n]$), preserving the following invariants for all $t$:

We now state our method as Algorithm 6; note that the update to $w^{(t)}$ is of the form in Lemma 2.

Under Assumption 2, for any iteration $t$ of Algorithm 6, suppose (35) held for all iterations $t^{\prime}\leq t-1$. Then, (35) holds at iteration $t$.

On the other hand, by (33) we know that the total quadratic form over $G$ is bounded as

Here, we applied the observation that the normalized $w_{i}$ restricted to $G$ are in $\mathfrak{S}^{n}_{1-3\epsilon}$ (e.g. using Lemma 14 inductively). However, since we did not terminate (Line 5), we must have by $u_{t}$ being a top eigenvector and Corollary 5 (we defer discussions of inexactness to Theorem 1) that

To obtain the last conclusion, we used (38). Finally, note that for all $i\in B\setminus I_{B}$,

Thus, the desired inequality (36) follows from combining the above derivations, e.g. using (37) and

Lemma 13 yields for all $t$ that $\sum_{i\in B}w^{(0)}_{i}-w^{(t)}_{i}\geq\sum_{i\in G}w^{(0)}_{i}-w^{(t)}_{i}$ by telescoping. Note that we can only remove at most $2\epsilon$ mass from $w$ total, as $\sum_{i\in B}w^{(0)}_{i}-w^{(t)}_{i}\leq\epsilon$. Denote for shorthand normalized weights $v^{(t)}:=\tfrac{w^{(t)}}{\left\lVert w^{(t)}\right\rVert_{1}}$. Then, the following is immediate by $\left\lVert w^{(t)}\right\rVert_{1}\geq 1-2\epsilon$.

Under Assumption 2, in all iterations $t$ of Algorithm 6, $v^{(t)}\in\mathfrak{S}_{2\epsilon}^{n}$.

Using Lemma 14, we show that the output has the desired quality of being a large eigenvector.

Under Assumption 2, let the output of Algorithm 6 be $u_{T}$. Then for a universal constant $C^{\star}$, $u_{T}^{\top}\bm{\Sigma}u_{T}\geq(1-C^{\star}\epsilon\log\epsilon^{-1})\|\bm{\Sigma}\|_{\infty}$.

We assume for now that $u_{T}$ is an exact top eigenvector, and discuss inexactness while proving Theorem 1. By (33) and Lemma 14, as then the normalized restriction of $w^{(T)}$ to $G$ is in $\mathfrak{S}_{3\epsilon}^{n}$,

We used the Courant-Fischer characterization of eigenvalues, and that $u_{T}$ is a top eigenvector of $\mathbf{M}(w^{(T)})$. Moreover, by termination conditions and Corollary 5 (correctness of $\mathsf{1DRobustVariance}$),

Combining these two bounds and rescaling yields the conclusion. ∎

Finally, we prove our main guarantee about Algorithm 6. See 1

First, we will operate under Assumption 3, which holds with probability at least $1-\delta$. It is clear that the analyses of Lemma 13 and 15 hold with $1-\Theta(\epsilon\log\epsilon^{-1})$ multiplicative approximations of top eigenvector computation, which the power method approximates with high probability. Thus, each iteration takes time $O\left(\frac{nd}{\epsilon}\log\frac{n}{\delta\epsilon}\right)$, where we will union bound over the number of iterations. We now give an iteration bound: in any iteration where we do not terminate, Lemma 2 implies

Appendix C Deferred proofs from Section 4

Since our notion of approximation is multiplicative, we can assume without more than constant loss that $\mathbf{A}$ has bounded entries. This observation is standard, and formalized in the following lemma.

Feasibility of Problem 2 is unaffected (up to constants in $\epsilon$) by removing columns of $\mathbf{A}$ with entries larger than $n\epsilon^{-1}$.

If $\mathbf{A}_{ji}>n\epsilon^{-1}$ for any entry, then $x_{i}\leq\tfrac{\epsilon(1+\epsilon)}{n}$, else $\left\lVert\mathbf{A}x\right\rVert_{p}$ is already larger than $1+\epsilon$. Ignoring all such entries of $x$ and rescaling can only change the objective by a $1+O(\epsilon)$ factor. ∎

As $g\leq\mathbf{1}\implies\delta\leq p^{-1}\mathbf{1}$, $\frac{\mathbf{A}(\delta\circ w_{t})}{\mathbf{A}w_{t}}\leq p^{-1}$ entrywise. Via $(1+x)^{p}\leq\exp(px)\leq 1+px+p^{2}x^{2}$ for $x\leq p^{-1}$, it follows that

By direct manipulation of the above quantity, and recalling we defined $v=\tfrac{\mathbf{A}w}{\left\lVert\mathbf{A}w\right\rVert_{p}}$,

Using $(1+x)^{p}>1+px$, i.e. $(1+px)^{1/p}<1+x$, we thus obtain

Cauchy-Schwarz yields that $[\mathbf{A}(\delta\circ w)]_{j}^{2}\leq[\mathbf{A}(\delta^{2}\circ w)]_{j}[\mathbf{A}w]_{j}$, $\forall j\in[d]$. Substituting into the above,

Finally, to bound this latter quantity, since $\delta=\eta g$, we observe that for all $j$ either $\delta_{j}=0$ or $1+p\delta_{j}=1+g_{j}=2-[\mathbf{A}^{\top}v^{p-1}]_{j}$, in which case

Thus, plugging this bound into (39) entrywise,

C.2 Proofs from Section 4.3

Our analysis of Algorithm 3 will use the following helper fact.

Feasibility of Problem 3 is unaffected (up to constants in $\epsilon$) by removing matrices $\mathbf{A}_{i}$ with an eigenvalue larger than $n\epsilon^{-1}$.

The proof is identical to Lemma 16; we also require the additional fact that the Schatten norm $\left\lVert\cdot\right\rVert_{p}$ is monotone in the Loewner order, forcing the constraint $x_{i}\leq\tfrac{\epsilon(1+\epsilon)}{n}$. ∎

We remark that we can perform this preprocessing procedure via power iteration on each $\mathbf{A}_{i}$.

Drop $t$ and define $\delta=\eta g$. For simplicity, define the matrices

We recall the Lieb-Thirring inequality $\textup{Tr}((\mathbf{A}\mathbf{B}\mathbf{A})^{p})\leq\textup{Tr}(\mathbf{A}^{2p}\mathbf{B}^{p})$. Applying this, we have

As $g\leq\mathbf{1}$, we have $\mathbf{M}_{0}^{-\frac{1}{2}}\mathbf{M}_{1}\mathbf{M}_{0}^{-\frac{1}{2}}\preceq p^{-1}\mathbf{I}$. Applying the bounds $(\mathbf{I}+\mathbf{M})^{p}\preceq\exp(p\mathbf{M})\preceq\mathbf{I}+p\mathbf{M}+p^{2}\mathbf{M}^{2}$ for $\mathbf{M}=\mathbf{M}_{0}^{-\frac{1}{2}}\mathbf{M}_{1}\mathbf{M}_{0}^{-\frac{1}{2}}$, where we use that $\mathbf{I}$ commutes with all $\mathbf{M}$, it follows that

Definitions of $\mathbf{M}_{0}$, $\mathbf{M}_{1}$, $\mathbf{M}_{2}$, and preservation of positiveness under Schur complements imply

Thus, $\mathbf{M}_{1}\mathbf{M}_{0}^{-1}\mathbf{M}_{1}\preceq\mathbf{M}_{2}$. Applying this and recalling $\mathbf{V}=\tfrac{\mathbf{M}_{0}}{\left\lVert\mathbf{M}_{0}\right\rVert_{p}}$,

By $(1+px)^{1/p}<1+x$, taking $p^{th}$ roots we thus have

Finally, the conclusion follows as in Lemma 3; by linearity of trace and $g=p\delta$,

Here, we used the inequality for all nonzero $g_{i}$,

The proof is analogous to that of Theorem 4; we sketch the main differences here. By applying Lemma 17 and monotonicity of Schatten norms in the Loewner order, we again have $\Phi_{0}\leq 1$, implying correctness whenever the algorithm terminates on Line 4. Correctness of dual certification again follows from lack of termination and the choice of $T$, as well as setting $u$ to indicate each coordinate. Finally, the returned matrix in Line 8 is correct by convexity of the Schatten-$q$ norm, and the fact that all $\mathbf{V}_{t}^{p-1}$ have unit Schatten-$q$ norm.

We now discuss issues regearding computing $g_{t}$ in Line 5 of the algorithm, the bottleneck step; these techniques are standard in the approximate SDP literature, and we defer a more formal discussion to e.g. [JLL+20]. First, note that each coordinate of $g_{t}$ requires us to compute

We estimate the two quantities in the above expression each to $1+\epsilon$ multiplicative error with high probability. Union bounding over iterations, and modifying Lemma 4 to use the potential $\left\lVert\sum_{i\in[n]}[w_{t}]_{i}\mathbf{A}_{i}\right\rVert_{p}-(1+O(\epsilon))\left\lVert w_{t}\right\rVert_{1}$, the analysis remains valid up to constants in $\epsilon$ with this multiplicative approximation quality. We now discuss our approximation strategies.

For shorthand, denote $\mathbf{M}=\sum_{i\in[n]}[w_{t}]_{i}\mathbf{A}_{i}$. To estimate the denominator of (40), it suffices to multiplicatively approximate $\left\lVert\mathbf{M}\right\rVert_{p}^{p}=\textup{Tr}[\mathbf{M}^{p}]$ within a $1+\epsilon$ factor, as raising to the $\tfrac{p-1}{p}$ power can only improve this. To do so, we use the well-known fact (e.g. [DG03]) that letting $\mathbf{Q}$ be a $k\times d$ matrix with independent entries $\sim\mathcal{N}(0,\tfrac{1}{k})$, for $k=O(\tfrac{\log(\frac{nd}{\epsilon})}{\epsilon^{2}})$, with probability $1-\textup{poly}((\tfrac{nd}{\epsilon})^{-1})$,

We can simultaneously compute all such quantities by first applying $O(p)$ matrix-vector multiplications through $\mathbf{M}$ to each row of $\mathbf{Q}$, and then computing all quadratic forms. In total, the computational cost per iteration of all approximations is $O(\textup{nnz}\cdot\tfrac{p\log(\frac{nd}{\epsilon})}{\epsilon^{2}})$ as desired. ∎

C.3 Proof of Proposition 2

In this section, following our prior developments, we prove the following claim.

Given access to an oracle for the following approximate decision problem, we can implement an efficient binary search for estimating OPT. Specifically, letting the range of OPT be $(\mu_{\text{lower}},\mu_{\text{upper}})$, we can subdivide the range into $O(\tfrac{1}{\epsilon}\log\tfrac{\mu_{\text{upper}}}{\mu_{\text{lower}}})$ multiplicative intervals of range $1+\epsilon$, and then compute a binary search using our decision oracle. This incurs a multiplicative $\log(\tfrac{nd}{\epsilon})$ overhead in the setting of Proposition 2 (see Appendix A, [JLL+20], for a more formal treatment).

C.3.2 Preliminaries

up to multiplicative $1+\epsilon$ tolerance on either side. Consider the potential function

It is clear that the first term of $\Phi(w)$ approximates the left hand side of (41) up to a $\log 2$ additive factor, so if any of $\left\lVert\mathcal{A}(w)\right\rVert_{p}$, $\left\lVert\mathcal{A}(w)\right\rVert_{p^{\prime}}$, or $\left\lVert w\right\rVert_{1}$ reaches the scale $3\epsilon^{-1}$ and $\Phi(w)$ is bounded by $1$, we can safely terminate. and conclude primal feasibility for Problem 4. Next, we compute

The following helper lemma will be useful in concluding dual infeasibility of Problem 4.

In the setting of Problem 4, suppose there exists $x^{*}\in\Delta^{n}$ with

From the definitions in (43), it is clear that $\left\lVert\mathbf{Y}(w)\right\rVert_{q}=\left\lVert z(w)\right\rVert_{q^{\prime}}=1$, where $q$, $q^{\prime}$ are the dual norms of $p$, $p^{\prime}$ respectively. Moreover, by the definition of $x^{*}$, we have for all $\left\lVert\mathbf{Y}\right\rVert_{q}=\left\lVert z\right\rVert_{q^{\prime}}=1$,

as desired (here, we used positivity of all relevant quantities). ∎

C.3.3 Potential monotonicity

We prove a monotonicity property regarding the potential $\Phi$ in (42).

Denote for simplicity the threshold $K=3\epsilon^{-1}$ and the step vector $\delta=\eta g$. First, by prior calculations in Lemma 3 and Lemma 4, it follows that

Next, note that by $\delta\leq\eta$ entrywise and lack of termination (i.e. the threshold $K$),

Therefore, by $\exp(x)\leq 1+x+x^{2}$ for $x\leq 1$,

Moreover, by applying Cauchy-Schwarz and the threshold $\left\lVert\mathcal{A}(w)\right\rVert_{p}\leq K$ once more,

Combining (44) and (45) (and applying similar reasoning to the term $\Delta_{\mathbf{S}}$), we conclude

Recall the inequality $\log(1+x)\leq x$ for nonnegative $x$. Expanding the definition of $\Phi$ and $\nabla\Phi$ (cf. (42)), and plugging in the above bounds, we conclude that

As before, we show that this sum is entrywise nonpositive. For any $i\in[n]$ with $\delta_{i}\neq 0$, we have

as desired, where we used that $\eta^{-1}\geq p^{\prime}+2K$. This yields the conclusion $\Phi(w^{\prime})\leq\Phi(w)$. ∎

C.3.4 Algorithm and analysis

Finally, we state Algorithm 7 and prove Proposition 2.

Correctness of the reduction to deciding Problem 4 follows from the discussion in Section C.3.1. Moreover, by the given Algorithm 7, it is clear (following e.g. the preprocessing of Lemma 17) that $\Phi(w_{t})\leq 1$ throughout the algorithm, so whenever the algorithm terminates we have primal feasibility. It suffices to prove that whenever the problem admits $x^{*}$ with

then the algorithm terminates on Line 5 in $T$ iterations. Analogously to Theorem 4, we have

Next, since $g_{t}$ is an upwards truncation of $\nabla\Phi(w_{t})$, applying Lemma 18 implies that

The conclusion follows by the definition of $T$, as desired. Finally, the iteration complexity follows analogously to the discussion in Theorem 5’s proof, where the only expensive cost is estimating coordinates of the $\mathcal{A}$ component of $\nabla\Phi(w_{t})$ every iteration. ∎

Finally, we remark that by opening up the dual certificates $\mathbf{Y}(w)$, $\mathbf{Z}(w)$ of our mirror descent analysis, we can in fact implement a stronger version of the decision Problem 4 which returns a feasible dual certificate whenever the primal problem is infeasible. We omit this extension for brevity, as it is unnecessary for our applications, but it is analogous to the analysis of Theorem 5.

Appendix D Deferred proofs from Section 5

We claim that Algorithm 1 in [MM15] applied to the matrix $\mathbf{A}^{p}$ with a careful choice of exponent $q$ in their Algorithm 1 yields this guarantee. Specifically, we choose $q_{1},q_{2}$, both of which satisfy the criteria in their main theorem, such that the iterates produced by simultaneous power iteration $\mathbf{M}^{p}$ with exponent $q_{1}$ and $\mathbf{M}^{p-1}$ with exponent $q_{2}$ are identical; it suffices to choose $q$ a multiple of $p(p-1)$. Thus, we can also apply their guarantees to $\mathbf{A}^{p-1}$ and apply a union bound. Notice that their Algorithm 1 also contains some postprocessing to ensure that they obtain singular values in the right space, which is unnecessary for us, as our matrices are Hermitian. ∎

D.2 Proof of Lemma 5

D.3 Proof of Lemma 6

We follow the notation of (10). First, by the guarantees of Corollary 1,

Therefore, again applying Corollary 1, for all $i\in G$,

We conclude that the set of weights $\{\tfrac{w_{i}}{w_{G}}\}_{i\in G}$ belong to $\mathfrak{S}_{3\epsilon}^{(1-\epsilon)n}$. By applying Corollary 4 to these weights and adjusting the definition of $C_{3}$ by a constant, we conclude with probability at least $1-\tfrac{\delta}{2}$