Much Faster Algorithms for Matrix Scaling

Abstract

We develop several efficient algorithms for the classical Matrix Scaling problem, which is used in many diverse areas, from preconditioning linear systems to approximation of the permanent. On an input n × n matrix A, this problem asks to find diagonal (scaling) matrices X and Y (if they exist), so that XAY ε-approximates a doubly stochastic matrix, or more generally a matrix with prescribed row and column sums.We address the general scaling problem as well as some important special cases. In particular, if A has m nonzero entries, and if there exist X and Y with polynomially large entries such that XAY is doubly stochastic, then we can solve the problem in total complexity O(m + n 4/3 ). This greatly improves on the best known previous results, which were either O(n 4 ) or O(mn 1/2 /ε).Our algorithms are based on tailor-made first and second order techniques, combined with other recent advances in continuous optimization, which may be of independent interest for solving similar problems. arXiv:1704.02315v1 [cs.DS] 7 Apr 2017 problems. We say that A is asymptotically (r, c)-scalable if the row and column sums can reach r and c asymptotically: that is, for every > 0, there exist positive diagonal matrices X, Y such that, letting B = XAY , we have B1 − r ≤ ε and 1 B − c ≤ ε. 2 The combinatorial essence of asymptotic scaling follows from a well-known characterization (see Proposition 2.2). A matrix A is asymptotically (1, 1)-scalable if and only if the permanent of A is positive, namely if the bipartite graph defined by the positive entries in A has a perfect matching. A matrix A is asymptotically (r, c)-scalable if and only if a natural flow on the same bipartite graph 3 has a solution. Duality (Hall's theorem and max-flow-min-cut theorem) gives simple certificates of non-scalability in terms of the patterns of 0's in the matrix A.The main computational problem we study is: given a matrix A, vectors r, c and ε > 0, determine if A is ε-approximately (r, c) scalable, and if so, find the scaling matrices X, Y .Before diving into the history of matrix scaling, we explain one of its most basic applications, which also demonstrates its algorithmic importance.Preconditioning Linear Systems. When solving a linear system Az = b, it is often desirable -for numerical stability and efficiency purposes-to have matrix A be well-conditioned. When this is not the case, one tries to transform A into a "better conditioned" matrix A . Matrix scaling provides a natural and efficient reduction to do so. For instance, one would hope that a scaled matrix A , in which e.g. all row and column p-norms are (say) 1, is better conditioned. 4 For this reason, we can use a matrix scaling algorithm to obtain diagonal matrices X, Y , and define A = XAY . Now, the solution to Az = b can be obtained by solving the (hopefully more numerically stable) linear system A z = Xb and setting z = Y −1 z . We stress here that A and A have the same sparsity. History and Prior WorkThe matrix (r, c)-scaling problem is so natural and important that it was discove...

Spectral Analysis of Matrix Scaling and Operator Scaling

“…Our results imply that the continuous operator scaling algorithm can be used to find a fractional perfect matching in an almost regular bipartite expander graph. We remark that our results also imply that the second-order methods for matrix scaling in [13,2] are near linear time algorithms for the instances in Corollary 4.6. This is because the condition number κ of the scaling solution for those instances is a constant by Theorem 1.7 and the algorithms in [13,2] have time complexity O(|E| log κ).…”

“…The dependency on n in these algorithms is at least Ω(n 7/2 ) even for sparse matrices. Recently, two independent groups [13,2] developed a fast second order method for matrix scaling, and this method is extended to geodesic convex optimization in [1] for the operator scaling problem. Theorem 1.3 ([13, 2, 1]).…”

Towards a Theory of Non-Commutative Optimization: Geodesic 1st and 2nd Order Methods for Moment Maps and Polytopes

“…Here we discuss our second order algorithm for Problem 1.10, the approximate norm minimization problem. As mentioned in Section 1.4, the paper [AZGL + 18] (following the algorithms developed in [AZLOW17,CMTV17] for the commutative Euclidean case) developed a second order polynomialtime algorithm for approximating the capacity for the simultaneous left-right action (Example 1.5) with running time polynomial in the bit description of the approximation parameter ε. In Section 5, we generalize this algorithm to arbitrary groups and actions (Algorithm 5.1).…”

“…Our second order method greatly generalizes the one used for the particular group action corresponding to operator scaling in [AZGL + 18]. It may be thought of as a geodesic analog of the "trust region method" [CGT00] or the "box-constrained Newton method" [CMTV17,AZLOW17] applied to a regularized function. For both methods, in this non-commutative setting, we recover the familiar convergence behavior of the classical commutative case: to achieve "proximity" ε to the optimum, our first order method converges in O(1/ε) iterations and our second order method in O(poly log(1/ε)) iterations.…”

Faster p-norm minimizing flows, via smoothed q-norm problems

“…The seminal work of Spielman and Teng [ST04] gave the first nearly-linear time algorithm for solving the weighted ℓ 2 version to highaccuracy (a (1+ε)-approximate solution in time O(m·log 1 ε ) 1 ). The work of Spielman-Teng and the several followup-works have led to the fastest algorithms for maximum matching [Mad13], shortest paths with negative weights [Coh+17b], graph partitioning [OSV12], sampling random spanning trees [KM09; MST15;Sch18], matrix scaling [Coh+17a; All+17], and resulted in dramatic progress on the problem of computing maximum flows.…”