Sparsified Cholesky and Multigrid Solvers for Connection Laplacians

Introduction

We introduce the sparsified Cholesky and sparsified multigrid algorithms for solving systems of linear equations. Two advantages of these algorithms over other recently introduced nearly-linear time algorithms for solving systems of equations in Laplacian matrices [Vai90, ST14, KMP10, KMP11, KOSZ13, CKM+14] are:

They give nearly-linear time algorithms for solving systems of equations in a much broader class of matrices—the connection Laplacians and Hermitian block diagonally dominant matrices. Connection Laplacians [SW12, BSS13] are a generalization of graph Laplacians that arise in many applications, including celebrated work on cryo-electron microscopy [SS11a, SS12, ZS14], phase retrieval [ABFM14, MTW14], and many image processing problems (e.g. [OSB15, ANKKS+12]). Previous algorithms for solving systems of equations in graph Laplacians cannot be extended to solve equations in connection Laplacians because the previous algorithms relied on some form of low stretch spanning trees, a concept that has no analog for the more general connection Laplacians.

They provide linear-sized approximate inverses of connection Laplacian matrices. That is, for every $n$-dimensional connection Laplacian $\boldsymbol{\mathit{M}}$, the sparsified Cholesky factorization algorithm produces an block-upper-triangular matrix $\boldsymbol{\mathit{U}}$ and block diagonal $\boldsymbol{\mathit{D}}$ with $O(n)$ nonzero entries so that $\boldsymbol{\mathit{U}}^{T}\boldsymbol{\mathit{D}}\boldsymbol{\mathit{U}}$ is a constant-factor approximation of $\boldsymbol{\mathit{M}}$. Such matrices $\boldsymbol{\mathit{U}}$ and $\boldsymbol{\mathit{D}}$ allow one to solve systems of equations in $\boldsymbol{\mathit{M}}$ to accuracy $\epsilon$ in time $O(m\log\epsilon^{-1})$, where $m$ is the number of nonzero entries of $\boldsymbol{\mathit{M}}$. Even for ordinary Laplacian matrices, the existence of such approximate inverses is entirely new. The one caveat of this result is that we do not yet know how to compute those approximate inverses in nearly linear time.

The sparsified Cholesky and sparsified multigrid algorithms work by sparsifying the matrices produced during Gaussian elimination. Recall that Cholesky factorization is the version of Gaussian elimination for symmetric matrices, and that the high running time of Gaussian elimination comes from fill—new nonzero matrix entries that are created by row operations. If Gaussian elimination never produced rows with a super-constant number of entries, then it would run in linear time. The sparsified Cholesky algorithm accelerates Gaussian elimination by sparsifying the rows that are produced by elimination, thereby guaranteeing that the elimination will be fast. Sparsified Cholesky is inspired by one of the major advances in algorithms for solving linear equations in Laplacian matrices—the Incomplete Cholesky factorization (ICC) [MV77]. However, ICC merely drops entries produced by elimination, whereas sparsification also carefully increases ones that remain. The difference is crucial, and is why ICC does not provide a nearly-linear time solver.

To control the error introduced by sparsification, we have to be careful not to do it too often. This means that our algorithm actually chooses a large set of rows and columns to eliminate at once, and then sparsifies the result. This is the basis of our first algorithms, which establish the existence of linear time solvers and linear sized approximate inverses, after precomputation. This precomputation is analogous to the procedure of computing a matrix inverse: the approximate inverse takes much less time to apply than to compute.

To produce entire algorithms that run in nearly linear time (both to compute and apply the approximate inverse) requires a little more work. To avoid the work of computing the matrix obtained by eliminating the large set of rows and columns, we design a fast algorithm for approximating it quickly. This resulting algorithm produces a solver routine that does a little more work at each level, and so resembles a multigrid V-cycle [TOS00]. We call the resulting algorithm the sparsified multigrid. We note that Krishnan, Fattal, and Szeliski [KFS13] present experimental results from the use of a sparsification heuristic in a multigrid algorithm for solving problems in computer vision.

Our new algorithms are most closely related to the Laplacian solver recently introduced by Peng and Spielman [PS14]: unlike the other Laplacian solvers, they rely only on sparsification and do not directly rely on “support theory” preconditioners or any form of low stretch spanning trees. This is why our algorithms can solve a much broader family of linear systems. To sparsify without using graph theoretic algorithms, we employ a recently developed algorithm of Cohen et. al. [CLM+14] that allows us to sparsify a matrix by solving systems of equations in a subsampled matrix.

In this section, we define block diagonally dominant (bDD) matrices—the most general family of matrices such that the associated systems of linear equations can be solved by our algorithms. We begin by defining our motivating case: the connection Laplacians

Connection Laplacians may be thought of as a generalization of graph Laplacians where every vertex is associated with a vector, instead of a real number, and every edge is associated with a unitary matrix. Like graph Laplacians, they describe a natural quadratic form. Let ${\boldsymbol{\mathit{M}}}_{[i,j]}\in\mathds{C}^{r\times r}$ be the unitary matrix associated with edge $(i,j)$, and let $w_{i,j}$ be the (nonnegative) weight of edge $(i,j)$. We require that ${\boldsymbol{\mathit{M}}}_{[i,j]}={\boldsymbol{\mathit{M}}}_{[j,i]}^{*}$, where $*$ denotes the conjugate transpose. The quadratic form associated with this connection Laplacian is a function of vectors $\boldsymbol{\mathit{v}}^{(i)}\in\mathds{C}^{r}$, one for each vertex $i$, that equals

The matrix corresponding to this quadratic form is a block matrix with blocks ${\boldsymbol{\mathit{M}}}_{[i,j]}$. Most applications of the connection Laplacian require one to either solve systems of linear equations in this matrix, or to compute approximations of its smallest eigenvalues and eigenvectors. By applying the inverse power method (or inverse Lanczos), we can perform these eigenvector calculations by solving a logarithmic number of linear systems in the matrix (see [ST14, Section 7]).

The matrices obtained from the connection Laplacian are a special case of block diagonally dominant (bDD) matrices, which we now define. Throughout this paper, we consider block matrices having entries in $\mathds{C}^{r\times r},$ where $r>0$ is a fixed integer. We say that $\boldsymbol{\mathit{M}}\in\left({\mathds{C}^{r\times r}}\right)^{m\times n}$ has $m$ block-rows, and $n$ block-columns. For $i\in[m],j\in[n],$ we let ${\boldsymbol{\mathit{M}}}_{[i,j]}\in\mathds{C}^{r\times r}$ denote the $i,j$-block in $\boldsymbol{\mathit{M}},$ and ${\boldsymbol{\mathit{M}}}_{[j]}$ denote the $j^{\textrm{th}}$ block-column, i.e., ${\boldsymbol{\mathit{M}}}_{[j]}=[({\boldsymbol{\mathit{M}}}_{[1,j]})^{*},({\boldsymbol{\mathit{M}}}_{[2,j]})^{*},\ldots,({\boldsymbol{\mathit{M}}}_{[m,j]})^{*}]^{*}.$ For sets $F\subseteq[m],C\subseteq[n],$ we let ${\boldsymbol{\mathit{M}}}_{[F,C]}\in\left({\mathds{C}^{r\times r}}\right)^{|F|\times|C|}$ denote the block-submatrix with blocks ${\boldsymbol{\mathit{M}}}_{[i,j]}$ for $i\in F,j\in C.$ $\boldsymbol{\mathit{M}}$ is block-diagonal, if ${\boldsymbol{\mathit{M}}}_{[i,j]}=0$ for $i\neq j.$ We emphasize that all computations are done over $\mathds{C}$, not over a matrix group.

To define bSDD matrices, let $\left\|\boldsymbol{\mathit{A}}\right\|$ denote the operator normRecall that the operator norm is the largest singular value of $\boldsymbol{\mathit{A}}$ and square root of the largest eigenvalue of $\boldsymbol{\mathit{A}}^{*}\boldsymbol{\mathit{A}}$. of a matrix $\boldsymbol{\mathit{A}}$. For Hermitian matrices $\boldsymbol{\mathit{A}},\boldsymbol{\mathit{B}},$ write $\boldsymbol{\mathit{A}}\succeq\boldsymbol{\mathit{B}}$ iff $\boldsymbol{\mathit{A}}-\boldsymbol{\mathit{B}}$ is positive semidefinite.

A Hermitian block-matrix $\boldsymbol{\mathit{M}}\in\left({\mathds{C}^{r\times r}}\right)^{n\times n}$ is block diagonally dominant (or bDD) if

Throughout the paper we treat $r$ as a constant. The hidden dependence on $r,$ of the running times of the algorithms we present, is polynomial. The major results of this paper are the following.

There is an algorithm that, when given a bDD matrix $\boldsymbol{\mathit{M}}$ with $n$ block rows and $m$ nonzero blocks, produces a solver for $\boldsymbol{\mathit{M}}$ in $O(m\log n+n\log^{2+o(1)}n)$ work and $O(n^{o(1)})$ depth so that the solver finds $\epsilon$-approximate solutions to systems of equations in $\boldsymbol{\mathit{M}}$ in $O((m+n\log^{1+o(1)}n)\log(1/\epsilon))$ work and $O(\log^{2}n\log\log n\log(1/\epsilon))$ depth.

Previously, the existence of nearly-linear time solvers was even unknown for the 1-dimensional case ($r=1$) when the off-diagonal entries were allowed to be complex numbers.

We can use the above algorithm to find approximations of the smallest eigenvalues and eigenvectors of such matrices at an additional logarithmic cost.

For every bDD matrix $\boldsymbol{\mathit{M}}$ with $n$ block-rows there exists a diagonal matrix $\boldsymbol{\mathit{D}}$ and an upper triangular matrix $\boldsymbol{\mathit{U}}$ with $O(n)$ nonzero blocks so that

Moreover, linear equations in $\boldsymbol{\mathit{U}}$, $\boldsymbol{\mathit{U}}^{T}$, and $\boldsymbol{\mathit{D}}$ can be solved with linear work in depth $O(\log^{2}n)$, and these matrices can be computed in polynomial time.

These matrices allow one to solve systems of linear equations in $\boldsymbol{\mathit{M}}$ to $\epsilon$-accuracy in parallel time $O(\log^{2}n\log^{-1}\epsilon)$ and work $O(m\log^{-1}\epsilon)$. Results of this form were previously unknown even for graph Laplacians.

In the next two sections we explain the ideas used to prove these theorems. Proofs may be found in the appendices that follow.

Background

We require some standard facts about the $\preccurlyeq$ order on matrices.

For $\boldsymbol{\mathit{A}}$ and $\boldsymbol{\mathit{B}}$ positive definite, $\boldsymbol{\mathit{A}}\succcurlyeq\boldsymbol{\mathit{B}}$ if and only if $\boldsymbol{\mathit{B}}^{-1}\succcurlyeq\boldsymbol{\mathit{A}}^{-1}$.

If $\boldsymbol{\mathit{A}}\succcurlyeq\boldsymbol{\mathit{B}}$ and $\boldsymbol{\mathit{C}}$ is any matrix of compatible dimension, then $\boldsymbol{\mathit{C}}\boldsymbol{\mathit{A}}\boldsymbol{\mathit{C}}^{T}\succcurlyeq\boldsymbol{\mathit{C}}\boldsymbol{\mathit{B}}\boldsymbol{\mathit{C}}^{T}$.

If $\boldsymbol{\mathit{A}}\approx_{\epsilon}\boldsymbol{\mathit{B}}$ and $\boldsymbol{\mathit{B}}\approx_{\delta}\boldsymbol{\mathit{C}}$, then $\boldsymbol{\mathit{A}}\approx_{\epsilon+\delta}\boldsymbol{\mathit{C}}$.

We now address one technicality of dealing with ${bDD}$ matrices: it is not immediate whether or not a ${bDD}$ matrix is singular. Moreover, if it is singular, the structure of its null space is not immediately clear either. Throughout the rest of this paper, we will consider ${bDD}$ matrices to which a small multiple of the identity have been added. These matrices will be nonsingular. To reduce the problem of solving equations in a general ${bDD}$ matrix $\boldsymbol{\mathit{M}}$ to that of solving equations in a nonsingular matrix, we require an estimate of the smallest nonzero eigenvalue of $\boldsymbol{\mathit{M}}$.

Hence, we can solve systems in $\boldsymbol{\mathit{M}}$ by approximately solving systems in $\boldsymbol{\mathit{M}}+\epsilon\mu\boldsymbol{\mathit{I}}$. Any lower bound $\mu$ on the smallest nonzero eigenvalue of $\boldsymbol{\mathit{M}}$ will suffice. It only impacts the running times of our algorithms in the numerical precision with which one must carry out the computations.

The above claim allows us to solve systems in $\boldsymbol{\mathit{M}}$ that have a solution. If $\boldsymbol{\mathit{M}}$ is singular and we want to apply its pseudoinverse (that is, to find the closest solution in the range of $\boldsymbol{\mathit{M}}$), we can do so by pre and post multiplying by $\boldsymbol{\mathit{M}}$ to project onto its range. The resulting algorithm, which is implicit in the following claim, requires applying a solver for $\boldsymbol{\mathit{M}}$ three times. It also takes $O(\log\kappa)$ times as long to run, where $\kappa$ is the finite condition numberThe finite condition number is the ratio of the largest singular value to the smallest nonzero singular value. of $\boldsymbol{\mathit{M}}$.

Let $\kappa$ be an upper bound on the finite condition number of $\boldsymbol{\mathit{M}}$. Given an error parameter $0<\epsilon<1$, let $\delta=\nicefrac{{\epsilon}}{{56\kappa^{3}}}$. If all nonzero eigenvalues of $\boldsymbol{\mathit{M}}$ are at least $\mu$, and if $\boldsymbol{\mathit{Z}}\approx_{\delta}(\boldsymbol{\mathit{M}}+\epsilon\mu\boldsymbol{\mathit{I}})^{-1}$, then

Overview of the algorithms

We index the block rows and columns of a block matrix by a set of vertices (or indices) $V$. When we perform an elimination, we eliminate a set $F\subset V$, and let $C=V-F$. Here, $F$ stands for “fine” and $C$ stands for “coarse”. In contrast with standard multigrid methods, we will have $\left|F\right|\leq\left|C\right|$.

To describe block-Cholesky factorization, we write the matrix $\boldsymbol{\mathit{M}}$ with the rows and columns in $F$ first:

Cholesky factorization writes the inverse of this matrix as

is the Schur complement of $\boldsymbol{\mathit{M}}$ with respect to $F$.

Our algorithms rely on two fundamental facts about ${bDD}$ matrices: that the Schur complement of a ${bDD}$ matrix is a ${bDD}$ matrix (Lemma B.3) and that one can sparsify ${bDD}$ matrices. The following theorem is implicit in [dCSHS11].

For every $\epsilon\leq 1$, every bDD matrix $\boldsymbol{\mathit{M}}$ with $n$ $r$-dimensional block rows and columns can be $\epsilon$-approximated by a bDD matrix having at most $10nr/\epsilon^{2}$ nonzero blocks.

We use the identity (1) to reduce the problem of solving a system of equations in $\boldsymbol{\mathit{M}}$ to that of solving equations in its Schur complement. The easiest part of this is multiplication by ${\boldsymbol{\mathit{M}}}_{[C,F]}$: the time is proportional to the number of nonzero entries in the submatrix, and we can sparsify $\boldsymbol{\mathit{M}}$ to guarantee that this is small.

The costlier part of the reduction is the application of the inverse of ${\boldsymbol{\mathit{M}}}_{[F,F]}$ three times. This would be fast if ${\boldsymbol{\mathit{M}}}_{[F,F]}$ were block-diagonal, which corresponds to $F$ being an independent set. We cannot find a sufficiently large independent set $F$, but we can find a large set that is almost independent. This results in a matrix ${\boldsymbol{\mathit{M}}}_{[F,F]}$ that is well approximated by its diagonal, and thus linear equations in this matrix can be quickly solved to high accuracy by a few Jacobi iterations (see Theorem 3.7 and Section E).

We can prove the existence of linear sized approximate inverses by explicitly computing $Sc\left(\boldsymbol{\mathit{M}},F\right)$, sparsifying it, and then recursively applying the algorithm just described. To make this algorithm efficient, we must compute a sparse approximation to $Sc\left(\boldsymbol{\mathit{M}},F\right)$ without constructing $Sc\left(\boldsymbol{\mathit{M}},F\right)$. This is the problem of spectral vertex sparsification, and we provide a fast algorithm for this task in Sections 3.4 and G.

We encode this recursive algorithm by a Schur complement chain (SCC). An SCC defines a linear operator that can be used to approximately solve equations in an initial matrix $\boldsymbol{\mathit{M}}^{(0)}$. If the ${bDD}$ matrix $\boldsymbol{\mathit{M}}$ is sparse, then it is the same as $\boldsymbol{\mathit{M}}^{(0)}$; if not, then $\boldsymbol{\mathit{M}}^{(0)}$ is a sparse approximation to $\boldsymbol{\mathit{M}}$. Let $F_{1}$ be the first set of vertices eliminated, $\boldsymbol{\mathit{M}}^{(1)}$ a sparse approximation to the Schur complement of $\boldsymbol{\mathit{M}}^{(0)}$ with respect to $F_{1},$ and $\boldsymbol{\mathit{Z}}^{(1)}$ an operators that approximates the inverse of ${\boldsymbol{\mathit{M}}}_{[F_{1},F_{1}]}^{(0)}.$

An $\boldsymbol{\epsilon}$-Schur complement chain ($\boldsymbol{\epsilon}$-SCC) for a matrix $\boldsymbol{\mathit{M}}^{{0}}$ indexed by vertex set $V$ is a sequence of operators and subsets,

so that for $C_{0}=V$ and $C_{i+1}=C_{i}\setminus F_{i+1}$, $\left|C_{d}\right|\leq 1000$ and for $1\leq i\leq d$,

The algorithm ApplyChain, described in Section C, applies an SCC to solve equations in $\boldsymbol{\mathit{M}}^{(0)}$ in the natural way and satisfies the following guarantee.

Consider an $\boldsymbol{\epsilon}$-SCC for $\boldsymbol{\mathit{M}}^{(0)}$ where $\boldsymbol{\mathit{M}}^{(i)}$ and $\boldsymbol{\mathit{Z}}^{(i)}$ can be applied to a vector in work $W_{\boldsymbol{\mathit{M}}^{(i)}},W_{\boldsymbol{\mathit{D}}^{(i)}}$ and depth $D_{\boldsymbol{\mathit{M}}^{(i)}},D_{\boldsymbol{\mathit{Z}}^{(i)}}$ respectively.

The algorithm $\textsc{ApplyChain}\left(\small(\boldsymbol{\mathit{M}}^{(1)},\boldsymbol{\mathit{Z}}^{(1)}),\ldots,(\boldsymbol{\mathit{M}}^{(d)},\boldsymbol{\mathit{Z}}^{(d)});F_{1},\ldots,F_{d};\boldsymbol{\mathit{b}}\right)$ corresponds to a linear operator $\boldsymbol{\mathit{W}}$ acting on $\boldsymbol{\mathit{b}}$ such that

$\small\boldsymbol{\mathit{W}}^{-1}\approx_{\sum_{i=1}^{d}2\epsilon_{i}}\boldsymbol{\mathit{M}}^{(0)},$ and

for any vector $\boldsymbol{\mathit{b}}$, $\textsc{ApplyChain}\left(\small(\boldsymbol{\mathit{M}}^{(1)},\boldsymbol{\mathit{Z}}^{(1)}),\ldots,(\boldsymbol{\mathit{M}}^{(d)},\boldsymbol{\mathit{Z}}^{(d)});F_{1},\ldots,F_{d};\boldsymbol{\mathit{b}}\right)$ runs in $O\left(\small\sum_{i=1}^{d}\left(D_{\boldsymbol{\mathit{M}}^{(i)}}+D_{\boldsymbol{\mathit{Z}}^{(i)}}\right)\right)$ depth and $O\left(\small\sum_{i=1}^{d}\left(W_{\boldsymbol{\mathit{M}}^{(i)}}+W_{\boldsymbol{\mathit{Z}}^{(i)}}\right)\right)$ work.

For an $\boldsymbol{\epsilon}$-SCC chain for $\boldsymbol{\mathit{M}}^{(0)},$ we define the depth and work of the chain to be the work and depth required by ApplyChain to apply the SCC.

We must choose the set of vertices $F_{i}$ so that we can approximate the inverse of ${\boldsymbol{\mathit{M}}}_{[F_{i},F_{i}]}^{(i)}$ by an operator $\boldsymbol{\mathit{Z}}^{(i)}$ that is efficiently computable. We do this by requiring that the matrix ${\boldsymbol{\mathit{M}}}_{[F_{i},F_{i}]}^{(i)}$ be $\alpha$-block diagonally dominant ($\alpha$-bDD), a term that we now define.

A Hermitian block-matrix $\boldsymbol{\mathit{M}}$ is $\alpha$-bDD if

We remark that a -bDD matrix is simply a bDD matrix. In particular, for $r=1,$ Laplacian matrices are -bDD.

By picking a subset of rows at random and discarding those that violate condition (2), the algorithm bDDSubset (described in Section D) finds a linear sized subset $F$ of the block-rows of a bDD matrix $\boldsymbol{\mathit{M}}$ so that ${\boldsymbol{\mathit{M}}}_{[F,F]}$ is $\alpha$-bDD.

Given a bDD matrix $\boldsymbol{\mathit{M}}$ with $n$ block-rows, and an $\alpha\geq 0$, bDDSubset computes a subset $F$ of size least $n/(8(1+\alpha))$ such that ${\boldsymbol{\mathit{M}}}_{[F,F]}$ is $\alpha$-bDD. It runs in runs in $O(m)$ expected work and $O(\log{n})$ expected depth, where $m$ is the number of nonzero blocks in $\boldsymbol{\mathit{M}}$.

We can express an $\alpha$-bDD matrix as a sum of a block diagonal matrix and a bDD matrix so that it is well-approximated by the diagonal.

Every $\alpha$-bDD matrix $\boldsymbol{\mathit{M}}$ can be written in the form $\boldsymbol{\mathit{X}}+\boldsymbol{\mathit{L}}$ where $\boldsymbol{\mathit{X}}$ is block-diagonal, $\boldsymbol{\mathit{L}}$ is bDD, and $\boldsymbol{\mathit{X}}\succcurlyeq\frac{\alpha}{2}\boldsymbol{\mathit{L}}$.

As $\boldsymbol{\mathit{L}}$ is positive semidefinite, $\left(1+\nicefrac{{2}}{{\alpha}}\right)\boldsymbol{\mathit{X}}\succcurlyeq\boldsymbol{\mathit{M}}\succcurlyeq\boldsymbol{\mathit{X}},$ which means that $\boldsymbol{\mathit{X}}$ is a good approximation of $\boldsymbol{\mathit{M}}$ when $\alpha$ is reasonably big. As block-diagonal matrices like $\boldsymbol{\mathit{X}}$ are easy to invert, systems in these well-conditioned matrices can be solved rapidly using preconditioned iterative methods. In Section E, we show that a variant of Jacobi iteration provides an operator that satisfies the requirements of Definition 3.2.

Let $\boldsymbol{\mathit{M}}$ be a bDD matrix with index set $V$, and let $F\subseteq V$ such that ${\boldsymbol{\mathit{M}}}_{[F,F]}$ is $\alpha$-bDD for some $\alpha\geq 4,$ and has $m_{FF}$ nonzero blocks. The algorithm $\textsc{Jacobi}(\epsilon,{\boldsymbol{\mathit{M}}}_{[F,F]},\boldsymbol{\mathit{b}})$ acts as a linear operator $\boldsymbol{\mathit{Z}}$ on $\boldsymbol{\mathit{b}}$ that satisfies

The algorithm takes $O(m_{FF}\log(\frac{1}{\epsilon}))$ work and $O(\log{n}\log(\frac{1}{\epsilon}))$ depth.

We now explain how Theorems 3.7 and 3.1 allow us to construct Schur complement chains that can be applied in nearly linear time. We optimize the construction in the next section.

Theorem 3.1 tells us that there is a ${bDD}$ matrix $\boldsymbol{\mathit{M}}^{(0)}$ with $O(n/\epsilon^{2})$ nonzero blocks that that $\epsilon$-approximates $\boldsymbol{\mathit{M}}$, and that for every $i$ there is a ${bDD}$ matrix $\boldsymbol{\mathit{M}}^{(i+1)}$ with $O(\left|C_{i}\right|/\epsilon^{2})$ nonzero blocks that is an $\epsilon$-approximation of $Sc\left(\boldsymbol{\mathit{M}}^{(i)},F_{i}\right)$. We will pick $\epsilon$ later. Lemma 3.5 provides a set $F_{i}$ containing a constant fraction of the block-rows of $\boldsymbol{\mathit{M}}^{(i)}$ so that ${\boldsymbol{\mathit{M}}}_{[F_{i},F_{i}]}^{(i)}$ is $4$-bDD. Theorem 3.7 then provides an operator that solves systems in ${\boldsymbol{\mathit{M}}}_{[F_{i},F_{i}]}^{(i)}$ to $\epsilon$ accuracy in time $O(\log\epsilon^{-1})$ times the number of nonzero entries in ${\boldsymbol{\mathit{M}}}_{[F_{i},F_{i}]}^{(i)}$. This is at most the number of nonzero entries in $\boldsymbol{\mathit{M}}^{(i)}$, and thus at most $O(\left|C_{i}\right|/\epsilon^{2})$. As each $F_{i}$ contains at least a constant fraction of the rows of $\boldsymbol{\mathit{M}}^{(i)}$, the depth of the recursion, $d$, is $O(\log n)$. Thus, we can obtain constant accuracy by setting $\epsilon=\Theta(1/\log n)$. The time required to apply the resulting SCC would thus be $O(n\log^{2}n\log\log n)$.

We can reduce the running time by setting $\epsilon_{1}$ to a constant and allowing it to shrink as $i$ increases. For example, setting $\epsilon_{i}=1/2(i+1)^{2}$ results in a linear time algorithm that produces a constant-factor approximation of the inverse of $\boldsymbol{\mathit{M}}$. We refine this idea in the next section.

3 Linear Sized Approximate Inverses

In this section we sketch the proof of Theorem 1.3, which tells us that every ${bDD}$ matrix has a linear-sized approximate inverse. The rest of the details appear in Section F.

The linear-sized approximate inverse of a ${bDD}$ matrix $\boldsymbol{\mathit{M}}$ with $n$ block rows and columns is provided by a $3/4$-approximate of the form $\boldsymbol{\mathit{U}}^{T}\boldsymbol{\mathit{D}}\boldsymbol{\mathit{U}}$ where $\boldsymbol{\mathit{D}}$ is block diagonal and $\boldsymbol{\mathit{U}}$ is block upper-triangular and has $O(n)$ nonzero blocks. As systems of linear equations in block-triangular matrices like $\boldsymbol{\mathit{U}}$ and $\boldsymbol{\mathit{U}}^{T}$ can be solved in time proportional to their number of nonzero blocks, this provides a linear time algorithm for computing a $3/4$ approximation of the inverse of $\boldsymbol{\mathit{M}}$. By iteratively refining the solutions provided by the approximate inverse, this allows one to find $\epsilon$-accurate solutions to systems of equations in $\boldsymbol{\mathit{M}}$ in time $O(m\log\epsilon^{-1})$.

The matrix $\boldsymbol{\mathit{U}}$ that we construct has meta-block structure that allows solves in $\boldsymbol{\mathit{U}}$ and $\boldsymbol{\mathit{U}}^{T}$ to be performed with linear work and depth $O(\log^{2}n)$. This results in a parallel algorithm for solving equations in $\boldsymbol{\mathit{M}}$ in work $O(m\log\epsilon^{-1})$ and depth $O(\log^{2}n\log\epsilon^{-1})$. We remark that with some additional work (analogous to Section 7 of [LPS15]), one could reduce this depth to $O(\log n\log\log n\log\epsilon^{-1})$.

The key to constructing $\boldsymbol{\mathit{U}}$ is realizing that the algorithm Jacobi corresponds to multiplying by the matrix $\boldsymbol{\mathit{Z}}^{(k)}$ defined in equation (7). Moreover, this matrix is a polynomial in $\boldsymbol{\mathit{X}}$ and $\boldsymbol{\mathit{L}}$ of degree $\log(3/\epsilon)$, where ${\boldsymbol{\mathit{M}}}_{[F_{i},F_{i}]}^{(i)}=\boldsymbol{\mathit{X}}+\boldsymbol{\mathit{L}}$ where $\boldsymbol{\mathit{L}}$ is bDD, and $\boldsymbol{\mathit{X}}$ block diagonal such that $\boldsymbol{\mathit{X}}\succcurlyeq 2\boldsymbol{\mathit{L}}.$ To force $\boldsymbol{\mathit{Z}}^{k}$ to be a sparse matrix, we require that $\boldsymbol{\mathit{L}}$ be sparse.

If we use the algorithm bSDDSubset to choose $F_{i}$, then $\boldsymbol{\mathit{L}}$ need not be sparse. However, this problem is easily remedied by forbidding algorithm bSDDSubset from choosing any vertex of more than twice average degree. Thus, we can ensure that all $\boldsymbol{\mathit{Z}}^{(i)}$ are sparse.

For every bDD matrix $\boldsymbol{\mathit{M}}$ and every $\alpha\geq 0$, there is a subset $F$ of size at least $\frac{n}{16(1+\alpha)}$ such that ${M}_{[F,F]}$ is $\alpha$-bDD and the number of nonzeros blocks in each block-row of $F$ at most twice the average number of nonzero blocks in a block-row of $\boldsymbol{\mathit{M}}$.

Discard every block-row of $\boldsymbol{\mathit{M}}$ that has more than twice the average number of nonzeros blocks per row-block. Then remove the corresponding row blocks. The remaining matrix has dimension at least $n/2$. We can now use Lemma 3.5 to find an $\alpha$-bDD submatrix. ∎

We obtain the $\boldsymbol{\mathit{U}}^{T}\boldsymbol{\mathit{D}}\boldsymbol{\mathit{U}}$ factorization by applying the inverse of the factorization (1):

In the left and right triangular matrices we replace ${\boldsymbol{\mathit{M}}}_{[F,F]}^{-1}$ with the polynomial we obtain from Jacobi. In the middle matrix, it suffices to approximate ${\boldsymbol{\mathit{M}}}_{[F,F]}$ by a block diagonal matrix, and $Sc\left(\boldsymbol{\mathit{M}},F\right)$ by a factorization of its sparse approximation given by Theorem 3.1. The details, along with a careful setting of $\epsilon_{i}$, are carried out in Section F.

4 Spectral Vertex Sparsification Algorithm

In this section, we outline a procedure ApproxSchur that efficiently approximates the Schur complement of a bDD matrix $\boldsymbol{\mathit{M}}$ w.r.t. a set of indices $F$ s.t. ${\boldsymbol{\mathit{M}}}_{[F,F]}$ is $\alpha$-bDD. The following lemma summarizes the guarantees of ApproxSchur.

Let $\boldsymbol{\mathit{M}}$ be a bDD matrix with index set $V$, and $m$ nonzero blocks. Let $F\subseteq V$ be such that ${\boldsymbol{\mathit{M}}}_{[F,F]}$ is $\alpha$-bDD for some $\alpha\geq 4.$ The algorithm $\textsc{ApproxSchur}(\boldsymbol{\mathit{M}},F,\epsilon)$, returns a matrix $\boldsymbol{\widetilde{\mathit{M}}}_{SC}$ s.t.

$\boldsymbol{\widetilde{\mathit{M}}}_{SC}$ has $O(m(\epsilon^{-1}\log\log\epsilon^{-1})^{O(\log\log\epsilon^{-1})})$ nonzero blocks, and

$\boldsymbol{\widetilde{\mathit{M}}}_{SC}\approx_{\epsilon}Sc\left(\boldsymbol{\mathit{M}},F\right)$,

in $O(m(\epsilon^{-1}\log\log\epsilon^{-1})^{O(\log\log\epsilon^{-1})})$ work and $O(\log n(\log\log\epsilon^{-1}))$ depth.

We sketch a proof of the above lemma in this section. A complete proof and pseudocode for ApproxSchur are given in Section G.

First consider a very simple special case: where $F$ is a singleton, $F=\{i\}$. Let $C=V\setminus F$ be the remaining indices.

Thus, if ${\boldsymbol{\mathit{M}}}_{[C,i]}$ has $k$ nonzero blocks, $Sc\left(\boldsymbol{\mathit{M}},i\right)$ could have $k^{2}$ additional nonzero blocks compared to ${\boldsymbol{\mathit{M}}}_{[C,C]}$, potentially making it dense. If $\boldsymbol{\mathit{M}}$ were a graph Laplacian, then ${\boldsymbol{\mathit{M}}}_{[C,i]}{\boldsymbol{\mathit{M}}}_{[i,i]}^{-1}{\boldsymbol{\mathit{M}}}_{[i,C]}$ would represent the adjacency matrix of a weighted clique. In Section J, we construct weighted expanders that allow us to $\epsilon$-approximate $Sc\left(\boldsymbol{\mathit{M}},i\right)$ in this case using $m+O(k\epsilon^{-4})$ edges. In Section G.4, we show how to use such weighted expanders to sparsify $Sc\left(\boldsymbol{\mathit{M}},i\right)$ when $\boldsymbol{\mathit{M}}$ is a bDD matrix.

This reduction can also be performed in parallel. If $F$ is such that ${\boldsymbol{\mathit{M}}}_{[F,F]}$ is block diagonal, we can approximate $Sc\left(\boldsymbol{\mathit{M}},F\right)$ by expressing ${\boldsymbol{\mathit{M}}}_{[C,i]}{\boldsymbol{\mathit{M}}}_{[i,i]}^{-1}{\boldsymbol{\mathit{M}}}_{[i,C]}$ as $\sum_{i\in F}{\boldsymbol{\mathit{M}}}_{[C,i]}{\boldsymbol{\mathit{M}}}_{[i,i]}^{-1}{\boldsymbol{\mathit{M}}}_{[i,C]},$ and using $|F|$ weighted expanders. However, ${\boldsymbol{\mathit{M}}}_{[F,F]}$ may not be diagonal. Instead, we give a procedure SchurSquare that generates $\boldsymbol{\mathit{M}}^{\prime}$ that is better approximated by its diagonal.

Invoking SchurSquare a few times leads a sequence of matrices ${\boldsymbol{\mathit{M}}}_{[F,F]}^{(0)},{\boldsymbol{\mathit{M}}}_{[F,F]}^{(1)},\ldots,{\boldsymbol{\mathit{M}}}_{[F,F]}^{(i)}$. We will show that ${\boldsymbol{\mathit{M}}}_{[F,F]}^{(i)}$ is $\epsilon$-approximated by its diagonal and we call the procedure LastStep to approximate $Sc\left(\boldsymbol{\mathit{M}}^{(i)},F\right)$. An additional caveat is that replacing ${\boldsymbol{\mathit{M}}}_{[F,F]}^{(i)}$ by its diagonal at this step gives errors that are difficult to bound. We discuss the correct approximation below

SchurSquare is based on a squaring identity for matrix inverse developed in [PS14]. Given a splitting of ${\boldsymbol{\mathit{M}}}_{[F,F]}$ into $\boldsymbol{\mathit{D}}-\boldsymbol{\mathit{A}},$ where $\boldsymbol{\mathit{D}}$ is block-diagonal, and $\boldsymbol{\mathit{A}}$ has its diagonal blocks as zero, it relies on the fact that the matrix

satisfies $Sc\left(\boldsymbol{\mathit{M}},F\right)=Sc\left(\boldsymbol{\mathit{M}}_{2},F\right)$. Furthermore, we can show that if $\boldsymbol{\mathit{M}}$ is $\alpha$-bDD, $\boldsymbol{\mathit{D}}-\boldsymbol{\mathit{A}}\boldsymbol{\mathit{D}}^{-1}\boldsymbol{\mathit{A}}$ is $\alpha^{2}$-bDD, which indicates that the block on $[F,F]$ rapidly approaches being diagonal.

As $\boldsymbol{\mathit{M}}_{2}$ may be dense, we construct a sparse approximation to it. Since $\boldsymbol{\mathit{D}}$ is diagonal, we can construct sparse approximations to the blocks ${(\boldsymbol{\mathit{M}}_{2})}_{[F,F]}$ and ${(\boldsymbol{\mathit{M}}_{2})}_{[C,C]}$ in a manner analogous to the case of diagonal ${\boldsymbol{\mathit{M}}}_{[F,F]}$. Similarly, we use bipartite expanders to construct sparse approximations to ${(\boldsymbol{\mathit{M}}_{2})}_{[C,F]}$ and ${(\boldsymbol{\mathit{M}}_{2})}_{[F,C]}$.

Our sequence of calls to SchurSquare terminates with ${\boldsymbol{\mathit{M}}}_{[F,F]}^{(i)}$ being roughly $\epsilon^{-1}$-bDD. We then return $\textsc{LastStep}(\boldsymbol{\mathit{M}}^{(i)},F,\epsilon^{-1},\epsilon)$. As mentioned above, we cannot just replace ${\boldsymbol{\mathit{M}}}_{[F,F]}^{(i)}$ by its diagonal. Instead, LastStep performs one step of squaring similar to SchurSquare with a key difference: Rather than expressing ${\boldsymbol{\mathit{M}}}_{[F,F]}^{(i)}$ as $\boldsymbol{\mathit{D}}-\boldsymbol{\mathit{A}},$ it expresses it as $\boldsymbol{\mathit{X}}+\boldsymbol{\mathit{L}},$ where $\boldsymbol{\mathit{X}}$ is block-diagonal, and $\boldsymbol{\mathit{L}}$ is just barely bDD. With this splitting, it constructs $\boldsymbol{\mathit{M}}^{\left(last\right)}$ after performing one iteration similar to Eq. (3). After this step, it replaces the ${\boldsymbol{\mathit{M}}}_{[F,F]}^{\left(last\right)}$ block with the block-diagonal matrix $\boldsymbol{\mathit{X}}$. Again, we directly produce sparse approximations to $\boldsymbol{\mathit{M}}^{\left(last\right)}$ and its Schur complement via weighted (bipartite) expanders. A precise description and proofs are given in Section G.2.

5 Sparsifying bDD matrices

The main technical hurdle left to address is how we sparsify bDD matrices. We to do this both to approximate $\boldsymbol{\mathit{M}}$ by $\boldsymbol{\mathit{M}}^{(0)}$, if $\boldsymbol{\mathit{M}}$ is not already sparse, and to ensure that all the matrices $\boldsymbol{\mathit{M}}^{(i)}$ remain sparse. While the spectral vertex sparsification algorithm described in the previous section allows us to compute an approximation to a Schur complement $Sc\left(\boldsymbol{\mathit{M}}^{(i)},F_{i}\right)$, it is sparse only when $\boldsymbol{\mathit{M}}^{(i)}$ is already sparse. As we iteratively apply this procedure, the density of the matrices produced will grow unacceptably. We overcome this problem by occasionally applying another sparsification routine that substantially decreases the number of nonzero blocks. The cost of this sparsification routine is that it requires solving systems of equations in sparse ${{bDD}}$ matrices. We, of course, do this recursively.

Our sparsification procedure begins by generalizing the observation that graph Laplacians can be sparsified by sampling edges with probabilities determined by their effective resistances [SS11b]. There is a block analog of leverage scores (37) that provides probabilities of sampling blocks so that the resulting sampled matrix approximates the original and has $O(n\log n)$ nonzero blocks with high probability. To compute this block analog of leverage scores we employ a recently introduced procedure of Cohen et. al. [CLM+14].

Once we generalize their results to block matrices, they show that we can obtain sufficiently good estimates of the block leverage scores by computing leverage scores in a bDD matrix obtained by randomly subsampling blocks of the original. The block leverage scores in this matrix are obtained by solving a logarithmic number of linear equations in this subsampled matrix. Thus, our sparsification procedure requires constructing a solver for a subsampled matrix and then applying that solver a logarithmic number of times. We compute this solver recursively.

There is a tradeoff between the number of nonzero blocks in the subsampled system and in the resulting approximation of the original matrix. If the original matrix has $m$ nonzero blocks and we subsample to a system of $m/K$ nonzero blocks, then we obtain an $\epsilon$-approximation of the original matrix with $O(K\epsilon^{-2}n\log n)$ nonzero blocks.

The details of the analysis of the undersampling procedure appear in Section H.

6 The main algorithm

We now explain how we prove Theorem 1.2. The details supporting this exposition appear in Section I. Our main goal is to control the density of the Schur complement chain as we repeatedly invoke Lemma 3.9.

Starting from some $\boldsymbol{\mathit{M}}^{(0)}$, we compute sets $F_{i}$ (via calls to bDDSubset), approximate solvers $\boldsymbol{\mathit{Z}}^{(i)}$ (via Jacobi), and approximations of Schur complements $\boldsymbol{\mathit{M}}^{(i)}$ (via ApproxSchur), until we obtain a matrix $\boldsymbol{\mathit{M}}^{(i)}$ such that its dimension is a smaller than that of $\boldsymbol{\mathit{M}}^{(0)}$ by a large constant factor (like 4). While the dimension of $\boldsymbol{\mathit{M}}^{(i)}$ is much smaller, its number of nonzero blocks is potentially larger by an even larger factor. We use the procedure described in the previous section to sparsify it. This sparsification procedure produces a sparse approximation of the matrix at the cost of solving systems of equations in a subsampled version of that matrix. Some care is required to balance the cost of the resulting recursion.

We now sketch an analysis of a nearly linear time algorithm that results from a simple choice of parameters. We optimize the parameter choice and analysis in Section I. Let $n$ be the dimension of $\boldsymbol{\mathit{M}}^{(0)}$ and let $m$ be its number of nonzero blocks. To begin, assume that $m\leq n\Delta\log^{3}n$, for a $\Delta$ to be specified later. We call this the sparse case, and address the case of dense $\boldsymbol{\mathit{M}}$ later.

We consider fixing $\epsilon_{i}=c/\log n$ for all $i$, for some constant $c$. As the depth of the Schur complement chain is $O(\log n)$, this results in a solver with constant accuracy. A constant number of iterations of the procedure described above are required to produce an $\boldsymbol{\mathit{M}}^{(i)}$ whose dimension is a factor of 4 smaller than $\boldsymbol{\mathit{M}}^{(0)}$. Lemma 3.9 tells us that the edge density of this $\boldsymbol{\mathit{M}}^{(i)}$ is potentially higher than that of $\boldsymbol{\mathit{M}}^{(0)}$ by a factor of

Set $\Delta$ to be this factor. We use the sparsification procedure from the previous section to guarantee that no matrix in the chain has density higher than that of this matrix, which is upper bounded by $(m/n)\Delta=\Delta^{2}\log^{3}n$.

Setting $K=2\Delta$, the subsampling produces a matrix of density half that of $\boldsymbol{\mathit{M}}^{(0)}$, and it produces a sparse approximation of $\boldsymbol{\mathit{M}}^{(i)}$ of density $O(K\epsilon_{i}^{-2}\log n)\leq O(\Delta\log^{3}n)$, which, by setting constants appropriately, we can also force to be half that of $\boldsymbol{\mathit{M}}^{(0)}$. In order to perform the sparsification procedure, we need to construct a Schur complement chain for the subsampled matrix, and then use it to solve $O(\log n)$ systems of linear equations. The cost of using this chain to solve equations in the subsampled system is at most $O(n\Delta^{2}\log^{4}n)$, and the cost of using the solutions to these equations to sparsify $\boldsymbol{\mathit{M}}^{(i)}$ is $O(m\log n)$. The cost of the calls to bDDSubset and ApproxSchur are proportional to the number of edges in the matrices, which is $O(n\Delta\log^{4}{n})$.

We repeat this procedure all the way down the chain, only using sparsification when the dimension of $\boldsymbol{\mathit{M}}^{(i)}$ shrinks by a factor of $4$. Since none of the matrices that we generate have density higher than $\Delta\log^{3}n$, we remain in the sparse case. Let $T_{sparse}(n)$ be the time required to construct a solver chain on systems of size $n$ with $m\leq n\Delta\log^{3}{n}$. We obtain the following recurrence

To handle the case of dense $\boldsymbol{\mathit{M}}$, we repeatedly sparsify while keeping $n$ fixed until we obtain a matrix with fewer than $n\Delta\log^{3}n$ edges, at which point we switch to the algorithm described above. The running time of this algorithm on a graph with $m$ edges, $T_{dense}(m)$, satisfies the recurrence

Thus $T_{dense}(m)$ is upper bounded bounded by $O(mn^{o(1)}+n^{1+o(1)})$.

We tighten this bound in Section I to prove Theorem 1.2 by carefully choosing the parameters to accompany a sequence $\boldsymbol{\epsilon}$ that starts constant and decreases slowly.

Summary

We introduce a new approach to solving systems of linear equations that gives the first nearly linear time algorithms for solving systems in connection Laplacians and the first proof that connection Laplacians have linear-sized approximate inverses. This was unknown even for graph Laplacians.

Our algorithms build on ideas introduced in [PS14] and are a break from those used in the previous work on solving systems of equations in graph Laplacians [Vai90, ST14, KMP10, KMP11, KOSZ13, CKM+14]. Those algorithms all rest on support theory [BGH+06], originally introduced by Vaidya [Vai90], and rely on the fact that the Laplacian of one edge is approximated by the Laplacian of a path between its endpoints. No analogous fact is true for connection Laplacians, even those with complex entries for $r=1$.

Instead, our algorithms rely on many new ideas, the first being that of sparsifying the matrices that appear during elimination. Other critical ideas are finding $\alpha$-bDD subsets of vertices to eliminate in bulk, approximating Schur complements without computing them explicitly, and the use of sub-sampling to sparsify in a recursive fashion. To efficiently compute approximations of the Schur complements, we introduce a new operation that transforms a matrix into one with the same Schur complement but a much better conditioned upper block (3). To obtain the sharp bounds in Theorem 1.2, we exploit a new linear-time algorithm for constructing linear-sized sparse approximations to implicitly represented weighted cliques whose edge weights are products of weights at vertices (Section J), and extend this to the analog for bDD matrices (Section G.4).

References

Appendix A Background

Since all nonzero eigenvalues of $\boldsymbol{\mathit{M}}$ are at least $\mu$, the eigenvalues of $\boldsymbol{\mathit{M}}^{1/2}(\boldsymbol{\mathit{M}}+\epsilon\mu\boldsymbol{\mathit{I}})^{-1}\boldsymbol{\mathit{M}}^{1/2}$ lie between $1/(1+\epsilon)$ and 1. Using $\boldsymbol{\mathit{Z}}\approx_{\epsilon}(\boldsymbol{\mathit{M}}+\epsilon\mu\boldsymbol{\mathit{I}})^{-1}$, we see that the eigenvalues of $\boldsymbol{\mathit{M}}^{1/2}\boldsymbol{\mathit{Z}}\boldsymbol{\mathit{M}}^{1/2}$ lie between $e^{-2\epsilon}$ and $e^{\epsilon}$. Using $0<\epsilon<1/2$, we have

Let $\boldsymbol{\mathit{A}}$ be a matrix of condition number $\kappa$ and let $\boldsymbol{\mathit{A}}\approx_{\epsilon}\boldsymbol{\mathit{B}}$ for $\epsilon\leq(56\kappa^{3})^{-1}$. Then, $\boldsymbol{\mathit{A}}^{3}\approx_{28\kappa^{3}\epsilon}\boldsymbol{\mathit{B}}^{3}$.

First, observe that $\boldsymbol{\mathit{A}}\approx_{\epsilon}\boldsymbol{\mathit{B}}$ implies that $\left\|\boldsymbol{\mathit{B}}\right\|\leq(1+e^{\epsilon})\left\|\boldsymbol{\mathit{A}}\right\|\leq 2\left\|\boldsymbol{\mathit{A}}\right\|$. It also implies that $\left\|\boldsymbol{\mathit{A}}-\boldsymbol{\mathit{B}}\right\|\leq 2\epsilon\left\|\boldsymbol{\mathit{A}}\right\|$. As

$\left\|\boldsymbol{\mathit{A}}^{2}-\boldsymbol{\mathit{B}}^{2}\right\|\leq 6\epsilon\left\|\boldsymbol{\mathit{A}}\right\|^{2}$. Similarly, as

$\left\|\boldsymbol{\mathit{A}}^{3}-\boldsymbol{\mathit{B}}^{3}\right\|\leq 28\epsilon\left\|\boldsymbol{\mathit{A}}\right\|^{3}$.

Let $\kappa=\left\|\boldsymbol{\mathit{A}}\right\|/\lambda_{min}(\boldsymbol{\mathit{A}})$. The above relation implies that

By Claim A.1, using $\delta=\frac{\epsilon}{56\kappa^{3}}$,

$\boldsymbol{\mathit{M}}$ has an eigendecomposition in the same basis as $(\boldsymbol{\mathit{M}}+\epsilon\mu I)^{-3}$, and so it follows that $\boldsymbol{\mathit{M}}(\boldsymbol{\mathit{M}}+\epsilon\mu I)^{-3}\boldsymbol{\mathit{M}}$ has the same eigenbasis, and the same null space as $\boldsymbol{\mathit{M}}$.

So $\boldsymbol{\mathit{M}}(\boldsymbol{\mathit{M}}+\epsilon\mu I)^{-3}\boldsymbol{\mathit{M}}\approx_{3\epsilon}\boldsymbol{\mathit{M}}^{+}$, and $\boldsymbol{\mathit{M}}\boldsymbol{\mathit{Z}}^{3}\boldsymbol{\mathit{M}}\approx_{4\epsilon}\boldsymbol{\mathit{M}}^{+}$. ∎

For every $\boldsymbol{\mathit{d}}\in\mathds{C}^{r\times r},$ we can find $\boldsymbol{\mathit{Q}}^{(1)},\boldsymbol{\mathit{Q}}^{(2)}\in\mathds{C}^{r\times r}$ where

Letting $\boldsymbol{\mathit{Q}}^{(k)}=\boldsymbol{\mathit{U}}\boldsymbol{\mathit{D}}^{(k)}\boldsymbol{\mathit{V}}^{*}$ for $k=1,2,$ we get $\hat{\boldsymbol{\mathit{d}}}=\frac{1}{2}(\boldsymbol{\mathit{Q}}^{(1)}+\boldsymbol{\mathit{Q}}^{(2)})$ and hence $\boldsymbol{\mathit{d}}=\frac{1}{2}\left\|\boldsymbol{\mathit{d}}\right\|(\boldsymbol{\mathit{Q}}^{(1)}+\boldsymbol{\mathit{Q}}^{(2)}).$ Moreover,

If $\boldsymbol{\mathit{M}}$ and $\boldsymbol{\widetilde{\mathit{M}}}$ are positive semidefinite matrices satisfying $\boldsymbol{\mathit{M}}\preceq\boldsymbol{\widetilde{\mathit{M}}}$, then

This fact can be proven via an energy minimization definition of Schur complement. More details on this formulation can be found in [MP13].

Appendix B Block Diagonally Dominant Matrices

In this section, we prove a few basic facts about bDD matrices. The following lemma gives an equivalent definition of bDD matrices.

The only if direction is easy. For bDD matrix $\boldsymbol{\mathit{M}},$ if we let $\boldsymbol{\mathit{D}}$ be the block diagonal matrix such that ${\boldsymbol{\mathit{D}}}_{[i,i]}={\boldsymbol{\mathit{M}}}_{[i,i]},$ and $\boldsymbol{\mathit{A}}$ be the matrix such that

where we used that since $\boldsymbol{\mathit{A}},\boldsymbol{\mathit{M}}$ are Hermitian, $\boldsymbol{\mathit{D}}$ is also Hermitian, and thus ${\boldsymbol{\mathit{D}}}_{[i,i]},{\boldsymbol{\mathit{A}}}_{[i,i]}$ are Hermitian. ∎

This immediately implies the corollary that flipping the sign of off-diagonal blocks preserves bDD-ness.

Given a bDD matrix $\boldsymbol{\mathit{M}},$ write it as $\boldsymbol{\mathit{D}}-\boldsymbol{\mathit{A}},$ where $\boldsymbol{\mathit{D}}$ is a block-diagonal, and $\boldsymbol{\mathit{A}}$ has its diagonal blocks as zero. Then, $\boldsymbol{\mathit{D}}+\boldsymbol{\mathit{A}}$ is also PSD.

First observe that for all $i,$ ${(\boldsymbol{\mathit{D}}+\boldsymbol{\mathit{A}})}_{[i,i]}={(\boldsymbol{\mathit{D}}-\boldsymbol{\mathit{A}})}_{[i,i]}={\boldsymbol{\mathit{M}}}_{[i,i]},$ i.e., their diagonal blocks are identical. Moreover, for all $i\neq j,$ we have $\left\|{(\boldsymbol{\mathit{D}}+\boldsymbol{\mathit{A}})}_{[i,j]}\right\|=\left\|{(\boldsymbol{\mathit{D}}-\boldsymbol{\mathit{A}})}_{[i,j]}\right\|=\left\|{\boldsymbol{\mathit{M}}}_{[i,j]}\right\|.$ Thus $\boldsymbol{\mathit{D}}+\boldsymbol{\mathit{A}}$ is also bDD, and hence PSD. ∎

Next, we show that the class of bDD matrices is closed under Schur complement.

The class of bDD matrices is closed under Schur complement.

Since Schur complementation does not depend on the order of indices eliminated, it suffices to prove that for any bDD matrix $\boldsymbol{\mathit{M}}\in(\mathds{C}^{r\times r})^{n\times n},$ $Sc\left(\boldsymbol{\mathit{M}},1\right)$ is a bDD matrix. Let $C=V\setminus\{1\}.$

We have $Sc\left(\boldsymbol{\mathit{M}},1\right)={\boldsymbol{\mathit{M}}}_{[C,C]}-{\boldsymbol{\mathit{M}}}_{[C,1]}{\boldsymbol{\mathit{M}}}_{}^{-1}{\boldsymbol{\mathit{M}}}_{[1,C]}.$ Let $\boldsymbol{\mathit{D}}\in(\mathds{C}^{r\times r})^{C\times C}$ be the block diagonal matrix such that ${\boldsymbol{\mathit{D}}}_{[i,i]}={\boldsymbol{\mathit{M}}}_{[i,i]}$ for $i\in C.$ Expressing $Sc\left(\boldsymbol{\mathit{M}},1\right)$ as $\boldsymbol{\mathit{D}}+({\boldsymbol{\mathit{M}}}_{[C,C]}-\boldsymbol{\mathit{D}}+{\boldsymbol{\mathit{M}}}_{[C,1]}{\boldsymbol{\mathit{M}}}_{}^{-1}{\boldsymbol{\mathit{M}}}_{[1,C]}),$ we have for any $i\in C,$

The next definition describes a special form that we can express any bDD matrix in, which will occasionally be useful.

Every bDD matrix $\boldsymbol{\mathit{M}}\in(\mathds{C}^{r\times r})^{n\times n}$ with $m$ nonzero off-diagonal blocks can be written as $\boldsymbol{\mathit{X}}+\boldsymbol{\mathit{B}}\boldsymbol{\mathit{B}}^{*}$ where $\boldsymbol{\mathit{B}}\in(\mathds{C}^{r\times r})^{n\times 2m}$ is a unitary edge-vertex transfer matrix and $\boldsymbol{\mathit{X}}$ is a block diagonal PSD matrix. This implies that every bDD matrix is PSD. Furthermore, for every block diagonal matrix $\boldsymbol{\mathit{Y}}$ s.t. $\boldsymbol{\mathit{M}}-\boldsymbol{\mathit{Y}}$ is bDD, we have $\boldsymbol{\mathit{X}}\succcurlyeq\boldsymbol{\mathit{Y}}$. This decomposition can be found in $O(m)$ time and $O(\log n)$ depth.

We let $\boldsymbol{\mathit{X}}=\boldsymbol{\mathit{M}}-\sum_{\{i,j\}:i\neq j}\left(\boldsymbol{\mathit{B}}^{(1)}_{\{i,j\}}(\boldsymbol{\mathit{B}}^{(1)}_{\{i,j\}})^{*}+\boldsymbol{\mathit{B}}^{(2)}_{\{i,j\}}(\boldsymbol{\mathit{B}}^{(2)}_{\{i,j\}})^{*}\right),$ which must be block-diagonal. We now show that for all $i\in V,$ the block ${\boldsymbol{\mathit{X}}}_{[i,i]}$ is PSD. We have for all $i\in V,$

where the last inequality holds since $\boldsymbol{\mathit{M}}$ is bDD.

Thus, if we define $\boldsymbol{\mathit{B}}\in(\mathds{C}^{r\times r})^{n\times 2m}$ such that its columns are all the vectors $\boldsymbol{\mathit{B}}^{(1)}_{\{i,j\}},\boldsymbol{\mathit{B}}^{(2)}_{\{i,j\}}$ defined above, we have $\boldsymbol{\mathit{M}}=\boldsymbol{\mathit{X}}+\boldsymbol{\mathit{B}}\boldsymbol{\mathit{B}}^{*},$ and every column of $\boldsymbol{\mathit{B}}$ has exactly 2 nonzero blocks.

To show that for every block diagonal $\boldsymbol{\mathit{Y}}$ s.t. $\boldsymbol{\mathit{M}}-\boldsymbol{\mathit{Y}}$ is bDD, $\boldsymbol{\mathit{X}}\succcurlyeq\boldsymbol{\mathit{Y}}$, first consider applying the decomposition described above to $\boldsymbol{\mathit{M}}-\boldsymbol{\mathit{Y}}$ instead of $\boldsymbol{\mathit{M}}$. Since the construction of $\boldsymbol{\mathit{B}}$ only depends on the off-diagonal blocks, we get $\boldsymbol{\mathit{M}}-\boldsymbol{\mathit{Y}}=\boldsymbol{\Lambda}+\boldsymbol{\mathit{B}}{\boldsymbol{\mathit{B}}^{*}}$, where $\boldsymbol{\Lambda}$ is block diagonal and PSD. So, $\boldsymbol{\mathit{X}}-\boldsymbol{\mathit{Y}}=\boldsymbol{\mathit{M}}-\boldsymbol{\mathit{B}}{\boldsymbol{\mathit{B}}^{*}}-\boldsymbol{\mathit{Y}}=\boldsymbol{\Lambda}\succcurlyeq 0$.

It is immediate that the decomposition can be found in $O(m)$ time and $O(\log n)$ depth.

Appendix C Schur Complement Chains

In this section, we give a proof of Lemma 3.3. We restate the lemma here for convenience. See 3.3

The pseudocode for procedure ApplyChain that uses an $\boldsymbol{\epsilon}$-vertex sparsifier chain to approximately solve a system of equations in $\boldsymbol{\mathit{M}}^{(0)}$ is given in Figure 1.

We begin by observing that the output vector $\boldsymbol{\mathit{x}}^{(0)}$ is a linear transformation of the input vector $\boldsymbol{\mathit{b}}^{(0)}$. Let $\boldsymbol{\mathit{W}}^{(0)}$ be the matrix that realizes this transformation. Similarly, for $1\leq i\leq d$, define $\boldsymbol{\mathit{W}}^{(i)}$ to be the matrix so that

An examination of the algorithm reveals that

We will now prove by backwards induction on $i$ that

The base case of $i=d$ follows from (4). Using the definition of an $\boldsymbol{\epsilon}$-SCC, we know that $0\preceq(\boldsymbol{\mathit{Z}}^{(i)})^{-1}-{\boldsymbol{\mathit{M}}}_{[F_{i},F_{i}]}^{(i-1)}\preceq\epsilon_{i}\cdot Sc\left(\boldsymbol{\mathit{M}}^{(i-1)},C_{i}\right).$ We show in Lemma C.1 that this implies

As $\boldsymbol{\mathit{M}}^{(i)}\approx_{\epsilon_{i}}Sc\left(\boldsymbol{\mathit{M}}^{(i-1)},F_{i}\right)$,

By combining this identity with (5) and our inductive hypothesis, we obtain

The whole algorithm involves a constant number of applications of $\boldsymbol{\mathit{Z}}^{(i)},{\boldsymbol{\mathit{M}}}_{[F_{i},C_{i}]}^{(i-1)},$ and ${\boldsymbol{\mathit{M}}}_{[C_{i},F_{i}]}^{(i-1)}.$ We observe that in order to compute ${\boldsymbol{\mathit{M}}}_{[F_{i},C_{i}]}^{(i-1)}x_{C_{i}},$ using a multiplication procedure for $\boldsymbol{\mathit{M}}^{(i-1)},$ we can pad $x_{C_{i}}$ with zeros, multiply by $\boldsymbol{\mathit{M}}^{(i-1)},$ and read off the answer on the indices in $C_{i}.$ Similarly, we can multiply vectors with ${\boldsymbol{\mathit{M}}}_{[C_{i},F_{i}]}^{(i-1)}.$ This immediately gives the claimed bounds on work and depth. ∎

We now prove the deferred claims from the above proof.

Let $\boldsymbol{\mathit{M}}\in(\mathds{C}^{r\times r})^{n\times n}$ be a bDD matrix, $F\subseteq V$ be a subset of the indices, and $\boldsymbol{\mathit{Z}}\in(\mathds{C}^{r\times r})^{|F|\times|F|}$ be an hermitian operator satisfying $0\preceq\boldsymbol{\mathit{Z}}^{-1}-{\boldsymbol{\mathit{M}}}_{[F,F]}\preceq\epsilon\cdot Sc\left(\boldsymbol{\mathit{M}},C\right).$ Then,

Using Lemma C.2, we know that the assumption on $\boldsymbol{\mathit{Z}}$ is equivalent to

When we use Fact 2.2 to substitute this inequality into the one above, we obtain

Given a bDD matrix $\boldsymbol{\mathit{M}},$ a partition of its indices $(F,C)$ such that ${\boldsymbol{\mathit{M}}}_{[F,F]},{\boldsymbol{\mathit{M}}}_{[C,C]}$ are invertible, and an invertible hermitian operator $\boldsymbol{\mathit{Z}}\in(\mathds{C}^{r\times r})^{|F|\times|F|},$ the following two conditions are equivalent:

$0\preceq(\boldsymbol{\mathit{Z}})^{-1}-{\boldsymbol{\mathit{M}}}_{[F,F]}\preceq\epsilon\cdot Sc\left(\boldsymbol{\mathit{M}},C\right).$

The left inequality in this statement is equivalent to the left inequality in condition 1. Thus, it suffices to prove the right sides are equivalent.

To this end, using ${\boldsymbol{\mathit{M}}}_{[F,C]}={\boldsymbol{\mathit{M}}}_{[C,F]}^{*},$ it suffices to prove that $\forall x\in(\mathds{C}^{r})^{|F|},y\in(\mathds{C}^{r})^{|C|},$

This is equivalent to proving $\forall x\in(\mathds{C}^{r})^{|F|},$

Since ${\boldsymbol{\mathit{M}}}_{[C,C]}\succcurlyeq 0,$ the rhs is a convex function of $y.$ The minimum is achieved at $y=-{\boldsymbol{\mathit{M}}}_{[C,C]}^{-1}{\boldsymbol{\mathit{M}}}_{[F,C]}^{*}x,$ and we obtain the equivalent condition

Since $Sc\left(\boldsymbol{\mathit{M}},C\right)={\boldsymbol{\mathit{M}}}_{[F,F]}-{\boldsymbol{\mathit{M}}}_{[F,C]}{\boldsymbol{\mathit{M}}}_{[C,C]}^{-1}{\boldsymbol{\mathit{M}}}_{[C,F]}={\boldsymbol{\mathit{M}}}_{[F,F]}-{\boldsymbol{\mathit{M}}}_{[F,C]}{\boldsymbol{\mathit{M}}}_{[C,C]}^{-1}{\boldsymbol{\mathit{M}}}_{[F,C]}^{*},$ we obtain our claim. ∎

Appendix D Finding α𝛼\alpha-bDD Subsets

In this section we check that a simple randomized sampling procedure leads to $\alpha$-bDD subsets. Specifically we will prove Lemma 3.5:

Pseudocode for this routine is given in Figure 2.

We first show that the set returned is guaranteed to be $\alpha$-strongly block diagonally dominant.

If bSDDSubset terminates, it returns $F$ such that ${\boldsymbol{\mathit{M}}}_{[F,F]}$ is $\alpha$-bDD.