We present the first polynomial-time algorithm for computing a near-optimal flow-expander decomposition. Given a graph G and a parameter ϕ, our algorithm removes at most a ϕlog^{1+o(1)}n fraction of edges so that every remaining connected component is a ϕ-flow-expander (a stronger guarantee than being a ϕ-cut-expander). This achieves overhead log^{1+o(1)}n, nearly matching the Ω(log n) graph-theoretic lower bound that already holds for cut-expander decompositions, up to a log^{o(1)}n factor. Prior polynomial-time algorithms required removing O(ϕlog^{1.5}n) and O(ϕlog²n) fractions of edges to guarantee ϕ-cut-expander and ϕ-flow-expander components, respectively.

Placing a dataset A = {ai}i∈[n]⊂ ℝdin radial isotropic position, i.e., finding an invertible R ∈ ℝd×dsuch that the unit vectors ${\left\{ {\left( {{\text{R}}{{\text{a}}_i}} \right)\left\| {{\text{R}}{{\text{a}}_i}} \right\|_2^{ - 1}} \right\}_{i \in [n]}}$ are in isotropic position, is a powerful tool with applications in functional analysis, communication complexity, coding theory, and the design of learning algorithms. When the transformed dataset has a second moment matrix within a exp(±ϵ) factor of a multiple of Id, we call R an ϵ-approximate Forster transform.We give a faster algorithm for computing approximate Forster transforms, based on optimizing an objective defined by Barthe [1]. When the transform has a polynomially-bounded aspect ratio, our algorithm uses $O\left( {n{d^{\omega - 1}}{{\left( {\frac{n}{\varepsilon }} \right)}^{o(1)}}} \right)$ time to output an ϵ-approximate Forster transform with high probability, when one exists. This is almost the natural limit of this approach, as even evaluating Barthe’s objective takes O(ndω−1) time. Previously, the state-of-the-art runtime in this regime was based on cutting-plane methods, and scaled at least as ≈ n3+n2dω−1. We also provide explicit estimates on the aspect ratio in the smoothed analysis setting, and show that our algorithm similarly improves upon those in the literature.To obtain our results, we develop a subroutine of potential broader interest: a reduction from almost-linear time sparsification of graph Laplacians to the ability to support almost-linear time matrix-vector products. We combine this tool with new stability bounds on Barthe’s objective to implicitly implement a box-constrained Newton’s method [2], [3].

We provide a general method to convert a "primal" black-box algorithm for solving regularized convex-concave minimax optimization problems into an algorithm for solving the associated dual maximin optimization problem. Our method adds recursive regularization over a logarithmic number of rounds where each round consists of an approximate regularized primal optimization followed by the computation of a dual best response. We apply this result to obtain new state-of-the-art runtimes for solving matrix games in specific parameter regimes, obtain improved query complexity for solving the dual of the CVaR distributionally robust optimization (DRO) problem, and recover the optimal query complexity for finding a stationary point of a convex function.

We provide an algorithm that, given an n-vertex m-edge Eulerian graph with polynomially bounded weights, computes an (n log2 n · ∈-2)-edge ε-approximate Eulerian sparsifier with high probability in (m log3 n ) time (where (·) hides polyloglog(n ) factors). Due to a reduction from [Peng-Song, STOC ’22], this yields an (m log3 n + n log6 n )-time algorithm for solving n-vertex m-edge Eulerian Laplacian systems with polynomially-bounded weights with high probability, improving upon the previous state- of-the-art runtime of Ω(m log8 n + n log23 n ). We also give a polynomial-time algorithm that computes O (min(n log n · ε-2 + n log5/3 n log3/2 n · ε-2))-edge sparsifiers, improving the best such sparsity bound of O (n log2 n · ε-2 + n log8/3 n · ε-4/3) [Sachdeva-Thudi-Zhao, ICALP ’24].

In this article, we provide an O ( m loglog O (1) n log (1/ε))-expected time algorithm for solving Laplacian systems on n -node m -edge graphs, improving upon the previous best expected runtime of \(O(m \sqrt {\log n} \mathrm{log log}^{O(1)} n \log (1/\epsilon))\) achieved by Cohen et al. [ 12 ]. To obtain this result, we provide efficient constructions of low spectral stretch graph approximations with improved stretch and sparsity bounds. As motivation for this work, we show that, for every set of vectors in \(\mathbb {R}^d\) (not just those induced by graphs) and all integer k > 1, there exist an ultra-sparsifier with d -1 + O ( d/k ) re-weighted vectors of relative condition number at most k 2 . For small k , this improves upon the previous best known multiplicative factor of \(k \cdot \tilde{O}(\log d)\) , which is only known for the graph case. Additionally, in the graph case, we employ our low-stretch subgraph construction to obtain n - 1 + O ( n/k )-edge ultrasparsifiers of relative condition number k 1+o(1) for k = ω (log δ n ) for any δ > 0: this improves upon the previous work for k = o (exp (log 1/2-δ n )).

The $k$-principal component analysis ($k$-PCA) problem is a fundamental algorithmic primitive that is widely-used in data analysis and dimensionality reduction applications. In statistical settings, the goal of $k$-PCA is to identify a top eigenspace of the covariance matrix of a distribution, which we only have black-box access to via samples. Motivated by these settings, we analyze black-box deflation methods as a framework for designing $k$-PCA algorithms, where we model access to the unknown target matrix via a black-box $1$-PCA oracle which returns an approximate top eigenvector, under two popular notions of approximation. Despite being arguably the most natural reduction-based approach to $k$-PCA algorithm design, such black-box methods, which recursively call a $1$-PCA oracle $k$ times, were previously poorly-understood. Our main contribution is significantly sharper bounds on the approximation parameter degradation of deflation methods for $k$-PCA. For a quadratic form notion of approximation we term ePCA (energy PCA), we show deflation methods suffer no parameter loss. For an alternative well-studied approximation notion we term cPCA (correlation PCA), we tightly characterize the parameter regimes where deflation methods are feasible. Moreover, we show that in all feasible regimes, $k$-cPCA deflation algorithms suffer no asymptotic parameter loss for any constant $k$. We apply our framework to obtain state-of-the-art $k$-PCA algorithms robust to dataset contamination, improving prior work in sample complexity by a $\mathsf{poly}(k)$ factor.

We develop a new parallel algorithm for minimizing Lipschitz, convex functions with a stochastic subgradient oracle. The total number of queries made and the query depth, i.e., the number of parallel rounds of queries, match the prior state-of-the-art, [CJJLLST23], while improving upon the computational depth by a polynomial factor for sufficiently small accuracy. When combined with previous state-of-the-art methods our result closes a gap between the best-known query depth and the best-known computational depth of parallel algorithms. Our method starts with a \emph{ball acceleration} framework of previous parallel methods, i.e., [CJJJLST20, ACJJS21], which reduce the problem to minimizing a regularized Gaussian convolution of the function constrained to Euclidean balls. By developing and leveraging new stability properties of the Hessian of this induced function, we depart from prior parallel algorithms and reduce these ball-constrained optimization problems to stochastic unconstrained quadratic minimization problems. Although we are unable to prove concentration of the asymmetric matrices that we use to approximate this Hessian, we nevertheless develop an efficient parallel method for solving these quadratics. Interestingly, our algorithms can be improved using fast matrix multiplication and run in nearly-linear time if the matrix multiplication exponent is 2.

In the recent literature on machine learning and decision making, calibration has emerged as a desirable and widely-studied statistical property of the outputs of binary prediction models. However, the algorithmic aspects of measuring model calibration have remained relatively less well-explored. Motivated by Blasiok et al '23, which proposed a rigorous framework for measuring distances to calibration, we initiate the algorithmic study of calibration through the lens of property testing. We define the problem of calibration testing from samples where given $n$ draws from a distribution $\mathcal{D}$ on $(\text{predictions}, \text{binary outcomes})$, our goal is to distinguish between the cases where $\mathcal{D}$ is perfectly calibrated or $\epsilon$-far from calibration. We make the simple observation that the empirical smooth calibration linear program can be reformulated as an instance of minimum-cost flow on a highly-structured graph, and design an exact dynamic programming-based solver for it which runs in time $O(n\log^2(n))$, and solves the calibration testing problem information-theoretically optimally in the same time. This improves upon state-of-the-art black-box linear program solvers requiring $\Omega(n^\omega)$ time, where $\omega > 2$ is the exponent of matrix multiplication. We also develop algorithms for tolerant variants of our testing problem improving upon black-box linear program solvers, and give sample complexity lower bounds for alternative calibration measures to the one considered in this work. Finally, we present experiments showing the testing problem we define faithfully captures standard notions of calibration, and that our algorithms scale efficiently to accommodate large sample sizes.

We design algorithms for minimizing maxi∈[n] fi(x) over a d-dimensional Euclidean or simplex domain. When each fi is 1-Lipschitz and 1-smooth, our method computes an ɛ-approximate solution using Õ(nɛ-1/3 + ɛ-2) gradient and function evaluations, and Õ(nɛ-4/3) additional runtime. For large n, our evaluation complexity is optimal up to polylogarithmic factors. In the special case where each fi is linear—which corresponds to finding a near-optimal primal strategy in a matrix game—our method finds an ɛ-approximate solution in runtime Õ(n(d/ɛ)2/3 +nd+dɛ-2). For n > d and this improves over all existing first-order methods. When additionally d = ω(n8/11) our runtime also improves over all known interior point methods.

Discrepancy theory has provided powerful tools for producing higher-quality objects which “beat the union bound” in fundamental settings throughout combinatorics and computer science. However, this quality has often come at the price of more computationally-expensive algorithms. We introduce a new framework for bridging this gap, by allowing for the efficient implementation of discrepancy-theoretic primitives. Our framework repeatedly solves regularized optimization problems to low accuracy to approximate the partial coloring method of [Rot17], and simplifies and generalizes recent work of [JSS23] on fast algorithms for Spencer's theorem. In particular, our framework only requires that the discrepancy body of interest has exponentially large Gaussian measure and is expressible as a sublevel set of a symmetric, convex function. We combine this framework with new tools for proving Gaussian measure lower bounds to give improved algorithms for a variety of sparsification and coloring problems.

We consider the sparsification of sums F : ℝn → ℝ+ where F(x) = f1(⟨ a1,x⟩) + ⋯ + fm(⟨ am,x⟩) for vectors a1,…,am ∈ ℝn and functions f1,…,fm : ℝ → ℝ+. We show that (1+ε)-approximate sparsifiers of F with support size n/ε2 (logn/ε)O(1) exist whenever the functions f1,…,fm are symmetric, monotone, and satisfy natural growth bounds. Additionally, we give efficient algorithms to compute such a sparsifier assuming each fi can be evaluated efficiently. Our results generalize the classical case of ℓp sparsification, where fi(z) = |z|p, for p ∈ (0, 2], and give the first near-linear size sparsifiers in the well-studied setting of the Huber loss function and its generalizations, e.g., fi(z) = min{|z|p, |z|2} for 0 < p ≤ 2. Our sparsification algorithm can be applied to give near-optimal reductions for optimizing a variety of generalized linear models including ℓp regression for p ∈ (1, 2] to high accuracy, via solving (logn)O(1) sparse regression instances with m ≤ n(logn)O(1), plus runtime proportional to the number of nonzero entries in the vectors a1, …, am.

We introduce a new tool for stochastic convex optimization (SCO): a Reweighted Stochastic Query (ReSQue) estimator for the gradient of a function convolved with a (Gaussian) probability density. Combining ReSQue with recent advances in ball oracle acceleration [CJJ+20], [ACJ+21], we develop algorithms achieving state-of-the-art complexities for SCO in parallel and private settings. For a SCO objective constrained to the unit ball in $\mathbb{R}^{d}$, we obtain the following results (up to polylogarithmic factors).1)We give a parallel algorithm obtaining optimization error $\epsilon_{\text {opt} }$ with $d^{1 / 3} \epsilon_{\text {opt} }^{-2 / 3}$ gradient oracle query depth and $d^{1 / 3} \epsilon_{\text {opt} }^{-2 / 3}+\epsilon_{\text {opt} }^{-2}$ gradient queries in total, assuming access to a bounded-variance stochastic gradient estimator. For $\epsilon_{\text {opt} } \in\left[d^{-1}, d^{-1 / 4}\right]$, our algorithm matches the state-of-the-art oracle depth of [BJL+19] while maintaining the optimal total work of stochastic gradient descent.2)Given n samples of Lipschitz loss functions, prior works [BFTT19], [BFGT20], [AFKT21], [KLL21] established that if $n \gt rsim d \epsilon_{\mathrm{dp}}^{-2},\left(\epsilon_{\mathrm{dp}}, \delta\right)$-differential privacy is attained at no asymptotic cost to the SCO utility. However, these prior works all required a superlinear number of gradient queries. We close this gap for sufficiently large $n \gt rsim d^{2} \epsilon_{{\bf d p}}^{-3}$, by using ReSQue to design an algorithm with near-linear gradient query complexity in this regime.

Abstract-For any norms $N_{1}, \ldots, N_{m}$ on $\mathbb{R}^{n}$ and $N(x):= N_{1}(x)+\cdots+N_{m}(x)$, we show there is a sparsified norm $\tilde{N}(x)= w_{1} N_{1}(x)+\cdots+w_{m} N_{m}(x)$ such that $|N(x)-\tilde{N}(x)| \leqslant \varepsilon N(x)$ for all $x \in \mathbb{R}^{n}$, where $w_{1}, \ldots, w_{m}$ are non-negative weights, of which only $O\left(\varepsilon^{-2} n \log (n / \varepsilon)(\log n)^{2.5}\right)$ are non-zero. Additionally, we show that such weights can be found with high probability in time $O\left(m(\log n)^{O(1)}+\right.$ poly $\left.(n)\right) T$, where T is the time required to evaluate a norm $N_{i}(x)$, assuming that $N(x)$ is poly $(n)$ equivalent to the Euclidean norm. This immediately yields analogous statements for sparsifying sums of symmetric submodular functions. More generally, we show how to sparsify sums of p th powers of norms when the sum is p-uniformly smooth.1

We develop a general framework for finding approximately-optimal preconditioners for solving linear systems. Leveraging this framework we obtain improved runtimes for fundamental preconditioning and linear system solving problems including: Diagonal preconditioning. We give an algorithm which, given positive definite $\mathbf{K} \in \mathbb{R}^{d \times d}$ with $\mathrm{nnz}(\mathbf{K})$ nonzero entries, computes an $\epsilon$-optimal diagonal preconditioner in time $\widetilde{O}(\mathrm{nnz}(\mathbf{K}) \cdot \mathrm{poly}(\kappa^\star,\epsilon^{-1}))$, where $\kappa^\star$ is the optimal condition number of the rescaled matrix. Structured linear systems. We give an algorithm which, given $\mathbf{M} \in \mathbb{R}^{d \times d}$ that is either the pseudoinverse of a graph Laplacian matrix or a constant spectral approximation of one, solves linear systems in $\mathbf{M}$ in $\widetilde{O}(d^2)$ time. Our diagonal preconditioning results improve state-of-the-art runtimes of $\Omega(d^{3.5})$ attained by general-purpose semidefinite programming, and our solvers improve state-of-the-art runtimes of $\Omega(d^{\omega})$ where $\omega > 2.3$ is the current matrix multiplication constant. We attain our results via new algorithms for a class of semidefinite programs (SDPs) we call matrix-dictionary approximation SDPs, which we leverage to solve an associated problem we call matrix-dictionary recovery.

We investigate area convexity [Sherman17], a mysterious tool introduced to tackle optimization problems under the challenging $\ell_\infty$ geometry. We develop a deeper understanding of its relationship with conventional analyses of extragradient methods [Nemirovski04, Nesterov07]. We also give improved solvers for the subproblems required by variants of the [Sherman17] algorithm, designed through the lens of relative smoothness [BBT17, LFN18}. Leveraging these new tools, we give a state-of-the-art first-order algorithm for solving box-simplex games (a primal-dual formulation of $\ell_\infty$ regression) in a $d \times n$ matrix with bounded rows, using $O(\log d \cdot \epsilon^{-1})$ matrix-vector queries. As a consequence, we obtain improved complexities for approximate maximum flow, optimal transport, min-mean-cycle, and other basic combinatorial optimization problems. We also develop a near-linear time algorithm for a matrix generalization of box-simplex games, capturing a family of problems closely related to semidefinite programs recently used as subroutines in robust statistics and numerical linear algebra.

We present an algorithm that given any n-vertex, m-edge, rank r hypergraph constructs a spectral sparsifier with O(n ε−2 logn logr) hyperedges in nearly-linear O(mr) time. This improves in both size and efficiency over a line of work [Bansal-Svensson-Trevisan 2019, Kapralov-Krauthgamer-Tardos-Yoshida 2021] for which the previous best size was O(min{n ε−4 log3 n,nr3 ε−2 logn}) and runtime was O(mr + nO(1)).

Box-simplex games are a family of bilinear minimax objectives which encapsulate graph-structured problems such as maximum flow [She17], optimal transport [JST19], and bipartite matching [AJJ+22]. We develop efficient near-linear time, high-accuracy solvers for regularized variants of these games. Beyond the immediate applications of such solvers for computing Sinkhorn distances, a prominent tool in machine learning, we show that these solvers can be used to obtain improved running times for maintaining a (fractional) $ε$-approximate maximum matching in a dynamic decremental bipartite graph against an adaptive adversary. We give a generic framework which reduces this dynamic matching problem to solving regularized graph-structured optimization problems to high accuracy. Through our reduction framework, our regularized box-simplex game solver implies a new algorithm for dynamic decremental bipartite matching in total time $\tilde{O}(m \cdot ε^{-3})$, from an initial graph with $m$ edges and $n$ nodes. We further show how to use recent advances in flow optimization [CKL+22] to improve our runtime to $m^{1 + o(1)} \cdot ε^{-2}$, thereby demonstrating the versatility of our reduction-based approach. These results improve upon the previous best runtime of $\tilde{O}(m \cdot ε^{-4})$ [BGS20] and illustrate the utility of using regularized optimization problem solvers for designing dynamic algorithms.

The accelerated proximal point method (APPA), also known as "Catalyst", is a well-established reduction from convex optimization to approximate proximal point computation (i.e., regularized minimization). This reduction is conceptually elegant and yields strong convergence rate guarantees. However, these rates feature an extraneous logarithmic term arising from the need to compute each proximal point to high accuracy. In this work, we propose a novel Relaxed Error Criterion for Accelerated Proximal Point (RECAPP) that eliminates the need for high accuracy subproblem solutions. We apply RECAPP to two canonical problems: finite-sum and max-structured minimization. For finite-sum problems, we match the best known complexity, previously obtained by carefully-designed problem-specific algorithms. For minimizing max_y f(x,y) where f is convex in x and strongly-concave in y, we improve on the best known (Catalyst-based) bound by a logarithmic factor.

We develop a variant of the Monteiro-Svaiter (MS) acceleration framework that removes the need to solve an expensive implicit equation at every iteration. Consequently, for any $p\ge 2$ we improve the complexity of convex optimization with Lipschitz $p$th derivative by a logarithmic factor, matching a lower bound. We also introduce an MS subproblem solver that requires no knowledge of problem parameters, and implement it as either a second- or first-order method by solving linear systems or applying MinRes, respectively. On logistic regression problems our method outperforms previous accelerated second-order methods, but under-performs Newton's method; simply iterating our first-order adaptive subproblem solver is competitive with L-BFGS.

We provide Õ(∊–1)-pass semi-streaming algorithms for computing (1–∊)-approximate maximum cardinality matchings in bipartite graphs. Our most efficient methods are deterministic and use optimal, O(n), space, improving upon the space complexity of the previous state-of-the-art Õ(∊–1)-pass algorithm of [AG18]. To obtain our results we provide semi-streaming adaptations of more general continuous optimization tools. Further, we leverage these techniques to obtain improvements for streaming variants of approximate linear programming, optimal transport, exact matching, transshipment, and shortest path problems.

We make several advances broadly related to the maintenance of electrical flows in weighted graphs undergoing dynamic resistance updates, including:

In this paper we obtain improved iteration complexities for solving ℓp regression. We provide methods which given any full-rank A ∈ ℝn × d with n ≥ d, b ∈ ℝn, and p ≥ 2 solve minx ∈ ℝd ||A x − b||p to high precision in time dominated by that of solving Op(dp−2/3p−2) linear systems in A⊤D A for positive diagonal matrices D. This improves upon the previous best iteration complexity of Op(np−2/3p−2) (Adil, Kyng, Peng, Sachdeva 2019). As a corollary, we obtain an O(d1/3є−2/3) iteration complexity for approximate ℓ∞ regression. Further, for q ∈ (1, 2] and dual norm q = p/(p−1) we provide an algorithm that solves ℓq regression in O(dp−2/2p−2) iterations.

We characterize the complexity of minimizing $\max_{i\in[N]} f_i(x)$ for convex, Lipschitz functions $f_1,\ldots, f_N$. For non-smooth functions, existing methods require $O(N\epsilon^{-2})$ queries to a first-order oracle to compute an $\epsilon$-suboptimal point and $\widetilde{O}(N\epsilon^{-1})$ queries if the $f_i$ are $O(1/\epsilon)$-smooth. We develop methods with improved complexity bounds of $\widetilde{O}(N\epsilon^{-2/3} + \epsilon^{-8/3})$ in the non-smooth case and $\widetilde{O}(N\epsilon^{-2/3} + \sqrt{N}\epsilon^{-1})$ in the $O(1/\epsilon)$-smooth case. Our methods consist of a recently proposed ball optimization oracle acceleration algorithm (which we refine) and a careful implementation of said oracle for the softmax function. We also prove an oracle complexity lower bound scaling as $\Omega(N\epsilon^{-2/3})$, showing that our dependence on $N$ is optimal up to polylogarithmic factors.

We develop a new primitive for stochastic optimization: a low-bias, low-cost estimator of the minimizer $x_\star$ of any Lipschitz strongly-convex function $f$. In particular, we use a multilevel Monte-Carlo approach due to Blanchet and Glynn to turn any optimal stochastic gradient method into an estimator of $x_\star$ with bias $\delta$, variance $O(\log(1/\delta))$, and an expected sampling cost of $O(\log(1/\delta))$ stochastic gradient evaluations. As an immediate consequence, we obtain cheap and nearly unbiased gradient estimators for the Moreau envelope of any Lipschitz convex function. We demonstrate the potential of our estimator through four applications. First, we develop a method for minimizing the maximum of $N$ functions, improving on recent results and matching a lower bound up to logarithmic factors. Second and third, we recover state-of-the-art rates for projection-efficient and gradient-efficient optimization using simple algorithms with a transparent analysis. Finally, we show that an improved version of our estimator would yield a nearly linear-time, optimal-utility, differentially-private non-smooth stochastic optimization method.

We study fast algorithms for statistical regression problems under the strong contamination model, where the goal is to approximately optimize a generalized linear model (GLM) given adversarially corrupted samples. Prior works in this line of research were based on the \emph{robust gradient descent} framework of \cite{PrasadSBR20}, a first-order method using biased gradient queries, or the \emph{Sever} framework of \cite{DiakonikolasKK019}, an iterative outlier-removal method calling a stationary point finder. We present nearly-linear time algorithms for robust regression problems with improved runtime or estimation guarantees compared to the state-of-the-art. For the general case of smooth GLMs (e.g.\ logistic regression), we show that the robust gradient descent framework of \cite{PrasadSBR20} can be \emph{accelerated}, and show our algorithm extends to optimizing the Moreau envelopes of Lipschitz GLMs (e.g.\ support vector machines), answering several open questions in the literature. For the well-studied case of robust linear regression, we present an alternative approach obtaining improved estimation rates over prior nearly-linear time algorithms. Interestingly, our algorithm starts with an identifiability proof introduced in the context of the sum-of-squares algorithm of \cite{BakshiP21}, which achieved optimal error rates while requiring large polynomial runtime and sample complexity. We reinterpret their proof within the Sever framework and obtain a dramatically faster and more sample-efficient algorithm under fewer distributional assumptions.

In this paper we provide an O(mloglogO(1) n log(1/∊))-expected time algorithm for solving Laplacian systems on n-node m-edge graphs, improving improving upon the previous best expected runtime of achieved by (Cohen, Kyng, Miller, Pachocki, Peng, Rao, Xu 2014). To obtain this result we provide efficient constructions of ℓp-stretch graph approximations with improved stretch and sparsity bounds. Additionally, as motivation for this work, we show that for every set of vectors in ℝd (not just those induced by graphs) and all k > 0 there exist an ultra-sparsifiers with d – 1 + O(d/k) re-weighted vectors of relative condition number at most k2. For small k, this improves upon the previous best known multiplicative factor of k · Õ(log d), which is only known for the graph case.

Consider an oracle which takes a point x and returns the minimizer of a convex function f in an l2 ball of radius r around x. It is straightforward to show that roughly r^{-1}\log(1/epsilon) calls to the oracle suffice to find an \epsilon-approximate minimizer of f in an l2 unit ball. Perhaps surprisingly, this is not optimal: we design an accelerated algorithm which attains an epsilon-approximate minimizer with roughly r^{-2/3} \log(1/epsilon) oracle queries, and give a matching lower bound. Further, we implement ball optimization oracles for functions with a locally stable Hessian using a variant of Newton's method and, in certain cases, stochastic first-order methods. The resulting algorithms apply to a number of problems of practical and theoretical import, improving upon previous results for logistic and linfinity regression and achieving guarantees comparable to the state-of-the-art for lp regression.

We develop two methods for the following fundamental statistical task: given an $\eps$-corrupted set of $n$ samples from a $d$-dimensional sub-Gaussian distribution, return an approximate top eigenvector of the covariance matrix. Our first robust PCA algorithm runs in polynomial time, returns a $1 - O(\eps\log\eps^{-1})$-approximate top eigenvector, and is based on a simple iterative filtering approach. Our second, which attains a slightly worse approximation factor, runs in nearly-linear time and sample complexity under a mild spectral gap assumption. These are the first polynomial-time algorithms yielding non-trivial information about the covariance of a corrupted sub-Gaussian distribution without requiring additional algebraic structure of moments. As a key technical tool, we develop the first width-independent solvers for Schatten-$p$ norm packing semidefinite programs, giving a $(1 + \eps)$-approximate solution in $O(p\log(\tfrac{nd}{\eps})\eps^{-1})$ input-sparsity time iterations (where $n$, $d$ are problem dimensions).

We give the first approximation algorithm for mixed packing and covering semidefinite programs (SDPs) with polylogarithmic dependence on width. Mixed packing and covering SDPs constitute a fundamental algorithmic primitive with applications in combinatorial optimization, robust learning, and quantum complexity. The current approximate solvers for positive semidefinite programming can handle only pure packing instances, and technical hurdles prevent their generalization to a wider class of positive instances. For a given multiplicative accuracy of є, our algorithm takes O(log3(ndρ) · є−3) parallelizable iterations, where n, d are the problem dimensions and ρ is a width parameter of the instance, generalizing or improving all previous parallel algorithms in the positive linear and semidefinite programming literature. When specialized to pure packing SDPs, our algorithm’s iteration complexity is O(log2 (nd) · є−2), a slight improvement and derandomization of the state-of-the-art due to Allen-Zhu et. al. ’16, Peng et. al. ’16, and Wang et. al. ’15. For several common structured instances, the iterations of our algorithm run in nearly-linear time. In doing so, we give matrix analytic techniques to overcome obstacles that have stymied prior approaches to this problem, as stated in Peng et. al. ’16 and Mahoney et. al. ’16. Crucial to our analysis are a simplification of existing algorithms for mixed positive linear programs, achieved by removing an asymmetry from modifying covering constraints, and a suite of matrix inequalities with proofs based on analyzing the Schur complements of matrices in a higher dimension. We hope that our algorithm and techniques open the door to improved solvers for positive semidefinite programming and its applications.

In this paper we provide a parallel algorithm that given any n-node m-edge directed graph and source vertex s computes all vertices reachable from s with Õ(m) work and n{1/2+o(1)} depth with high probability in n. This algorithm also computes a set of Õ(n) edges which when added to the graph preserves reachability and ensures that the diameter of the resulting graph is at most n{1/2+o(1)}. Our result improves upon the previous best known almost linear work reachability algorithm due to Fineman [1] which had depth Õ(n2/3). Further, we show how to leverage this algorithm to achieve improved distributed algorithms for single source reachability in the CONGEST model. In particular, we provide a distributed algorithm that given a n-node digraph of undirected hop-diameter D solves the single source reachability problem with Õ(n1/2+ n1/3+o(1)D2/3) rounds of the communication in the CONGEST model with high probability in n. Our algorithm is nearly optimal whenever D = O(n1/4-ε) for any constant ε > 0 and is the first nearly optimal algorithm for general graphs whose diameter is Ω(nδ) for any constant δ.

Optimal transportation, or computing the Wasserstein or ``earth mover's'' distance between two $n$-dimensional distributions, is a fundamental primitive which arises in many learning and statistical settings. We give an algorithm which solves the problem to additive $\epsilon$ accuracy with $\tilde{O}(1/\epsilon)$ parallel depth and $\tilde{O}\left(n^2/\epsilon\right)$ work. \cite{BlanchetJKS18, Quanrud19} obtained this runtime through reductions to positive linear programming and matrix scaling. However, these reduction-based algorithms use subroutines which may be impractical due to requiring solvers for second-order iterations (matrix scaling) or non-parallelizability (positive LP). Our methods match the previous-best work bounds by \cite{BlanchetJKS18, Quanrud19} while either improving parallelization or removing the need for linear system solves, and improve upon the previous best first-order methods running in time $\tilde{O}(\min(n^2 / \epsilon^2, n^{2.5} / \epsilon))$ \cite{DvurechenskyGK18, LinHJ19}. We obtain our results by a primal-dual extragradient method, motivated by recent theoretical improvements to maximum flow \cite{Sherman17}.

In this paper we provide nearly linear time algorithms for several problems closely associated with the classic Perron-Frobenius theorem, including computing Perron vectors, i.e. entrywise non-negative eigenvectors of non-negative matrices, and solving linear systems in asymmetric M-matrices, a generalization of Laplacian systems. The running times of our algorithms depend nearly linearly on the input size and polylogarithmically on the desired accuracy and problem condition number. Leveraging these results we also provide improved running times for a broader range of problems including computing random walk-based graph kernels, computing Katz centrality, and more. The running times of our algorithms improve upon previously known results which either depended polynomially on the condition number of the problem, required quadratic time, or only applied to special cases. We obtain these results by providing new iterative methods for reducing these problems to solving linear systems in Row-Column Diagonally Dominant (RCDD) matrices. Our methods are related to the classic shift-and-invert preconditioning technique for eigenvector computation and constitute the first alternative to the result in Cohen et al. (2016) for reducing stationary distribution computation and solving directed Laplacian systems to solving RCDD systems.

In this paper we consider the problem of efficiently computing ∊-sketches for the Laplacian and its pseudoinverse. Given a Laplacian and an error tolerance ∊, we seek to construct a function f such that for any vector x (chosen obliviously from f), with high probability (1 – ∊)x⊤ Ax ≤ f(x) ≤ (1 + ∊)x⊤ Ax where A is either the Laplacian or its pseudoinverse. Our goal is to construct such a sketch f efficiently and to store it in the least space possible. We provide nearly-linear time algorithms that, when given a Laplacian matrix ℒ ∊ ℝn×n and an error tolerance ∊, produce Õ(n/∊)-size sketches of both ℒ and its pseudoinverse. Our algorithms improve upon the previous best sketch size of Õ(n/∊1.6) for sketching the Laplacian form by [1] and O(n/∊2) for sketching the Laplacian pseudoinverse by [2]. Furthermore we show how to compute all-pairs effective resistances from our Õ(n/∊) size sketch in Õ(n2/∊) time. This improves upon the previous best running time of Õ(n2/∊2) by [3].