List Decodable Mean Estimation in Nearly Linear Time

Introduction

Estimating the mean of data is a cardinal scientific task. The population mean can be shifted arbitrarily by a single outlier, a problem which is compounded in high dimensions where outliers can conspire to destroy the performance of even sophisticated estimators of central tendency. Robust statistics, beginning with the works of Tukey and Huber [Tuk60, Hub64], endeavors to design, model, and mitigate the effect of data deviating from statistical assumptions [Hub11].

A canonical model of data corruption is the Huber contamination model [Hub64]. Let $I(\mu)$ be a probability distribution parameterized by $\mu$. We say a dataset $X_{1},X_{2},...,X_{N}$ is $\alpha$-Huber contaminated for some constant $\alpha\in$ if it is drawn i.i.d from

where $\mathcal{O}$ is an arbitrary outlier distribution which can be adversarial and dependent on $\mathcal{I}(\mu)$. The goal is to estimate $\mu$ with an estimator $\hat{\mu}$ such that the two are close with respect to a meaningful metric. The Huber contamination model captures the setting where only an $\alpha$ fraction of the dataset is subject to statistical assumptions. One would hope to design estimators $\hat{\mu}$ for which $\alpha$ is as small as possible thereby tolerating the largest fraction of outliers–a quantity known as the breakdown point . The study of estimators with large breakdown points is the focus of a long and extensive body of work, which we do not attempt to survey here. For review see [Hub11, HRRS86].

A first observation, is that the breakdown point of a single estimator must be smaller than $\frac{1}{2}$. For concreteness, consider the problem of estimating the mean of a standard normal. The adversary can set up a mixture of $\frac{1}{\alpha}$ standard normals for which the means of the mixture components are far apart. This intrinsic difficulty also gives rise to a natural notion of recovery in the presence of overwhelming outliers. Instead of outputting a single estimator, consider outputting a list of candidate estimators $\mathcal{L}=\{\hat{\mu}_{1},\hat{\mu}_{2},...,\hat{\mu}_{\frac{1}{\alpha}}\}$ with the guarantee that the true $\mu$ is amongst the elements of the list. This is the setting of ’List Decodable Learning’ [BBV08, CSV17], analogous to list decoding in the theory of error correcting codes.

In their influential work [CSV17] introduces list decodable learning in the context of robust statistics. They consider the problem of estimating the mean $\mu$ of a $d$-dimensional distribution $\mathcal{I}(\mu)$ with a bounded covariance $\mathbf{Cov}(\mathcal{I}(\mu))\preceq\sigma^{2}I$ for a constant $\sigma$ from $N=\frac{d}{\alpha}$ samples. Their algorithm recovers a list $\mathcal{L}$ of $O(\frac{1}{\alpha})$ candidate means with the guarantee that there exists a $\hat{\mu}^{*}\in\mathcal{L}$ achieving the recovery guarantee $\lVert\hat{\mu}^{*}-\mu\rVert\leqslant O\left(\sigma\sqrt{\frac{\log\left(\frac{1}{\alpha}\right)}{\alpha}}\right)$ with high probability $1-\frac{1}{\text{\rm poly}(d)}$. Furthermore, their algorithm is ’efficient’, running in time $\text{\rm poly}(N,d,\frac{1}{\alpha})$ via the polynomial time solvability of ellipsoidal convex programming.

Our first contribution is an algorithm for list decodable mean estimation of covariance bounded distributions, which outputs a list $\mathcal{L}$ of length $O(\frac{1}{\alpha})$, achieving (up to constants) the information theoretically optimal recovery $O(\frac{\sigma}{\sqrt{\alpha}})$, with linear sample complexity $N=\frac{d}{\alpha}$, and running in nearly linear time $\widetilde{O}(Nd\text{\rm poly}(\frac{1}{\alpha}))$ where $\widetilde{O}$ omits logarithmic factors in $d$. For the matching minimax $\Omega(\frac{\sigma}{\sqrt{\alpha}})$ lower bound see [DKS18]. Formally, we state our main theorem.

For precise constants and failure probability see Section 4. At a high level, we define a nonconvex cost function for which $\mu$ is an approximate minimizer and build a ’descent style’ algorithm to find $\mu$. As with most nonconvex algorithms, our approach is susceptible to falling in suboptimal minima. Our key algorithmic insight is that our algorithm fails to descend the cost function exactly when a corresponding dual procedure succeeds in ”sanitizing” the dataset by removing a large fraction of outliers — a win-win.

First observed in [CSV17], the list decoding problem lends itself to applications for which our algorithm offers immediate improvements. Firstly, it is perhaps surprising that a succinct list of estimators can be procured from a dataset overwhelmed by outliers. Perhaps more surprising is that the optimal candidate mean can be isolated from the list $\mathcal{L}$ with additional access to a mere $\log(\frac{1}{\alpha})$ clean samples drawn from $\mathcal{I}(\mu)$. This ”semi-supervised” learning is compelling in settings where large quantities of data are collected from unreliable providers (crowdsourcing, multiple sensors, etc.). Although it is resource intensive to ensure the cleanliness of a large dataset, it is easier to audit a small, in our case $\log(\frac{1}{\alpha})$, set of samples for cleanliness. Given access to this small set of samples as side information, our algorithm returns estimators for mean estimation with breakdown points higher than $\frac{1}{2}$ in nearly linear time.

Faster list decodable mean estimation also accelerates finding planted partitions in semirandom graphs. In particular, consider the problem where $G$ is a directed graph where the (outgoing) neighborhoods of an $\alpha$ fraction of vertices $S$ are random while the neighborhoods of the remaining vertices are arbitrary, and the goal is to output $O(1/\alpha)$ lists such that one of them is “close” to $S$. Our algorithm for list decodable mean estimation implies a faster algorithm for this problem as well.

Lastly, list decodable mean estimation is a superset of learning mixture models of bounded covariance distributions with minimum mixture weight $\alpha$. By treating a single cluster as the inliers, one can recover the list of means comprising the mixture model. Notably, this can be done without any separation assumptions between the mixture components and is robust to outliers.

Fast Semidefinite Programming:

Rapidly computing our cost function necessitates the design of new packing/covering solvers for Positive Semidefinite Programs (SDP) over general Fantopes (the convex hull of the projection matrices). Positive SDP’s have seen remarkable success in areas spanning quantum computing, spectral graph theory, and approximation algorithms (See [AHK12, ALO16, JLL+20] and the references therein). Informally, a packing SDP computes the fractional number of ellipses that can be packed into a spectral norm ball which involves optimization over the spectrahedron. A natural question is whether the packing concept can be extended to balls equipped with general norms, say the sum of the top $k$ eigenvalues (the Ky Fan norm), where for $k=1$ we recover the oft studied spectral norm packing. We use results from Loewner’s theory of operator monotonicity and operator algebras to design fast, and as far as we know the first solvers for packing/covering positive SDP’s under Ky Fan norms (See Theorem 5.13).

2 Related Work

Robust statistics has a long history [Tuk60, Tuk75, Hub64, Ham71]. This extensive body of work develops the theory of estimators with high breakdown points, of influence functions and sensitivity curves, and of designing robust M-estimators. See [Hub11, HRRS86]. However, little was understood about the computational aspects of robustness which features prominently in high dimensional settings.

Recent work in theoretical computer science [DKK+16, LRV16] designed the first algorithms for estimating the mean and covariance of high dimensional gaussians tolerating a constant fraction of outliers in polynomial time $poly(N,d,\frac{1}{1-\alpha})$. Since then, a flurry of work has emerged studying robust regression [KKM18, DKS19], sparse robust regression [BDLS17, DKK+19], fast algorithms for robustly estimating mean/covariance [CDG19, DHL19, CDGW19], statistical query hardness of robustness [DKS17], worst case hardness [HL19], robust graphical models [CDKS18], and applications of the sum of squares algorithm to robust statistics [KSS18]. See survey [DK19] for an overview.

List Decodable Learning

Despite the remarkable progress in robust statistics for large $\alpha$ contamination, progress on the list decoding problem has been slower. This is partially owed to the intrinsic computational hardness of the problem. Even for the natural question of list decoding the mean of a high dimensional gaussian, [DKS18] exhibits a quasipolynomial time lower bound against Statistical Query algorithms for achieving the information theoretically optimal recovery of $\Theta\big(\sqrt{\log(\frac{1}{\alpha})}\big)$. This stands in contrast to large $\alpha$ robust mean estimation where nearly linear time algorithms [CDG19] achieve optimal recovery.

In light of this hardness, a natural question is to determine whether polynomial time algorithms can at least approach the optimal recovery for list decoding the mean of a gaussian. In a series of concurrent works [KS17] [DKS18], develop the first algorithms approaching the $\Theta\big(\sqrt{\log(\frac{1}{\alpha})}\big)$ recovery guarantee. At a high level, both papers achieve recovery $O(\frac{\sigma}{\alpha^{c/k}})$ for different fixed constants $c>1$ in time $\text{\rm poly}(\frac{d}{\alpha})^{O(k)}$ for $k$ a positive integer greater than $2$. The [DKS18] algorithm, known as the ”multi-filter”, is a spectral approach reasoning about high degree polynomials of the moments of data. Furthermore, the ”low degree” multi-filter achieves a suboptimal $O\big(\sqrt{\frac{\log(\frac{1}{\alpha})}{\alpha}}\big)$ recovery guarantee for list decoding the mean of subgaussian distributions, which is fast and may be of practical value. [KS17] develop a convex hierarchy (sum of squares) style approach, which achieve similar guarantees for more general distributional families satisfying a poincare inequality. In particular for list decoding the mean of bounded covariance distributions they achieve the optimal $O(\frac{1}{\sqrt{\alpha}})$ guarantee via the polynomial time solvability of convex concave optimization. Finally, [DKS18, KS17] and a concurrent work [HL18] develop tools for reasoning about the high degree moments of data to break the longstanding ”single-linkage” barrier in clustering mixtures of spherical gaussians.

In other statistical settings a series of concurrent works [RY20a, KKK19] demonstrate information theoretic impossibility for list decoding regression even under subguassian design. Similar barriers arise in the context of list decodable subspace recovery [RY20b, BK20] where it is information theoretically impossible to list decode a dataset for which an $\alpha$ fraction is drawn from a subgaussian distribution in a subspace. Indeed, since list decoding is a superset of learning mixture models, these hardness considerations stem from barriers in learning mixtures of linear regressions and subspace clustering. On the other hand, the above works also construct polynomial time, $d^{\text{\rm poly}(\frac{1}{\alpha})}$, algorithms for regression and subspace recovery for Gaussian design and Gaussian subspaces respectively, which holds true for a larger class of ”certifiably anticoncentrated” distributions.

In this backdrop of computational and statistical hardness, and given the practical value of robust statistics, it is a natural challenge to design list decoding algorithms that are both fast and statistically optimal. The current work is a step in this direction.

SDP Solvers

There has been much recent interest in designing fast algorithms for positive SDP solvers due to the ubiquity of their application in approximation algorithms. We do not attempt to survey the full breadth of these results and their applications in this section. We refer the interested reader to [JLL+20, ALO16, PTZ12, AHK12] for more context on these developments. We will restrict ourselves to the following class of SDPs relevant to our work:

Semirandom Graph Inference

The study of problems that are typically computationally hard in the worst case in semirandom graph models was initiated by [BS95] and perpetuated by [FK01]. A specific problem of interest to us studied by [FK01] for which nearly optimal algorithms were given by [MMT20] is the semirandom independent set problem where the set of edges between a planted independent set and the remaining (adversarially chosen) graph come from a randomized model. In a similar vein [CSV17] studies a planted partition where instead of an independent set the given graph is some other sparse random graph (albeit directed). Our results improve upon the statistical guarantees of [CSV17] as well as give faster algorithms, however both [CSV17] and our work fall short of capturing the results of [MMT20] due to the directed model we work in. However, we believe the hurdle is a technical point rather than an inherent shortcoming of our approach.

Sample Complexity:

The following lemma of [CSV17] achieves linear sample complexity which suffices for our algorithm.

Taking $N=O\left(\frac{d}{\alpha}\right)$, for the rest of the paper we will adjust $\sigma$ by a constant and assume the inlier set $I$ satisfies $\lVert\frac{1}{|I|}\sum_{i\in\mathcal{I}}(x_{i}-\mu)(x_{i}-\mu)^{T}\rVert\preceq\sigma^{2}I$.

Notation:

Organization:

Our paper is organized as follows: In Section 2, we outline the key ideas underlying the design of our algorithm for list-decodable mean estimation, our solver for the generalized class of Packing/Covering SDPs considered in this paper and the technical challenges involved in doing so. Then, in Sections 4 and 8, we formally describe and analyze our algorithm for list-decodable mean estimation and its application to the semirandom graph model considered in [CSV17]. Sections 6, 5 and 7 contain our refined power method analysis, a formal description and analysis of our solver and the hard thresholding based operator required to implement the solver in nearly-linear time. Finally, Appendices B, A and C contain supporting results required by the previous sections.

Techniques

First we present an inefficient algorithm for list decodable mean estimation. Although it is inefficient, it captures the core ideas and foreshadows the difficulties encountered by our efficient algorithm. At a high level, the inefficient algorithm greedily searches through the dataset for subsets of points with small covariance with the goal of finding the subset of inliers.

Second, append $\hat{\mu}=\sum_{i=1}^{N}\hat{w}_{i}x_{i}$ to $\mathcal{L}$

Third, update $b$ such that $b_{i}=b_{i}-\hat{w}_{i}$

We claim the algorithm outputs a list $\mathcal{L}$ of length $\frac{2}{\alpha}$ and that there exists a $\hat{\mu}^{*}\in L$ satisfying $\lVert\hat{\mu}^{*}-\mu\rVert\leqslant O(\frac{\sigma}{\sqrt{\alpha}})$. Next we outline the proof of correctness.

Proof Outline:

Sanitizing the Dataset:

Abstracting the guarantees of our inefficient algorithm, we say that an algorithm ”sanitizes” a dataset if it outputs a tuple $(\hat{\mu},\hat{w})$ where $\sum_{i=1}^{N}\hat{w}\geqslant\Omega(1)$ satisfying the following conditions. If $\lVert\hat{\mu}-\mu\rVert\geqslant O(\frac{\sigma}{\sqrt{\alpha}})$ then $\sum_{i\in O}\hat{w}_{i}\geqslant\frac{2}{\alpha}\sum_{i\in I}\hat{w}_{i}$. Any algorithm that sanitizes the dataset iteratively, is guaranteed to succeed as a list decoding algorithm. This is made formal in Section 4.

Descent Style Formulation:

The optimization problem Eq. 1 is nonconvex and hard to solve directly. A novel approach to minimizing Eq. 1 is to replace $\mu(w)$ with a parameter $\nu$ and define a cost function $f(\nu)$. First introduced in [CDG19] in the context of robust mean estimation and later in robust covariance estimation [CDGW19] consider the function $f(\nu)$ defined as follows:

where $b_{i}=\frac{1}{\alpha N}$ for all $i\in[N]$. This formulation has two appealing aspects. Firstly, the cost function can be computed efficiently via convex concave optimization. Indeed, the operator norm can be replaced by the maximization over its associated fantope $\mathcal{F}_{1}$

Secondly, for $\alpha>\frac{2}{3}$ (robust mean estimation), a crucial insight of [CDG19] is that $f(\nu)$ approximates the squared distance from $\nu$ to the mean $\mu$. Then a good estimate of the mean is the minimizer of the cost.

In their setting the minimization in Eq. 2 can be performed by a descent style algorithm.

Substantial challenges arise when designing such a cost function for list decodable mean estimation. Chiefly, the inliers are unidentifiable from the dataset so there is no function of the data that approximates the distance to the true mean. Our solution is to design a function that either approximates the distance to the true mean, or when the approximation is poor, prove there exists a corresponding dual procedure that sanitizes the dataset. This win-win observation can be made algorithmic and is the subject of Section 4

1 Our Approach

We call the above min-max formulation the dual and the associated minimizer $w^{*}$ the dual minimizer or dual weights. By Von Neumann’s min max theorem we have

An Easier Problem:

List Decoding Main Lemma:

In analogy to clustering, one should hope that for any $\nu\in R^{d}$ further than $O(\frac{\sigma}{\sqrt{\alpha}})$ from $\mu$, that $Cost(\nu)\approx\lVert\nu-\mu\rVert^{2}$. Although this is impossible, it turns out that when it is false, there exists a corresponding ”dual procedure” for outputting a sanitizing tuple. More precisely, we claim that either $0.4\lVert\nu-\mu\rVert^{2}\leqslant Cost(\nu)\leqslant 1.1\lVert\nu-\mu\rVert^{2}$, or a simple procedure outputs a set of weights $\hat{w}$ identifying vastly more outliers than inliers i.e $\sum_{i\in\mathcal{O}}\hat{w}_{i}\geqslant\frac{\alpha}{2}\sum_{i\in\mathcal{I}}\hat{w}_{i}$, or both.

The dual procedure is as follows. Let $\widehat{\Sigma}:=\sum_{i=1}^{n}w^{*}_{i}(x_{i}-\nu)(x_{i}-\nu)^{T}$ be the weighted second moment matrix centered at $\nu$. Let $V$ be the top $O(\frac{1}{\alpha})$ eigenspace of $\widehat{\Sigma}$. We project the dataset onto the affine subspace $V$ with offset $\nu$. We then sort the points $\{\Pi_{V}(x_{i}-\nu)\}_{i=1}^{N}$ by Euclidean lengths. This sorting determines an ordering of the weights $w^{*}_{1},...,w^{*}_{N}$. We pass through the sorted list, and find the smallest $m\in[N]$ such that $\sum_{i=1}^{m}w^{*}_{i}\geqslant 0.5$. We set $\hat{w}_{i}=w^{*}_{i}$ for $i=1,...,m$ and $\hat{w}_{i}=0$ for $i>m$. The following lemma guarantees $\sum_{i\in\mathcal{O}}\hat{w}_{i}\geqslant\frac{\alpha}{2}\sum_{i\in\mathcal{I}}\hat{w}_{i}$.

2 Generalized Packing/Covering Solvers and Improved Power Method Analysis

We start by considering the simpler problem of computing $Cost_{X,b,1}(\nu)$. The approach taken in [CDG19] is to reduce the problem to a packing SDP via the introduction of an additional parameter $\lambda$; specifically, they solve the following packing SDP:

for which there exist fast linear-time solvers [PTZ12, ALO16]. It can be shown that the value of the above program when viewed as a function of $\lambda$ is monotonic, continuous and attains the value $1$ precisely when $\lambda=Cost_{X,b,1}(\nu)$. Therefore, by performing a binary search over $\lambda$, one obtains accurate estimates of $w^{*}$ and $Cost_{X,b,1}(\nu)$.

which does not fall into the standard class of packing SDPs. We extend and generalize fast linear time solvers for packing/covering SDPs from [PTZ12] to this broader class of problems. However, this generalization is not straightforward.

To demonstrate the main difficulties, we will delve more deeply into the solver from [PTZ12] and state the packing/covering primal dual pairs they consider:

The algorithm then proceeds to increment the weights of all $i$ such that $\langle P_{1},A_{i}\rangle\leqslant(1+\varepsilon)$ for a user defined accuracy parameter, $\varepsilon$, by a multiplicative factor. Intuitively, these indices correspond to “directions”, $A_{i}$, along which $\sum_{i=1}^{N}w_{i}A_{i}$ is small and therefore, their weights can be increased in the dual formulation. By incorporating a standard regret analysis from [AK16] for the matrices, $P_{1}$, they show that one either outputs a primal feasible, $M$, with $\Tr M\leqslant 1$ or a dual feasible $w$ with $\sum_{i=1}^{N}w_{i}\geqslant(1-\varepsilon)$.

Preliminaries

The following can be found in [Cha15, Example 13(iii)]:

Suppose $A$ and $B$ are positive semidefinite matrices such that $A\succcurlyeq B$, then $\log A\succcurlyeq\log B$.

2 Optimization

Similarly, we call $f$ $\alpha$-strongly concave with respect to $\|\cdot\|_{\square}$ if

A key property of von Neumann entropy we use is:

$\mathsf{vNE}$ is $1$-strongly concave with respect to the trace norm.

The following is an immediate consequence of Fact 3.5:

We emphasize that we slightly deviate from the convention that von Neumann entropy and quantum relative entropy are defined only on PSD matrices of trace exactly $1$.

Algorithm for List-Decodable Mean Estimation

If $\delta=0$ we would compute cost exactly. This is computationally expensive so we take $\delta$ to be $0.01$. For ease of reading, one can first set $\delta=0$ with the understanding that the algorithmic lemmas succeed for small $\delta$.

For a dataset $X=\{x_{1},\dots,x_{N}\}$, an inlier set $I\subseteq[N]$ with $|I|=\alpha N$, budgets $b_{1},\dots,b_{N}\in\left[0,\frac{2}{\alpha N}\right]$, we say that $(\hat{\mu},\hat{w})$ is a sanitizing tuple for $(X,I,b)$ if the it satisfies:

If $\lVert\hat{\mu}-\mu\rVert\geqslant 2\cdot 10^{3}\frac{\sigma}{\sqrt{\alpha}}$, then $\sum_{i\in I}\hat{w}_{i}\leqslant\frac{\alpha}{4}$.

For all $i$: $0\leqslant\hat{w}_{i}\leqslant b_{i}$, $\|\hat{w}\|_{1}\geqslant 0.5$.

Recall that when the algorithm terminates $\theta^{(t)}\leqslant\sigma^{2}$ or $\theta^{(t+1)}\geqslant 0.5\cdot\theta^{(t)}$. We prove that $(\hat{\mu},\hat{w})$ is a sanitizing tuple when $\theta^{(t)}\leqslant\sigma^{2}$ in Lemma A.1, and when $\theta^{(t+1)}\geqslant 0.5\cdot\theta^{(t)}$ in Lemma 4.6. ∎

Let $X=\{x_{1},...,x_{N}\}$ be a dataset with an inlier set $I$ of size $|I|=\alpha N$ satisfying $Cov_{I}(x)\preceq\sigma^{2}I$. Let $\mu=\sum_{i\in I}x_{i}$. $OutputList(X,\alpha)$ returns a list $\mathcal{L}$ of length $O(\frac{1}{\alpha})$ such that there exists $\hat{\mu}^{*}\in\mathcal{L}$ satisfying $\lVert\hat{\mu}^{*}-\mu\rVert\leqslant r\frac{\sigma}{\sqrt{\alpha}}$ with high probability $1-\frac{1}{d^{10}}$

The proof of the corollary is elementary and similar to the proof of correctness for the inefficient algorithm. For example, see Appendix A.

A brute force search through $V$ is inefficient. Instead, we project the dataset $X^{(t)}$ onto $V$. We evaluate the cost at $p=O(\frac{\log(d)}{\log(\frac{1}{1-\alpha})})$ randomly chosen projected datapoints. We choose $\nu^{(t+1)}$ to be the projected datapoint with the smallest cost. Although the projected datapoints are by no means an exhaustive search of the subspace $V$, if this procedure fails to make sufficient progress i.e $\theta^{(t+1)}\geqslant 0.5\theta^{(t)}$ then a corresponding procedure uses the dual weights to output a sanitizing tuple.

The weight removal procedure is as follows. We sort the vectors $\{\Pi_{V}(x_{i}-\nu^{(T)})\}_{i=1}^{N}$ by euclidean norm. This sorting determines an ordering of the weights $\bar{w}^{(T)}_{1},...,\bar{w}^{(T)}_{N}$. We pass through the sorted list, and find the smallest $m\in[N]$ such that $\sum_{i=1}^{m}\bar{w}^{(T)}_{i}\geqslant 0.5$. We set $\hat{w}_{i}=\bar{w}^{(T)}_{i}$ for $i=1,...,m$ and $\hat{w}_{i}=0$ for $i>m$. Finally we output $(\nu^{(T)},\hat{w})$ as the sanitizing tuple.

DescendCost(X,b) is given dataset $X$ and weight upper bound $b$ satisfying the assumptions of Theorem 4.4. If at iteration $t$, $\theta^{(t+1)}\geqslant 0.5\theta^{(t)}$ then $DescendCost(X,b)$ outputs $(\hat{\mu},\hat{w})$ satisfying $\|\hat{\mu}-\mu\|\leqslant r\frac{1}{\sqrt{\alpha}}$ or $\sum_{i\in I}\hat{w}_{i}\leqslant\frac{\alpha}{4}$.

Proving Lemma 4.6 is the primary objective of the remainder of this section. We will make use of the following two lemmas and defer their proofs to the appendix. The first Lemma 4.7 states that if $(\nu^{(t)},\bar{w}^{(t)})$ is not a sanitizing tuple then the descent procedure succeeds.

Let $(X,b)$ satisfying the assumptions of Theorem 4.4. If at iteration $t$ of DescendCost(X,b), the tuple $(\nu^{(t)},\bar{w}^{(t)})$ is not a sanitizing tuple, then $\theta^{(t+1)}\leqslant 0.04\lVert\mu-\nu^{(t)}\rVert^{2}+ck\sigma^{2}$ for $c=10^{5}$ with high probability $1-\frac{1}{d^{10}}$.

We will also need Lemma 4.8 which states that if the cost is a constant factor smaller than $\lVert\mu-\nu^{(t)}\rVert$ then the weight removal procedure outputs a sanitizing tuple.

Using the above two lemmas we prove Lemma 4.6.

(Proof of Lemma 4.6) Firstly, we observe that at any given iteration, if $\lVert\nu^{(t)}-\mu\rVert\leqslant r\frac{\sigma}{\sqrt{\alpha}}$ then either the descent makes progress or the weight removal procedure outputs a sanitizing tuple. Likewise, if $\sum_{i\in I}\bar{w}^{(t)}_{i}\leqslant\frac{\alpha}{4}$ then either descent makes progress or weight removal outputs a sanitizing tuple. So without loss of generality, we assume $\lVert\nu^{(t)}-\mu\rVert\geqslant r\frac{\sigma}{\sqrt{\alpha}}$ and $\sum_{i\in I}\hat{w}^{(t)}_{i}\geqslant\frac{\alpha}{4}$.

1 Analysis I: Descending Cost

Our first step is to prove that $\mu$ has a large component in $V$ i.e

For some fixed constant $c_{1}$. By definition of projection we have

Here, the inequality follows by dropping squared terms. Now using the inequality $2ab\leqslant\frac{a^{2}}{c^{2}}+c^{2}b^{2}$ for $a=\langle\varphi,\mu-\nu\rangle$, $b=\langle\varphi,\mu-x_{i}\rangle$, $c=10$ we obtain

We use the fact that $\sum_{i\in I}w_{i}\geqslant\sum_{i\in I}\bar{w}_{i}\geqslant\frac{1}{4k}$ to lower bound the first term. And we use the fact $w_{i}\leqslant\frac{2}{(1-\delta)\alpha N}$ to lower bound the second term to obtain

Consider the second term $\sum_{i\in I}\frac{1}{\alpha N}\langle\varphi,x_{i}-\mu\rangle^{2}$. We can upper bound it by the fact that the inliers are covariance bounded $\sum_{i\in I}\frac{1}{\alpha N}\langle\varphi,x_{i}-\mu\rangle^{2}\leqslant\sigma^{2}$. Plugging this bound into (7) we obtain

For a fixed constant $c_{1}$. Rearranging the LHS and RHS we upper bound

with probability greater than $1-\frac{1}{d^{10}}$. Here we aim for a $\frac{1}{d^{10}}$ failure probability so that by union bound over the $O(\frac{1}{\alpha}\log(d))$ iterations of the algorithm we continue to succeed with high probability. Thus we have

Where the first equality is by definition of $\theta^{(t+1)}$, and the inequality follows by Lemma A.2. The above is then:

2 Analysis II: Removing Weights

In this section we prove Lemma 4.8 See 4.8

We need to prove two facts. Firstly, $\mathbf{Pr}_{w}[\lVert\Pi_{V}(x_{i}-\nu)\rVert^{2}<0.8\lVert\nu-\mu\rVert^{2}]\geqslant 0.5$, and secondly $\mathbf{Pr}_{i\in I}[\lVert\Pi_{V}(x_{i}-\nu)\rVert^{2}<0.8\lVert\nu-\mu\rVert^{2}]\leqslant\frac{\alpha}{4}$. Taken together, this implies that a sort of the list $\{\lVert\Pi_{V}(x_{i}-\nu)\rVert\}_{i=1}^{N}$ succeeds in isolating at least $0.5$ weight where $\sum_{i\in I}\hat{w}_{i}\leqslant\frac{\alpha}{4}$. We use Markov’s inequality to prove the first statement.

Now we prove the second statement $\mathbf{Pr}_{i\in I}[\lVert\Pi_{V}(x_{i}-\nu)\rVert^{2}<0.8\lVert\nu-\mu\rVert^{2}]\leqslant\frac{\alpha}{4}$. Let $x\in I$ be an inlier. Let $\rho:=\frac{\Pi_{V}(\mu-\nu)}{\lVert\Pi_{V}(\mu-\nu)\rVert}$. We have that $\|\Pi_{V}(x-\nu)\|^{2}$ is lower bounded by

Where the first inequality follows because $\rho$ is a unit vector in $V$. The second inequality follows by the fact that $2ab\leqslant\frac{a^{2}}{c^{2}}+c^{2}b^{2}$ for $a=\langle\mu-\nu,\rho\rangle$, $b=\langle\mu-x,\rho\rangle$, $c=10$. We further lower bound by

Plugging this lower bound for $\|\Pi_{V}(x-\nu)\|^{2}$ into $\mathbf{Pr}_{i\in I}[\lVert\Pi_{V}(x_{i}-\nu)\rVert^{2}<0.8\lVert\nu-\mu\rVert^{2}]$ we obtain

Where the second inequality follows from applying the bounded covariance of the inliers. The last inequality follows from the assumption that $\lVert\nu-\mu\rVert\geqslant r\frac{\sigma}{\sqrt{\alpha}}$ for $r=2\cdot 10^{3}$. ∎

Fantope optimization in nearly linear time

In this section, we will design a solver for solving the following class of generalized packing/covering SDPs that we will need to solve in the course of our algorithm.

Find either a dual feasible, $w$, with $\sum_{i=1}^{n}w_{i}\geqslant(1-\epsilon)$ or a dual feasible $(M,W)$, satisfying $\Tr M+\Tr W\leqslant 1+\epsilon$.

In this subsection, we will establish a regret guarantee useful for designing fast solvers for our class of SDPs. First, let $\mathcal{S}$ defined as:

The game takes place over $T$ rounds where for each round $t\in 1,\dots,T$:

The player plays two psd matrices $(M_{t},W_{t})\in\mathcal{S}$.

The environment then reveals two gain matrices $(F_{t},G_{t})$ with $\lVert F_{t}\rVert\leqslant 1$ and $\lVert G_{t}\rVert\leqslant 1$ and the player achieves a gain of $\langle F_{t},M_{t}\rangle+\langle G_{t},W_{t}\rangle$.

The goal of the player is to minimize their total regret:

We will first provide a regret guarantee for the following strategy where in each iteration $(M_{t},W_{t})$ are defined for $\eta>0$ by:

Before we move on to the regret bound, we will require 3.5:

The function $f(M,W)=\mathsf{vNE}(M)+\mathsf{vNE}(W)$ is $1$-strongly concave with respect to the following norm:

We will now state a standard regret guarantee (See, for example, Theorem 5.2 from [Haz19]) for the update rule defined in Equation 13:

For a sequence of gain matrices, $\{(F_{t},G_{t})\}_{t=1}^{T}$ satisfying $\lVert F_{t}\rVert\leqslant 1$ and $\lVert G_{t}\rVert\leqslant 1$, the update rule defined in Equation 13 satisfies:

The lemma follows immediately from Theorem 5.2 in [Haz19]. ∎

We will use the following corollary in the analysis of our solver:

Let $\{(F_{t},G_{t})\}_{t=1}^{T}$ be any sequence of gain matrices satisfying $\lVert F_{t}\rVert\leqslant 1$ and $\lVert G_{t}\rVert\leqslant 1$ and let $(M_{t},W_{t})$ be defined as in Equation 13 and suppose that $(M^{(t)}{},W^{(t)}{})$ satisfy:

The corollary follows from the fact that for each $1\leqslant t\leqslant T$, we have:

where the first inequality follows from Matrix-Hölders inequality. ∎

2 Analysis of the Solver

In this subsection, we formally introduce our solver and incorporate the regret analysis from the previous subsection into its analysis. We first introduce the following notation:

Our algorithm and its subsequent analysis follow along the lines of [PTZ12]:

For the rest of the proof, we will assume that the algorithm terminates at the end of the $T^{th}$ loop for some $T\leqslant R$. For ease of exposition, we now define the following variables:

Note that $M^{(t)}{}$ and $W^{(t)}{}$ are meant to be approximations to $\widetilde{M}^{(t)}{}$ and $\widetilde{W}^{(t)}{}$ respectively and the correctness of these projections is guaranteed by Theorem 7.20. Also, observe that $\omega^{(t)}{}=\varepsilon^{\dagger}\sum_{i=1}^{t-1}F^{(i)}$ and $\theta^{(t)}{}=\varepsilon^{\dagger}\sum_{i=1}^{t-1}G^{(i)}$. In the next few lemmas proving the correctness of Algorithm 3, we will simplify presentation by making the following assumptions. We will prove in the main theorem of the section that these assumptions hold with the desired probability.

We assume the following about the running of Algorithm 3 for all $t\in[T]$:

The projections $(M^{(t)}{},W^{(t)}{})$ satisfy $\lVert W^{(t)}{}\rVert\leqslant\Tr(W^{(t)}{})/k$ and satisfy:

The estimates, $y^{(t)}{}$ and $z^{(t)}{}$ satisfy for all $i\in[n]$:

We also make the following non-probabilistic assumptions about the problem:

We assume that the problem instance satisfies:

The following three claims are analogues of Claims 3.3-3.5 from [PTZ12]:

Assume 5.5 and 5.6. Then, for $t=1,\dots,T$:

We have from the definition of $S^{(t)}$ and 5.5:

Assume 5.5 and 5.6. Then, for $t=0,\dots,T$:

It suffices to prove the claim for $t=T$ as for $t<T$, the claim is true from the fact that the while loop continued till the next iteration. Now, we have from the fact that $\xi^{(t)}_{i}\leqslant\alpha w^{(t-1)}_{i}$:

We start with the following decomposition of $\psi^{(t)}{}$:

Under 5.5 and 5.6, we have for every $t=0,\dots,T$:

As in the proof of Lemma 3.2 in [PTZ12], we will prove the claim via strong induction on $t$. We have from the definitions of $\alpha$ and $\varepsilon^{\dagger}$:

We can now apply the results of Corollary 5.4 and the definition of $S^{(t)}{}$ along with 5.5 to obtain:

From the previous two inequalities, we get for $\psi^{(t)}{}$ from Equation 14:

Finally, we get for $\phi^{(t)}{}$ from Equation 15:

Under 5.5 and 5.6, Algorithm 3 terminates with $\lVert w^{(R)}\rVert_{1}\leqslant K$, we have for all $i\in[n]$:

Suppose for the sake of contradiction, that there exists $i\in[n]$ such that:

Now, let $U$ denote the steps in algorithm where the dual variable, $w^{(t)}_{i}$ was incremented. From 5.5, we get that $w^{(t)}_{i}$ is at least incremented for every iteration in the set $Y$ defined as:

By Markov’s inequality and the definition of $\varepsilon^{\dagger}$, we must have $\lvert Y\rvert\geqslant\frac{\varepsilon}{2(1+\varepsilon)}\cdot R$. We must have as $w^{(t)}_{i}$ is incremented by a factor of $(1+\alpha)$ each time:

which is a contradiction. This concludes the proof of the lemma. ∎

Assume 5.6. Then, 5.5 holds in the running of Algorithm 3 with probability at least $1-\delta$. Furthermore, the total runtime of Algorithm 3 is at most:

where $t_{C_{i}}$ and $t_{D_{i}}$ denote the time taken to compute one matrix vector multiplication with $C_{i}$ and $D_{i}$ respectively, $t_{C}=\sum_{i=1}^{n}t_{C_{i}}$ and $t_{D}=\sum_{i=1}^{n}t_{D_{i}}$.

We will prove that 5.5 hold by induction on the number of steps of the Algorithm. Our induction hypothesis will be that 5.5 hold with probability $\delta^{\dagger}(3t)$ up to iteration $t$. The hypothesis is trivially true at $t=0$. Now, we will inductively prove that the assumptions hold true when $t=q+1$ given that they hold at $t=1\dots q$. We start by computing a bound on the matrices $\omega^{(t)}{}$ and $\theta^{(t)}{}$. We have by the application of Lemma 5.10 up to iteration $q$ that:

Therefore, the upper bounds computed on $\lVert\omega^{(t)}\rVert$ and $\lVert\omega^{(t)}\rVert$ remain valid even in iteration $q+1$. Therefore, conditioned on 5.5 holding true for iteration $q$, the conclusions of Theorems 7.20 and B.4 hold for Algorithms $\mathsf{FantopeProjection}$ and $\mathsf{InnerProductEstimation}$ for iteration $q+1$ with probability at least $1-3\delta^{\dagger}$. Hence, 5.6 hold for iteration $q+1$ with probability at least $(1-3\delta^{\dagger})(1-3\delta^{\dagger}q)\geqslant 1-(3\delta^{\dagger})(q+1)$.

The runtime guarantees follow from the runtime guarantees in Lemmas B.4 and 7.20 along with the fact that $2Kk$ is $O(\text{\rm poly}(\frac{1}{\varepsilon},\log(l+m+n),k))$ and matrix-vector multiplies with $\omega^{(t)}{}$ and $\theta^{(t)}{}$ can be implemented in time $O(t_{C})$ and $O(t_{D})$ respectively. And furthermore, a matrix vector product for all the $C_{i}$ and $\mathcal{P}_{V^{(t)}{}}^{\perp}D_{i}$ required by Lemma B.4 can be implemented in time $O(t_{C})$ and $O(mk+t_{D})$ respectively as for any vector $v$, computing $v^{\top}\mathcal{P}_{V^{(t)}{}}^{\perp}$ takes $O(mk)$ time and subsequently, the resultant is multiplied with each of the $D_{i}$. ∎

We now conclude with the main theorem of the section.

There exists an Algorithm, $\mathsf{PackingCoveringDecision}$, which when given an instance of 5.1, with $A_{i}=C_{i}C_{i}^{\top}$, $B_{i}=D_{i}D_{i}^{\top}$, error tolerance $\varepsilon\geqslant\frac{1}{n^{2}}$ and failure probability $\delta$, runs in time:

where $t_{C_{i}}$ and $t_{D_{i}}$ are the time taken to perform a matrix-vector product with $C_{i}$ and $D_{i}$ respectively and $t_{C}=\sum_{i=1}^{n}t_{C_{i}}$ and $t_{D}=\sum_{i=1}^{n}t_{D_{i}}$, and outputs a correct answer to 5.1 with probability at least $1-\delta$.

We first discard $A_{i}$ and $B_{i}$ for those indices $i$ satisfying,

We will now run Algorithm 3 instantiated with error parameter set to $\varepsilon/20$ and failure probability $\delta$. We first quickly address the case where $T=0$. In this case, it must be that $\lVert w^{(0)}\rVert_{1}>\frac{20(1+\log(n+m+l))}{\varepsilon}$. In this case, $w^{*}$ returned by the algorithm satisfies by definition $\sum_{i=1}^{n}w_{i}\geqslant(1-\varepsilon/2)$ and furthermore, by 5.7, is a valid dual solution. In this case, we can simply output $\hat{w}=w^{*}$ as a valid answer to 5.1.

Now, after discarding the above two cases, we have that 5.6 hold for the input passed to Algorithm 3. We have from Lemma 5.12 that Algorithm 3 runs in time:

and that 5.5 hold in the running of Algorithm 3 with probability at least $1-\delta$. Conditioned on this event, we consider two possible cases:

The algorithm returns a dual solution, $w^{\ast}$.

The algorithm returns a primal solution, $(M^{\ast},W^{\ast})$.

In the first case, we have by Lemma 5.10 and the definition of $w^{\ast}$ that $w^{\ast}$ is a feasible dual solution and furthermore, that $\sum_{i=1}^{n}w^{\ast}\geqslant 1-\varepsilon/2$ from our setting of the arguments to Algorithm 3. In this case, we simply define $\hat{w}=w^{\ast}$ for indices that are included in the input to Algorithm 3 and $0$ for the discarded indices. Clearly, $\hat{w}$ is feasible dual solution to the original $\varepsilon$-decision problem.

In the second case, we construct a new primal solution, $(\widehat{M},\widehat{W})=(M^{*}+I/(n+l+m)^{5},W^{*}/(n+l+m)^{5})$. Note that for our bounds on $\varepsilon$ and $n$, the trace of $(\widehat{M},\widehat{W})$ from 5.5 is at most:

and furthermore, from Lemma 5.11, $(\widehat{M},\widehat{W})$ satisfies all the primal constraints for the indices passed to Algorithm 3 and finally for any discarded index, $i$, we have:

Furthermore, from 5.5 since $W^{*}$ satisfied $\lVert W^{*}\rVert\leqslant\Tr W^{*}/k$, we have

Therefore, $(\hat{M},\hat{W})$ is a valid primal solution to the original $\varepsilon$-decision problem. Now, the run time guarantees follow from the fact that the run-time is dominated by the running of Algorithm 3 and the probabilistic guarantees follow from the fact that Algorithm 3 runs correctly with probability at least $1-\delta$ as established previously. ∎

Power Method Analysis

Thus, we analyze the left hand side of the above expression. A short calculation reveals that

So far we have not used the fact that $\widetilde{\phi}_{1}$ is the output of Algorithm 5. In particular, our progress so far which is recorded in Eq. 17 holds true for arbitrary $\widetilde{\phi}_{1}$. We now discuss and prove the relevant properties of $\widetilde{\phi}_{1}$ we use for showing Lemma 6.2 (more specifically, for showing Section 6).

Recall that $\phi_{1},\dots,\phi_{d}$ is the basis of eigenvectors of $A$.

$\langle\bm{g},\phi_{1}\rangle$ is distributed as a scalar standard Gaussian random variable and hence

The following can be found in [Tao, Theorem 2.1.12]:

From Algorithm 4 $\widetilde{\phi}_{1}$ is equal to $\frac{A^{t}g}{\|A^{t}g\|}$ for $t=\Theta\left(\frac{\log d+\log\frac{1}{\delta}+\log\frac{1}{\varepsilon}}{\varepsilon}\right)$. Additionally, from Proposition 6.4, Lemma 6.5 $g$ is a $\delta$-tempered vector and

$\|\Pi_{[0,(1-\varepsilon)\lambda_{1}]}\widetilde{\phi}_{1}\|\leqslant\varepsilon$.

We start by expressing $g$ in the basis $\{\phi_{1},\dots,\phi_{d}\}$ as

Since $t=L\frac{\log d+\log\frac{1}{\delta}+\log\frac{1}{\varepsilon}}{\varepsilon}$, we can choose constant $L$ large enough so that the above is bounded by $\varepsilon^{2}$. ∎

Let $S_{\varepsilon}=\{i:\lambda_{i}\geqslant(1-\varepsilon)\lambda_{1}\}$ and $T=\{i:\lambda_{i}>0\}$. We now establish the following result:

First fix one particular $i\in T\setminus S_{\varepsilon}$ and a in the proof of Proposition 6.6:

where the first inequality follows from the fact that $\lvert c_{i}\rvert\leqslant\left\lvert\frac{\lambda_{i}^{t}\hat{g}_{i}}{\lambda_{1}^{t}\hat{g}_{1}}\right\rvert$, the second inequality follows from the assumption that $\hat{g}$ is $g$-tempered and has a bounded norm and the final inequality from our definition of $t$. By summing up over all the terms, the statement of the proposition follows. ∎

$(1-2\varepsilon)\lambda_{1}\leqslant\widetilde{\phi}_{1}^{\top}A\widetilde{\phi}_{1}\leqslant\lambda_{1}$.

The upper bound follows from $\lambda_{1}$ being the maximum eigenvalue of $A$. As a consequence of Proposition 6.6 and the fact that $\|\widetilde{\phi}_{1}\|=1$,

The same proof as Proposition 6.8 also shows that:

2 Wrapup and proof of Theorem 6.1

In this section we will first prove Lemma 6.2 and then prove Theorem 6.1.

As in Proposition 6.7, we will define the sets $S_{\varepsilon}\coloneqq\{i:\lambda_{i}\geqslant(1-\varepsilon)\lambda_{1}\}$ and $T=\{i:\lambda_{i}>0\}$. Recall that it suffices to prove:

Note that we may write $\widetilde{\phi}_{1}=\sum_{i\in T}c_{i}\phi_{i}$ as we run at least one iteration of the power method which ensures that $\widetilde{\phi}_{1}$ is in the row/column space of $A$. Therefore, we can assume that the sums in Eq. 17 only go over the elements in $T$. Using this as our starting point, we have:

We start by bounding the first term in Eq. 18:

where the first inequality follows from Proposition 6.8 and the definition of the set $S_{\varepsilon}$, the second inequality follows from Cauchy-Schwarz, the third follows again from the definition of the set $S_{\varepsilon}$ and the final inequality from the fact that $\sum_{i}c_{i}^{2}=1$.

where the second inequality follows from Cauchy-Schwarz and the final inequality follows from the definition of $S_{\varepsilon}$ and Proposition 6.7.

and the proof proceeds as before. For the final term, we have from Cauchy-Schwarz and Proposition 6.7:

Putting the bounds on the four terms on Eq. 18, we get the desired result. ∎

to obtain $\widetilde{\phi}_{m}$. Let us define:

since $\widetilde{\phi}_{m}$ is orthogonal to the space spanned by $\widetilde{\phi}_{1},\dots,\widetilde{\phi}_{m-1}$, and further note that

to Section 6.2 along with the definitions of $\widetilde{A}$ and $F$ gives us

Hence our induction is complete and our goal statement is proved. ∎

Fast Projection on Fantopes

In this section we define $\mathcal{S}$ as follows

We will be concerned with solving the following optimization problem:

In this subsection, we will prove that Algorithm 6 correctly computes the optimizer to (20). To do this, we will first analyze the following simpler problem:

Henceforth, we use $p(G)$ to denote $W^{*}$.

We first prove that the optimizer, $W^{*}$, of (21) has the same eigenvectors as that of $G$.

Given $G\succcurlyeq 0$, the optimizer, $W^{*}$, of (21) has the same eigenvectors as $G$.

Let $\sigma_{1}\geqslant\dots\geqslant\sigma_{m}\geqslant 0$ and $\lambda_{1}\geqslant\dots\lambda_{m}\geqslant 0$ denote the eigenvalues of $W^{*}$ and $G$ respectively. Now, we have the von Neumann’s trace inequality:

with equality when the eigenvectors of $W^{*}$ corresponding to the eigenvalue $\sigma_{i}$ coincide with the eigenvectors of $G$ for the eigenvalue $\lambda_{i}$. Therefore, the optimizer $W^{*}$ must share the same set of eigenvectors as $G$. ∎

Given, $G\succcurlyeq 0$ and let $H=\exp G$ with eigenvalue decomposition $H=\sum_{i=1}^{m}\lambda_{i}u_{i}u_{i}^{\top}$. Let $\nu^{*}$ be defined as follows:

Then, the optimizer, $W^{*}$, of (21) is given by:

From Lemma 7.2, we know that the eigenvectors for $W^{*}$ and $G$ and hence, $H$, coincide. Let $\sigma_{1}^{*},\dots,\sigma_{m}^{*}$ denote the eigenvalues of $W^{*}$ corresponding to the eigenvectors $u_{1},\dots,u_{m}$. Then, we see from (21) that:

Since, the above optimization problem is convex, we compute its Lagrangian (Note that we must set $\alpha_{i}\geqslant 0$):

Now, picking $\beta^{\prime}<-\max_{i\in[m]}\lvert\log\lambda_{i}\rvert$. Note that $\log\lambda_{i}$ are the eigenvalues of $G$ and hence $\beta$ is finite. We now have:

by noting that $-x\log x$ is maximized at $x=1/e$. Note that the above conclusion holds true for $\alpha$ satisfying $\max_{i}\alpha_{i}\leqslant-\max_{i\in[m]}\lvert\log\lambda_{i}\rvert-\beta^{\prime}$. Therefore, Slaters’ condition holds for both the primal problem, Prog, and its dual. Furthermore, the optimal value of Prog is bounded as both Prog and its dual have a feasible point with finite objective value. Therefore, strong duality holds for Prog and its dual and their optimal value is attained. Let $\{\sigma_{i}^{*}\}$ and $(\{\alpha_{i}^{*}\geqslant 0\},\beta^{*})$ denote the primal and dual optimal points respectively. Note that by a simple exchange argument $\sigma^{*}_{i}\neq 0$. Therefore, the KKT conditions apply to Prog and we get:

From the condition of primal feasibility, we get that $e^{\beta^{*}}=(\sum_{i=1}^{m}e^{-\alpha^{*}_{i}}\lambda_{i})^{-1}$. Also, note that we get from complementary slackness that $\alpha^{*}_{i}>0$ implies that $\sigma^{*}_{i}=1/k$. Additionally, from complementary slackness, we obtain that $\sigma^{*}_{i}\geqslant\sigma^{*}_{j}$ for $i\geqslant j$.

Let $l=\#\{i:\alpha^{*}_{i}>0\}$. We first tackle the case where $l=0$. In this case, the optimizer is simply $W^{*}=H/\Tr(H)$ and the statement of the lemma is true.

Now assume that $l>0$. Let us now consider the function, $f$, defined as:

When $\lambda_{m}>0$ which holds in this case, $f(\nu)$ is a strictly increasing, continuous function of $\nu$ in the interval $[\lambda_{m},\infty)$ and its value increases from $1/m$ to $\infty$. For $i,j\in[l]$, we have $\sigma^{*}_{i}=\sigma^{*}_{j}=1/k$ by complementary slackness and therefore $e^{-\alpha^{*}_{i}}\lambda_{i}=e^{-\alpha^{*}_{j}}\lambda_{j}=\widehat{\nu}$. For $i\notin[l]$, we have $\sigma^{*}_{i}=e^{\beta^{*}}\min(\lambda_{i},\widehat{\nu})$ as we have $\sigma^{*}_{i}=e^{\beta^{*}}\lambda_{i}\leqslant e^{\beta^{*}}\widehat{\nu}=1/k$. From the previous two statements, we have $\sigma^{*}_{i}=e^{\beta^{*}}\min(\widehat{\nu},\lambda_{i})$ for all $i\in[m]$. Finally, we have $f(\widehat{\nu})=1/k$ from complementary slackness which implies that $\widehat{\nu}=\nu^{*}$ as $f$ is strictly increasing and continuous. Which implies that the optimal value of $\sigma^{*}_{i}$ is given by $\sigma^{*}_{i}=\min(\lambda_{i},\nu^{*})/(\sum_{j=1}^{m}\min(\lambda_{j},\nu^{*}))$, thus proving the lemma.

Finally, we will now show how to use solutions to (21) to obtain solutions to the following:

The result is detailed in the following lemma:

Let $F,G\succcurlyeq 0$ and let $Q=\exp(F)$, $Z_{1}=\Tr(\exp(F))$, $H=\exp(G)$ with eigenvalue decomposition $H=\sum_{i=1}^{m}\lambda_{i}u_{i}u_{i}^{\top}$ and $\nu^{*}$ and $Z_{2}$ be defined as:

Then, the optimizers, $(M^{*},W^{*})$, of Equation (22) are given by:

where $\gamma=\log Z_{1}$, $\zeta=\log(Z_{2})+k^{-1}\sum_{i=1}^{k}(\log(\lambda_{i})-\log\min(\lambda_{i},\nu^{*}))$ and $\widehat{W}$ and $\widehat{M}$ are defined as:

Let $(M^{*},W^{*})$ denote the solutions of (22) and let $\alpha=\Tr{M^{*}}$. Then, we must have:

Now consider the case where $\alpha>0$. For the first equation, we have:

When $\alpha=0$, the conclusion of the previous manipulation is trivially true. By a similar manipulation, from Lemma 7.3 we have $W^{*}=(1-\alpha)\widehat{W}$. We now have:

We now proceed for a similar computation for $\widehat{W}$:

where the second-to-last equality follows because at most $k$ of the $\lambda_{i}$ are greater than $\nu^{*}$ and the final inequality follows from the fact that $\lambda_{i}\geqslant\nu^{*}$ implies that $\min(\lambda_{i},\nu^{*})=\nu^{*}/Z_{2}=1/k$. By putting the previous two results together, we get that:

whose optimal value is given by $\alpha=e^{\gamma}/(e^{\gamma}+e^{\zeta})$ which concludes the proof the lemma. ∎

It is unclear how to exactly solve the optimization problem (20) fast, so we give an algorithm that outputs a solution close to the exact optimizer in trace norm. As a first step, we give an algorithm to approximately solve the optimization problem (21). In the algorithm below, not all matrices are explicitly computed and we obtain an implicit representation of $W^{*}$ rather than an explicit $m\times m$ matrix. For simplicity of exposition, we defer the details of this implicit representation to later subsections. In the algorithm below, when we say $\Tr_{\varepsilon}(H)$, we mean running the trace estimation algorithm from Corollary B.5 from $H$, which with high probability produces a $(1\pm\varepsilon)$-approximation of the trace.

Our first goal is to show that the output of Algorithm 7 on input $G$ is close in trace norm to $p(G)$ where $p$ is as defined in Remark 7.1. Concretely, we prove:

Let $G$ be a positive semidefinite matrix, and let $W^{*}$ be the output of Algorithm 7 on input $G$. Then:

The full statement of the above, which states some more technical properties of $W^{*}$ can be found in Theorem 7.15.

2 Closeness in trace norm I: projections of spectrally similar matrices

In this section, let $A$ and $\widetilde{A}$ be positive semidefinite matrices such that

for some $0<\varepsilon<1/2$. $\widetilde{A}$ will be a matrix obtained via the power iteration based PCA algorithm

$\|\log\widetilde{A}-\log A\|\leqslant 4\varepsilon$.

where the last inequality follows from $\varepsilon<1/2$. ∎

For the rest of this section, let $\displaystyle M\coloneqq\arg\min_{\begin{subarray}{c}X\succcurlyeq 0\\ \Tr(X)=1\\ \|X\|\leqslant 1/k\end{subarray}}\mathsf{QRE}(X,A)$ and let $\displaystyle\widetilde{M}\coloneqq\arg\min_{\begin{subarray}{c}X\succcurlyeq 0\\ \Tr(X)=1\\ \|X\|\leqslant 1/k\end{subarray}}\mathsf{QRE}(X,\widetilde{A})$. In the language of Remark 7.1, $M=p(\log A)$.