We present an algorithm that computes exact maximum flows and minimum-cost flows on directed graphs with m edges and polynomially bounded integral demands, costs, and capacities in \(m^{1+o(1)}\) time. Our algorithm builds the flow through a sequence of \(m^{1+o(1)}\) approximate undirected minimum-ratio cycles, each of which is computed and processed in amortized \(m^{o(1)}\) time using a new dynamic graph data structure. Our framework extends to algorithms running in \(m^{1+o(1)}\) time for computing flows that minimize general edge-separable convex functions to high accuracy. This gives almost-linear time algorithms for several problems including entropy-regularized optimal transport, matrix scaling, p -norm flows, and p -norm isotonic regression on arbitrary directed acyclic graphs.

We present a nearly-linear time algorithm for finding a minimum-cost flow in planar graphs with polynomially-bounded integer costs and capacities. The previous fastest algorithm for this problem is based on interior point methods (IPMs) and works for general sparse graphs in O ( n 1.5 ⋅ poly (log n )) time [Daitch-Spielman, STOC’08]. Intuitively, Ω ( n 1.5 ) is a natural runtime barrier for IPM-based methods, since they require \(\sqrt {n}\) iterations, each routing a possibly-dense electrical flow. To break this barrier, we develop a new implicit representation for flows based on generalized nested dissection [Lipton-Rose-Tarjan, SINUM’79] and approximate Schur complements [Kyng-Sachdeva, FOCS’16]. This implicit representation permits us to design a data structure to route an electrical flow with sparse demands in roughly \(\sqrt {n}\) update time, resulting in a total runtime of O ( n ⋅ poly (log n )). Our results immediately extend to all families of separable graphs.

Can linear systems be solved faster than matrix multiplication? Although there has been much progress on systems with additional structures, in the general setting, the complexity of solving an n × n linear system Ax = b is at least Ω ( n ω ) bit operations, where ω < 2.372 is the matrix multiplication exponent. Improving on this bound of n < has been an open problem even for sparse linear systems with poly( n ) condition number. In this article, we present an algorithm that solves linear systems in sparse matrices asymptotically faster than matrix multiplication or any ω > 2. This speedup holds for any input matrix A with o ( n ω-1 /log (κ (A))) nonzeros, where κ ( A ) is the condition number of A . For poly( n )-conditioned matrices with Õ( n ) nonzeros, and the current value of ω, the bit complexity of our algorithm to solve to within any 1/poly( n ) error is O ( n 2.331 ). Our algorithm can be viewed as an efficient, randomized implementation of the block Krylov method via recursive low displacement rank factorizations. It is inspired by the algorithm of [Eberly et al. ISSAC ‘06 ‘07] for inverting matrices over finite fields. In our analysis of numerical stability, we use matrix anti-concentration techniques to bound the smallest eigenvalue and the smallest gap in eigenvalues of semi-random matrices. Incorporating the matrix anti-concentration bounds by [Nie STOC‘22] gives an improved bit complexity of O ( n 2.271 ).

Bipartite graphs are ubiquitous in many domains, e.g., e-commerce platforms, social networks, and academia, by modeling interactions between distinct entity sets. Within these graphs, the butterfly motif, a complete 2×2 biclique, represents the simplest yet significant subgraph structure, crucial for analyzing complex network patterns. Counting the butterflies offers significant benefits across various applications, including community analysis and recommender systems. Additionally, the temporal dimension of bipartite graphs, where edges activate within specific time frames, introduces the concept of historical butterfly counting, i.e., counting butterflies within a given time interval. This temporal analysis sheds light on the dynamics and evolution of network interactions, offering new insights into their mechanisms. Despite its importance, no existing algorithm can efficiently solve the historical butterfly counting task. To address this, we design two novel indices whose memory footprints are dependent on #butterflies and #wedges, respectively. Combining these indices, we propose a graph structure-aware indexing approach that significantly reduces memory usage while preserving exceptional query speed. To further reduce the index size and boost the query efficiency, we design an index compression strategy, enabling the fast, high-quality, and unbiased approximation of historical butterfly counts. We theoretically prove that our approach is particularly advantageous on power-law graphs, a common characteristic of real-world bipartite graphs, by surpassing traditional complexity barriers for general graphs. Extensive experiments reveal that our query algorithms outperform existing methods by up to five magnitudes, effectively balancing speed with manageable memory requirements.

Signed networks, characterized by edges labeled as either positive or negative, offer nuanced insights into interaction dynamics beyond the capabilities of unsigned graphs. Central to this is the task of identifying the maximum balanced subgraph, crucial for applications like polarized community detection in social networks and portfolio analysis in finance. Traditional models, however, are limited by an assumption of perfect partitioning, which fails to mirror the complexities of real-world data. Addressing this gap, we introduce an innovative generalized balanced subgraph model that incorporates tolerance for imbalance. Our proposed region-based heuristic algorithm, tailored for this NP -hard problem, strikes a balance between low time complexity and high-quality outcomes. Comparative experiments validate its superior performance against leading solutions, delivering enhanced effectiveness (notably larger subgraph sizes) and efficiency (achieving up to 100× speedup) in both traditional and generalized contexts.

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 the first dynamic distance oracles for interval graphs, which are intersection graphs of a set of intervals on the real line, and for proper interval graphs, which are intersection graphs of a set of intervals in which no interval is properly contained in another. For proper interval graphs, we design a linear space data structure which supports distance queries (computing the distance between two query vertices) and vertex insertion or deletion in O(lg⁡n) worst-case time, where n is the number of vertices currently in G. Under incremental (insertion only) or decremental (deletion only) settings in general interval graphs, we design linear space data structures that support distance queries in O(lg⁡n) worst-case time and vertex insertion or deletion in O(lg⁡n) amortized time, where n is the maximum number of vertices in the graph. Under fully dynamic settings in general interval graphs, we design a data structure that represents an interval graph G in O(n) words of space to support distance queries in O(nlg⁡n/S(n)) worst-case time and vertex insertion or deletion in O(S(n)+lg⁡n) worst-case time, where n is the number of vertices currently in G and S(n) is an arbitrary function that satisfies S(n)=Ω(1) and S(n)=O(n). This implies an O(n)-word solution with O(nlg⁡n)-time support for both distance queries and updates. All four data structures can answer shortest path queries by reporting the vertices in the shortest path between two query vertices in O(lg⁡n) worst-case time per vertex. We also study the hardness of supporting distance queries under updates over an intersection graph of 3D axis-aligned line segments, which generalizes our problem to 3D. Finally, we solve the problem of computing the diameter of a dynamic connected interval graph.

The \(\ell _p\) -norm regression problem is a classic problem in optimization with wide ranging applications in machine learning and theoretical computer science. The goal is to compute \(\boldsymbol {\mathit {x}}^{\star } =\arg \min _{\boldsymbol {\mathit {A}}\boldsymbol {\mathit {x}}=\boldsymbol {\mathit {b}}}\Vert \boldsymbol {\mathit {x}}\Vert _p^p\) , where \(\boldsymbol {\mathit {x}}^{\star }\in \mathbb {R}^n,\boldsymbol {\mathit {A}}\in \mathbb {R}^{d\times n},\boldsymbol {\mathit {b}}\in \mathbb {R}^d\) and \(d\le n\) . Efficient high-accuracy algorithms for the problem have been challenging both in theory and practice and the state-of-the-art algorithms require \(poly(p)\cdot n^{\frac{1}{2}-\frac{1}{p}}\) linear system solves for \(p\ge 2\) . In this article, we provide new algorithms for \(\ell _p\) -regression (and a more general formulation of the problem) that obtain a high-accuracy solution in \(O(p n^{ {(p-2)}{(3p-2)}})\) linear system solves. We further propose a new inverse maintenance procedure that speeds-up our algorithm to \(\widetilde{O}(n^{\omega })\) total runtime, where \(O(n^{\omega })\) denotes the running time for multiplying \(n \times n\) matrices. Additionally, we give the first Iteratively Reweighted Least Squares (IRLS) algorithm that is guaranteed to converge to an optimum in a few iterations. Our IRLS algorithm has shown exceptional practical performance, beating the currently available implementations in MATLAB/CVX by 10–50×.

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 analyze the bit complexity of efficient algorithms for fundamental optimization problems, such as linear regression, p-norm regression, and linear programming (LP). State-of-the-art algorithms are iterative, and in terms of the number of arithmetic operations, they match the current time complexity of multiplying two n-by-n matrices (up to polylogarithmic factors). However, previous work has typically assumed infinite precision arithmetic, and due to complicated inverse maintenance techniques, the actual running times of these algorithms are unknown. To settle the running time and bit complexity of these algorithms, we demonstrate that a core common subroutine, known as inverse maintenance, is backward-stable. Additionally, we show that iterative approaches for solving constrained weighted regression problems can be accomplished with bounded-error preconditioners. Specifically, we prove that linear programs can be solved approximately in matrix multiplication time multiplied by polylog factors that depend on the condition number $\kappa$ of the matrix and the inner and outer radius of the LP problem. p-norm regression can be solved approximately in matrix multiplication time multiplied by polylog factors in $\kappa$. Lastly, linear regression can be solved approximately in input-sparsity time multiplied by polylog factors in $\kappa$. Furthermore, we present results for achieving lower than matrix multiplication time for p-norm regression by utilizing faster solvers for sparse linear systems.

Over the last two decades, a significant line of work in theoretical algorithms has made progress in solving linear systems of the form 𝐋𝐱 = 𝐛, where 𝐋 is the Laplacian matrix of a weighted graph with weights w(i,j) > 0 on the edges. The solution 𝐱 of the linear system can be interpreted as the potentials of an electrical flow in which the resistance on edge (i,j) is 1/w(i,j). Kelner, Orrechia, Sidford, and Zhu [Kelner et al., 2013] give a combinatorial, near-linear time algorithm that maintains the Kirchoff Current Law, and gradually enforces the Kirchoff Potential Law by updating flows around cycles (cycle toggling). In this paper, we consider a dual version of the algorithm that maintains the Kirchoff Potential Law, and gradually enforces the Kirchoff Current Law by cut toggling: each iteration updates all potentials on one side of a fundamental cut of a spanning tree by the same amount. We prove that this dual algorithm also runs in a near-linear number of iterations. We show, however, that if we abstract cut toggling as a natural data structure problem, this problem can be reduced to the online vector-matrix-vector problem (OMv), which has been conjectured to be difficult for dynamic algorithms [Henzinger et al., 2015]. The conjecture implies that the data structure does not have an O(n^{1-ε}) time algorithm for any ε > 0, and thus a straightforward implementation of the cut-toggling algorithm requires essentially linear time per iteration. To circumvent the lower bound, we batch update steps, and perform them simultaneously instead of sequentially. An appropriate choice of batching leads to an Õ(m^{1.5}) time cut-toggling algorithm for solving Laplacian systems. Furthermore, we show that if we sparsify the graph and call our algorithm recursively on the Laplacian system implied by batching and sparsifying, we can reduce the running time to O(m^{1 + ε}) for any ε > 0. Thus, the dual cut-toggling algorithm can achieve (almost) the same running time as its primal cycle-toggling counterpart.

We develop a framework for graph sparsification and sketching, based on a new tool, short cycle decomposition, which is a decomposition of an unweighted graph into an edge-disjoint collection of short cycles, plus a small number of extra edges. A simple observation shows that every graph $G$ on $n$ vertices with $m$ edges can be decomposed in $O(mn)$ time into cycles of length at most $2 \log n,$ and at most $2n$ extra edges. We give an $(m^{1+o(1)})$-time algorithm for constructing a short cycle decomposition, with cycles of length $n^{o(1)},$ and $n^{1+o(1)}$ extra edges. Both the existential and algorithmic variants of this decomposition enable us to make the following progress on several open problems in randomized graph algorithms: (1) We present an algorithm that runs in time $m^{1+o(1)}\varepsilon^{-1.5}$ and returns $(1\pm\varepsilon)$-approximations to effective resistances of all edges, improving over the previous best runtime of $\widetilde{{O}}(\min\{m\varepsilon^{-2}, n^{2} \varepsilon^{-1}\})$. This routine in turn gives an algorithm for approximating the determinant of a graph Laplacian up to a factor of $(1\pm \varepsilon)$ in $m^{1 + o(1)} + n^{\nicefrac{15}{8}+o(1)}\varepsilon^{-\nicefrac{7}{4}}$ time. (2) We show the existence of graphical spectral sketches with about $n\varepsilon^{-1}$ edges, and also give efficient algorithms to construct them. A graphical spectral sketch is a distribution over sparse graphs $H$ such that for a fixed vector ${\mathit{x}}$, we have ${{x}}^{\top} {L}_H {{x}} = (1\pm\varepsilon) {{x}}^{\top} {L}_G {{x}}$ and ${{x}}^{\top} {L}^{+}_H {{x}} = (1\pm\varepsilon) {{x}}^{\top} {L}^{+}_G {{x}}$ with high probability, where ${L}$ is the graph Laplacian and ${L}^{+}$ is its pseudoinverse. This implies the existence of resistance sparsifiers with about $n \varepsilon^{-1}$ edges that preserve the effective resistance between every pair of vertices up to $(1\pm\varepsilon)$. (3) By combining short cycle decompositions with known tools in graph sparsification, we show the existence of nearly linear sized degree-preserving spectral sparsifiers, as well as significantly sparser approximations of Eulerian directed graphs. The latter is critical to recent breakthroughs on faster algorithms for solving linear systems in directed Laplacians. The running time and output qualities of our spectral sketch and degree-preserving (directed) sparsification algorithms are limited by the efficiency of our routines for constructing short cycle decompositions. Improved algorithms for short cycle decompositions will lead to improvement in each of these algorithms.

We proposeLightNE 2.0, a cost-effective, scalable, automated, and high-quality network embedding system that scales to graphs with hundreds of billions of edges on a single machine. In contrast to the mainstream belief that distributed architecture and GPUs are needed for large-scale network embedding with good quality, we prove that we can achieve higher quality, better scalability, lower cost, and faster runtime with shared-memory, CPU-only architecture.LightNE 2.0combines two theoretically grounded embedding methods NetSMF and ProNE. We introduce the following techniques to network embedding for the first time: (1) a newly proposed downsampling method to reduce the sample complexity of NetSMF while preserving its theoretical advantages; (2) a high-performance parallel graph processing stack GBBS to achieve high memory efficiency and scalability; (3) sparse parallel hash table to aggregate and maintain the matrix sparsifier in memory; (4) a fast randomized singular value decomposition (SVD) enhanced by power iteration and fast orthonormalization to improve vanilla randomized SVD in terms of both efficiency and effectiveness; (5) Intel MKL for proposed fast randomized SVD and spectral propagation; and (6) a fast and lightweight AutoML library FLAML for automated hyperparameter tuning. Experimental results show thatLightNE 2.0can be up to 84× faster than GraphVite, 30× faster than PBG and 9× faster than NetSMF while delivering better performance.LightNE 2.0can embed very large graph with 1.7 billion nodes and 124 billion edges in half an hour on a CPU server, while other baselines cannot handle very large graphs of this scale.

We give an algorithm that computes exact maximum flows and minimum-cost flows on directed graphs with m edges and polynomially bounded integral demands, costs, and capacities in $m^{1+o(1)}$ time. Our algorithm builds the flow through a sequence of $m^{1+o(1)}$ approximate undirected minimum-ratio cycles, each of which is computed and processed in amortized $m^{o(1)}$ time using a new dynamic graph data structure. Our framework extends to algorithms running in $m^{1+o(1)}$ time for computing flows that minimize general edge-separable convex functions to high accuracy. This gives almost-linear time algorithms for several problems including entropy-regularized optimal transport, matrix scaling, p-norm flows, and p-norm isotonic regression on arbitrary directed acyclic graphs.

We present a nearly-linear time algorithm for finding a minimum-cost flow in planar graphs with polynomially bounded integer costs and capacities. The previous fastest algorithm for this problem was based on interior point methods (IPMs) and worked for general sparse graphs in O(n1.5 poly(log n)) time [Daitch-Spielman, STOC'08]. Intuitively, Ω(n1.5) is a natural runtime barrier for IPM based methods, since they require iterations, each routing a possibly-dense electrical flow. To break this barrier, we develop a new implicit representation for flows based on generalized nested-dissection [Lipton-Rose-Tarjan, JSTOR'79] and approximate Schur complements [Kyng-Sachdeva, FOCS'16]. This implicit representation permits us to design a data structure to route an electrical flow with sparse demands in roughly update time, resulting in a total running time of O(n · poly(log n)). Our results immediately extend to all families of separable graphs.

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

We show that the sparsified block elimination algorithm for solving undirected Laplacian linear systems from [Kyng-Lee-Peng-Sachdeva-Spielman STOC’16] directly works for directed Laplacians. Given access to a sparsification algorithm that, on graphs with n vertices and m edges, takes time TS(m) to output a sparsifier with NS(n) edges, our algorithm solves a directed Eulerian system on n vertices and m edges to є relative accuracy in time O(TS(m) + NS(n)lognlog(n/є)) + Õ(TS(NS(n)) logn), where the Õ(·) notation hides loglog(n) factors. By previous results, this implies improved runtimes for linear systems in strongly connected directed graphs, PageRank matrices, and asymmetric M-matrices. When combined with slower constructions of smaller Eulerian sparsifiers based on short cycle decompositions, it also gives a solver algorithm that, after pre-processing the matrix in O(n2 logO(1) n) time, takes O(n log5n log(n / є)) time per solve. At the core of our analyses are constructions of augmented matrices whose Schur complements encode error matrices.

Diffusion is a fundamental graph procedure and has been a basic building block in a wide range of theoretical and empirical applications such as graph partitioning and semi-supervised learning on graphs. In this paper, we study computationally efficient diffusion primitives beyond random walk. We design an near-linear time randomized algorithm for the 2-norm flow diffusion problem, a recently proposed diffusion model based on network flow with demonstrated graph clustering related applications both in theory and in practice. Examples include finding locally-biased low conductance cuts. Using a known connection between the optimal dual solution of the flow diffusion problem and the local cut structure, our algorithm gives an alternative approach for finding such cuts in nearly linear time. From a technical point of view, our algorithm contributes a novel way of dealing with inequality constraints in graph optimization problems. It adapts the high-level algorithmic framework of nearly linear time Laplacian system solvers, but requires several new tools: vertex elimination under constraints, a new family of graph ultra-sparsifiers, and ac-celerated proximal gradient methods with inexact proximal mapping computation. See https://arxiv.org/abs/2105.14629 for the full version of this paper.

We study distributed algorithms built around minor-based vertex sparsifiers, and give the first algorithm in the CONGEST model for solving linear systems in graph Laplacian matrices to high accuracy. Our Laplacian solver has a round complexity of$O(n^{o(1)}(\sqrt{n}+D))$, and thus almost matches the lower bound of$\widetilde{\Omega}(\sqrt{n}+D)$, where$n$is the number of nodes in the network and$D$is its diameter. We show that our distributed solver yields new sublinear round algorithms for several cornerstone problems in combinatorial optimization. This is achieved by leveraging the powerful algorithmic framework of Interior Point Methods (IPMs) and the Laplacian paradigm in the context of distributed graph algorithms, which entails numerically solving optimization problems on graphs via a series of Laplacian systems. Problems that benefit from our distributed algorithmic paradigm include exact mincost flow, negative weight shortest paths, maxflow, and bipartite matching on sparse directed graphs. For the maxflow problem, this is the first exact distributed algorithm that applies to directed graphs, while the previous work by [Ghaffari et al. SICOMP'18] considered the approximate setting and works only for undirected graphs. For the mincost flow and the negative weight shortest path problems, our results constitute the first exact distributed algorithms running in a sublinear number of rounds. Given that the hybrid between IPMs and the Laplacian paradigm has proven useful for tackling numerous optimization problems in the centralized setting, we believe that our distributed solver will find future applications. At the heart of our distributed Laplacian solver is the notion of spectral subspace sparsifiers of [Li, Schild FOCS'18]. We present a nontrivial distributed implementation of their construction by (i) giving a parallel variant of their algorithm that avoids the sampling of random spanning trees and uses approximate leverage scores instead, and (ii) showing that the algorithm still produces a high-quality subspace spectral sparsifier by carefully setting up and analyzing matrix martingales. Combining this vertex reduction recursively with both tree and elimination-based preconditioners leads to our algorithm for solving Laplacian systems. The construction of the elimination-based preconditioners is based on computing short random walks, and we introduce a new technique for reducing the congestion incurred by the simulation of these walks on weighted graphs.

We give an algorithm for computing exact maximum flows on graphs with$m$edges and integer capacities in the range [$1,U$] in$\tilde{O}(m^{\frac{3}{2}-\frac{1}{328}}\log U)$time.11We use$\tilde{O}(\cdot)$to suppress logarithmic factors in$m$. For sparse graphs with polynomially bounded integer capacities, this is the first improvement over the$\tilde{O}(m^{1.5}\log U)$time bound from [Goldberg-Rao JACM '98]. Our algorithm revolves around dynamically maintaining the augmenting electrical flows at the core of the interior point method based algorithm from [Mądry JACM '16]. This entails designing data structures that, in limited settings, return edges with large electric energy in a graph undergoing resistance updates.

We propose LightNE, a cost-effective, scalable, and high quality network embedding system that scales to graphs with hundreds of billions of edges on a single machine. In contrast to the mainstream belief that distributed architecture and GPUs are needed for large-scale network embedding with good quality, we prove that we can achieve higher quality, better scalability, lower cost and faster runtime with shared-memory, CPU-only architecture. LightNE combines two theoretically grounded embedding methods NetSMF and ProNE. We introduce the following techniques to network embedding for the first time: (1) a newly proposed downsampling method to reduce the sample complexity of NetSMF while preserving its theoretical advantages; (2) a high-performance parallel graph processing stack GBBS to achieve high memory efficiency and scalability; (3) sparse parallel hash table to aggregate and maintain the matrix sparsifier in memory; and (4) Intel MKL for efficient randomized SVD and spectral propagation.

Graph compression or sparsification is a basic information-theoretic and computational question. A major open problem in this research area is whether (1 + ∊)-approximate cut-preserving vertex sparsifiers with size close to the number of terminals exist. As a step towards this goal, we study a thresholded version of the problem: for a given parameter c, find a smaller graph, which we call connectivity-c mimicking network, which preserves connectivity among k terminals exactly up to the value of c. We show that connectivity-c mimicking networks with O(kc4) edges exist and can be found in time m(c log n)O(c). We also give a separate algorithm that constructs such graphs with k · O(c)2c edges in time mcO(c) logO(1) n. These results lead to the first data structures for answering fully dynamic offline c-edge-connectivity queries for c ≥ 4 in polylogarithmic time per query, as well as more efficient algorithms for survivable network design on bounded treewidth graphs.

Can linear systems be solved faster than matrix multiplication? While there has been remarkable progress for the special cases of graph structured linear systems, in the general setting, the bit complexity of solving an n × n linear system Ax = b is Õ(nω), where ω < 2.372864 is the matrix multiplication exponent. Improving on this has been an open problem even for sparse linear systems with poly(n) condition number. In this paper, we present an algorithm that solves linear systems in sparse matrices asymptotically faster than matrix multiplication for any ω > 2. This speedup holds for any input matrix A with o(nω–1/log(κ(A))) non-zeros, where κ(A) is the condition number of A. For poly(n)-conditioned matrices with Õ(n) nonzeros, and the current value of ω, the bit complexity of our algorithm to solve to within any 1/poly(n) error is O(n2.331645). Our algorithm can be viewed as an efficient, randomized implementation of the block Krylov method via recursive low displacement rank factorizations. It is inspired by the algorithm of [Eberly et al. ISSAC ‘06 ‘07] for inverting matrices over finite fields. In our analysis of numerical stability, we develop matrix anti-concentration techniques to bound the smallest eigenvalue and the smallest gap in eigenvalues of semi-random matrices.

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 present a general framework of designing efficient dynamic approximate algorithms for optimization problems on undirected graphs. In particular, we develop a technique that, given any problem that admits a certain notion of vertex sparsifiers, gives data structures that maintain approximate solutions in sub-linear update and query time. We illustrate the applicability of our paradigm to the following problems. (1)A fully-dynamic algorithm that approximates all-pair maximum-flows/minimum-cuts up to a nearly logarithmic factor in ~O(n2/3)11The ~O(·) notation is used in this paper to hide poly-logarithmic factors. amortized time against an oblivious adversary, and ~O(m3/4) time against an adaptive adversary. (2)An incremental data structure that maintains O(1) - approximate shortest path in no(1)time per operation, as well as fully dynamic approximate all-pair shortest path and transshipment in ~O(n2/3+o(1)) amortized time per operation. (3)A fully-dynamic algorithm that approximates all-pair effective resistance up to an ( 1+ε) factor in ~O(n2/3+o(1)ε-O(1)) amortized update time per operation. The key tool behind result (1) is the dynamic maintenance of an algorithmic construction due to Madry [FOCS' 10], which partitions a graph into a collection of simpler graph structures (known as j-trees) and approximately captures the cut-flow and metric structure of the graph. The O(1)-approximation guarantee of (2) is by adapting the distance oracles by [Thorup-Zwick JACM '05]. Result (3) is obtained by invoking the random-walk based spectral vertex sparsifier by [Durfee et al. STOC '19] in a hierarchical manner, while carefully keeping track of the recourse among levels in the hierarchy. See https://arxiv.org/pdf/2005.02368.pdf for the full version of this paper.

We consider the classical Minimum Balanced Cut problem: given a graph G, compute a partition of its vertices into two subsets of roughly equal volume, while minimizing the number of edges connecting the subsets. We present the first deterministic, almost-linear time approximation algorithm for this problem. Specifically, our algorithm, given an n-vertex m-edge graph G and any parameter 1 ≤ r ≤ O(logn), computes a (logm)r2-approximation for Minimum Balanced Cut in G, in time O(m1+O(1/r)+o(1)·(logm)O(r2)). In particular, we obtain a (logm)1/ε-approximation in time m1+O(√{ε})for any constant , and a (logm)f(m)-approximation in time m1+o(1), for any slowly growing function f(m). We obtain deterministic algorithms with similar guarantees for the Sparsest Cut and the Lowest-Conductance Cut problems. Our algorithm for the Minimum Balanced Cut problem in fact provides a stronger guarantee: it either returns a balanced cut whose value is close to a given target value, or it certifies that such a cut does not exist by exhibiting a large subgraph of G that has high conductance. We use this algorithm to obtain deterministic algorithms for dynamic connectivity and minimum spanning forest, whose worst-case update time on an n-vertex graph is no(1), thus resolving a major open problem in the area of dynamic graph algorithms. Our work also implies deterministic algorithms for a host of additional problems, whose time complexities match, up to subpolynomial in n factors, those of known randomized algorithms. The implications include almost-linear time deterministic algorithms for solving Laplacian systems and for approximating maximum flows in undirected graphs.

Graph embeddings are a ubiquitous tool for machine learning tasks, such as node classification and link prediction, on graph-structured data. However, computing the embeddings for large-scale graphs is prohibitively inefficient even if we are interested only in a small subset of relevant vertices. To address this, we present an efficient graph coarsening approach, based on Schur complements, for computing the embedding of the relevant vertices. We prove that these embeddings are preserved exactly by the Schur complement graph that is obtained via Gaussian elimination on the non-relevant vertices. As computing Schur complements is expensive, we give a nearly-linear time algorithm that generates a coarsened graph on the relevant vertices that provably matches the Schur complement in expectation in each iteration. Our experiments involving prediction tasks on graphs demonstrate that computing embeddings on the coarsened graph, rather than the entire graph, leads to significant time savings without sacrificing accuracy.

We prove a Chernoff-type bound for sums of matrix-valued random variables sampled via a regular (aperiodic and irreducible) finite Markov chain. Specially, consider a random walk on a regular Markov chain and a Hermitian matrix-valued function on its state space. Our result gives exponentially decreasing bounds on the tail distributions of the extreme eigenvalues of the sample mean matrix. Our proof is based on the matrix expander (regular undirected graph) Chernoff bound [Garg et al. STOC '18] and scalar Chernoff-Hoeffding bounds for Markov chains [Chung et al. STACS '12]. Our matrix Chernoff bound for Markov chains can be applied to analyze the behavior of co-occurrence statistics for sequential data, which have been common and important data signals in machine learning. We show that given a regular Markov chain with n states and mixing time t, we need a trajectory of length O(t(log(n) + log(t))/e^2) to achieve an estimator of the co-occurrence matrix with error bound e. We conduct several experiments and the experimental results are consistent with the exponentially fast convergence rate from theoretical analysis. Our result gives the first bound on the convergence rate of the co-occurrence matrix and the first sample complexity analysis in graph representation learning.

In this paper we study the problem of dynamically maintaining graph properties under batches of edge insertions and deletions in the massively parallel model of computation. In this setting, the graph is stored on a number of machines, each having space strongly sublinear with respect to the number of vertices, that is, nϵ for some constant 0 < ϵ < 1. Our goal is to handle batches of updates and queries where the data for each batch fits onto one machine in constant rounds of parallel computation, as well as to reduce the total communication between the machines. This objective corresponds to the gradual buildup of databases over time, while the goal of obtaining constant rounds of communication for problems in the static setting has been elusive for problems as simple as undirected graph connectivity. We give an algorithm for dynamic graph connectivity in this setting with constant communication rounds and communication cost almost linear in terms of the batch size. Our techniques combine a new graph contraction technique, an independent random sample extractor from correlated samples, as well as distributed data structures supporting parallel updates and queries in batches. We also illustrate the power of dynamic algorithms in the MPC model by showing that the batched version of the adaptive connectivity problem is P-complete in the centralized setting, but sub-linear sized batches can be handled in a constant number of rounds. Due to the wide applicability of our approaches, we believe it represents a practically-motivated workaround to the current difficulties in designing more efficient massively parallel static graph algorithms.

The problem of finding dense components of a graph is a major primitive in graph mining and data analysis. The densest subgraph problem (DSP) that asks to find a subgraph with maximum average degree forms a basic primitive in dense subgraph discovery with applications ranging from community detection to unsupervised discovery of biological network modules [16]. The DSP is exactly solvable in polynomial time using maximum flows [14, 17, 22]. Due to the high computational cost of maximum flows, Charikar’s greedy approximation algorithm is usually preferred in practice due to its linear time and linear space complexity [3, 8]. It constitutes a key algorithmic idea in scalable solutions for large-scale dynamic graphs [5, 7]. However, its output density can be a factor 2 off the optimal solution.

We show that variants of spectral sparsification routines can preserve the total spanning tree counts of graphs. By Kirchhoff's matrix-tree theorem, this is equivalent to preserving the determinant of a graph Laplacian minor or, equivalently, of any symmetric diagonally dominant matrix (SDDM). Our analyses utilize this combinatorial connection to bridge the gap between statistical leverage scores/effective resistances and the analysis of random graphs by Janson [ Combin. Probab. Comput., 3 (1994), pp. 97--126]. This leads to a routine that, in quadratic time, sparsifies a graph down to about $n^{1.5}$ edges in a way that preserves both the determinant and the distribution of spanning trees (provided the sparsified graph is viewed as a random object). Extending this algorithm to work with Schur complements and approximate Choleksy factorizations leads to algorithms for counting and sampling spanning trees which are nearly optimal for dense graphs. We give an algorithm that computes a $(1 \pm \delta)$ approximation to the determinant of any SDDM matrix with constant probability in about $n^2 \delta^{-2}$ time. This is the first routine for graphs that outperforms general-purpose routines for computing determinants of arbitrary matrices. We also give an algorithm that generates, in about $n^2 \delta^{-2}$ time, a spanning tree of a weighted undirected graph from a distribution with a total variation distance of $\delta$ from the $\boldsymbol{\mathit{w}}$-uniform distribution.

Linear regression in Lp-norm is a canonical optimization problem that arises in several applications, including sparse recovery, semi-supervised learning, and signal processing. Generic convex optimization algorithms for solving Lp-regression are slow in practice. Iteratively Reweighted Least Squares (IRLS) is an easy to implement family of algorithms for solving these problems that has been studied for over 50 years. However, these algorithms often diverge for p > 3, and since the work of Osborne (1985), it has been an open problem whether there is an IRLS algorithm that converges for p > 3. We propose p-IRLS, the first IRLS algorithm that provably converges geometrically for any p \in [2,\infty). Our algorithm is simple to implement and is guaranteed to find a high accuracy solution in a sub-linear number of iterations. Our experiments demonstrate that it performs even better than our theoretical bounds, beats the standard Matlab/CVX implementation for solving these problems by 10–50x, and is the fastest among available implementations in the high-accuracy regime.

We give improved algorithms for the ℓp-regression problem, minx ‖x‖p such that Ax = b, for all p ∊ (1, 2) ∪ (2, ∞). Our algorithms obtain a high accuracy solution in iterations, where each iteration requires solving an m × m linear system, with m being the dimension of the ambient space. Incorporating a procedure for maintaining an approximate inverse of the linear systems that we need to solve at each iteration, we give algorithms for solving ℓp-regression to 1/poly(n) accuracy that runs in time Õp(mmax{ω, 7/3}), where ω is the matrix multiplication constant. For the current best value of ω > 2.37, this means that we can solve ℓp regression as fast as ℓ2 regression, for all constant p bounded away from 1. Our algorithms can be combined with nearly-linear time solvers for linear systems in graph Laplacians to give minimum ℓp-norm flow / voltage solutions to 1/poly(n) accuracy on an undirected graph with m edges in time. For sparse graphs and for matrices with similar dimensions, our iteration counts and running times improve upon the p-norm regression algorithm by [Bubeck-Cohen-Lee-Li STOC'18], as well as general purpose convex optimization algorithms. At the core of our algorithms is an iterative refinement scheme for ℓp-norms, using the quadratically-smoothed ℓp-norms introduced in the work of Bubeck et al. Formally, given an initial solution, we construct a problem that seeks to minimize a quadratically-smoothed ℓp norm over a subspace, such that a crude solution to this problem allows us to improve the initial solution by a constant factor, leading to algorithms with fast convergence.

We study dynamic algorithms for maintaining spectral vertex sparsifiers of graphs with respect to a set of terminals T of our choice. Such objects preserve pairwise resistances, solutions to systems of linear equations, and energy of electrical flows between the terminals in T. We give a data structure that supports insertions and deletions of edges, and terminal additions, all in sublinear time. We then show the applicability of our result to the following problems.

We present algorithms for solving a large class of flow and regression problems on unit weighted graphs to (1 + 1 / poly(n)) accuracy in almost-linear time. These problems include ℓp-norm minimizing flow for p large (p ∈ [ω(1), o(log2/3 n) ]), and their duals, ℓp-norm semi-supervised learning for p close to 1.

The problem of selecting a group of vertices under certain constraints that maximize their joint centrality arises in many practical scenarios. In this paper, we extend the notion of current flow closeness centrality (CFCC) to a set of vertices in a graph, and investigate the problem of selecting a subset S to maximizes its CFCC C(S), with the cardinality constraint |S| = k. We show the NP-hardness of the problem, but propose two greedy algorithms to minimize the reciprocal of C(S). We prove the approximation ratios by showing the monotonicity and supermodularity. A proposed deterministic greedy algorithm has an approximation factor and cubic running time. To compare with, a proposed randomized algorithm gives -approximation in nearly-linear time, for any ? > 0. Extensive experiments on model and real networks demonstrate the effectiveness and efficiency of the proposed algorithms, with the randomized algorithm being applied to massive networks with more than a million vertices.

We develop a framework for graph sparsification based on a new tool, short cycle decomposition for graphs - a decomposition of a graph into a collection of short cycles, plus a small number of extra edges. A simple observation gives that every graph G on n vertices with m edges can be decomposed in O(mn) time into cycles of length at most 2 log n, and at most 2n extra edges. We give an m1+o(1)time algorithm for constructing a short cycle decomposition of the graph, with cycles of length no(1), and n1+o(1)extra edges. Both the existential and algorithmic variants of this decomposition enable us to make progress on several open problems in randomized graph algorithms. 1. We present an algorithm that runs in time m1+o(1)ε-1.5and returns (1 ± ε)-approximations to effective resistances of all edges, improving over the previous best of Õ(min{mε-2, n2ε-1}) This gives an algorithm to approximate the determinant of a graph Laplacian up to a factor of (1 ± ε) in roughly m + n15/8ε-7/4. 2. We show existence and efficient algorithms for constructing graphical spectral sketches - a distribution over sparse graphs H with about nε-1edges such that for a fixed vector x, we have xTLHx = (1 ± eps) xTLGx and xTL+Hx = (1 ± ε) xTL+Gx with high probability, where L is the graph Laplacian and L+ is its pseudoinverse. This implies resistance-sparsifiers with about nε edges that preserve the effective resistances between every pair of vertices up to (1 + eps). 3. By combining short cycle decomposition with importance sampling, we show the existence of nearly-linear sized degree-preserving spectral sparsifiers, as well as significantly sparser approximations of directed graphs. The latter is critical to recent breakthroughs on faster algorithms for directed random walks and linear systems in directed Laplacian. The running time and output qualities of our spectral sketch and degree-preserving (directed) sparsification algorithms are limited by the efficiency of our routines for producing short cycle decompositions. Improved algorithms for short cycle decompositions will lead to improvements for each of these algorithms.

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.

Motivated by the study of matrix elimination orderings in combinatorial scientific computing, we utilize graph sketching and local sampling to give a data structure that provides access to approximate fill degrees of a matrix undergoing elimination in polylogarithmic time per elimination and query. We then study the problem of using this data structure in the minimum degree algorithm, which is a widely-used heuristic for producing elimination orderings for sparse matrices by repeatedly eliminating the vertex with (approximate) minimum fill degree. This leads to a nearly-linear time algorithm for generating approximate greedy minimum degree orderings. Despite extensive studies of algorithms for elimination orderings in combinatorial scientific computing, our result is the first rigorous incorporation of randomized tools in this setting, as well as the first nearly-linear time algorithm for producing elimination orderings with provable approximation guarantees. While our sketching data structure readily works in the oblivious adversary model, by repeatedly querying and greedily updating itself, it enters the adaptive adversarial model where the underlying sketches become prone to failure due to dependency issues with their internal randomness. We show how to use an additional sampling procedure to circumvent this problem and to create an independent access sequence. Our technique for decorrelating interleaved queries and updates to this randomized data structure may be of independent interest.

We present an asymptotically faster algorithm for solving linear systems in well-structured 3-dimensional truss stiffness matrices. These linear systems arise from linear elasticity problems, and can be viewed as extensions of graph Laplacians into higher dimensions. Faster solvers for the 2-D variants of such systems have been studied using generalizations of tools for solving graph Laplacians [Daitch-Spielman CSC’07, Shklarski-Toledo SIMAX’08].

We give faster algorithms for producing sparse approximations of the transition matrices of k-step random walks on undirected and weighted graphs. These transition matrices also form graphs, and arise as intermediate objects in a variety of graph algorithms. Our improvements are based on a better understanding of processes that sample such walks, as well as tighter bounds on key weights underlying these sampling processes. On a graph with n vertices and m edges, our algorithm produces a graph with about nlog(n) edges that approximates the k-step random walk graph in about m + k^2 nlog^4(n) time. In order to obtain this runtime bound, we also revisit "density independent" algorithms for sparsifying graphs whose runtime overhead is expressed only in terms of the number of vertices.

We show variants of spectral sparsification routines can preserve the total spanning tree counts of graphs, which by Kirchhoff's matrix-tree theorem, is equivalent to determinant of a graph Laplacian minor, or equivalently, of any SDDM matrix. Our analyses utilizes this combinatorial connection to bridge between statistical leverage scores/effective resistances and the analysis of random graphs by [Janson, Combinatorics, Probability and Computing `94]. This leads to a routine that in quadratic time, sparsifies a graph down to about n1.5 edges in ways that preserve both the determinant and the distribution of spanning trees (provided the sparsified graph is viewed as a random object). Extending this algorithm to work with Schur complements and approximate Cholesky factorizations leads to algorithms for counting and sampling spanning trees which are nearly optimal for dense graphs. We give an algorithm that computes a (1±δ) approximation to the determinant of any SDDM matrix with constant probability in about n2δ-2time. This is the first routine for graphs that outperforms general-purpose routines for computing determinants of arbitrary matrices. We also give an algorithm that generates in about n2δ-2time a spanning tree of a weighted undirected graph from a distribution with total variation distance of δ from the w-uniform distribution.

We develop a framework for graph sparsification and sketching, based on a new tool, short cycle decomposition, which is a decomposition of an unweighted graph into an edge-disjoint collection of short cycles, plus a small number of extra edges. A simple observation shows that every graph $G$ on $n$ vertices with $m$ edges can be decomposed in $O(mn)$ time into cycles of length at most $2 \log n,$ and at most $2n$ extra edges. We give an $(m^{1+o(1)})$-time algorithm for constructing a short cycle decomposition, with cycles of length $n^{o(1)},$ and $n^{1+o(1)}$ extra edges. Both the existential and algorithmic variants of this decomposition enable us to make the following progress on several open problems in randomized graph algorithms: (1) We present an algorithm that runs in time $m^{1+o(1)}\varepsilon^{-1.5}$ and returns $(1\pm\varepsilon)$-approximations to effective resistances of all edges, improving over the previous best runtime of $\widetilde{{O}}(\min\{m\varepsilon^{-2}, n^{2} \varepsilon^{-1}\})$. This routine in turn gives an algorithm for approximating the determinant of a graph Laplacian up to a factor of $(1\pm \varepsilon)$ in $m^{1 + o(1)} + n^{\nicefrac{15}{8}+o(1)}\varepsilon^{-\nicefrac{7}{4}}$ time. (2) We show the existence of graphical spectral sketches with about $n\varepsilon^{-1}$ edges, and also give efficient algorithms to construct them. A graphical spectral sketch is a distribution over sparse graphs $H$ such that for a fixed vector ${\mathit{x}}$, we have ${{x}}^{\top} {L}_H {{x}} = (1\pm\varepsilon) {{x}}^{\top} {L}_G {{x}}$ and ${{x}}^{\top} {L}^{+}_H {{x}} = (1\pm\varepsilon) {{x}}^{\top} {L}^{+}_G {{x}}$ with high probability, where ${L}$ is the graph Laplacian and ${L}^{+}$ is its pseudoinverse. This implies the existence of resistance sparsifiers with about $n \varepsilon^{-1}$ edges that preserve the effective resistance between every pair of vertices up to $(1\pm\varepsilon)$. (3) By combining short cycle decompositions with known tools in graph sparsification, we show the existence of nearly linear sized degree-preserving spectral sparsifiers, as well as significantly sparser approximations of Eulerian directed graphs. The latter is critical to recent breakthroughs on faster algorithms for solving linear systems in directed Laplacians. The running time and output qualities of our spectral sketch and degree-preserving (directed) sparsification algorithms are limited by the efficiency of our routines for constructing short cycle decompositions. Improved algorithms for short cycle decompositions will lead to improvement in each of these algorithms.

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.

In this paper we study variants of the widely used spectral clustering that partitions a graph into $k$ clusters by (1) embedding the vertices of a graph into a low-dimensional space using the bottom eigenvectors of the Laplacian matrix and (2) grouping the embedded points into $k$ clusters via $k$-means algorithms. We show that, for a wide class of graphs, spectral clustering gives a good approximation of the optimal clustering. While this approach was proposed in the early 1990s and has comprehensive applications, prior to our work similar results were known only for graphs generated from stochastic models. We also give a nearly linear time algorithm for partitioning well-clustered graphs based on computing a matrix exponential and approximate nearest neighbor data structures.

We present a simple spectral approach to the well-studied constrained clustering problem. It captures constrained clustering as a generalized eigenvalue problem with graph Laplacians. The algorithm works in nearly-linear time and provides concrete guarantees for the quality of the clusters, at least for the case of 2-way partitioning. In practice this translates to a very fast implementation that consistently outperforms existing spectral approaches both in speed and quality.

We initiate the study of fast dynamic algorithms for graph sparsification problems and obtain fully dynamic algorithms, allowing both edge insertions and edge deletions, that take polylogarithmic time after each update in the graph. Our three main results are as follows. First, we give a fully dynamic algorithm for maintaining a (1 ± ϵ)-spectral sparsifier with amortized update time poly(log n, ϵ-1). Second, we give a fully dynamic algorithm for maintaining a (1 ± ϵ)-cut sparsifier with worst-case update time poly(log n, ϵ-1). Both sparsifiers have size n · poly(log n, ϵ-1). Third, we apply our dynamic sparsifier algorithm to obtain a fully dynamic algorithm for maintaining a (1 - ϵ)-approximation to the value of the maximum flow in an unweighted, undirected, bipartite graph with amortized update time poly(log n, ϵ-1).

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.

Tensor CANDECOMP/PARAFAC (CP) decomposition is a powerful but computationally challenging tool in modern data analytics. In this paper, we show ways of sampling intermediate steps of alternating minimization algorithms for computing low rank tensor CP decompositions, leading to the sparse alternating least squares (SPALS) method. Specifically, we sample the the Khatri-Rao product, which arises as an intermediate object during the iterations of alternating least squares. This product captures the interactions between different tensor modes, and form the main computational bottleneck for solving many tensor related tasks. By exploiting the spectral structures of the matrix Khatri-Rao product, we provide efficient access to its statistical leverage scores. When applied to the tensor CP decomposition, our method leads to the first algorithm that runs in sublinear time per-iteration and approximates the output of deterministic alternating least squares algorithms. Empirical evaluations of this approach show significantly speedups over existing randomized and deterministic routines for performing CP decomposition. On a tensor of the size 2.4m by 6.6m by 92k with over 2 billion nonzeros formed by Amazon product reviews, our routine converges in two minutes to the same error as deterministic ALS.

We give the first O(mpolylog(n)) time algorithms for approximating maximum flows in undirected graphs and constructing polylog(n)-quality cut-approximating hierarchical tree decompositions. Our algorithm invokes existing algorithms for these two problems recursively while gradually incorporating size reductions. These size reductions are in turn obtained via ultra-sparsifiers, which are key tools in solvers for symmetric diagonally dominant (SDD) linear systems.

We introduce the sparsified Cholesky and sparsified multigrid algorithms for solving systems of linear equations. These algorithms accelerate Gaussian elimination by sparsifying the nonzero matrix entries created by the elimination process. We use these new algorithms to derive the first nearly linear time algorithms for solving systems of equations in connection Laplacians---a generalization of Laplacian matrices that arise in many problems in image and signal processing. We also prove that every connection Laplacian has a linear sized approximate inverse. This is an LU factorization with a linear number of nonzero entries that is a strong approximation of the original matrix. Using such a factorization one can solve systems of equations in a connection Laplacian in linear time. Such a factorization was unknown even for ordinary graph Laplacians.

We study algorithms for spectral graph sparsification. The input is a graph G with n vertices and m edges, and the output is a sparse graph G˜ that approximates G in an algebraic sense. Concretely, for all vectors x and any ϵ > 0, the graph G˜ satisfies (1-ϵ ) x T L G x ≤ x T L G˜ x ≤ (1+ϵ) x T L G x , where L G and G˜ are the Laplacians of G and G˜ respectively. The first contribution of this article applies to all existing sparsification algorithms that rely on solving solving linear systems on graph Laplacians. These algorithms are the fastest known to date. Specifically, we show that less precision is required in the solution of the linear systems, leading to speedups by an O (log n ) factor. We also present faster sparsification algorithms for slightly dense graphs: — An O ( m log n ) time algorithm that generates a sparsifier with O ( n log 3 n /ϵ 2 ) edges. — An O ( m log log n ) time algorithm for graphs with more than n log 5 n log log n edges. — An O ( m ) algorithm for graphs with more than n log 10 n edges. — An O ( m ) algorithm for unweighted graphs with more than n log 8 n edges. These bounds hold up to factors that are in O ( poly (log log n )) and are conjectured to be removable.

Motivated by a sampling problem basic to computational statistical inference, we develop a toolset based on spectral sparsification for a family of fundamental problems involving Gaussian sampling, matrix functionals, and reversible Markov chains. Drawing on the connection between Gaussian graphical models and the recent breakthroughs in spectral graph theory, we give the first nearly linear time algorithm for the following basic matrix problem: Given an n\times n Laplacian matrix \mathbfM and a constant -1 ≤p ≤1, provide efficient access to a sparse n\times n linear operator \tilde\mathbfC such that $\mathbfM^p ≈\tilde\mathbfC \tilde\mathbfC^⊤, where ≈denotes spectral similarity. When p is set to -1, this gives the first parallel sampling algorithm that is essentially optimal both in total work and randomness for Gaussian random fields with symmetric diagonally dominant (SDD) precision matrices. It only requires \em nearly linear work and 2n \em i.i.d. random univariate Gaussian samples to generate an n-dimensional \em i.i.d. Gaussian random sample in polylogarithmic depth. The key ingredient of our approach is an integration of spectral sparsification with multilevel method: Our algorithms are based on factoring \mathbfM^p$ into a product of well-conditioned matrices, then introducing powers and replacing dense matrices with sparse approximations. We give two sparsification methods for this approach that may be of independent interest. The first invokes Maclaurin series on the factors, while the second builds on our new nearly linear time spectral sparsification algorithm for random-walk matrix polynomials. We expect these algorithmic advances will also help to strengthen the connection between machine learning and spectral graph theory, two of the most active fields in understanding large data and networks.

In this work we study the widely used \emphspectral clustering algorithms, i.e. partition a graph into k clusters via (1) embedding the vertices of a graph into a low-dimensional space using the bottom eigenvectors of the Laplacian matrix, and (2) partitioning embedded points via k-means algorithms. We show that, for a wide class of \emphwell-clustered graphs, spectral clustering algorithms can give a good approximation of the optimal clustering. To the best of our knowledge, it is the \emphfirst theoretical analysis of spectral clustering algorithms for a wide family of graphs, even though such approach was proposed in the early 1990s and has comprehensive applications. We also give a nearly-linear time algorithm for partitioning well-clustered graphs, which is based on heat kernel embeddings and approximate nearest neighbor data structures.

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.

Extracting dense subgraphs from large graphs is a key primitive in a variety of graph mining applications, ranging from mining social networks and the Web graph to bioinformatics [41]. In this paper we focus on a family of poly-time solvable formulations, known as the k-clique densest subgraph problem (k-Clique-DSP) [57]. When k=2, the problem becomes the well-known densest subgraph problem (DSP) [22, 31, 33, 39]. Our main contribution is a sampling scheme that gives densest subgraph sparsifier, yielding a randomized algorithm that produces high-quality approximations while providing significant speedups and improved space complexity. We also extend this family of formulations to bipartite graphs by introducing the (p,q)-biclique densest subgraph problem ((p,q)-Biclique-DSP), and devise an exact algorithm that can treat both clique and biclique densities in a unified way.

We use exponential start time clustering to design faster parallel graph algorithms involving distances. Previous algorithms usually rely on graph decomposition routines with strict restrictions on the diameters of the decomposed pieces. We weaken these bounds in favor of stronger local probabilistic guarantees. This allows more direct analyses of the overall process, giving:

We give a simple algorithm to efficiently sample the rows of a matrix while preserving the p-norms of its product with vectors. Given an n * d matrix A, we find with high probability and in input sparsity time an A' consisting of about d log d rescaled rows of A such that |Ax|1 is close to |A'x|1 for all vectors x. We also show similar results for all Lp that give nearly optimal sample bounds in input sparsity time. Our results are based on sampling by "Lewis weights", which can be viewed as statistical leverage scores of a reweighted matrix. We also give an elementary proof of the guarantees of this sampling process for L1.

We present an efficient algorithm for solving a linear system arising from the 1-Laplacian corresponding to a collapsible simplicial complex with a known collapsing sequence. When combined with a result of Chillingworth, our algorithm is applicable to convex simplicial complexes embedded in ℝ3. The running time of our algorithm is nearly-linear in the size of the complex and is logarithmic on its numerical properties. Our algorithm is based on projection operators and combinatorial steps for transferring between them. The former relies on decomposing flows into circulations and potential flows using fast solvers for graph Laplacians, and the latter relates Gaussian elimination to topological properties of simplicial complexes.

We show an algorithm for solving symmetric diagonally dominant (SDD) linear systems with m non-zero entries to a relative error of ε in O(m log1/2 n logc n log(1/ε)) time. Our approach follows the recursive preconditioning framework, which aims to reduce graphs to trees using iterative methods. We improve two key components of this framework: random sampling and tree embeddings. Both of these components are used in a variety of other algorithms, and our approach also extends to the dual problem of computing electrical flows.

We present the first parallel algorithm for solving systems of linear equations in symmetric, diagonally dominant (SDD) matrices that runs in polylogarithmic time and nearly-linear work. The heart of our algorithm is a construction of a sparse approximate inverse chain for the input matrix: a sequence of sparse matrices whose product approximates its inverse. Whereas other fast algorithms for solving systems of equations in SDD matrices exploit low-stretch spanning trees, our algorithm only requires spectral graph sparsifiers.

We present an algorithm that on input of an $n$-vertex $m$-edge weighted graph $G$ and a value $k$ produces an incremental sparsifier $\hat{G}$ with $n-1 + m/k$ edges, such that the relative condition number of $G$ with $\hat{G}$ is bounded above by $\tilde{O}(k\log^2 n)$, with probability $1-p$ (we use the $\tilde{O}()$ notation to hide a factor of at most $(\log\log n)^4$). The algorithm runs in time $\tilde{O}((m \log{n} + n\log^2{n})\log(1/p)).$ As a result, we obtain an algorithm that on input of an $n\times n$ symmetric diagonally dominant matrix $A$ with $m$ nonzero entries and a vector $b$ computes a vector ${x}$ satisfying $||{x}-A^{+}b||_A<\epsilon ||A^{+}b||_A $, in expected time $\tilde{O}(m\log^2{n}\log(1/\epsilon)).$ The solver is based on repeated applications of the incremental sparsifier that produces a chain of graphs which is then used as input to the recursive preconditioned Chebyshev iteration.

We present the first data structures that maintain near optimal maximum cardinality and maximum weighted matchings on sparse graphs in sub linear time per update. Our main result is a data structure that maintains a (1+ε) approximation of maximum matching under edge insertions/deletions in worst case Õ(?mε-2) time per update. This improves the 3/2 approximation given by Neiman and Solomon [20] which runs in similar time. The result is based on two ideas. The first is to re-run a static algorithm after a chosen number of updates to ensure approximation guarantees. The second is to judiciously trim the graph to a smaller equivalent one whenever possible. We also study extensions of our approach to the weighted setting, and combine it with known frameworks to obtain arbitrary approximation ratios. For a constant ε and for graphs with edge weights between 1 and N, we design an algorithm that maintains an (1+ε) approximate maximum weighted matching in Õ(?m log N) time per update. The only previous result for maintaining weighted matchings on dynamic graphs has an approximation ratio of 4.9108, and was shown by An and et al. [2], [3].

There has been significant interest and progress recently in algorithms that solve regression problems involving tall and thin matrices in input sparsity time. Given a n * d matrix where n ≥ d, these algorithms find an approximation with fewer rows, allowing one to solve a poly(d) sized problem instead. In practice, the best performances are often obtained by invoking these routines in an iterative fashion. We show these iterative methods can be adapted to give theoretical guarantees comparable to and better than the current state of the art. Our approaches are based on computing the importances of the rows, known as leverage scores, in an iterative manner. We show that alternating between computing a short matrix estimate and finding more accurate approximate leverage scores leads to a series of geometrically smaller instances. This gives an algorithm whose runtime is input sparsity plus an overhead comparable to the cost of solving a regression problem on the smaller approximation. Our results build upon the close connection between randomized matrix algorithms, iterative methods, and graph sparsification.

We study theoretical runtime guarantees for a class of optimization problems that occur in a wide variety of inference problems. These problems are motivated by the LASSO framework and have applications in machine learning and computer vision. Our work shows a close connection between these problems and core questions in algorithmic graph theory. While this connection demonstrates the difficulties of obtaining runtime guarantees, it also suggests an approach of using techniques originally developed for graph algorithms.

We present faster algorithms for approximate maximum flow in undirected graphs with good separator structures, such as bounded genus, minor free, and geometric graphs. Given such a graph with n vertices, m edges along with a recursive -vertex separator structure, our algorithm finds an 1 − ∊ approximate maximum flow in time Õ(m6/5poly(∊−1)), ignoring poly-logarithmic terms. Similar speedups are also achieved for separable graphs with larger size separators albeit with larger run times. These bounds also apply to image problems in two and three dimensions. Key to our algorithm is an intermediate problem that we term grouped L2 flow, which exists between maximum flows and electrical flows. Our algorithm also makes use of spectral vertex sparsifiers in order to remove vertices while preserving the energy dissipation of electrical flows. We also give faster spectral vertex sparsification algorithms on well separated graphs, which may be of independent interest.

We show an improved parallel algorithm for decomposing an undirected unweighted graph into small diameter pieces with a small fraction of the edges in between. These decompositions form critical subroutines in a number of graph algorithms. Our algorithm builds upon the shifted shortest path approach introduced in [Blelloch, Gupta, Koutis, Miller, Peng, Tangwongsan, SPAA 2011]. By combining various stages of the previous algorithm, we obtain a significantly simpler algorithm with the same asymptotic guarantees as the best sequential algorithm.

This paper studies the problem of finding a (1+ε)-approximate solution to positive semidefinite programs. These are semidefinite programs in which all matrices in the constraints and objective are positive semidefinite and all scalars are non-negative. At FOCS'11, Jain and Yao gave an NC algorithm that requires O(t 1/ε13 log13 m log n) iterations on input n constraint matrices of dimension m-by-m, where each iteration performs at least Ω(mω) work since it involves computing the spectral decomposition. We present a simpler NC parallel algorithm that on input with n constraint matrices, requires O(1/ε4 log4 n log(1/ε)) iterations, each of which involves only simple matrix operations and computing the trace of the product of a matrix exponential and a positive semidefinite matrix. Further, given a positive SDP in a factorized form, the total work of our algorithm is nearly-linear in the number of non-zero entries in the factorization. Our algorithm can be viewed as a generalization of Young's algorithm and analysis techniques for positive linear programs (Young, FOCS'01) to the semidefinite programming setting.

We present three spectral sparsification algorithms that, on input a graph G with n vertices and m edges, return a graph H with n vertices and O(n log n/epsilon^2) edges that provides a strong approximation of G. Namely, for all vectors x and any epsilon>0, we have (1-epsilon) x^T L_G x <= x^T L_H x <= (1+epsilon) x^T L_G x, where L_G and L_H are the Laplacians of the two graphs. The first algorithm is a simple modification of the fastest known algorithm and runs in tilde{O}(m log^2 n) time, an O(log n) factor faster than before. The second algorithm runs in tilde{O}(m log n) time and generates a sparsifier with tilde{O}(n log^3 n) edges. The third algorithm applies to graphs where m>n log^5 n and runs in tilde{O}(m log_{m/ n log^5 n} n time. In the range where m>n^{1+r} for some constant r this becomes softO(m). The improved sparsification algorithms are employed to accelerate linear system solvers and algorithms for computing fundamental eigenvectors of dense SDD matrices.

The maximum multicommodity flow problem is a natural generalization of the maximum flow problem to route multiple distinct flows. Obtaining a 1-ε approximation to the multicommodity flow problem on graphs is a well-studied problem. In this paper we present an adaptation of recent advances in single-commodity flow algorithms to this problem. As the underlying linear systems in the electrical problems of multicommodity flow problems are no longer Laplacians, our approach is tailored to generate specialized systems which can be preconditioned and solved efficiently using Laplacians. Given an undirected graph with m edges and k commodities, we give algorithms that find 1-ε approximate solutions to the maximum concurrent flow problem and maximum weighted multicommodity flow problem in time O(m4/3poly(k,ε-1)).

We present an improved algorithm for solving symmetrically diagonally dominant linear systems. On input of an n×n symmetric diagonally dominant matrix A with m non-zero entries and a vector b such that Ax̅ = b for some (unknown) vector x̅, our algorithm computes a vector x such that ∥x-x̅∥A≤ϵ∥x̅∥A1in time Õ (m log n log (1/ϵ))2. The solver utilizes in a standard way a 'preconditioning' chain of progressively sparser graphs. To claim the faster running time we make a two-fold improvement in the algorithm for constructing the chain. The new chain exploits previously unknown properties of the graph sparsification algorithm given in [Koutis,Miller,Peng, FOCS 2010], allowing for stronger preconditioning properties.We also present an algorithm of independent interest that constructs nearly-tight low-stretch spanning trees in time Õ (m log n), a factor of O (log n) faster than the algorithm in [Abraham,Bartal,Neiman, FOCS 2008]. This speedup directly reflects on the construction time of the preconditioning chain.

We consider the problem of approximating a signal P with another signal F consisting of a few piecewise constant segments. This problem arises naturally in applications including databases (e.g., histogram construction), speech recognition, computational biology (e.g., denoising aCGH data) and many more. Specifically, let P = (P1, P2, …, Pn), Pi ∊ ℝ for all i, be a signal and let C be a constant. Our goal is to find a function F : [n] → ℝ which optimizes the following objective function: The above optimization problem reduces to solving the following recurrence, which can be done using dynamic programming in O(n2) time: This recurrence arises naturally in several applications where one wants to approximate a given signal P with a signal F which ideally consists of few piecewise constant segments. Such applications include histogram construction in databases, determining DNA copy numbers in cancer cells from micro-array data, speech recognition, data mining and many others. In this work we present two new techniques for optimizing dynamic programming that can handle cost functions not treated by other standard methods. The basis of our first algorithm is the definition of a constant-shifted variant of the objective function that can be efficiently approximated using state of the art methods for range searching. Our technique approximates the optimal value of our objective function within additive ∊ error and runs in time, where δ is an arbitrarily small positive constant and . The second algorithm we provide solves a similar recurrence that's within a multiplicative factor of (1+∊) and runs in O(n log n/∊). The new technique introduced by our algorithm is the decomposition of the initial problem into a small (logarithmic) number of Monge optimization subproblems which we can speed up using existing techniques.

This paper presents the design and analysis of a near linear-work parallel algorithm for solving symmetric diagonally dominant (SDD) linear systems. On input an SDD n-by-n matrix A with m non-zero entries and a vector b, our algorithm computes a vector x such that Ax - A+b ≤ ε • A+b in O(m logO(1) n log 1/ε) work and O(m1/3+θ log 1/ε) depth for any fixed θ > 0.

We present parallel greedy approximation algorithms for set cover and related problems. These algorithms build on an algorithm for solving a graph problem we formulate and study called Maximal Nearly Independent Set (MaNIS)---a graph abstraction of a key component in existing work on parallel set cover.

We present an algorithm that on input of an n-vertex m-edge weighted graph G and a value k, produces an incremental sparsifier G with n-1+m/k edges, such that the condition number of G with G is bounded above by Õ(k log2n), with probability 1-p. The algorithm runs in time Õ((m log n + n log n) log(1/p)). As a result, we obtain an algorithm that on input of an n × n symmetric diagonally dominant matrix A with m non-zero entries and a vector b, computes a vector x satisfying ||x-A+b||A+b||A, in expected time Õ(m log2n log(1/ϵ)). The solver is based on repeated applications of the incremental sparsifier that produces a chain of graphs which is then used as input to a recursive preconditioned Chebyshev iteration.