We provide faster running times for exactly solving discounted Markov Decision Processes (DMDPs) in strongly polynomial time. We obtain our results by efficiently reducing computing optimal values and policies in DMDPs to the easier tasks of policy evaluation and computing approximately optimal values. We provide both a straightforward deterministic reduction and a more efficient randomized variant that, together with advances in approximately solving DMDPs, yield our results.

We provide algorithms that compute $\epsilon$-estimates of the $\ell_p$-Lewis weights of a matrix $A \in \mathbb{R}^{m \times n}$ for $p \geq 4$ using $O(p^2 \log(m/\epsilon))$ rounds of leverage score computation, where $\ell_p$-Lewis weights and leverage scores are both standard measures of row importance. This improves upon the state-of-the-art round complexity of $O(p^3 \log(m/\epsilon))$ due to Fazel, Lee, Padmanabha, and Sidford (2022). We obtain our results by carefully applying a local variant of relatively smooth gradient descent to primal and dual forms of the $\ell_p$-Lewis weight optimization problem and providing tools to convert between different notions of approximate $\ell_p$-Lewis weights. This work subsumes the note “On computing approximate Lewis weights” by Apers, Gribling, and Sidford.

We study the fundamental classification problem of computing a separating hyperplane for a binary-labeled dataset of size $n$ with normalized $d$-dimensional features. Letting $\Phi \in \mathbb{R}^{n \times d}$ denote the feature matrix and $\gamma$ the margin of the maximum-margin separating hyperplane, we present a randomized algorithm that solves this problem in $\tilde{O}(\gamma^{-2/3}\, \operatorname{nnz}(\Phi) + \gamma^{-2(\omega+1)/3})$-sequential running time (work), $\tilde{O}(\gamma^{-2/3})$-parallel (computational) depth, and accesses $\Phi$ only through $\tilde{O}(\gamma^{-2/3})$-matrix-vector queries (matvecs). We also present a second, faster randomized algorithm with a $\tilde{O}(\gamma^{-2/3}\, \operatorname{nnz}(\Phi) + \gamma^{-2})$-sequential running time that uses $\tilde{O}(\gamma^{-2/3})$-matvecs to $\Phi$, but achieves only $\tilde{O}(\gamma^{-4/3})$-parallel depth. Both algorithms match the near-optimal deterministic matvec complexity recently established by Kornowski and Shamir [2025], Karmarkar et al. [2026] and achieve improved sequential runtime and parallel depth, albeit at the expense of using randomness.

We provide faster strongly polynomial time algorithms solving maximum flow in structured \(n\)-node \(m\)-arc networks. Our results imply an \(n^{\omega+o(1)}\)-time strongly polynomial time algorithms for computing a maximum bipartite \(b\)-matching where \(\omega\) is the matrix multiplication constant. Additionally, they imply an \(m^{1+o(1)}W\)-time algorithm for solving the problem on graphs with a given tree decomposition of width \(W\).

We provide new high-accuracy randomized algorithms for solving linear systems and regression problems that are well-conditioned except for \(k\) large singular values. For solving such \(d \times d\) positive definite systems our algorithms succeed whp. and run in time \(\tilde{O}(d^{2} + k^{\omega})\). For solving such regression problems in a matrix \(\textbf A \in \mathbb{R}^{n \times d}\) our methods succeed whp. and run in time \(\tilde{O}(\mathrm{nnz}(\textbf A) + d^{2} + k^{\omega})\) where \(\omega\) is the matrix multiplication exponent and \(\mathrm{nnz}(\textbf A)\) is the number of non-zeros in \(\textbf A\). Our methods nearly-match a natural complexity limit under dense inputs for these problems and improve upon a trade-off in prior approaches that obtain running times of either \(\tilde{O}(d^{2.065} + k^{\omega})\) or \(\tilde{O}(d^{2} + d\,k^{\omega-1})\) for \(d \times d\) systems. Moreover, we show how to obtain these running times even under the weaker assumption that all but \(k\) of the singular values have a suitably bounded generalized mean. Consequently, we give the first nearly-linear time algorithm for computing a multiplicative approximation to the nuclear norm of an arbitrary dense matrix. Our algorithms are built on three general recursive preconditioning frameworks, where matrix sketching and low-rank update formulas are carefully tailored to the problems’ structure.

We study the problem of computing an є-approximate Nash equilibrium of a two-player, bilinear game with a bounded payoff matrix A ∈ ℝm × n, when the players’ strategies are constrained to lie in simple sets. We provide algorithms which solve this problem in Õ(є−2/3) matrix-vector multiplies (matvecs) in two well-studied cases: ℓ1-ℓ1 (or zero-sum) games, where the players’ strategies are both in the probability simplex, and ℓ2-ℓ1 games (encompassing hard-margin SVMs), where the players’ strategies are in the unit Euclidean ball and probability simplex respectively. These results improve upon the previous state-of-the-art complexities of Õ(є−8/9) for ℓ1-ℓ1 and Õ(є−7/9) for ℓ2-ℓ1 due to [KOS ’25]. In both settings our results are nearly-optimal as they match lower bounds of [KS ’25] up to polylogarithmic factors.

In this paper we consider generalized flow problems where there is an m-edge n-node directed graph $G=(V, E)$ and each edge $e \in E$ has a loss factor $\gamma_{e}\gt 0$ governing whether the flow is increased or decreased as it crosses edge e. We provide a randomized $\widetilde{O}\left(\left(m+n^{1.5}\right) \cdot \operatorname{polylog}\left(\frac{W}{\delta}\right)\right)$ time algorithm for solving the generalized maximum flow and generalized minimum cost flow problems in this setting where $\delta$ is the target accuracy and W is the maximum of all costs, capacities, and loss factors and their inverses. This improves upon the previous state-of-the-art $\widetilde{O}\left(m \sqrt{n} \cdot \log ^{2}\left(\frac{W}{\delta}\right)\right)$ time algorithm, obtained by combining the algorithm of [17] with techniques from [29]. To obtain this result we provide new dynamic data structures and spectral results regarding the matrices associated to generalized flows and apply them through the interior point method framework of [39].1.1The full version of this paper is available at https://arxiv.org/abs/2510.17740.

In this paper we consider the problem of computing an $\epsilon$-approximate Nash Equilibrium of a zerosum game in a payoff matrix $A \in \mathbb{R}^{m \times n}$ with $O(1)$-bounded entries given access to a matrix-vector product oracle for A and its transpose $A^{\top}$. We provide a deterministic algorithm that solves the problem using $\tilde{O}\left(\epsilon^{-8 / 9}\right)$-oracle queries, where $\tilde{O}(\cdot)$ hides factors polylogarithmic in m, n, and $\epsilon^{-1}$. Our result improves upon the state-of-the-art query complexity of $\tilde{O}\left(\epsilon^{-1}\right)$ established by [Nemirovski, 2004] and [Nesterov, 2005]. We obtain this result through a general framework that yields improved deterministic query complexities for solving a broader class of minimax optimization problems which includes computing a linear classifier (hard-margin support vector machine) as well as linear regression.11This paper is an extended abstract. The full paper can be accessed at https://arxiv.org/abs/2509.04426.

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 develop optimization methods which offer new trade-offs between the number of gradient and Hessian computations needed to compute the critical point of a non-convex function. We provide a method that for a twice-differentiable $f\colon \mathbb{R}^d \rightarrow \mathbb{R}$ with $L_2$-Lipschitz Hessian, and input initial point with $\Delta$-bounded sub-optimality and sufficiently small $\epsilon > 0$ outputs an $\epsilon$-critical point, i.e., a point $x$ such that $\|\nabla f(x)\| \leq \epsilon$, using $\tilde{O}(\Delta L_2^{1/4} n_H^{-1/2}\epsilon^{-9/4})$ queries to a gradient oracle and $n_H$ queries to a Hessian oracle. As a consequence, we obtain an improved gradient query complexity of $\tilde{O}(d^{1/3}L_2^{1/2}\Delta\epsilon^{-3/2})$ in the case of bounded dimension and of $\tilde{O}(\Delta^{3/2} L_2^{3/4}\epsilon^{-9/4})$ in the case where we are allowed only a single Hessian query. We obtain these results through a more general algorithm which can handle approximate Hessian computations and recovers known prior state-of-the-art bounds of computing an $\epsilon$-critical point, under the additional assumption that $f$ has an $L_1$-Lipschitz gradient, with $O(\Delta L_2^{1/4}\epsilon^{-7/4})$-gradient queries.

We consider the problem of minimizing a $d$-dimensional Lipschitz convex function using a stochastic gradient oracle. We introduce and motivate a setting where the noise of the stochastic gradient is isotropic in that it is bounded in every direction with high probability. We then develop an algorithm for this setting which improves upon prior results by a factor of $d$ in certain regimes, and as a corollary, achieves a new state-of-the-art complexity for sub-exponential noise. We give matching lower bounds (up to polylogarithmic factors) for both results. Additionally, we develop an efficient quantum isotropifier, a quantum algorithm which converts a variance-bounded quantum sampling oracle into one that outputs an unbiased estimate with isotropic error. Combining our results, we obtain improved dimension-dependent rates for quantum stochastic convex optimization.

We provide an algorithm that maintains, against an adaptive adversary, a (1 — ε )-approximate maximum matching in n-node m-edge general (not necessarily bipartite) undirected graph undergoing edge deletions with high probability with (amortized) O (poly(ε-1, log n )) time per update. We also obtain the same update time for maintaining a fractional approximate weighted matching (and hence an approximation to the value of the maximum weight matching) and an integral approximate weighted matching in dense graphs.1 Our unweighted result improves upon the prior state-of-the-art which includes a poly(log n ) · 2O (1/ɛ2) update time [Assadi-Bernstein-Dudeja 2022] and an update time [Gupta-Peng 2013], and our weighted result improves upon the log n ) update time due to [Gupta-Peng 2013].

We consider the foundational problem of maintaining a (1 — ε )-approximate maximum weight matching (MWM) in an n-node dynamic graph undergoing edge insertions and deletions. We provide a general reduction that reduces the problem on graphs with a weight range of poly(n ) to poly(1/ε ) at the cost of just an additive poly(1/ε ) in update time. This improves upon the prior reduction of Gupta-Peng (FOCS 2013) which reduces the problem to a weight range of ε-O (1/ε) with a multiplicative cost of O (log n ).

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 )).

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.

We consider the problem of estimating the spectral density of a normalized graph adjacency matrix. Concretely, given an undirected graph $G = (V, E, w)$ with $n$ nodes and positive edge weights $w \in \mathbb{R}^{E}_{> 0}$, the goal is to return eigenvalue estimates $\widehat{\lambda}_1 \le \cdots\le \widehat{\lambda}_n$ such that \begin{align*} \frac{1}{n} \sum_{i\in\{1,\ldots, n\}}|\widehat{\lambda}_i-\lambda_i(N_G)|\le \varepsilon, \end{align*} where ${\lambda}_1(N_G)\le \cdots\le{\lambda}_n(N_G)$ are the eigenvalues of $G$’s normalized adjacency matrix, $N_G$. This goal is equivalent to requiring that the Wasserstein-1 distance between the uniform distribution on $\lambda_1, \ldots, \lambda_n$ and the uniform distribution on $\widehat{\lambda}_1, \ldots, \widehat{\lambda}_n$ is less than $\varepsilon$. We provide a randomized algorithm that achieves the guarantee above with $O(n\varepsilon^{-2})$ queries to a degree and neighbor oracle and in $O(n\varepsilon^{-3})$ time. This improves on previous state-of-the-art methods, including an $O(n\varepsilon^{-7})$ time algorithm from [Braverman et al., STOC 2022] and, for sufficiently small $\varepsilon$, a $2^{O(\varepsilon^{-1})}$ time method from [Cohen-Steiner et al., KDD 2018]. To achieve this result, we introduce a new notion of graph sparsification, which we call \emph{nuclear sparsification}. We provide an $O(n\varepsilon^{-2})$-query and $O(n\varepsilon^{-2})$-time algorithm for computing $O(n\varepsilon^{-2})$-sparse nuclear sparsifiers. We show that this bound is optimal in both its sparsity and query complexity, and we separate our results from the related notion of additive spectral sparsification. Of independent interest, we show that our sparsification method also yields the first \emph{deterministic} algorithm for spectral density estimation that scales linearly with $n$ (sublinear in the representation size of the graph).

We provide faster randomized algorithms for computing an $\epsilon$-optimal policy in a discounted Markov decision process with $A_{\text{tot}}$-state-action pairs, bounded rewards, and discount factor $\gamma$. We provide an $\tilde{O}(A_{\text{tot}}[(1 - \gamma)^{-3}\epsilon^{-2} + (1 - \gamma)^{-2}])$-time algorithm in the sampling setting, where the probability transition matrix is unknown but accessible through a generative model which can be queried in $\tilde{O}(1)$-time, and an $\tilde{O}(s + (1-\gamma)^{-2})$-time algorithm in the offline setting where the probability transition matrix is known and $s$-sparse. These results improve upon the prior state-of-the-art which either ran in $\tilde{O}(A_{\text{tot}}[(1 - \gamma)^{-3}\epsilon^{-2} + (1 - \gamma)^{-3}])$ time [Sidford, Wang, Wu, Ye 2018] in the sampling setting, $\tilde{O}(s + A_{\text{tot}} (1-\gamma)^{-3})$ time [Sidford, Wang, Wu, Yang, Ye 2018] in the offline setting, or time at least quadratic in the number of states using interior point methods for linear programming. We achieve our results by building upon prior stochastic variance-reduced value iteration methods [Sidford, Wang, Wu, Yang, Ye 2018]. We provide a variant that carefully truncates the progress of its iterates to improve the variance of new variance-reduced sampling procedures that we introduce to implement the steps. Our method is essentially model-free and can be implemented in $\tilde{O}(A_{\text{tot}})$-space when given generative model access. Consequently, our results take a step in closing the sample-complexity gap between model-free and model-based methods.

We consider the well-studied problem of completing a rank-$r$, $\mu$-incoherent matrix $\mathbf{M} \in \mathbb{R}^{d \times d}$ from incomplete observations. We focus on this problem in the semi-random setting where each entry is independently revealed with probability at least $p = \frac{\textup{poly}(r, \mu, \log d)}{d}$. Whereas multiple nearly-linear time algorithms have been established in the more specialized fully-random setting where each entry is revealed with probablity exactly $p$, the only known nearly-linear time algorithm in the semi-random setting is due to [CG18], whose sample complexity has a polynomial dependence on the inverse accuracy and condition number and thus cannot achieve high-accuracy recovery. Our main result is the first high-accuracy nearly-linear time algorithm for solving semi-random matrix completion, and an extension to the noisy observation setting. Our result builds upon the recent short-flat decomposition framework of [KLLST23a, KLLST23b] and leverages fast algorithms for flow problems on graphs to solve adaptive reweighting subproblems efficiently.

We provide an algorithm which, with high probability, maintains a (1 — ɛ)-approximate maximum flow on an undirected graph undergoing m-edge additions in amortized mo(1)ɛ-3 time per update. To obtain this result, we provide a more general algorithm that solves what we call the incremental, thresholded, p-norm flow problem that asks to determine the first edge-insertion in an undirected graph that causes the minimum ℓp-norm flow to decrease below a given threshold in value. Since we solve this thresholded problem, our data structure succeeds against an adaptive adversary that can only see the data structure's output. Furthermore, since our algorithm holds for p = 2, we obtain improved algorithms for dynamically maintaining the effective resistance between a pair of vertices in an undirected graph undergoing edge insertions.

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.

We study dynamic (1−є)-approximate rounding of fractional matchings—a key ingredient in numerous breakthroughs in the dynamic graph algorithms literature. Our first contribution is a surprisingly simple deterministic rounding algorithm in bipartite graphs with amortized update time O(є−1 log2 (є−1 · n)), matching an (unconditional) recourse lower bound of Ω(є−1) up to logarithmic factors. Moreover, this algorithm’s update time improves provided the minimum (non-zero) weight in the fractional matching is lower bounded throughout. Combining this algorithm with novel dynamic partial rounding algorithms to increase this minimum weight, we obtain a number of algorithms that improve this dependence on n. For example, we give a high-probability randomized algorithm with Õ(є−1 · (loglogn)2)-update time against adaptive adversaries. Using our rounding algorithms, we also round known (1−є)-decremental fractional bipartite matching algorithms with no asymptotic overhead, thus improving on state-of-the-art algorithms for the decremental bipartite matching problem. Further, we provide extensions of our results to general graphs and to maintaining almost-maximal matchings.

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 show that any memory-constrained, first-order algorithm which minimizesd-dimensional, 1-Lipschitz convex functions over the unit ball to 1/poly(d) accuracy using at most d1.25 - δbits of memory must make at least \(\tilde{\Omega }(d^{1 + (4/3)\delta })\) first-order queries (for any constant \(\delta \in [0, 1/4]\) ). Consequently, the performance of such memory-constrained algorithms are at least a polynomial factor worse than the optimal Õ(d) query bound for this problem obtained by cutting plane methods that use Õ(d2) memory. This resolves one of the open problems in the COLT 2019 open problem publication of Woodworth and Srebro.

We present an algorithm which given any $m$-edge directed graph with positive integer capacities at most $U$, vertices $a$ and $b$, and an approximation parameter $\epsilon \in (0, 1)$ computes an additive $\epsilon mU$-approximate $a$-$b$ maximum flow in time $m^{1+o(1)}/\sqrt{\epsilon}$. By applying the algorithm for $\epsilon = (mU)^{-2/3}$, rounding to an integral flow, and using augmenting paths, we obtain an algorithm which computes an exact $a$-$b$ maximum flow in time $m^{4/3+o(1)}U^{1/3}$ and an algorithm which given an $m$-edge bipartite graph computes an exact maximum cardinality matching in time $m^{4/3+o(1)}$.

We study lower bounds for the problem of approximating a one dimensional distribution given (noisy) measurements of its moments. We show that there are distributions on $[-1,1]$ that cannot be approximated to accuracy $\epsilon$ in Wasserstein-1 distance even if we know \emph{all} of their moments to multiplicative accuracy $(1\pm2^{-\Omega(1/\epsilon)})$; this result matches an upper bound of Kong and Valiant [Annals of Statistics, 2017]. To obtain our result, we provide a hard instance involving distributions induced by the eigenvalue spectra of carefully constructed graph adjacency matrices. Efficiently approximating such spectra in Wasserstein-1 distance is a well-studied algorithmic problem, and a recent result of Cohen-Steiner et al. [KDD 2018] gives a method based on accurately approximating spectral moments using $2^{O(1/\epsilon)}$ random walks initiated at uniformly random nodes in the graph.As a strengthening of our main result, we show that improving the dependence on $1/\epsilon$ in this result would require a new algorithmic approach. Specifically, no algorithm can compute an $\epsilon$-accurate approximation to the spectrum of a normalized graph adjacency matrix with constant probability, even when given the transcript of $2^{\Omega(1/\epsilon)}$ random walks of length $2^{\Omega(1/\epsilon)}$ started at random nodes.

Sparse recovery is one of the most fundamental and well-studied inverse problems.Standard statistical formulations of the problem are provably solved by general convex programming techniques and more practical, fast (nearly-linear time) iterative methods. However, these latter “fast algorithms” have previously been observed to be brittle in various real-world settings.We investigate the brittleness of fast sparse recovery algorithms to generative model changes through the lens of studying their robustness to a “helpful” semi-random adversary, a framework for testing overfitting to input assumptions. We consider the following basic model: let $\mathbf{A} \in \mathbb{R}^{n \times d}$ be a measurement matrix containing an unknown subset of rows $\mathbf{G} \in \mathb{R}^{m \times d}$ which are bounded and satisfy the restricted isometry property (RIP), but is otherwise arbitrary. Letting $x^\star \in \mathbb{R}^d$ be $s$-sparse, and given either exact or noisy measurements, $b = \mathbf{A} x^\star$ or $b = \mathbf{A} x^\star + \xi$, we design algorithms recovering $x^\star$ information-theoretically optimally in nearly-linear time. We extend our algorithm to hold for weaker generative models relaxing our planted RIP row subset assumption to a natural weighted variant, and show that our method’s guarantees naturally interpolate the quality of the measurement matrix to, in some parameter regimes, run in sublinear time.Our approach differs from that of prior fast iterative methods with provable guarantees under semi-random generative models [CG18, LSTZ20], which typically separate the problem of learning the planted instance from the estimation problem, i.e. they attempt to first learn the planted “good” instance (in our case, the matrix $\mathbf{G}$). However, natural conditions on a submatrix which make sparse recovery tractable, such as RIP, are NP-hard to verify and hence first learning a sufficient row reweighting appears challenging. We eschew this approach and design a new iterative method, tailored to the geometry of sparse recovery, which is provably robust to our semi-random model. Our hope is that our approach opens the door to new robust, efficient algorithms for other natural statistical inverse problems.

In this paper, we introduce a new, spectral notion of approximation between directed graphs, which we call singular value (SV) approximation. SV-approximation is stronger than previous notions of spectral approximation considered in the literature, including spectral approximation of Laplacians for undirected graphs [ST04], standard approximation for directed graphs [CKP+17], and unit-circle (UC) approximation for directed graphs [AKM+20]. Further, SV approximation enjoys several useful properties not possessed by previous notions of approximation, e.g., it is preserved under products of randomwalk matrices and bounded matrices. We provide a nearly linear-time algorithm for SV-sparsifying (and hence UC-sparsifying) Eulerian directed graphs, as well as $\ell$-step random walks on such graphs, for any $\ell \leq \operatorname{poly}(n)$. Combined with the Eulerian scaling algorithms of [CKK+18], given an arbitrary (not necessarily Eulerian) directed graph and a set S of vertices, we can approximate the stationary probability mass of the $\left(S, S^{c}\right)$ cut in an $\ell$-step random walk to within a multiplicative error of $1 / \operatorname{polylog}(n)$ and an additive error of $1 / \operatorname{poly}(n)$ in nearly linear time. As a starting point for these results, we provide a simple black-box reduction from SV-sparsifying Eulerian directed graphs to SV-sparsifying undirected graphs; such a directed-to-undirected reduction was not known for previous notions of spectral approximation.

We give a deterministic $m^{1+o(1)}$ time algorithm that computes exact maximum flows and minimum-cost flows on directed graphs with m edges and polynomially bounded integral demands, costs, and capacities. As a consequence, we obtain the first running time improvement for deterministic algorithms that compute maximum-flow in graphs with polynomial bounded capacities since the work of Goldberg-Rao [J.ACM ’98].Our algorithm builds on the framework of Chen-Kyng-Liu-Peng-Gutenberg-Sachdeva [FOCS ’22] that computes an optimal flow by computing a sequence of $m^{1+o(1)}$-approximate undirected minimum-ratio cycles. We develop a deterministic dynamic graph data-structure to compute such a sequence of minimum-ratio cycles in an amortized $m^{o(1)}$ time per edge update. Our key technical contributions are deterministic analogues of the vertex sparsification and edge sparsification components of the data-structure from Chen et al. For the vertex sparsification component, we give a method to avoid the randomness in Chen et al. which involved sampling random trees to recurse on. For the edge sparsification component, we design a deterministic algorithm that maintains an embedding of a dynamic graph into a sparse spanner. We also show how our dynamic spanner can be applied to give a deterministic data structure that maintains a fully dynamic low-stretch spanning tree on graphs with polynomially bounded edge lengths, with subpolynomial average stretch and subpolynomial amortized time per edge update.

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.

In this paper we study the problem of minimizing a submodular function $f: 2^{V} \rightarrow \mathbb{R}$ that is guaranteed to have a k-sparse minimizer. We give a deterministic algorithm that computes an additive $\epsilon$-approximate minimizer of such f in $\widetilde{O}(\operatorname{poly}(k) \log (|f| / \epsilon))$ parallel depth using a polynomial number of queries to an evaluation oracle of f, where $|f|=\max_{S \subseteq V}|f(S)|$. Further, we give a randomized algorithm that computes an exact minimizer of f with high probability using $\widetilde{O}(|V| \cdot \operatorname{poly}(k))$ queries and polynomial time. When $k=\widetilde{O}(1)$, our algorithms use either nearly-constant parallel depth or a nearly-linear number of evaluation oracle queries. All previous algorithms for this problem either use $\Omega(|V|)$ parallel depth or $\Omega\left(|V|^{2}\right)$ queries. In contrast to state-of-the-art weakly-polynomial and strongly-polynomial time algorithms for SFM, our algorithms use first-order optimization methods, e.g., mirror descent and follow the regularized leader. We introduce what we call sparse dual certificates, which encode information on the structure of sparse minimizers, and both our parallel and sequential algorithms provide new algorithmic tools for allowing first-order optimization methods to efficiently compute them. Correspondingly, our algorithm does not invoke fast matrix multiplication or general linear system solvers and in this sense is more combinatorial than previous state-of-the-art methods.

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 give a new framework for solving the fundamental problem of low-rank matrix completion, i.e., approximating a rank- r matrix $\mathbf{M} \in \mathbb{R}^{m \times n}$ (where $m \geq n$) from random observations. First, we provide an algorithm which completes M on $99 \%$ of rows and columns under no further assumptions on M from $\approx m r$ samples and using $\approx m r^{2}$ time. Then, assuming the row and column spans of M satisfy additional regularity properties, we show how to boost this partial completion guarantee to a full matrix completion algorithm by aggregating solutions to regression problems involving the observations. In the well-studied setting where M has incoherent row and column spans, our algorithms complete M to high precision from $m r^{2+o(1)}$ observations in $m r^{3+o(1)}$ time (omitting logarithmic factors in problem parameters), improving upon the prior state-of-the-art [JN15] which used $\approx m r^{5}$ samples and $\approx m r^{7}$ time. Under an assumption on the row and column spans of M we introduce (which is satisfied by random subspaces with high probability), our sample complexity improves to an almost information-theoretically optimal $m r^{1+o(1)}$, and our runtime improves to $m r^{2+o(1)}$. Our runtimes have the appealing property of matching the best known runtime to verify that a rankr decomposition $\mathrm{UV}^{\top}$ agrees with the sampled observations. We also provide robust variants of our algorithms that, given random observations from $\mathrm{M}+\mathrm{N}$ with $\|\mathrm{N}\|_{\mathrm{F}} \leq \Delta$, complete M to Frobenius norm distance $\approx r^{1.5} \Delta$ in the same runtimes as the noiseless setting. Prior noisy matrix completion algorithms [CP10] only guaranteed a distance of $\approx \sqrt{n} \Delta$.

We give a quantum algorithm for computing an $\epsilon$-approximate Nash equilibrium of a zero-sum game in a $m \times n$ payoff matrix with bounded entries. Given a standard quantum oracle for accessing the payoff matrix our algorithm runs in time $\widetilde{O}(\sqrt{m + n}\cdot \epsilon^{-2.5} + \epsilon^{-3})$ and outputs a classical representation of the $\epsilon$-approximate Nash equilibrium. This improves upon the best prior quantum runtime of $\widetilde{O}(\sqrt{m + n} \cdot \epsilon^{-3})$ obtained by [van Apeldoorn, Gilyen ’19] and the classical $\widetilde{O}((m + n) \cdot \epsilon^{-2})$ runtime due to [Grigoradis, Khachiyan ’95] whenever $\epsilon = \Omega((m +n)^{-1})$. We obtain this result by designing new quantum data structures for efficiently sampling from a slowly-changing Gibbs distribution.

Minimizing a convex function with access to a first order oracle---that returns the function evaluation and (sub)gradient at a query point---is a canonical optimization problem and a fundamental primitive in machine learning. Gradient-based methods are the most popular approaches used for solving the problem, owing to their simplicity and computational efficiency. These methods, however, do not achieve the information-theoretically optimal query complexity for minimizing the underlying function to small error, which are achieved by more expensive techniques based on cutting-plane methods. Is it possible to achieve the information-theoretically query complexity without using these more complex and computationally expensive methods? In this work, we use memory as a lens to understand this, and show that is is not possible to achieve optimal query complexity without using significantly more memory than that used by gradient descent.

We study the complexity of computing stationary Nash equilibrium (NE) in n-player infinite-horizon general-sum stochastic games. We focus on the problem of computing NE in such stochastic games when each player is restricted to choosing a stationary policy and rewards are discounted. First, we prove that computing such NE is in PPAD (in addition to clearly being PPAD-hard). Second, we consider turn-based specializations of such games where at each state there is at most a single player that can take actions and show that these (seemingly-simpler) games remain PPAD-hard. Third, we show that under further structural assumptions on the rewards computing NE in such turn-based games is possible in polynomial time. Towards achieving these results we establish structural facts about stochastic games of broader utility, including monotonicity of utilities under single-state single-action changes and reductions to settings where each player controls a single state.

We consider the parallel complexity of submodular function minimization (SFM). We provide a pair of methods which obtain two new query versus depth trade-offs a submodular function defined on subsets of $n$ elements that has integer values between $-M$ and $M$. The first method has depth $2$ and query complexity $n^{O(M)}$ and the second method has depth $\widetilde{O}(n^{1/3} M^{2/3})$ and query complexity $O(\mathrm{poly}(n, M))$. Despite a line of work on improved parallel lower bounds for SFM, prior to our work the only known algorithms for parallel SFM either followed from more general methods for sequential SFM or highly-parallel minimization of convex $\ell_2$-Lipschitz functions. Interestingly, to obtain our second result we provide the first highly-parallel algorithm for minimizing $\ell_\infty$-Lipschitz function over the hypercube which obtains near-optimal depth for obtaining constant accuracy.

We provide new algorithms and conditional hardness for the problem of estimating effective resistances in $n$-node $m$-edge undirected, expander graphs. We provide an $\widetilde{O}(m\epsilon^{-1})$-time algorithm that produces with high probability, an $\widetilde{O}(n\epsilon^{-1})$-bit sketch from which the effective resistance between any pair of nodes can be estimated, to $(1 \pm \epsilon)$-multiplicative accuracy, in $\widetilde{O}(1)$-time. Consequently, we obtain an $\widetilde{O}(m\epsilon^{-1})$-time algorithm for estimating the effective resistance of all edges in such graphs, improving (for sparse graphs) on the previous fastest runtimes of $\widetilde{O}(m\epsilon^{-3/2})$ [Chu et. al. 2018] and $\widetilde{O}(n^2\epsilon^{-1})$ [Jambulapati, Sidford, 2018] for general graphs and $\widetilde{O}(m + n\epsilon^{-2})$ for expanders [Li, Sachdeva 2022]. We complement this result by showing a conditional lower bound that a broad set of algorithms for computing such estimates of the effective resistances between all pairs of nodes require $\widetilde{\Omega}(n^2 \epsilon^{-1/2})$-time, improving upon the previous best such lower bound of $\widetilde{\Omega}(n^2 \epsilon^{-1/13})$ [Musco et. al. 2017]. Further, we leverage the tools underlying these results to obtain improved algorithms and conditional hardness for more general problems of sketching the pseudoinverse of positive semidefinite matrices and estimating functions of their eigenvalues.

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 consider the problem of minimizing a continuous function given given access to a natural quantum generalization of a stochastic gradient oracle. We provide two new methods for the special case of minimizing a Lipschitz convex function. Each method obtains a dimension versus accuracy trade-off which is provably unachievable classically and we prove that one method is asymptotically optimal in low-dimensional settings. Additionally, we provide quantum algorithms for computing a critical point of a smooth non-convex function at rates not known to be achievable classically. To obtain these results we build upon the quantum multivariate mean estimation result of Cornelissen et al. and provide a general quantum variance reduction technique of independent interest.

Abstract. Let [Formula: see text] be an [Formula: see text]-node [Formula: see text]-edge weighted undirected graph, where [Formula: see text] is a real length function defined on its edges, and let [Formula: see text] denote the girth of [Formula: see text], i.e., the length of a shortest cycle. We present an algorithm that, for any input, integer [Formula: see text], in [Formula: see text] expected time finds a cycle of length at most [Formula: see text]. This algorithm nearly matches an [Formula: see text]-time algorithm of Kadria et al. [ Algorithmic trade-offs for girth approximation in undirected graphs, in Proceedings of the 2022 Annual ACM-SIAM Symposium on Discrete Algorithms (SODA), SIAM, 2022, pp. 1471–1492] which applied to unweighted graphs of girth 3. For weighted graphs, this result also improves upon the previous state-of-the-art algorithm that in [Formula: see text] time, where [Formula: see text] is an integral length function, finds a cycle of length at most [Formula: see text] of Kadria et al. [ Algorithmic trade-offs for girth approximation in undirected graphs, in Proceedings of the 2022 Annual ACM-SIAM Symposium on Discrete Algorithms (SODA), SIAM, 2022, pp. 1471–1492]. For [Formula: see text], this result improves upon the result of Roditty and Tov [ ACM Trans. Algorithms, 9 (2013), pp. 15:1–15:13].

In this paper we provide an algorithm for maintaining a (1−є)-approximate maximum flow in a dynamic, capacitated graph undergoing edge insertions. Over a sequence of m insertions to an n-node graph where every edge has capacity O(poly(m)) our algorithm runs in time O(m √n · є−1). To obtain this result we design dynamic data structures for the more general problem of detecting when the value of the minimum cost circulation in a dynamic graph undergoing edge insertions achieves value at most F (exactly) for a given threshold F. Over a sequence m insertions to an n-node graph where every edge has capacity O(poly(m)) and cost O(poly(m)) we solve this thresholded minimum cost flow problem in O(m √n). Both of our algorithms succeed with high probability against an adaptive adversary. We obtain these results by dynamizing the recent interior point method by [Chen et al. ‍FOCS 2022] used to obtain an almost linear time algorithm for minimum cost flow, and introducing a new dynamic data structure for maintaining minimum ratio cycles in an undirected graph that succeeds with high probability against adaptive adversaries.

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)).

We design accelerated algorithms with improved rates for several fundamental classes of optimization problems. Our algorithms all build upon techniques related to the analysis of primal-dual extragradient methods via relative Lipschitzness proposed recently by Cohen, Sidford, and Tian ’21. (1) We study separable minimax optimization problems of the form $\min_x \max_y f(x) - g(y) + h(x, y)$, where $f$ and $g$ have smoothness and strong convexity parameters $(L^x, \mu^x)$, $(L^y, \mu^y)$, and h is convex-concave with a $(\Lambda^{xx}, \Lambda^{xy}, \Lambda^{yy})$-blockwise operator norm bounded Hessian. We provide an algorithm using $\tilde{O}(\sqrt{\frac{L^x}{\mu^x}} + \sqrt{\frac{L^y}{\mu^y}} + \frac{\Lambda^{xx}}{\mu^x} + \frac{\Lambda^{xy}}{\sqrt{\mu^x\mu^y}} + \frac{\Lambda^{yy}}{\mu^y})$ gradient queries. Notably, for convex-concave minimax problems with bilinear coupling (e.g. quadratics), where $\Lambda^{xx} = \Lambda^{yy} = 0$, our rate matches a lower bound of Zhang, Hong, and Zhang ’19. (2) We study finite sum optimization problems of the form $\min_x \frac 1 n \sum_{i \in [n]} f_i(x)$, where each $f_i$ is $L_i$-smooth and the overall problem is $\mu$-strongly convex. We provide an algorithm using $\tilde{O}(n + \sum_{i \in [n]} \sqrt{\frac{L_i}{n\mu}} )$ gradient queries. Notably, when the smoothness bounds $\{L_i\}_{i\in[n]}$ are non-uniform, our rate improves upon accelerated SVRG (Lin et al., Frostig et al. ’15) and Katyusha (Allen-Zhu ’17) by up to a $\sqrt{n}$ factor. (3) We generalize our algorithms for minimax and finite sum optimization to solve a natural family of minimax finite sum optimization problems at an accelerated rate, encapsulating both above results up to a logarithmic factor.

We provide new gradient-based methods for efficiently solving a broad class of ill-conditioned optimization problems. We consider the problem of minimizing a function $f : \mathbb{R}^d \rightarrow \mathbb{R}$ which is implicitly decomposable as the sum of $m$ unknown non-interacting smooth, strongly convex functions and provide a method which solves this problem with a number of gradient evaluations that scales (up to logarithmic factors) as the product of the square-root of the condition numbers of the components. This complexity bound (which we prove is nearly optimal) can improve almost exponentially on that of accelerated gradient methods, which grow as the square root of the condition number of $f$. Additionally, we provide efficient methods for solving stochastic, quadratic variants of this multiscale optimization problem. Rather than learn the decomposition of $f$ (which would be prohibitively expensive), our methods apply a clean recursive “Big-Step-Little-Step” interleaving of standard methods. The resulting algorithms use $\tilde{\mathcal{O}}(d m)$ space, are numerically stable, and open the door to a more fine-grained understanding of the complexity of convex optimization beyond condition number.

We show that any memory-constrained, first-order algorithm which minimizes $d$-dimensional, $1$-Lipschitz convex functions over the unit ball to $1/\mathrm{poly}(d)$ accuracy using at most $d^{1.25 - \delta}$ bits of memory must make at least $\Omega(d^{1 + (4/3)\delta})$ first-order queries (for any constant $\delta \in [0, 1/4]$). Consequently, the performance of such memory-constrained algorithms are a polynomial factor worse than the optimal $\tilde{O}(d)$ query bound for this problem obtained by cutting plane methods that use $\tilde{O}(d^2)$ memory. This resolves one of the open problems in the COLT 2019 open problem publication of Woodworth and Srebro.

We provide a generic technique for constructing families of submodular functions to obtain lower bounds for submodular function minimization (SFM). Applying this technique, we prove that any deterministic SFM algorithm on a ground set of n elements requires at least $\Omega(n\log n)$ queries to an evaluation oracle. This is the first super-linear query complexity lower bound for SFM and improves upon the previous best lower bound of 2n given by [Graur et al., ITCS 2020]. Using our construction, we also prove that any (possibly randomized) parallel SFM algorithm, which can make up to poly $(n)$ queries per round, requires at least $\Omega(n/\log n)$ rounds to minimize a submodular function. This improves upon the previous best lower bound of $\tilde{\Omega}(n^{1/3})$ rounds due to [Chakrabarty et al., FOCS 2021], and settles the parallel complexity of query-efficient SFM up to logarithmic factors due to a recent advance in [Jiang, SODA 2021].

Designing dynamic graph algorithms against an adaptive adversary is a major goal in the field of dynamic graph algorithms. While a few such algorithms are known for spanning trees, matchings, and single-source shortest paths, very little was known for an important primitive like graph sparsifiers. The challenge is how to approximately preserve so much information about the graph (e.g., all-pairs distances and all cuts) without revealing the algorithms' underlying randomness to the adaptive adversary. In this paper we present the first non-trivial efficient adaptive algorithms for maintaining spanners and cut sparisifers. These algorithms in turn imply improvements over existing algorithms for other problems. Our first algorithm maintains a polylog$(n)$-spanner of size $\tilde O(n)$ in polylog$(n)$ amortized update time. The second algorithm maintains an $O(k)$-approximate cut sparsifier of size $\tilde O(n)$ in $\tilde O(n^{1/k})$ amortized update time, for any $k\ge1$, which is polylog$(n)$ time when $k=\log(n)$. The third algorithm maintains a polylog$(n)$-approximate spectral sparsifier in polylog$(n)$ amortized update time. The amortized update time of both algorithms can be made worst-case by paying some sub-polynomial factors. Prior to our result, there were near-optimal algorithms against oblivious adversaries (e.g. Baswana et al. [TALG'12] and Abraham et al. [FOCS'16]), but the only non-trivial adaptive dynamic algorithm requires $O(n)$ amortized update time to maintain $3$- and $5$-spanner of size $O(n^{1+1/2})$ and $O(n^{1+1/3})$, respectively [Ausiello et al. ESA'05]. Our results are based on two novel techniques. The first technique, is a generic black-box reduction that allows us to assume that the graph undergoes only edge deletions and, more importantly, remains an expander with almost-uniform degree. The second technique we call proactive resampling. [...]

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 an efficient unified plug-in approach for estimating symmetric properties of distributions given $n$ independent samples. Our estimator is based on profile-maximum-likelihood (PML) and is sample optimal for estimating various symmetric properties when the estimation error $\epsilon \gg n^{-1/3}$. This result improves upon the previous best accuracy threshold of $\epsilon \gg n^{-1/4}$ achievable by polynomial time computable PML-based universal estimators \cite{ACSS20, ACSS20b}. Our estimator reaches a theoretical limit for universal symmetric property estimation as \cite{Han20} shows that a broad class of universal estimators (containing many well known approaches including ours) cannot be sample optimal for every $1$-Lipschitz property when $\epsilon \ll n^{-1/3}$.

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 present an algorithm for computing approximate ℓp Lewis weights to high precision. Given a full-rank A ∊ ℝm × n with m ≥ n and a scalar p > 2, our algorithm computes ∊-approximate ℓp Lewis weights of A in Õp(log(1/∊)) iterations; the cost of each iteration is linear in the input size plus the cost of computing the leverage scores of DA for diagonal D ∊ ℝm × m. Prior to our work, such a computational complexity was known only for p ∊ (0,4) [CP15], and combined with this result, our work yields the first polylogarithmic-depth polynomial-work algorithm for the problem of computing ℓp Lewis weights to high precision for all constant p > 0. An important consequence of this result is also the first polylogarithmic-depth polynomial-work algorithm for computing a nearly optimal self-concordant barrier for a polytope.

We present several new efficient algorithms for approximating the girth, g, of weighted and unweighted n-vertex, m-edge undirected graphs. For undirected graphs with polynomially bounded, integer, non-negative edge weights, we provide an algorithm that for every integer k ≥ 1, runs in Õ(m + n1 + 1/k log g) time and returns a cycle of length at most 2kg. For unweighted, undirected graphs we present an algorithm that for every k ≥ 1, runs in Õ(n1 + 1/k) time and returns a cycle of length at most 2k[g/2], an almost k-approximation. Both algorithms provide trade-offs between the running time and the quality of the approximation. We also obtain faster algorithms for approximation factors better than 2, and improved approximations when the girth is odd or small (e.g., 3 and 4).

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.

In this paper we consider the problem of computing the likelihood of the profile of a discrete distribution, i.e., the probability of observing the multiset of element frequencies, and computing a profile maximum likelihood (PML) distribution, i.e., a distribution with the maximum profile likelihood. For each problem we provide polynomial time algorithms that given $n$ i.i.d. samples from a discrete distribution, achieve an approximation factor of $\exp\left(-O(\sqrt{n} \log n) \right)$, improving upon the previous best-known bound achievable in polynomial time of $\exp(-O(n^{2/3} \log n))$ (Charikar, Shiragur and Sidford, 2019). Through the work of Acharya, Das, Orlitsky and Suresh (2016), this implies a polynomial time universal estimator for symmetric properties of discrete distributions in a broader range of error parameter. To obtain our results on PML we establish new connections between PML and the well-studied Bethe and Sinkhorn approximations to the permanent (Vontobel, 2012 and 2014). It is known that the PML objective is proportional to the permanent of a certain Vandermonde matrix (Vontobel, 2012) with $\sqrt{n}$ distinct columns, i.e. with non-negative rank at most $\sqrt{n}$. This allows us to show that the convex approximation to computing PML distributions studied in (Charikar, Shiragur and Sidford, 2019) is governed, in part, by the quality of Sinkhorn approximations to the permanent. We show that both Bethe and Sinkhorn permanents are $\exp(O(k \log(N/k)))$ approximations to the permanent of $N \times N$ matrices with non-negative rank at most $k$. This improves upon the previous known bounds of $\exp(O(N))$ and combining these insights with careful rounding of the convex relaxation yields our results.

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 prove new upper and lower bounds for sample complexity of finding an $\epsilon$-optimal policy of an infinite-horizon average-reward Markov decision process (MDP) given access to a generative model. When the mixing time of the probability transition matrix of all policies is at most $t_\mathrm{mix}$, we provide an algorithm that solves the problem using $\widetilde{O}(t_\mathrm{mix} \epsilon^{-3})$ (oblivious) samples per state-action pair. Further, we provide a lower bound showing that a linear dependence on $t_\mathrm{mix}$ is necessary in the worst case for any algorithm which computes oblivious samples. We obtain our results by establishing connections between infinite-horizon average-reward MDPs and discounted MDPs of possible further utility.

We show that standard extragradient methods (i.e. mirror prox and dual extrapolation) recover optimal accelerated rates for first-order minimization of smooth convex functions. To obtain this result we provide a fine-grained characterization of the convergence rates of extragradient methods for solving monotone variational inequalities in terms of a natural condition we call relative Lipschitzness. We further generalize this framework to handle local and randomized notions of relative Lipschitzness and thereby recover rates for box-constrained $\ell_\infty$ regression based on area convexity and complexity bounds achieved by accelerated (randomized) coordinate descent for smooth convex function minimization.

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.

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.

In this paper we provide new randomized algorithms with improved runtimes for solving linear programs with two-sided constraints. In the special case of the minimum cost flow problem on n-vertex m-edge graphs with integer polynomially-bounded costs and capacities we obtain a randomized method which solves the problem in Õ(m + n1.5) time. This improves upon the previous best runtime of Õ(m √n) [Lee-Sidford’14] and, in the special case of unit-capacity maximum flow, improves upon the previous best runtimes of m4/3 + o(1) [Liu-Sidford’20, Kathuria’20] and Õ(m √n) [Lee-Sidford’14] for sufficiently dense graphs.

We give a deterministic $O(\log n\cdot\log\log n)$-space algorithm for approximately solving linear systems given by Laplacians of undirected graphs, and consequently also approximating hitting times, commute times, and escape probabilities for undirected graphs. Previously, such systems were known to be solvable by randomized algorithms using $O(\log n)$ space [D. Doron, F. Le Gall, and A. Ta-Shma, Probabilistic logarithmic-space algorithms for Laplacian solvers, in APPROX/RANDOM 2017, LIPIcs. Leibniz Int. Proc. Inform. 81, Schloss Dagstuhl. Leibniz-Zent. Inform., Wadern, Germany, 2017, 41] and hence by deterministic algorithms using $O(\log^{3/2} n)$ space [M. Saks and S. Zhou, J. Comput. System Sci., 58 (1999), pp. 376--403]. Our algorithm combines ideas from time-efficient Laplacian solvers [D. A. Spielman and S.-H. Teng, Nearly-linear time algorithms for graph partitioning, graph sparsification, and solving linear systems, in STOC 2004, ACM, New York, 2004, pp. 81--90; R. Peng and D. A. Spielman, An efficient parallel solver for SDD linear systems, in STOC 2014, ACM, New York, 2014, pp. 333--342] with ideas used to show that Undirected S-T Connectivity is in deterministic logspace [O. Reingold, J. ACM, 55 (2008); E. Rozenman and S. Vadhan, Derandomized squaring of graphs, in RANDOM 2005, Lecture Notes in Comput. Sci. 3624, Springer, Berlin, 2005, pp. 436--447].

In this paper we settle the sampling complexity of solving discounted two-player turn-based zero-sum stochastic games up to polylogarithmic factors. Given a stochastic game with discount factor $\gamma\in(0,1)$ we provide an algorithm that computes an $\epsilon$-optimal strategy with high-probability given $\tilde{O}((1 - \gamma)^{-3} \epsilon^{-2})$ samples from the transition function for each state-action-pair. Our algorithm runs in time nearly linear in the number of samples and uses space nearly linear in the number of state-action pairs. As stochastic games generalize Markov decision processes (MDPs) our runtime and sample complexities are optimal due to \cite{azar2013minimax}. We achieve our results by showing how to generalize a near-optimal Q-learning based algorithms for MDP, in particular \cite{sidford2018near}, to two-player strategy computation algorithms. This overcomes limitations of standard Q-learning and strategy iteration or alternating minimization based approaches and we hope will pave the way for future reinforcement learning results by facilitating the extension of MDP results to multi-agent settings with little loss.

Given a matrix $\mathbf{A}\in\R^{n\times d}$ and a vector $b\in\R^{d}$, we show how to compute an $\epsilon$-approximate solution to the regression problem $ \min_{x\in\R^{d}}\frac{1}{2} \norm{\mathbf{A} x-b}_{2}^{2} $ in time $ \widetilde{O} ((n+\sqrt{d\cdot\kappa_{\text{sum}}}) s \log\epsilon^{-1}) $ where $\kappa_{\text{sum}}=\tr\left(\mathbf{A}^{\top}\mathbf{A}\right)/\lambda_{\min}(\mathbf{A}^{\top}\mathbf{A})$ and $s$ is the maximum number of non-zero entries in a row of $\mathbf{A}$. This improves upon the previous best running time of $ \widetilde{O} ((n+\sqrt{n \cdot\kappa_{\text{sum}}}) s \log\epsilon^{-1})$. We achieve our result through an interesting combination of leverage score sampling, proximal point methods, and accelerated coordinate descent methods. Further, we show that our method not only matches the performance of previous methods up to polylogarithmic factors, but further improves whenever leverage scores of rows are small. We also provide a non-linear generalization of these results that improves the running time for solving a broader class of ERM problems and expands the set of ERM problems provably solvable in nearly linear time.

In this paper, we provide near-optimal accelerated first-order methods for minimizing a broad class of smooth nonconvex functions that are unimodal on all lines through a minimizer. This function class, which we call the class of smooth quasar-convex functions, is parameterized by a constant $$\gamma \in (0,1]$$: $$\gamma = 1$$ encompasses the classes of smooth convex and star-convex functions, and smaller values of $$\gamma$$ indicate that the function can be "more nonconvex." We develop a variant of accelerated gradient descent that computes an $$\epsilon$$-approximate minimizer of a smooth $$\gamma$$-quasar-convex function with at most $$O(\gamma^{-1} \epsilon^{-1/2} \log(\gamma^{-1} \epsilon^{-1}))$$ total function and gradient evaluations. We also derive a lower bound of $$\Omega(\gamma^{-1} \epsilon^{-1/2})$$ on the worst-case number of gradient evaluations required by any deterministic first-order method, showing that, up to a logarithmic factor, no deterministic first-order method can improve upon ours.

In this paper, we provide a deterministic ~O(log N)-space algorithm for estimating random walk probabilities on undirected graphs, and more generally Eulerian directed graphs, to within inverse polynomial additive error (ε = 1/poly(N)) where N is the length of the input. Previously, this problem was known to be solvable by a randomized algorithm using space O(log N) (following Aleliunas et al., FOCS '79) and by a deterministic algorithm using space O(log3/2N) (Saks and Zhou, FOCS '95 and JCSS '99), both of which held for arbitrary directed graphs but had not been improved even for undirected graphs. We also give improvements on the space complexity of both of these previous algorithms for non-Eulerian directed graphs when the error is negligible (ε = 1/Nω(1)), generalizing what Hoza and Zuckerman (FOCS '18) recently showed for the special case of distinguishing whether a random walk probability is 0 or greater than ε. We achieve these results by giving new reductions between powering Eulerian random-walk matrices and inverting Eulerian Laplacian matrices, providing a new notion of spectral approximation for Eulerian graphs that is preserved under powering, and giving the first deterministic ~O(log N)-space algorithm for inverting Eulerian Laplacian matrices. The latter algorithm builds on the work of Murtagh et al. (FOCS '17) that gave a deterministic ~O(log N)-space algorithm for inverting undirected Laplacian matrices, and the work of Cohen et al. (FOCS '19) that gave a randomized ~O(N)-time algorithm for inverting Eulerian Laplacian matrices. A running theme throughout these contributions is an analysis of “cycle-lifted graphs,” where we take a graph and “lift” it to a new graph whose adjacency matrix is the tensor product of the original adjacency matrix and a directed cycle (or variants of one).

We present an ~O(m+n1.5)-time randomized algorithm for maximum cardinality bipartite matching and related problems (e.g. transshipment, negative-weight shortest paths, and optimal transport) on m-edge, n-node graphs. For maximum cardinality bipartite matching on moderately dense graphs, i.e. m=Ω(n1.5), our algorithm runs in time nearly linear in the input size and constitutes the first improvement over the classic O(m√n)-time [Dinic 1970; Hopcroft-Karp 1971; Karzanov 1973] and ~O(nω)-time algorithms [Ibarra-Moran 1981] (where currently ω ≈ 2.373). On sparser graphs, i.e. when m=n9/8+δfor any constant , our result improves upon the recent advances of [Madry 2013] and [Liu-Sidford 2020b, 2020a] which achieve an ~O(m4/3+o(1)) runtime. We obtain these results by combining and advancing recent lines of research in interior point methods (IPMs) and dynamic graph algorithms. First, we simplify and improve the IPM of [v.d.Brand-Lee-Sidford-Song 2020], providing a general primal-dual IPM framework and new sampling-based techniques for handling infeasibility induced by approximate linear system solvers. Second, we provide a simple sublinear-time algorithm for detecting and sampling high-energy edges in electric flows on expanders and show that when combined with recent advances in dynamic expander decompositions, this yields efficient data structures for maintaining the iterates of both [v.d.Brand et al.] and our new IPMs. Combining this general machinery yields a simpler ~O(n√m) time algorithm for matching based on the logarithmic barrier function, and our state-of-the-art ~O(m+n1.5) time algorithm for matching based on the [Lee-Sidford 2014] barrier (as regularized in [v.d.Brand et al.]).

We develop primal-dual coordinate methods for solving bilinear saddle-point problems of the form minx ∈ Xmaxy ∈ YyTAx which contain linear programming, classification, and regression as special cases. Our methods push existing fully stochastic sublinear methods and variance-reduced methods towards their limits in terms of per-iteration complexity and sample complexity. We obtain nearly-constant per-iteration complexity by designing efficient data structures leveraging Taylor approximations to the exponential and a binomial heap. We improve sample complexity via low-variance gradient estimators using dynamic sampling distributions that depend on both the iterates and the magnitude of the matrix entries. Our runtime bounds improve upon those of existing primal-dual methods by a factor depending on sparsity measures of the m by n matrix A. For example, when rows and columns have constant l1/l2 norm ratios, we offer improvements by a factor of m+n in the fully stochastic setting and √{m+n} in the variance-reduced setting. We apply our methods to computational geometry problems, i.e. minimum enclosing ball, maximum inscribed ball, and linear regression, and obtain improved complexity bounds. For linear regression with an elementwise nonnegative matrix, our guarantees improve on exact gradient methods by a factor of √{nnz(A)/(m+n)}.

We present an algorithm, which given any m-edge n-vertex directed graph with positive integer capacities at most U computes a maximum s-t flow for any vertices s and t in O(m4/3+o(1)U1/3) time. This improves upon the previous best running times of O(m11/8+o(1)U1/4) [1], Õ(m√nlogU) [2] and O(mn) [3] when the graph is not too dense and doesn't have large capacities. We build upon advances for sparse maxflow based on interior point methods [1], [4], [5]. Whereas these methods increase the energy of local ℓ2-norm minimizing electrical flows, we instead increase the Bregman divergence value of flows which minimize the Bregman divergence with respect to a weighted log barrier. This allows us to trace the central path with progress depending only on ℓ∞norm bounds on the congestion vector as opposed to the ℓ4norm, which arises in these prior works. Further, we show that smoothed ℓ2-ℓpflows [6], [7] which were used to maximize energy [1] can also be used to efficiently maximize divergence, thereby yielding our desired runtimes. We believe our approach towards Bregman divergences of barriers may be of further interest.

We present a unified framework based on primal-dual stochastic mirror descent for approximately solving infinite-horizon Markov decision processes (MDPs) given a generative model. When applied to an average-reward MDP with $A_{tot}$ total actions and mixing time bound $t_{mix}$ our method computes an $\epsilon$-optimal policy with an expected $\widetilde{O}(t_{mix}^2 A_{tot} \epsilon^{-2})$ samples from the state-transition matrix, removing the ergodicity dependence of prior art. When applied to a $\gamma$-discounted MDP with $A_{tot}$ total actions our method computes an $\epsilon$-optimal policy with an expected $\widetilde{O}((1-\gamma)^{-4} A_{tot} \epsilon^{-2})$ samples, improving over the best-known primal-dual methods while matching the state-of-the-art up to a $(1-\gamma)^{-1}$ factor. Both methods are model-free, update state values and policies simultaneously, and run in time linear in the number of samples taken. We achieve these results through a more general stochastic mirror descent framework for solving bilinear saddle-point problems with simplex and box domains and we demonstrate the flexibility of this framework by providing further applications to constrained MDPs.

In this paper we provide a new efficient algorithm for approximately computing the profile maximum likelihood (PML) distribution, a prominent quantity in symmetric property estimation. We provide an algorithm which matches the previous best known efficient algorithms for computing approximate PML distributions and improves when the number of distinct observed frequencies in the given instance is small. We achieve this result by exploiting new sparsity structure in approximate PML distributions and providing a new matrix rounding algorithm, of independent interest. Leveraging this result, we obtain the first provable computationally efficient implementation of PseudoPML, a general framework for estimating a broad class of symmetric properties. Additionally, we obtain efficient PML-based estimators for distributions with small profile entropy, a natural instance-based complexity measure. Further, we provide a simpler and more practical PseudoPML implementation that matches the best-known theoretical guarantees of such an estimator and evaluate this method empirically.

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 propose and analyze algorithms for distributionally robust optimization of convex losses with conditional value at risk (CVaR) and $\chi^2$ divergence uncertainty sets. We prove that our algorithms require a number of gradient evaluations independent of training set size and number of parameters, making them suitable for large-scale applications. For $\chi^2$ uncertainty sets these are the first such guarantees in the literature, and for CVaR our guarantees scale linearly in the uncertainty level rather than quadratically as in previous work. We also provide lower bounds proving the worst-case optimality of our algorithms for CVaR and a penalized version of the $\chi^2$ problem. Our primary technical contributions are novel bounds on the bias of batch robust risk estimation and the variance of a multilevel Monte Carlo gradient estimator due to [Blanchet & Glynn, 2015]. Experiments on MNIST and ImageNet confirm the theoretical scaling of our algorithms, which are 9-36 times more efficient than full-batch methods.

In this paper we provide improved running times and oracle complexities for approximately minimizing a submodular function. Our main result is a randomized algorithm, which given any submodular function defined on n-elements with range [–1, 1], computes an ε-additive approximate minimizer in Õ(n/ε2) oracle evaluations with high probability. This improves over the Õ(n5/3/ε2) oracle evaluation algorithm of Chakrabarty et al. (STOC 2017) and the Õ(n3/2/ε2) oracle evaluation algorithm of Hamoudi et al… Further, we leverage a generalization of this result to obtain efficient algorithms for minimizing a broad class of nonconvex functions. For any function f with domain [0, 1]n that satisfies for all i ≠ j and is L-Lipschitz with respect to the L∞-norm we give an algorithm that computes an ε-additive approximate minimizer with Õ(n · poly(L/ε) function evaluation with high probability.

In this paper, we resolve the complexity problem of spectral graph sparcification in dynamic streams up to polylogarithmic factors. Using a linear sketch we design a streaming algorithm that uses Õ(n) space, and with high probability, recovers a spectral sparsifier from the sketch in Õ(n) time.1 Prior results either achieved near optimal Õ(n) space, but Ω(n2) recovery time [Kapralov et al. ‘14], or ran in o(n2) time, but used polynomially suboptimal space [Ahn et al ‘13]. Our main technical contribution is a novel method for recovering graph edges with high effective resistance from a linear sketch. We show how to do so in nearly linear time by ‘bucketing’ vertices of the input graph into clusters using a coarse approximation to the graph's effective resistance metric. A second main contribution is a new pseudorandom generator (PRG) for linear sketching algorithms. Constructed from a locally computable randomness extractor, our PRG stretches a seed of Õ(n) random bits polynomially in length with just logO(1) n run-time cost per evaluation. This improves on Nisan's commonly used PRG, which in our setting would require Õ(n) time per evaluation. Our faster PRG is essential to simultaneously achieving near optimal space and time complexity.

In this paper we provide an O(nd+d 3) time randomized algorithm for solving linear programs with d variables and n constraints with high probability. To obtain this result we provide a robust, primal-dual O(√d)-iteration interior point method inspired by the methods of Lee and Sidford (2014, 2019) and show how to efficiently implement this method using new data-structures based on heavy-hitters, the Johnson–Lindenstrauss lemma, and inverse maintenance. Interestingly, we obtain this running time without using fast matrix multiplication and consequently, barring a major advance in linear system solving, our running time is near optimal for solving dense linear programs among algorithms that do not use fast matrix multiplication.

In this paper we provide a Õ(m√n) time algorithm that computes a 3-multiplicative approximation of the girth of a n-node m-edge directed graph with non-negative edge lengths. This is the first algorithm which approximates the girth of a directed graph up to a constant multiplicative factor faster than All-Pairs Shortest Paths (APSP) time, i.e. O(mn). Additionally, for any integer k ≥ 1, we provide a deterministic algorithm for a O(kloglogn)-multiplicative approximation to the girth in directed graphs in Õ(m 1+1/k ) time. Combining the techniques from these two results gives us an algorithm for a O(klogk)-multiplicative approximation to the girth in directed graphs in Õ(m 1+1/k ) time. Our results naturally also provide algorithms for improved constructions of roundtrip spanners, the analog of spanners in directed graphs.

In this paper we provide an algorithm which given any m-edge n-vertex directed graph with integer capacities at most U computes a maximum s-t flow for any vertices s and t in m 11/8+o(1) U 1/4 time with high probability. This running time improves upon the previous best of Õ(m 10/7 U 1/7) (Mądry 2016), Õ(m √n logU) (Lee Sidford 2014), and O(mn) (Orlin 2013) when the graph is not too dense or has large capacities.

We give a deterministic, nearly logarithmic-space algorithm that given an undirected graph G, a positive integer r, and a set S of vertices, approximates the conductance of S in the r-step random walk on G to within a factor of 1+epsilon, where epsilon>0 is an arbitrarily small constant. More generally, our algorithm computes an epsilon-spectral approximation to the normalized Laplacian of the r-step walk. Our algorithm combines the derandomized square graph operation [Eyal Rozenman and Salil Vadhan, 2005], which we recently used for solving Laplacian systems in nearly logarithmic space [Murtagh et al., 2017], with ideas from [Cheng et al., 2015], which gave an algorithm that is time-efficient (while ours is space-efficient) and randomized (while ours is deterministic) for the case of even r (while ours works for all r). Along the way, we provide some new results that generalize technical machinery and yield improvements over previous work. First, we obtain a nearly linear-time randomized algorithm for computing a spectral approximation to the normalized Laplacian for odd r. Second, we define and analyze a generalization of the derandomized square for irregular graphs and for sparsifying the product of two distinct graphs. As part of this generalization, we also give a strongly explicit construction of expander graphs of every size.

We propose a near-optimal method for highly smooth convex optimization. More precisely, in the oracle model where one obtains the $p^{th}$ order Taylor expansion of a function at the query point, we propose a method with rate of convergence $\tilde{O}(1/k^{\frac{ 3p +1}{2}})$ after $k$ queries to the oracle for any convex function whose $p^{th}$ order derivative is Lipschitz.

We show that a simple randomized sketch of the matrix multiplicative weight (MMW) update enjoys (in expectation) the same regret bounds as MMW, up to a small constant factor. Unlike MMW, where every step requires full matrix exponentiation, our steps require only a single product of the form $e^A b$, which the Lanczos method approximates efficiently. Our key technique is to view the sketch as a \emph{randomized mirror projection}, and perform mirror descent analysis on the \emph{expected projection}. Our sketch solves the online eigenvector problem, improving the best known complexity bounds by $\Omega(\log^5 n)$. We also apply this sketch to semidefinite programming in saddle-point form, yielding a simple primal-dual scheme with guarantees matching the best in the literature.

In this merged paper, we consider the problem of minimizing a convex function with Lipschitz-continuous $p$-th order derivatives. Given an oracle which when queried at a point returns the first $p$-derivatives of the function at that point we provide some methods which compute an $\e$ approximate minimizer in $O\left(\e^{-\frac{2}{3p+1}} \right)$ iterations. These methods match known lower bounds up to polylogarithmic factors for constant $p$.

In this paper we consider the classic matroid intersection problem: given two matroids M1= (V, I1) and M2= (V, I2) defined over a common ground set V , compute a set S ∈ I1∩ I2of largest possible cardinality, denoted by r. We consider this problem both in the setting where each Mi is accessed through an independence oracle, i.e. a routine which returns whether or not a set S ∈ Iiin Tindtime, and the setting where each Mi is accessed through a rank oracle, i.e. a routine which returns the size of the largest independent subset of S in Miin Tranktime. In each setting we provide faster exact and approximate algorithms. Given an independence oracle, we provide an exact O(nr log r · Tind) time algorithm. This improves upon previous best known running times of O(nr1.5·Tind) due to Cunningham O(n2·Tindin 1986 and + n3) due to Lee, Sidford, and Wong in 2015. We also provide two algorithms which compute a (1- ε-approximate solution to matroid intersection running in times O(n1.5/ε1.5· Tind) and O((n2r-1ε-2+ r1.5ε-4.5) · Tind), respectively. These results improve upon the O(nr/ε · Tind)time algorithm of Cunningham (noted recently by Chekuri and Quanrud). Given a rank oracle, we provide algorithms with even better dependence on n and r. We provide an O(n√r log n · Trank)time exact algorithm and an O(nε-1log n · Trank)-time algorithm which obtains a (1 - 0)-approximation to the matroid intersection problem. The former result improves over the O(nr · Trank+ n3)-time algorithm by Lee, Sidford, and Wong. The rank oracle is of particular interest as the matroid intersection problem with this oracle is a special case (via Edmond's minimax characterization of matroid intersection) of the submodular function minimization (SFM) problem with an evaluation oracle, and understanding SFM query complexity is an outstanding open question.

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 δ.

A landmark result of non-smooth convex optimization is that gradient descent is an optimal algorithm whenever the number of computed gradients is smaller than the dimension $d$. In this paper we study the extension of this result to the parallel optimization setting. Namely we consider optimization algorithms interacting with a highly parallel gradient oracle, that is one that can answer $\mathrm{poly}(d)$ gradient queries in parallel. We show that in this case gradient descent is optimal only up to $\tilde{O}(\sqrt{d})$ rounds of interactions with the oracle. The lower bound improves upon a decades old construction by Nemirovski which proves optimality only up to $d^{1/3}$ rounds (as recently observed by Balkanski and Singer), and the suboptimality of gradient descent after $\sqrt{d}$ rounds was already observed by Duchi, Bartlett and Wainwright. In the latter regime we propose a new method with improved complexity, which we conjecture to be optimal. The analysis of this new method is based upon a generalized version of the recent results on optimal acceleration for highly smooth convex optimization.

We present a randomized primal-dual algorithm that solves the problem minx maxy y^T A x to additive error epsilon in time nnz(A) + sqrt{nnz(A) n} / epsilon, for matrix A with larger dimension n and nnz(A) nonzero entries. This improves the best known exact gradient methods by a factor of sqrt{nnz(A) / n} and is faster than fully stochastic gradient methods in the accurate and/or sparse regime epsilon < sqrt{n / nnz(A)$. Our results hold for x,y in the simplex (matrix games, linear programming) and for x in an \ell_2 ball and y in the simplex (perceptron / SVM, minimum enclosing ball). Our algorithm combines the Nemirovski's "conceptual prox-method" and a novel reduced-variance gradient estimator based on "sampling from the difference" between the current iterate and a reference point.

In this paper we provide a general framework for estimating symmetric properties of distributions from i.i.d. samples. For a broad class of symmetric properties we identify the {\em easy} region where empirical estimation works and the {\em difficult} region where more complex estimators are required. We show that by approximately computing the profile maximum likelihood (PML) distribution \cite{ADOS16} in this difficult region we obtain a symmetric property estimation framework that is sample complexity optimal for many properties in a broader parameter regime than previous universal estimation approaches based on PML. The resulting algorithms based on these \emph{pseudo PML distributions} are also more practical.

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}.

Given a n-by-d data matrix A, principal component projection (PCP) and principal component regression (PCR), i.e. projection and regression restricted to the top-eigenspace of A, are fundamental problems in machine learning, optimization, and numerical analysis. In this paper we provide the first algorithms that solve these problems in nearly linear time for fixed eigenvalue distribution and large n. This improves upon previous methods which had superlinear running times when either the number of top eigenvalues or gap between the eigenspaces were large. We achieve our results by applying rational polynomial approximations to reduce the problem to solving asymmetric linear systems which we solve by a variant of SVRG. We corroborate these findings with preliminary empirical experiments.

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.

Estimating symmetric properties of a distribution, e.g. support size, coverage, entropy, distance to uniformity, are among the most fundamental problems in algorithmic statistics. While these properties have been studied extensively and separate optimal estimators have been produced, in striking recent work Acharya et al. provided a single estimator that is competitive for each. They showed that the value of the property on the distribution that approximately maximizes profile likelihood (PML), i.e. the probability of observed frequency of frequencies, is sample competitive with respect to a broad class of estimators. Unfortunately, prior to this work, there was no known polynomial time algorithm to compute such an approximation or use PML to obtain a universal plug-in estimator.

We consider the problem of performing linear regression over a stream of d-dimensional examples, and show that any algorithm that uses a subquadratic amount of memory exhibits a slower rate of convergence than can be achieved without memory constraints. Specifically, consider a sequence of labeled examples (a1,b1), (a2,b2)…, with ai drawn independently from a d-dimensional isotropic Gaussian, and where bi = ⟨ ai, x⟩ + ηi, for a fixed x ∈ ℝd with ||x||2 = 1 and with independent noise ηi drawn uniformly from the interval [−2−d/5,2−d/5]. We show that any algorithm with at most d2/4 bits of memory requires at least Ω(d loglog1/є) samples to approximate x to ℓ2 error є with probability of success at least 2/3, for є sufficiently small as a function of d. In contrast, for such є, x can be recovered to error є with probability 1−o(1) with memory O(d2 log(1/є)) using d examples. This represents the first nontrivial lower bounds for regression with super-linear memory, and may open the door for strong memory/sample tradeoffs for continuous optimization.

There is widespread sentiment that fast gradient methods (e.g. Nesterov’s acceleration, conjugate gradient, heavy ball) are not effective for the purposes of stochastic optimization due to their instability and error accumulation. Numerous works have attempted to quantify these instabilities in the face of either statistical or non-statistical. This work considers these issues for the special case of stochastic approximation for the least squares regression problem, and our main result refutes this conventional wisdom by showing that acceleration can be made robust to statistical errors. In particular, this work introduces an accelerated stochastic gradient method that provably achieves the minimax optimal statistical risk faster than stochastic gradient descent. Critical to the analysis is a sharp characterization of accelerated stochastic gradient descent as a stochastic process. We hope this characterization gives insights towards the broader question of designing simple and effective accelerated stochastic methods for more general convex and non-convex optimization problems.

We consider the problem of minimizing a convex function over a convex set given access only to an evaluation oracle for the function and a membership oracle for the set. We give a simple algorithm which solves this problem with $\tilde{O}(n^{2})$ oracle calls and $\tilde{O}(n^{3})$ additional arithmetic operations. Using this result, we obtain more efficient reductions among the five basic oracles for convex sets and functions defined by Gr{ö}tschel, Lov{á}sz, and Schrijver (1988).

In this paper, we show how to solve directed Laplacian systems in nearly-linear time. Given a linear system in an n × n Eulerian directed Laplacian with m nonzero entries, we show how to compute an ε-approximate solution in time O(m logO(1)(n) log (1/ε)). Through reductions from [Cohen et al. FOCS'16], this gives the first nearly-linear time algorithms for computing ε-approximate solutions to row or column diagonally dominant linear systems (including arbitrary directed Laplacians) and computing ε-approximations to various properties of random walks on directed graphs, including stationary distributions, personalized PageRank vectors, hitting times, and escape probabilities. These bounds improve upon the recent almost-linear algorithms of [Cohen et al. STOC'17], which gave an algorithm to solve Eulerian Laplacian systems in time O((m+n2O(√ log n log log n))logO(1)(n ε-1)). To achieve our results, we provide a structural result that we believe is of independent interest. We show that Eulerian Laplacians (and therefore the Laplacians of all strongly connected directed graphs) have sparse approximate LU-factorizations. That is, for every such directed Laplacian there are lower upper triangular matrices each with at most Õ(n) nonzero entries such that there product spectrally approximates the directed Laplacian in an appropriate norm. This claim can be viewed as an analog of recent work on sparse Cholesky factorizations of Laplacians of undirected graphs. We show how to construct such factorizations in nearly-linear time and prove that once constructed they yield nearly-linear time algorithms for solving directed Laplacian systems.

In this paper we provide faster algorithms for approximately solving ℓ∞regression, a fundamental problem prevalent in both combinatorial and continuous optimization. In particular we provide an accelerated coordinate descent method which converges in k iterations at a O(1/k) rate independent of the dimension of the problem, and whose iterations can be implemented cheaply for many structured matrices. Our algorithm can be viewed as an alternative approach to the recent breakthrough result of Sherman [She17] which achieves a similar running time improvement over classic algorithmic approaches, i.e. smoothing and gradient descent, which either converge at a O(1/√k) rate or have running times with a worse dependence on problem parameters. Our running times match those of [She17] across a broad range of parameters and in certain cases, improves upon it. We demonstrate the efficacy of our result by providing faster algorithms for the well-studied maximum flow problem. We show how to leverage our algorithm to achieve a runtime of Õ(m + √ns/ε) to compute an ε-approximate maximum flow, for an undirected graph with m edges, n vertices, and where s is the squared ℓ2norm of the congestion of any optimal flow. As s = O(m) this yields a running time of Õ(m + √nm/ε), generically improving upon the previous best known runtime of Õ(m/ε) in [She17] whenever the graph is slightly dense. Moreover, we show how to leverage this result to achieve improved exact algorithms for maximum flow on a variety of unit capacity graphs. We achieve these results by providing an accelerated coordinate descent method capable of provably exploiting dynamic measures of coordinate smoothness for smoothed versions of ℓ∞regression. Our analysis leverages the structure of the Hessian of the smoothed problem via a simple bound on its trace, as well as techniques for exploiting column sparsity of the constraint matrix for faster sampling and improved smoothness estimates. We hope that the work of this paper can serve as an important step towards achieving even faster maximum flow algorithms.

Understanding the singular value spectrum of a matrix $A \in \mathbb{R}^{n \times n}$ is a fundamental task in countless applications. In matrix multiplication time, it is possible to perform a full SVD and directly compute the singular values $σ_1,...,σ_n$. However, little is known about algorithms that break this runtime barrier. Using tools from stochastic trace estimation, polynomial approximation, and fast system solvers, we show how to efficiently isolate different ranges of $A$'s spectrum and approximate the number of singular values in these ranges. We thus effectively compute a histogram of the spectrum, which can stand in for the true singular values in many applications. We use this primitive to give the first algorithms for approximating a wide class of symmetric matrix norms in faster than matrix multiplication time. For example, we give a $(1 + ε)$ approximation algorithm for the Schatten-$1$ norm (the nuclear norm) running in just $\tilde O((nnz(A)n^{1/3} + n^2)ε^{-3})$ time for $A$ with uniform row sparsity or $\tilde O(n^{2.18} ε^{-3})$ time for dense matrices. The runtime scales smoothly for general Schatten-$p$ norms, notably becoming $\tilde O (p \cdot nnz(A) ε^{-3})$ for any $p \ge 2$. At the same time, we show that the complexity of spectrum approximation is inherently tied to fast matrix multiplication in the small $ε$ regime. We prove that achieving milder $ε$ dependencies in our algorithms would imply faster than matrix multiplication time triangle detection for general graphs. This further implies that highly accurate algorithms running in subcubic time yield subcubic time matrix multiplication. As an application of our bounds, we show that precisely computing all effective resistances in a graph in less than matrix multiplication time is likely difficult, barring a major algorithmic breakthrough.

In this paper, we obtain improved running times for regression and top eigenvector computation for numerically sparse matrices. Given a data matrix $\mat{A} \in \R^{n \times d}$ where every row $a \in \R^d$ has $\|a\|_2^2 \leq L$ and numerical sparsity $\leq s$, i.e. $\|a\|_1^2 / \|a\|_2^2 \leq s$, we provide faster algorithms for these problems for many parameter settings. For top eigenvector computation, when $\gap > 0$ is the relative gap between the top two eigenvectors of $\mat{A}^\top \mat{A}$ and $r$ is the stable rank of $\mat{A}$ we obtain a running time of $\otilde(nd + r(s + \sqrt{r s}) / \gap^2)$ improving upon the previous best unaccelerated running time of $O(nd + r d / \gap^2)$. As $r \leq d$ and $s \leq d$ our algorithm everywhere improves or matches the previous bounds for all parameter settings. For regression, when $\mu > 0$ is the smallest eigenvalue of $\mat{A}^\top \mat{A}$ we obtain a running time of $\otilde(nd + (nL / \mu) \sqrt{s nL / \mu})$ improving upon the previous best unaccelerated running time of $\otilde(nd + n L d / \mu)$. This result expands when regression can be solved in nearly linear time from when $L/\mu = \otilde(1)$ to when $L / \mu = \otilde(d^{2/3} / (sn)^{1/3})$. Furthermore, we obtain similar improvements even when row norms and numerical sparsities are non-uniform and we show how to achieve even faster running times by accelerating using approximate proximal point \cite{frostig2015regularizing} / catalyst \cite{lin2015universal}. Our running times depend only on the size of the input and natural numerical measures of the matrix, i.e. eigenvalues and $\ell_p$ norms, making progress on a key open problem regarding optimal running times for efficient large-scale learning.

In this paper we consider the problem of computing an $\epsilon$-optimal policy of a discounted Markov Decision Process (DMDP) provided we can only access its transition function through a generative sampling model that given any state-action pair samples from the transition function in $O(1)$ time. Given such a DMDP with states $\states$, actions $\actions$, discount factor $\gamma\in(0,1)$, and rewards in range $[0, 1]$ we provide an algorithm which computes an $\epsilon$-optimal policy with probability $1 - \delta$ where {\it both} the run time spent and number of sample taken is upper bounded by \[ O\left[\frac{|\cS||\cA|}{(1-\gamma)^3 \epsilon^2} \log \left(\frac{|\cS||\cA|}{(1-\gamma)\delta \epsilon} \right) \log\left(\frac{1}{(1-\gamma)\epsilon}\right)\right] ~. \] For fixed values of $\epsilon$, this improves upon the previous best known bounds by a factor of $(1 - \gamma)^{-1}$ and matches the sample complexity lower bounds proved in \cite{azar2013minimax} up to logarithmic factors. We also extend our method to computing $\epsilon$-optimal policies for finite-horizon MDP with a generative model and provide a nearly matching sample complexity lower bound.

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].

Theoretically elegant and ubiquitous in practice, the Lanczos method can approximate f (A)x for any symmetric matrix A ∊ ℝn × n, vector x ∊ ℝn, and function f. In exact arithmetic, the method's error after k iterations is bounded by the error of the best degree-k polynomial uniformly approximating the scalar function f(x) on the range [λmin(A), λmax(A)]. However, despite decades of work, it has been unclear if this powerful guarantee holds in finite precision. We resolve this problem, proving that when , Lanczos essentially matches the exact arithmetic guarantee if computations use roughly log(nC║A║) bits of precision. Our proof extends work of Druskin and Knizhnerman [11], leveraging the stability of the classic Chebyshev recurrence to bound the stability of any polynomial approximating f(x). We also study the special case of f(A) = A–1 for positive definite A, where stronger guarantees hold for Lanczos. In exact arithmetic the algorithm performs as well as the best polynomial approximating 1/x at each of A's eigenvalues, rather than on the full range [λmin (A), λmax (A)]. In seminal work, Greenbaum gives a natural approach to extending this bound to finite precision: she proves that finite precision Lanczos and the related conjugate gradient method match any polynomial approximating 1/x in a tiny range around each eigenvalue [17]. For A–1, Greenbaum's bound appears stronger than our result. However, we exhibit matrices with condition number κ where exact arithmetic Lanczos converges in polylog(k) iterations, but Greenbaum's bound predicts at best Ω(κ1/5) iterations in finite precision. It thus cannot offer more than a polynomial improvement over the O(κ1/2) bound achievable via our result for general f(A). Our analysis bounds the power of stable approximating polynomials and raises the question of if they fully characterize the behavior of finite precision Lanczos in solving linear systems. If they do, convergence in less than poly(K) iterations cannot be expected, even for matrices with clustered, skewed, or otherwise favorable eigenvalue distributions.

The girth of a graph, i.e. the length of its shortest cycle, is a fundamental graph parameter. Unfortunately all known algorithms for computing, even approximately, the girth and girth-related structures in directed weighted m-edge and n-node graphs require Ω(min{nω, mn}) time (for 2 ≤ ω < 2.373). In this paper, we drastically improve these runtimes as follows: • Multiplicative Approximations in Nearly Linear Time: We give an algorithm that in Õ(m) time computes an Õ(1)-multiplicative approximation of the girth as well as an Õ(1)-multiplicative roundtrip spanner with Õ(n) edges with high probability (w.h.p). • Nearly Tight Additive Approximations: For unweighted graphs and any a ∊ (0, 1) we give an algorithm that in Õ(mn1–a) time computes an O(na)-additive approximation of the girth, w.h.p. We show that the runtime of our algorithm cannot be significantly improved without a breakthrough in combinatorial boolean matrix multiplication. We also show that if the girth is O(na), then the same guarantee can be achieved via a deterministic algorithm. Our main technical contribution to achieve these results is the first nearly linear time algorithm for computing roundtrip covers, a directed graph decomposition concept key to previous roundtrip spanner constructions. Previously it was not known how to compute these significantly faster than Ω(mn) time. Given the traditional difficulty in efficiently processing directed graphs, we hope our techniques may find further applications.

In this paper we provide faster algorithms for approximately solving discounted Markov Decision Processes in multiple parameter regimes. Given a discounted Markov Decision Process (DMDP) with |S| states, |A| actions, discount factor γ ∊ (0, 1), and rewards in the range [–M, M], we show how to compute an ∊-optimal policy, with probability 1 – δ in time This contribution reflects the first nearly linear time, nearly linearly convergent algorithm for solving DMDP's for intermediate values of γ. We also show how to obtain improved sublinear time algorithms and provide an algorithm which computes an ∊-optimal policy with probability 1 – δ in time provided we can sample from the transition function in O(1) time. Interestingly, we obtain our results by a careful modification of approximate value iteration. We show how to combine classic approximate value iteration analysis with new techniques in variance reduction. Our fastest algorithms leverage further insights to ensure that our algorithms make monotonic progress towards the optimal value. This paper is one of few instances in using sampling to obtain a linearly convergent linear programming algorithm and we hope that the analysis may be useful more broadly.

We give a deterministic Õ(log n)-space algorithm for approximately solving linear systems given by Laplacians of undirected graphs, and consequently also approximating hitting times, commute times, and escape probabilities for undirected graphs. Previously, such systems were known to be solvable by randomized algorithms using O(log n) space (Doron, Le Gall, and Ta-Shma, 2017) and hence by deterministic algorithms using O(log3/2n) space (Saks and Zhou, FOCS 1995 and JCSS 1999). Our algorithm combines ideas from time-efficient Laplacian solvers (Spielman and Teng, STOC `04; Peng and Spielman, STOC `14) with ideas used to show that UNDIRECTED S-T CONNECTIVITY is in deterministic logspace (Reingold, STOC `05 and JACM `08; Rozenman and Vadhan, RANDOM `05).

This work provides a simplified proof of the statistical minimax optimality of (iterate averaged) stochastic gradient descent (SGD), for the special case of least squares. This result is obtained by analyzing SGD as a stochastic process and by sharply characterizing the stationary covariance matrix of this process. The finite rate optimality characterization captures the constant factors and addresses model mis-specification.

We develop and analyze a variant of Nesterov’s accelerated gradient descent (AGD) for minimization of smooth non-convex functions. We prove that one of two cases occurs: either our AGD variant converges quickly, as if the function was convex, or we produce a certificate that the function is “guilty” of being non-convex. This non-convexity certificate allows us to exploit negative curvature and obtain deterministic, dimension-free acceleration of convergence for non-convex functions. For a function $f$ with Lipschitz continuous gradient and Hessian, we compute a point $x$ with $\|\nabla f(x)\| \le \epsilon$ in $O(\epsilon^{-7/4} \log(1/ \epsilon) )$ gradient and function evaluations. Assuming additionally that the third derivative is Lipschitz, we require only $O(\epsilon^{-5/3} \log(1/ \epsilon) )$ evaluations.

Submodular function minimization (SFM) is a fundamental discrete optimization problem which generalizes many well known problems, has applications in various fields, and can be solved in polynomial time. Owing to applications in computer vision and machine learning, fast SFM algorithms are highly desirable. The current fastest algorithms [Lee, Sidford, Wong, 2015] run in O(n2lognM· EO + n3logO(1)nM) time and O(n3log2n· EO +n4logO(1)n)time respectively, where M is the largest absolute value of the function (assuming the range is integers) and is the time taken to evaluate the function on any set. Although the best known lower bound on the query complexity is only Ω(n) [Harvey, 2008], the current shortest non-deterministic proof [Cunningham, 1985] certifying the optimum value of a function requires Ω(n2) function evaluations.

In this paper, we begin to address the longstanding algorithmic gap between general and reversible Markov chains. We develop directed analogues of several spectral graph-theoretic tools that had previously been available only in the undirected setting, and for which it was not clear that directed versions even existed. In particular, we provide a notion of approximation for directed graphs, prove sparsifiers under this notion always exist, and show how to construct them in almost linear time. Using this notion of approximation, we design the first almost-linear-time directed Laplacian system solver, and, by leveraging the recent framework of [Cohen-Kelner-Peebles-Peng-Sidford-Vladu, FOCS '16], we also obtain almost-linear-time algorithms for computing the stationary distribution of a Markov chain, computing expected commute times in a directed graph, and more. For each problem, our algorithms improve the previous best running times of O((nm3/4 + n2/3 m) logO(1) (n κ ε-1)) to O((m + n2O(√lognloglogn)) logO(1) (n κε-1)) where n is the number of vertices in the graph, m is the number of edges, κ is a natural condition number associated with the problem, and ε is the desired accuracy. We hope these results open the door for further studies into directed spectral graph theory, and that they will serve as a stepping stone for designing a new generation of fast algorithms for directed graphs.

We present the first single pass algorithm for computing spectral sparsifiers for graphs in the dynamic semi-streaming model. Given a single pass over a stream containing insertions and deletions of edges to a graph $G$, our algorithm maintains a randomized linear sketch of the incidence matrix of $G$ into dimension $O (\frac{1}{\epsilon^2}n{polylog} (n))$. Using this sketch, at any point, the algorithm can output a $(1 \pm \epsilon)$ spectral sparsifier for $G$ with high probability. While $O (\frac{1}{\epsilon^2} n{polylog}(n))$ space algorithms are known for computing cut sparsifiers in dynamic streams [K. J. Ahn, S. Guha, and A. McGregor, in Proceedings of the 31st ACM Symposium on Principles of Database Systems, 2012, pp. 5--14; A. Goel, M. Kapralov, and I. Post, \hrefhttp://arXiv.org/abs/1203.4900 arXiv:1203.4900, 2002] and spectral sparsifiers in insertion-only streams [J. A. Kelner and A. Levin, Theory Comput. Syst., 53 (2013), pp. 243--262], prior to our work, the best known single pass algorithm for maintaining spectral sparsifiers in dynamic streams required sketches of dimension $\Omega (\frac{1}{\epsilon^2}n^{5/3})$ [K. J. Ahn, S. Guha, and A. McGregor, in Proceedings of the 16th International Workshop on Approximation Algorithms for Combinatorial Optimization Problems, 2013, pp. 1--10]. To achieve our result, we show that using a coarse sparsifier for $G$ and a linear sketch of $G$'s incidence matrix, it is possible to sample edges by effective resistance, obtaining a spectral sparsifier of arbitrary precision. Sampling from the sketch requires a novel application of $\ell_2/\ell_2$ sparse recovery, a natural extension of the $\ell_0$ methods used for cut sparsifiers in [K. J. Ahn, S. Guha, and A. McGregor, in Proceedings of the 31st ACM Symposium on Principles of Database Systems, 2012, pp. 5--14]. Recent work on row sampling for matrix approximation gives a recursive approach for obtaining the required coarse sparsifiers [G. L. Miller and R. Peng, ŭlhttp://arXiv.org/abs/1211.2713v1, 2012]. Under certain restrictions, our approach also extends to the problem of maintaining a spectral approximation for a general matrix $A^\top A$ given a stream of updates to rows in $A$.

In this paper we provide improved guarantees for streaming principal component analysis (PCA). Given A_1, \ldots, A_n∈\mathbbR^d\times d sampled independently from distributions satisfying \mathbbE[A_i] = Σfor Σ\succeq 0, we present an O(d)-space linear-time single-pass streaming algorithm for estimating the top eigenvector of Σ. The algorithm nearly matches (and in certain cases improves upon) the accuracy obtained by the standard batch method that computes top eigenvector of the empirical covariance \frac1n \sum_i ∈[n] A_i as analyzed by the matrix Bernstein inequality. Moreover, to achieve constant accuracy, our algorithm improves upon the best previous known sample complexities of streaming algorithms by either a multiplicative factor of O(d) or 1/\mathrmgap where \mathrmgap is the relative distance between the top two eigenvalues of Σ. We achieve these results through a novel analysis of the classic Oja’s algorithm, one of the oldest and perhaps, most popular algorithms for streaming PCA. We show that simply picking a random initial point w_0 and applying the natural update rule w_i + 1 = w_i + \eta_i A_i w_i suffices for suitable choice of \eta_i. We believe our result sheds light on how to efficiently perform streaming PCA both in theory and in practice and we hope that our analysis may serve as the basis for analyzing many variants and extensions of streaming PCA.

In this paper, we provide faster algorithms for computing variousfundamental quantities associated with random walks on a directedgraph, including the stationary distribution, personalized PageRankvectors, hitting times, and escape probabilities. In particular, ona directed graph with n vertices and m edges, we show how tocompute each quantity in time Õ(m3/4n + mn2/3), wherethe Õ notation suppresses polylog factors in n, the desired accuracy, and the appropriate condition number (i.e. themixing time or restart probability). Our result improves upon the previous fastest running times for these problems, previous results either invoke a general purpose linearsystem solver on a n × n matrix with m non-zero entries, or depend polynomially on the desired error or natural condition numberassociated with the problem (i.e. the mixing time or restart probability). For sparse graphs, we obtain a running time of Õ(n7/4), breaking the O(n2) barrier of the best running time one couldhope to achieve using fast matrix multiplication. We achieve our result by providing a similar running time improvementfor solving directed Laplacian systems, a natural directedor asymmetric analog of the well studied symmetric or undirected Laplaciansystems. We show how to solve such systems in time Õ(m3/4n + mn2/3), and efficiently reduce a broad range of problems to solving Õ(1) directed Laplacian systems on Eulerian graphs. We hope these resultsand our analysis open the door for further study into directedspectral graph theory.

We show how to efficiently project a vector onto the top principal components of a matrix, *without explicitly computing these components*. Specifically, we introduce an iterative algorithm that provably computes the projection using few calls to any black-box routine for ridge regression. By avoiding explicit principal component analysis (PCA), our algorithm is the first with no runtime dependence on the number of top principal components. We show that it can be used to give a fast iterative method for the popular principal component regression problem, giving the first major runtime improvement over the naive method of combining PCA with regression. To achieve our results, we first observe that ridge regression can be used to obtain a "smooth projection" onto the top principal components. We then sharpen this approximation to true projection using a low-degree polynomial approximation to the matrix step function. Step function approximation is a topic of long-term interest in scientific computing. We extend prior theory by constructing polynomials with simple iterative structure and rigorously analyzing their behavior under limited precision.

We give faster algorithms and improved sample complexities for the fundamental problem of estimating the top eigenvector. Given an explicit matrix $A \in \mathbb{R}^{n \times d}$, we show how to compute an $\epsilon$-approximate top eigenvector of $A^TA$ in time $\tilde O\left( \left[\text{nnz}(A) + \frac{d \text{sr}(A)}{\text{gap}^2} \right] \cdot \log 1/\epsilon\right)$. Here $\text{nnz}(A)$ is the number of nonzeros in $A$, $\text{sr}(A)$ is the stable rank, and gap is the relative eigengap. We also consider an online setting in which, given a stream of i.i.d. samples from a distribution D with covariance matrix $\Sigma$ and a vector $x_0$ which is an $O(\text{gap})$ approximate top eigenvector for $\Sigma$, we show how to refine $x_0$ to an $\epsilon$ approximation using $O \left( \frac{\text{var}(\mathcal{D})}{\text{gap}-\epsilon}\right)$ samples from $\mathcal{D}$. Here $\text{var}(\mathcal{D})$ is a natural notion of variance. Combining our algorithm with previous work to initialize $x_0$, we obtain improved sample complexities and runtimes under a variety of assumptions on D. We achieve our results via a robust analysis of the classic shift-and-invert preconditioning method. This technique lets us reduce eigenvector computation to approximately solving a series of linear systems with fast stochastic gradient methods.

This paper considers the problem of canonical-correlation analysis (CCA) and, more broadly, the generalized eigenvector problem for a pair of symmetric matrices. These are two fundamental problems in data analysis and scientific computing with numerous applications in machine learning and statistics. We provide simple iterative algorithms, with improved runtimes, for solving these problems that are globally linearly convergent with moderate dependencies on the condition numbers and eigenvalue gaps of the matrices involved. We obtain our results by reducing CCA to the top-k generalized eigenvector problem. We solve this problem through a general framework that simply requires black box access to an approximate linear system solver. Instantiating this framework with accelerated gradient descent we obtain a running time of \order\fracz k \sqrtκρ \log(1/ε) \log \left(kκ/ρ\right) where z is the total number of nonzero entries, κis the condition number and ρis the relative eigenvalue gap of the appropriate matrices. Our algorithm is linear in the input size and the number of components k up to a \log(k) factor. This is essential for handling large-scale matrices that appear in practice. To the best of our knowledge this is the first such algorithm with global linear convergence. We hope that our results prompt further research and ultimately improve the practical running time for performing these important data analysis procedures on large data sets.

In this paper we provide faster algorithms for solving the geometric median problem: given n points in d compute a point that minimizes the sum of Euclidean distances to the points. This is one of the oldest non-trivial problems in computational geometry yet despite a long history of research the previous fastest running times for computing a (1+є)-approximate geometric median were O(d· n4/3є−8/3) by Chin et. al, Õ(dexpє−4logє−1) by Badoiu et. al, O(nd+poly(d,є−1)) by Feldman and Langberg, and the polynomial running time of O((nd)O(1)log1/є) by Parrilo and Sturmfels and Xue and Ye.

We introduce the notion of balance for directed graphs: a weighted directed graph is α-balanced if for every cut S ⊆ V, the total weight of edges going from S to V∖ S is within factor α of the total weight of edges going from V∖ S to S. Several important families of graphs are nearly balanced, in particular, Eulerian graphs (with α = 1) and residual graphs of (1+є)-approximate undirected maximum flows (with α=O(1/є)).

In many estimation problems, e.g. linear and logistic regression, we wish to minimize an unknown objective given only unbiased samples of the objective function. Furthermore, we aim to achieve this using as few samples as possible. In the absence of computational constraints, the minimizer of a sample average of observed data – commonly referred to as either the empirical risk minimizer (ERM) or the M-estimator – is widely regarded as the estimation strategy of choice due to its desirable statistical convergence properties. Our goal in this work is to perform as well as the ERM, on \emphevery problem, while minimizing the use of computational resources such as running time and space usage. We provide a simple streaming algorithm which, under standard regularity assumptions on the underlying problem, enjoys the following properties: \beginenumerate \item The algorithm can be implemented in linear time with a single pass of the observed data, using space linear in the size of a single sample. \item The algorithm achieves the same statistical rate of convergence as the empirical risk minimizer on every problem, even considering constant factors. \item The algorithm’s performance depends on the initial error at a rate that decreases super-polynomially. \item The algorithm is easily parallelizable. \endenumerate Moreover, we quantify the (finite-sample) rate at which the algorithm becomes competitive with the ERM.

In this paper, we consider the following inverse maintenance problem: given A ∈ Rn×dand a number of rounds r, at round k, we receive a n x n diagonal matrix D(k)and we wish to maintain an efficient linear system solver for ATD(k)A under the assumption D(k)does not change too rapidly. This inverse maintenance problem is the computational bottleneck in solving multiple optimization problems. We show how to solve this problem with Õ (nnz(A) + dω) preprocessing time and amortized Õ(nnz(A) + d2) time per round, improving upon previous running times. Consequently, we obtain the fastest known running times for solving multiple problems including, linear programming and computing a rounding of a polytope. In particular given a feasible point in a linear program with n variables, d constraints, and constraint matrix A ∈ Rd×n, we show how to solve the linear program in time Õ((nnz(A) + d2)√ d log(∈-1)). We achieve our results through a novel combination of classic numerical techniques of low rank update, preconditioning, and fast matrix multiplication as well as recent work on subspace embeddings and spectral sparsification that we hope will be of independent interest.

In this paper we improve upon the running time for finding a point in a convex set given a separation oracle. In particular, given a separation oracle for a convex set K ⊂ Rnthat is contained in a box of radius R we show how to either compute a point in K or prove that K does not contain a ball of radius ϵ using an expected O(n log(nR/ϵ)) evaluations of the oracle and additional time O(n3logO(1)(nR/ϵ)). This matches the oracle complexity and improves upon the O(nω+1log(nR/ϵ)) additional time of the previous fastest algorithm achieved over 25 years ago by Vaidya [91] for the current value of the matrix multiplication constant w2log nM · EO + n3logO(1)nM) and O(n3log2n · EO + n4logO(1)n), improving upon the previous best of O((n4· EO + n5)logM) and O(n5· EO + n6) respectively. · Submodular Flow: n = |V|, m = |E|, C is the maximum edge cost in absolute value and U is maximum edge capacity in absolute value. We obtain a faster weakly polynomial running time of O(n2log nCU · EO + n3logO(1) nCU), improving upon the previous best of O(mn5log nU · EO) and O (n4h min {log C, log U}) from 15 years ago by a factor of Õ(n4). We also achieve faster strongly polynomial time algorithms as a consequence of our result on submodular minimization. · Matroid Intersection: n is the size of the ground set, r is the maximum size of independent sets, M is the maximum absolute value of element weight, Trankand Tindare the time for each rank and independence oracle query. We obtain a running time of O((nr log2nTrank+n3logO(1)n) log nM) and O((n2log nTind+n3logO(1)n) log nM), achieving the first quadratic bound on the query complexity for the independence and rank oracles. In the unweighted case, this is the first improvement since 1986 for independence oracle. · Semidefinite Programming: n is the number of constraints, m is the number of dimensions and S is the total number of non-zeros in the constraint matrices. We obtain a running time of O(n(n2+ mω+ S)), improving upon the previous best of Õ(n(nω+ mω+ S)) for the regime S is small.

We develop a family of accelerated stochastic algorithms that optimize sums of convex functions. Our algorithms improve upon the fastest running time for empirical risk minimization (ERM), and in particular linear least-squares regression, across a wide range of problem settings. To achieve this, we establish a framework, based on the classical proximal point algorithm, useful for accelerating recent fast stochastic algorithms in a black-box fashion. Empirically, we demonstrate that the resulting algorithms exhibit notions of stability that are advantageous in practice. Both in theory and in practice, the provided algorithms reap the computational benefits of adding a large strongly convex regularization term, without incurring a corresponding bias to the original ERM problem.

Random sampling has become a critical tool in solving massive matrix problems. For linear regression, a small, manageable set of data rows can be randomly selected to approximate a tall, skinny data matrix, improving processing time significantly. For theoretical performance guarantees, each row must be sampled with probability proportional to its statistical leverage score. Unfortunately, leverage scores are difficult to compute. A simple alternative is to sample rows uniformly at random. While this often works, uniform sampling will eliminate critical row information for many natural instances.

We present the first single pass algorithm for computing spectral sparsifiers of graphs in the dynamic semi-streaming model. Given a single pass over a stream containing insertions and deletions of edges to a graph, G, our algorithm maintains a randomized linear sketch of the incidence matrix into dimension O(1/∈2npolylog(n)). Using this sketch, the algorithm can output a (1±∈) spectral sparsifier for G with high probability. While O(1/∈2npolylog(n)) space algorithms are known for computing cut sparsifiers in dynamic streams [1], [2] and spectral sparsifiers in insertion-only streams [3], prior to our work, the best known single pass algorithm for maintaining spectral sparsifiers in dynamic streams required sketches of dimension Ω(1/∈2n5/3). To achieve our result, we show that, using a coarse sparsifier of G and a linear sketch of G's incidence matrix, it is possible to sample edges by effective resistance, obtaining a spectral sparsifier of arbitrary precision. Sampling from the sketch requires a novel application of ℓ2/ℓ2sparse recovery, a natural extension of the ℓ0methods used for cut sparsifiers in [1]. Recent work of [2] on row sampling for matrix approximation gives a recursive approach for obtaining the required coarse sparsifiers. Under certain restrictions, our approach also extends to the problem of maintaining a spectral approximation for a general matrix ATA given a stream of updates to rows in A.

In this paper, we present a new algorithm for '/ solving linear programs that requires only Õ(√rank(A)L) iterations where A is the constraint matrix of a linear program with m constraints, n variables, and bit complexity L. Each iteration of our method consists of solving Õ(1) linear systems and additional nearly linear time computation. Our method improves upon the previous best iteration bounds by factor Ω̃((m/rank (A)))1/4) of for methods with polynomial time computable iterations and by Ω̃((m/rank (A))1/2) for methods which solve at most Õ(1) linear systems in each iteration each achieved over 20 years ago. Applying our techniques to the linear program formulation of maximum flow yields an Õ(|E| √|V| log2U) time algorithm for solving the maximum flow problem on directed graphs with |E| edges, |V| vertices, and capacity ratio U. This improves upon the previous fastest running time of O(|E| min{|E|1/2, |V|2/3} log (|V|2/|E|) log(U)) achieved over 15 years ago by Goldberg and Rao and improves upon the previous best running times for solving dense directed unit capacity graphs of Õ(|E| min{|E|1/2, |V|2/3}) achieved by Even and Tarjan over 35 years ago and a running time of Õ(|E|10/7) achieved recently by Madry.

In this paper, we introduce a new framework for approximately solving flow problems in capacitated, undirected graphs and apply it to provide asymptotically faster algorithms for the maximum s-t flow and maximum concurrent multicommodity flow problems. For graphs with n vertices and m edges, it allows us to find an ∊-approximate maximum s-t flow in time O(m1+o(1)∊−2), improving on the previous best bound of Õ(mn1/3poly(∊−1)). Applying the same framework in the multicommodity setting solves a maximum concurrent multicommodity flow problem with k commodities in O(m1+o(1)∊−2k2) time, improving on the existing bound of Õ(m4/3poly(k, ∊−1)). Our algorithms utilize several new technical tools that we believe may be of independent interest: We give a non-Euclidean generalization of gradient descent and provide bounds on its performance. Using this, we show how to reduce approximate maximum flow and maximum concurrent flow to oblivious routing. We define and provide an efficient construction of a new type of flow sparsifier. Previous sparsifier constructions approximately preserved the size of cuts and, by duality, the value of the maximum flows as well. However, they did not provide any direct way to route flows in the sparsifier G′ back in the original graph G, leading to a longstanding gap between the efficacy of sparsification on flow and cut problems. We ameliorate this by constructing a sparsifier G' that can be embedded (very efficiently) into G with low congestion, allowing one to transfer flows from G′ back to G. We give the first almost-linear-time construction of an O(mo(1))-competitive oblivious routing scheme. No previous such algorithm ran in time better than . By reducing the running time to almost-linear, our work provides a powerful new primitive for constructing very fast graph algorithms. The interested reader is referred to the full version of the paper [8] for a more complete treatment of these results.

In this paper we show how to accelerate randomized coordinate descent methods and achieve faster convergence rates without paying per-iteration costs in asymptotic running time. In particular, we show how to generalize and efficiently implement a method proposed by Nesterov, giving faster asymptotic running times for various algorithms that use standard coordinate descent as a black box. In addition to providing a proof of convergence for this new general method, we show that it is numerically stable, efficiently implementable, and in certain regimes, asymptotically optimal. To highlight the power of this algorithm, we show how it can used to create faster linear system solvers in several regimes: - We show how this method achieves a faster asymptotic runtime than conjugate gradient for solving a broad class of symmetric positive definite systems of equations. - We improve the convergence guarantees for Kaczmarz methods, a popular technique for image reconstruction and solving over determined systems of equations, by accelerating an algorithm of Strohmer and Vershynin. - We achieve the best known running time for solving Symmetric Diagonally Dominant (SDD) system of equations in the unit-cost RAM model, obtaining a running time of O(m log3/2n (log log n)1/2 log((log n)/eps)) by accelerating a recent solver by Kelner et al. Beyond the independent interest of these solvers, we believe they highlight the versatility of the approach of this paper and we hope that they will open the door for further algorithmic improvements in the future.

In this paper, we present a simple combinatorial algorithm that solves symmetric diagonally dominant (SDD) linear systems in nearly-linear time. It uses little of the machinery that previously appeared to be necessary for a such an algorithm. It does not require recursive preconditioning, spectral sparsification, or even the Chebyshev Method or Conjugate Gradient. After constructing a "nice" spanning tree of a graph associated with the linear system, the entire algorithm consists of the repeated application of a simple update rule, which it implements using a lightweight data structure. The algorithm is numerically stable and can be implemented without the increased bit-precision required by previous solvers. As such, the algorithm has the fastest known running time under the standard unit-cost RAM model. We hope the simplicity of the algorithm and the insights yielded by its analysis will be useful in both theory and practice.