SODA Conference 2025 Conference Paper
Eulerian Graph Sparsification by Effective Resistance Decomposition
- Arun Jambulapati
- Sushant Sachdeva
- Aaron Sidford
- Kevin Tian
- Yibin Zhao 0003
Author name cluster
Possible papers associated with this exact author name in Arrow. This page groups case-insensitive exact name matches and is not a full identity disambiguation profile.
SODA Conference 2025 Conference Paper
FOCS Conference 2024 Conference Paper
We give the first almost-linear total time algorithm for deciding if a flow of cost at most $F$ still exists in a directed graph, with edge costs and capacities, undergoing decremental updates, i. e. , edge deletions, capacity decreases, and cost increases. This implies almost-linear time algorithms for approximating the minimum-cost flow value and s-t distance on such decremental graphs. Our framework additionally allows us to maintain decremental strongly connected components in almost-linear time deterministically. These algorithms also improve over the current best known runtimes for statically computing minimum-cost flow, in both the randomized and deterministic settings. We obtain our algorithms by taking the dual perspective, which yields cut-based algorithms. More precisely, our algorithm computes the flow via a sequence of $m^{1+o(1)}$ -dynamic min-ratio cut problems, the dual analog of the dynamic min-ratio cycle problem that underlies recent fast algorithms for minimum-cost flow. Our main technical contribution is a new data structure that returns an approximately optimal min-ratio cut in amortized $m^{o(1)}$ time by maintaining a tree-cut sparsifier. This is achieved by devising a new algorithm to maintain the dynamic expander hierarchy of [ $\text{Goranci-Racke-}$ SaranurakTan, SODA 2021] that also works in capacitated graphs. All our algorithms are deterministc, though they can be sped up further using randomized techniques while still working against an adaptive adversary.
SODA Conference 2024 Conference Paper
SODA Conference 2024 Conference Paper
We provide an algorithm which, with high probability, maintains a (1 — ɛ)-approximate maximum flow on an undirected graph undergoing m -edge additions in amortized m o (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. Our algorithm builds upon previous dynamic algorithms for approximately solving the minimum-ratio cycle problem that underlie previous advances on the maximum flow problem [Chen-Kyng-Liu-Peng-Probst Gutenberg-Sachdeva, FOCS ‘22] as well as recent dynamic maximum flow algorithms [v. d. Brand-Liu-Sidford, STOC ‘23]. Instead of using interior point methods, which were a key component of these recent advances, our algorithm uses an optimization method based on ℓ p -norm iterative refinement and the multiplicative weight update method. This ensures a monotonicity property in the minimum-ratio cycle subproblems that allows us to apply known data structures and bypass issues arising from adaptive queries. * The full version of the paper can be accessed at https: //arxiv. org/abs/2311. 03174
FOCS Conference 2023 Conference Paper
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.
SODA Conference 2023 Conference Paper
We demonstrate that for expander graphs, for all ε > 0, there exists a data structure of size Õ ( n ε -1 ) which can be used to return (1 + ε)-approximations to effective resistances in Õ (1) time per query. Short of storing all effective resistances, previous best approaches could achieve Õ(nε -2 ) size and Õ (ε -2 ) time per query by storing Johnson-Lindenstrauss vectors for each vertex, or Õ ( n ε -1 ) size and Õ ( n ε -1 ) time per query by storing a spectral sketch. Our construction is based on two key ideas: 1) ε -1 -sparse, ε-additive approximations to σ u for all u, vectors similar to DL + 1 u, can be used to recover (1 + ε)-approximations to the effective resistances, 2) In expander graphs, only Õ (ε -1 ) coordinates of σ u are larger than ε. We give an efficient construction for such a data structure in Õ ( m + n ε -2 ) time via random walks. This results in an algorithm on expander graphs for computing (1 + ε)-approximate effective resistances for s vertex pairs that runs in Õ ( m + n ε -2 + s ) time, improving over the previously best known running time of m 1+ o (1) + ( n + s)n o (1) ε -1. 5 for s = ω( n ε -0. 5 ). We employ the above algorithm to compute a (1 + δ)-approximation to the number of spanning trees in an expander graph, or equivalently, approximating the (pseudo)determinant of its Laplacian in Õ(m + n 1. 5 δ -1 ) time. This improves on the previously best known result of m 1+ o (1) + n 1. 875 + o (1) δ -1. 75 time, and matches the best known size of determinant sparsifiers. * This research is supported by an NSERC Discovery grant awarded to Sushant Sachdeva.
ICML Conference 2022 Conference Paper
We study a variant of a recently introduced min-max optimization framework where the max-player is constrained to update its parameters in a greedy manner until it reaches a first-order stationary point. Our equilibrium definition for this framework depends on a proposal distribution which the min-player uses to choose directions in which to update its parameters. We show that, given a smooth and bounded nonconvex-nonconcave objective function, access to any proposal distribution for the min-player’s updates, and stochastic gradient oracle for the max-player, our algorithm converges to the aforementioned approximate local equilibrium in a number of iterations that does not depend on the dimension. The equilibrium point found by our algorithm depends on the proposal distribution, and when applying our algorithm to train GANs we choose the proposal distribution to be a distribution of stochastic gradients. We empirically evaluate our algorithm on challenging nonconvex-nonconcave test-functions and loss functions arising in GAN training. Our algorithm converges on these test functions and, when used to train GANs, trains stably on synthetic and real-world datasets and avoids mode collapse.
FOCS Conference 2022 Conference Paper
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.
SODA Conference 2022 Conference Paper
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 ( 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 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.
NeurIPS Conference 2021 Conference Paper
We provide several algorithms for constrained optimization of a large class of convex problems, including softmax, $\ell_p$ regression, and logistic regression. Central to our approach is the notion of width reduction, a technique which has proven immensely useful in the context of maximum flow [Christiano et al. , STOC'11] and, more recently, $\ell_p$ regression [Adil et al. , SODA'19], in terms of improving the iteration complexity from $O(m^{1/2})$ to $\tilde{O}(m^{1/3})$, where $m$ is the number of rows of the design matrix, and where each iteration amounts to a linear system solve. However, a considerable drawback is that these methods require both problem-specific potentials and individually tailored analyses. As our main contribution, we initiate a new direction of study by presenting the first \emph{unified} approach to achieving $m^{1/3}$-type rates. Notably, our method goes beyond these previously considered problems to more broadly capture \emph{quasi-self-concordant} losses, a class which has recently generated much interest and includes the well-studied problem of logistic regression, among others. In order to do so, we develop a unified width reduction method for carefully handling these losses based on a more general set of potentials. Additionally, we directly achieve $m^{1/3}$-type rates in the constrained setting without the need for any explicit acceleration schemes, thus naturally complementing recent work based on a ball-oracle approach [Carmon et al. , NeurIPS'20].
ICML Conference 2020 Conference Paper
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.
SODA Conference 2020 Conference Paper
We present faster high-accuracy algorithms for computing ℓ p -norm minimizing flows. On a graph with m edges, our algorithm can compute a (1 + 1/poly( m ))-approximate unweighted ℓ p -norm minimizing flow with operations, for any p ≥ 2, giving the best bound for all p ≳ 5. 24. Combined with the algorithm from the work of Adil et al. (SODA '19), we can now compute such flows for any 2 ≤ p ≤ m o (1) in time at most O ( m 1. 24 ). In comparison, the previous best running time was Ω( m 1. 33 ) for large constant p. For p ∼ σ −1 log m, our algorithm computes a (1 + σ )-approximate maximum flow on undirected graphs using m 1+ o (1) σ −1 operations, matching the current best bound, albeit only for unit-capacity graphs. We also give an algorithm for solving general ℓ p -norm regression problems for large p. Our algorithm makes calls to a linear solver. This gives the first high-accuracy algorithm for computing weighted ℓ p -norm minimizing flows that runs in time o ( m 1. 5 ) for some p = m Ω(1). Our key technical contribution is to show that smoothed ℓ p -norm problems introduced by Adil et al. , are interreducible for different values of p. No such reduction is known for standard ℓ p -norm problems.
NeurIPS Conference 2020 Conference Paper
Our understanding of learning input-output relationships with neural nets has improved rapidly in recent years, but little is known about the convergence of the underlying representations, even in the simple case of linear autoencoders (LAEs). We show that when trained with proper regularization, LAEs can directly learn the optimal representation -- ordered, axis-aligned principal components. We analyze two such regularization schemes: non-uniform L2 regularization and a deterministic variant of nested dropout [Rippel et al, ICML' 2014]. Though both regularization schemes converge to the optimal representation, we show that this convergence is slow due to ill-conditioning that worsens with increasing latent dimension. We show that the inefficiency of learning the optimal representation is not inevitable -- we present a simple modification to the gradient descent update that greatly speeds up convergence empirically.
NeurIPS Conference 2019 Conference Paper
Linear regression in L p-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 L p-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.
STOC Conference 2019 Conference Paper
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 (log 2/3 n ) ]), and their duals, ℓ p -norm semi-supervised learning for p close to 1. As p tends to infinity, p -norm flow and its dual tend to max-flow and min-cut respectively. Using this connection and our algorithms, we give an alternate approach for approximating undirected max-flow, and the first almost-linear time approximations of discretizations of total variation minimization objectives. Our framework is inspired by the routing-based solver for Laplacian linear systems by Spielman and Teng (STOC ’04, SIMAX ’14), and is based on several new tools we develop, including adaptive non-linear preconditioning, tree-routings, and (ultra-)sparsification for mixed ℓ 2 and ℓ p norm objectives.
SODA Conference 2019 Conference Paper
SODA Conference 2019 Conference Paper
We present improved algorithms for short cycle decomposition of a graph – a decomposition of an undirected, unweighted graph into edge-disjoint cycles, plus a small number of additional edges. Short cycle decompositions were introduced in the recent work of Chu et al. (FOCS 2018), and were used to make progress on several questions in graph sparsification. For all constants δ ∊ (0, 1], we give an O ( mn δ ) time algorithm that, given a graph G, partitions its edges into cycles of length, with O ( n ) extra edges not in any cycle. This gives the first subquadratic, in fact almost linear time, algorithm achieving polylogarithmic cycle lengths. We also give an m · time algorithm that partitions the edges of a graph into cycles of length, with O ( n ) extra edges not in any cycle. This improves on the short cycle decomposition algorithms given by Chu et al. in terms of all parameters, and is significantly simpler. As a result, we obtain faster algorithms and improved guarantees for several problems in graph sparsification – construction of resistance sparsifiers, graphical spectral sketches, degree preserving sparsifiers, and approximating the effective resistances of all edges.
NeurIPS Conference 2019 Conference Paper
Increasing the batch size is a popular way to speed up neural network training, but beyond some critical batch size, larger batch sizes yield diminishing returns. In this work, we study how the critical batch size changes based on properties of the optimization algorithm, including acceleration and preconditioning, through two different lenses: large scale experiments and analysis using a simple noisy quadratic model (NQM). We experimentally demonstrate that optimization algorithms that employ preconditioning, specifically Adam and K-FAC, result in much larger critical batch sizes than stochastic gradient descent with momentum. We also demonstrate that the NQM captures many of the essential features of real neural network training, despite being drastically simpler to work with. The NQM predicts our results with preconditioned optimizers, previous results with accelerated gradient descent, and other results around optimal learning rates and large batch training, making it a useful tool to generate testable predictions about neural network optimization. We demonstrate empirically that the simple noisy quadratic model (NQM) displays many similarities to neural networks in terms of large-batch training. We prove analytical convergence results for the NQM model that predict such behavior and hence provide possible explanations and a better understanding for many large-batch training phenomena.
FOCS Conference 2018 Conference Paper
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 m 1+o(1) time algorithm for constructing a short cycle decomposition of the graph, 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 progress on several open problems in randomized graph algorithms. 1. We present an algorithm that runs in time m 1+o(1) ε -1. 5 and returns (1 ± ε)-approximations to effective resistances of all edges, improving over the previous best of Õ(min{mε -2, n 2 ε -1 }) This gives an algorithm to approximate the determinant of a graph Laplacian up to a factor of (1 ± ε) in roughly m + n 15/8 ε -7/4. 2. We show existence and efficient algorithms for constructing graphical spectral sketches - a distribution over sparse graphs H with about nε -1 edges such that for a fixed vector x, we have x T L H x = (1 ± eps) x T L G x and x T L+ H x = (1 ± ε) x T L+ G x 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.
SODA Conference 2018 Conference Paper
SODA Conference 2017 Conference Paper
A spectral sparsifier of a graph G is a sparser graph H that approximately preserves the quadratic form of G, i. e. , for all vectors x, x T L Gx ≈ x T L H x, where L G and L H denote the respective graph Laplacians. Spectral sparsifiers generalize cut sparsifiers, and have found many applications in designing graph algorithms. In recent years, there has been interest in computing spectral sparsifiers in semi-streaming and dynamic settings. Natural algorithms in these settings often involve repeated sparsification of a graph, and in turn accumulation of errors across these steps. We present a framework for analyzing algorithms that perform repeated sparsifications that only incur error corresponding to a single sparsification step, leading to better results for many of these reseparsification based algorithms. As an application, we show how to maintain a spectral sparsifier in the semi-streaming setting: We present a simple algorithm that, for a graph G on n vertices and m edges, computes a spectral sparsifier of G with O ( n log n ) edges in a single pass over G, using only O ( n log n ) space, and O ( m log 2 n ) total time. This improves on previous best semi-streaming algorithms for both spectral and cut sparsifiers by a factor of log n in both space and runtime. The algorithm also extends to semi-streaming row sampling for general PSD matrices. As another application, we use this framework to combine a spectral sparsification algorithm by Koutis with improved spanner constructions to give a parallel algorithm for constructing O ( n log 2 n log log n ) sized spectral sparsifiers in O ( m log 2 n log log n ) time. This is the best combinatorial graph sparsification algorithm to date, and the size of the sparsifiers produced is only a factor log n log log n more than ones produced by numerical routines.
STOC Conference 2017 Conference Paper
We present an algorithm that, with high probability, generates a random spanning tree from an edge-weighted undirected graph in ( n 5/3 m 1/3 ) time. The tree is sampled from a distribution where the probability of each tree is proportional to the product of its edge weights. This improves upon the previous best algorithm due to Colbourn et al. that runs in matrix multiplication time, O ( n ω ). For the special case of unweighted graphs, this improves upon the best previously known running time of Õ(min{ n ω , m √ n , m 4/3 }) for m ⪢ n 7/4 (Colbourn et al. '96, Kelner-Madry '09, Madry et al. '15). The effective resistance metric is essential to our algorithm, as in the work of Madry et al., but we eschew determinant-based and random walk-based techniques used by previous algorithms. Instead, our algorithm is based on Gaussian elimination, and the fact that effective resistance is preserved in the graph resulting from eliminating a subset of vertices (called a Schur complement). As part of our algorithm, we show how to compute -approximate effective resistances for a set S of vertex pairs via approximate Schur complements in Õ( m +( n + | S |)ε -2 ) time, without using the Johnson-Lindenstrauss lemma which requires Õ( min{( m + | S |) ε2 , m + n ε -4 +| S |ε 2 }) time. We combine this approximation procedure with an error correction procedure for handling edges where our estimate isn't sufficiently accurate.
FOCS Conference 2016 Conference Paper
We show how to perform sparse approximate Gaussian elimination for Laplacian matrices. We present a simple, nearly linear time algorithm that approximates a Laplacian by the product of a sparse lower triangular matrix with its transpose. This gives the first nearly linear time solver for Laplacian systems that is based purely on random sampling, and does not use any graph theoretic constructions such as low-stretch trees, sparsifiers, or expanders. Our algorithm performs a subsampled Cholesky factorization, which we analyze using matrix martingales. As part of the analysis, we give a proof of a concentration inequality for matrix martingales where the differences are sums of conditionally independent variables.
STOC Conference 2016 Conference Paper
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.
NeurIPS Conference 2015 Conference Paper
Given a directed acyclic graph $G, $ and a set of values $y$ on the vertices, the Isotonic Regression of $y$ is a vector $x$ that respects the partial order described by $G, $ and minimizes $\|x-y\|, $ for a specified norm. This paper gives improved algorithms for computing the Isotonic Regression for all weighted $\ell_{p}$-norms with rigorous performance guarantees. Our algorithms are quite practical, and their variants can be implemented to run fast in practice.
STOC Conference 2012 Conference Paper
We give a novel spectral approximation algorithm for the balanced (edge-)separator problem that, given a graph G, a constant balance b ∈ (0,1/2], and a parameter γ, either finds an Ω(b)-balanced cut of conductance O(√γ) in G, or outputs a certificate that all b-balanced cuts in G have conductance at least γ, and runs in time ~O(m). This settles the question of designing asymptotically optimal spectral algorithms for balanced separator. Our algorithm relies on a variant of the heat kernel random walk and requires, as a subroutine, an algorithm to compute exp(-L)v where L is the Laplacian of a graph related to G and v is a vector. Algorithms for computing the matrix-exponential-vector product efficiently comprise our next set of results. Our main result here is a new algorithm which computes a good approximation to exp(-A)v for a class of symmetric positive semidefinite (PSD) matrices A and a given vector v, in time roughly ~O(m A ), independent of the norm of A, where m A is the number of non-zero entries of A. This uses, in a non-trivial way, the result of Spielman and Teng on inverting symmetric and diagonally-dominant matrices in ~O(m A ) time. Finally, using old and new uniform approximations to e -x we show how to obtain, via the Lanczos method, a simple algorithm to compute exp(-A)v for symmetric PSD matrices that runs in time roughly O(t A ⋅ √norm(A)), where t A is the time required for the computation of the vector Aw for given vector w. As an application, we obtain a simple and practical algorithm, with output conductance O(√γ), for balanced separator that runs in time O(m/√γ). This latter algorithm matches the running time, but improves on the approximation guarantee of the Evolving-Sets-based algorithm by Andersen and Peres for balanced separator.
NeurIPS Conference 2012 Conference Paper
We present a new algorithm for Independent Component Analysis (ICA) which has provable performance guarantees. In particular, suppose we are given samples of the form $y = Ax + \eta$ where $A$ is an unknown $n \times n$ matrix and $x$ is chosen uniformly at random from $\{+1, -1\}^n$, $\eta$ is an $n$-dimensional Gaussian random variable with unknown covariance $\Sigma$: We give an algorithm that provable recovers $A$ and $\Sigma$ up to an additive $\epsilon$ whose running time and sample complexity are polynomial in $n$ and $1 / \epsilon$. To accomplish this, we introduce a novel ``quasi-whitening'' step that may be useful in other contexts in which the covariance of Gaussian noise is not known in advance. We also give a general framework for finding all local optima of a function (given an oracle for approximately finding just one) and this is a crucial step in our algorithm, one that has been overlooked in previous attempts, and allows us to control the accumulation of error when we find the columns of $A$ one by one via local search.