Approaching optimality for solving SDD systems

Introduction

Fast algorithms for solving linear systems and the related problem of finding a few fundamental eigenvectors is possibly one of the most important problems in algorithm design. It has motivated work on fast matrix multiplication methods, graph separators, and more recently graph sparsifiers. For most applications the matrix is sparse, and thus one would like algorithms whose run time is efficient in terms of the number of non-zero entries of the matrix. Little is known about how to efficiently solve general sparse systems, $Ax=b$. But substantial progress has been made in the case of symmetric and diagonally dominant (SDD) systems, where $A_{ii}\geq\sum_{j\not=i}|A_{ij}|$. In a seminal work, Spielman and Teng showed that SDD systems can be solved in nearly-linear time [ST04, EEST05, ST06].

Recent research, largely motivated by the Spielman and Teng solver (ST-solver), demonstrates the power of SDD solvers as an algorithmic primitive. The ST-solver is the key subroutine of the fastest known algorithms for a multitude of problems that include: (i) Computing the first non-trivial (Fiedler) eigenvector of the graph, or more generally the first few eigenvectors, with well known applications to the sparsest-cut problem [Fie73, ST96, Chu97]; (ii) Generating spectral sparsifiers that also act as cut-preserving sparsifiers [SS08]; (iii) Solving linear systems derived from elliptic finite elements discretizations of a significant class of partial differential equations [BHV04]. (iv) Generalized lossy flow problems [SD08]; (v) Generating random spanning trees [KM09]; and (vi) Several optimization problems in computer vision [KMT09, KMST09b] and graphics [MP08, JMD+07]; A more thorough discussion of applications of the solver can be found in [Spi10, Ten10].

The ST-solver is an iterative algorithm that produces a sequence of approximate solutions converging to the actual solution of the input system $Ax=b$. The performance of iterative methods is commonly measured in terms of the time required to reduce an appropriately defined approximation error by a constant factor. Even including recent improvements on some of its components, the time complexity of the ST-solver is at least $O(m\log^{15}n)$. The large exponent in the logarithm is indicative of the fact that the algorithm is quite complicated and lacks practicality. The design of a faster and simpler solver is a challenging open question.

Preliminaries

In this Section we briefly recall background facts about Laplacians of weighted graphs. For more details, we refer the reader to [RG97] and [BH03]. Throughout the paper, we discuss connected graphs with positive edge weights. We use $n$ and $m$ to denote $|V|$ and $|E|$.

A symmetric matrix $A$ is positive semi-definite if for any vector $x$, $x^{T}Ax\geq 0$. For such semi-definite matrices $A$, we can also define the $A$-norm of a vector $x$ by

Fix an arbitrary numbering of the vertices and edges of a graph $G$. Let $w_{i,j}$ denote the weight of the edge $(i,j)$. The Laplacian $L_{G}$ of $G$ is the matrix defined by: (i) $L_{G}(i,j)=-w_{i,j}$, (ii) $L_{G}(i,i)=\sum_{i\neq j}w_{i,j}$. For any vector $x$, one can check that

It follows that $L_{G}$ is positive semi-definite and $L_{G}$-norm is a valid norm.

We also define a partial order $\preceq$ on symmetric semi-definite matrices, where $A\preceq B$ if $B-A$ is positive semi-definite. This definition is equivalent to $x^{T}Ax\leq x^{T}Bx$ for all $x$. We say that a graph $H$ $\kappa$-approximates a graph $G$ if

By the definition of $\preceq$ from above, this relationship is equivalent to $x^{T}L_{H}x\leq x^{T}L_{G}x\leq\kappa x^{T}L_{H}x$ for all vectors $x$. This implies that the condition number of the pair $(L_{G},L_{H})$ is upper bounded by $\kappa$. The condition number is an algebraically motivated notion; upper bounds on it are used to predict the convergence rate of iterative numerical algorithms.

Prior work on SDD solvers and related graph theoretic problems

Symmetric diagonally dominant systems are linear-time reducible to linear systems whose matrix is the Laplacian of a weighted graph via a construction known as double cover that only doubles the number of non-zero entries in the system [GMZ95, Gre96]. The one-to-one correspondence between graphs and their Laplacians allows us to focus on weighted graphs, and interchangeably use the words graph and Laplacian.

In a ground-breaking approach, Vaidya [Vai91] proposed the use of spectral graph-theoretic properties for the design of provably good graph preconditioners, i.e. graphs that -in some sense- approximate the input graph, but yet are somehow easier to solve. Many authors built upon the ideas of Vaidya, to develop combinatorial preconditioning, an area on the border of numerical linear algebra and spectral graph theory [BGH+05]. The work in the present paper as well as the Spielman and Teng solver is based on this approach. It is worth noting that combinatorial preconditioning is only one of the rich connections between combinatorics and linear algebra [Chu97, RG97].

Vaidya originally proposed the construction of a preconditioner for a given graph, based on a maximum weight spanning tree of the graph and its subsequent augmentation with graph edges. This yielded the first non-trivial results, an $O((dn)^{1.75})$ time algorithm for maximum degree $d$ graphs, and an $O((dn)^{1.2})$ algorithm for maximum degree $d$ planar graphs [Jos97].

Later, Boman and Hendrickson [BH03] made the crucial observation that the notion of stretch (see Section 6 for a definition) is crucial for the construction of a good spanning tree preconditioner; they showed that if the non-tree edges have average stretch $s$ over a spanning tree, the spanning tree is an $O(sm)$-approximation of the graph. Armed with this observation and the low-stretch tree of Alon et al. [AKPW95], Spielman and Teng [ST03] presented a solver running in time $O(m^{1.31})$.

The major new notion introduced by Spielman and Teng [ST04] in their nearly-linear time algorithm was that of a spectral sparsifier, i.e. a graph with a nearly-linear number of edges that $\alpha$-approximates a given graph for a constant $\alpha$. Before the introduction of spectral sparsifiers, Benczúr and Karger [BK96] had presented an $O(m\log^{3}n)$ algorithm for the construction of a cut-preserving sparsifier with $O(n\log n)$ edges. A good spectral sparsifier is a also a good cut-preserving sparsifier, but the opposite is not necessarily true.

The ST-solver [ST04] consists of a number of major algorithmic components. The base routine is a local partitioning algorithm which is the main subroutine of a global nearly-linear time partitioning algorithm. The partitioning algorithm is used as a subroutine in the construction of the spectral sparsifier. Finally, the spectral sparsifier is combined with the ${O}(m\log^{2}n)$ total stretch spanning trees of [EEST05] to produce a $(k,O(k\log^{c}n))$ ultrasparsifier, i.e. a graph $\hat{G}$ with $n-1+(n/k)$ edges which $O(k\log^{c}n)$-approximates the given graph, for some $c>25$. The bottleneck in the complexity of the ST-solver lies in the running time of the ultra-sparsification algorithm and the approximation quality of the ultrasparsifier.

Our contribution

In an effort to design a faster sparsification algorithm, we ask: when and why the much simpler faster cut-preserving sparsifier of [BK96] fails to work as a spectral sparsifier? Perhaps the essential example is that of the cycle and the line graph; while the two graphs have roughly the same cuts, their condition number is $O(n)$. The missing edge has a stretch of $O(n)$ through the rest of the graph, and thus it has high effective resistance; the effective resistance-based algorithm of Spielman and Srivastava would have kept this edge. It is then natural to try to design a sparsification algorithm that avoids precisely to generate a graph whose “missing” edges have a high stretch over the rest of the original graph.

As we explain in Section 7 the incremental sparsifier is all we need to design a solver that runs in the claimed time. Precisely, we prove the following.

Sparsification by Oversampling

In this section we revisit a sampling scheme proposed by Spielman and Srivastava for sparsifying a graph [SS08]. Consider the following general sampling scheme:

Spielman and Srivastava pick $p^{\prime}_{e}=w_{e}R_{e}$ where $R_{e}$ is the effective resistance of $e$ in $G$, if $G$ is viewed as an electrical network with resistances $1/w_{e}$. This choice returns a spectral sparsifier. A key to bounding the number of required samples is the identity $\sum_{e}p_{e}^{\prime}=n-1$. Calculating good approximations to the effective resistances seems to be at least as hard as solving a system, but as we will see in Section 6, it is easier to compute numbers $p_{e}^{\prime}\geq(w_{e}R_{e})$, while still controlling the size of $t=\sum_{e}p_{e}^{\prime}$. The following Theorem considers a sampling scheme based on $p_{e}^{\prime}$’s with this property.

(Oversampling) Let $G=(V,E,w)$ be a graph. Assuming that $p^{\prime}_{e}\geq w_{e}R_{e}$ for each edge $e\in E$, and $\xi\in\Omega(1/n)$, the graph $G^{\prime}=\textsc{Sample}(G,p^{\prime},\xi)$ satisfies

The proof follows closely that Spielman and Srivastava [SS08], with only a minor difference in one calculation. Let us first review some necessary lemmas.

If we assign arbitrary orientations on the edges, then we can define the incidence matrix $\Gamma\in\Re^{m\times n}$ as follows:

If we let $W$ be the diagonal matrix containing edge weights, then $W^{1/2}$ is a real positive diagonal matrix as well since all edge weights are positive. The Laplacian $L$ can be written in terms of $W$ and $\Gamma$ as

Algorithm Sample forms a new graph by multiplying each edge $e$ by a nonnegative number $s_{e}$. If ${\mathbf{S}}$ is the diagonal matrix with $S(e,e)=s_{e}$, the Laplacian of the new graph can be seen to be equal to

Let $L^{+}$ denote the Moore-Penrose of $L$, i.e. the unique matrix sharing with $L$ its null space, and acting as the inverse of $L$ in its range. The key to the proofs of [SS08] is the $m\times m$ matrix

for which the following lemmas are proved.

(Lemma 3i in [SS08]) $\Pi$ is a projection matrix, i.e. $\Pi^{2}=\Pi$.

We also use Lemma 5.4 below, which is Theorem 3.1 from Rudelson and Vershynin [RV07]. The first part of the Lemma was also used as Lemma 5 in [SS08] in a similar way.

Proof (of Theorem 5.1) Following the pseudocode of Sample, let $t=\sum_{e}p_{e}^{\prime}$ and $q=C_{s}t\log t\log(1/\xi)$. It can be seen that

where the $y_{i}$ are drawn from the distribution

For the distribution $y$ we have $E(yy^{T})=\Pi\Pi=\Pi$. Since $\Pi$ is a projection matrix, we have $||\Pi||_{2}\leq 1$. So, the condition imposed by Lemma 5.4 on the distribution holds for $y$. The fact that $\Pi$ is a projection matrix also gives

The last inequality follows from the assumption about the $p_{e}^{\prime}$. Recall now that we have $\log(1/\xi)\leq\log{n}$ by assumption, $t\geq\sum_{e}w_{e}R_{e}$ by construction, and $\sum_{e}w_{e}R_{e}=n-1$ by Lemma 3 in [SS08]. Combining these facts and setting $q=c_{S}t\log t\log(1/\xi)$ for a proper constant $c_{S}$, part 1 of Lemma 5.4 gives

Now substituting $x=\frac{1}{2}$ into part 2 of Lemma 5.4, we get

It follows that with probability at least $1-\xi$ we have

which implies $||\Pi S\Pi-\Pi\Pi||_{2}\leq 1/2$. The theorem then follows by Lemma 5.3. $\blacksquare$

Note. The upper bound for $M$ in inequality 5.1 is in fact the only place where our proof differs from that of [SS08]. In their case the last inequality is replaced by an exact inequality, which is possible because the exact values for $w_{e}R_{e}$ are used. In our case, by using inexact values we get a weaker upper bound which reflects in the density (depending on $m$, not $n$) of the incremental sparsifier. It is however enough for the solver.

Incremental Sparsifier

Consider a spanning tree $T$ of $G=(V,E,w)$. Let $w^{\prime}(e)=1/{w(e)}$. If the unique path connecting the endpoints of $e$ consists of edges $e_{1}\dots e_{k}$, the stretch of $e$ by $T$ is defined to be

Let $R_{e}$ denote the effective resistance of $e$ in $G$ and $RT_{e}$ denote the effective resistance of $e$ in $T$. We have $RT_{e}=\sum_{i=1}^{k}1/w(e_{i})$. Thus $stretch_{T}(e)=w_{e}RT_{e}$. By Rayleigh’s monotonicity law [DS00], we have $RT_{e}\geq R_{e}$, so $stretch_{T}(e)\geq w_{e}R_{e}$. As the numbers $stretch_{T}(e)$ satisfy the condition of Theorem 5.1, we can use them for oversampling. But at the same time we want to control the total stretch, as it will directly affect the total number of samples required in SAMPLE. This leads to taking $T$ to be a low-stretch tree, with the guarantees provided by the following result of Abraham, Bartal, and Neiman [ABN08].

(Corollary 6 in [ABN08]) Given a graph $G=(V,E,w^{\prime})$, LowStretchTree(G) in time $O(m\log n+n\log^{2}n)$, outputs a spanning tree $T$ of $G$ satisfying $\sum_{e\in E}=O(m\log{n}\cdot\log{\log{n}}^{3}).$

Our key idea is to scale up the low-stretch tree by a factor of $\kappa$, incurring a condition number of $\kappa$ but allowing us to sample the non-tree edges aggressively using the upper bounds on their effective resistances given by the tree. The details are given in algorithm IncrementalSparsify.

Given a graph $G$ with $n$ vertices, $m$ edges and any values $\kappa<m$, $\xi\in\Omega(1/n)$, IncrementalSparsify computes a graph $H$ such that:

Proof We first bound the condition number. Since the weight of an edge is increased by at most a factor of $\kappa$, we have $G\preceq G^{\prime}\preceq\kappa G$. Furthermore, the effective resistance along the tree of each non-tree edge decreases by a factor of $\kappa$. Thus IncrementalSparsify sets $p^{\prime}_{e}=1$ if $e\in T$ and $stretch_{T}(e)/\kappa$ otherwise, and invokes Sample to compute a graph $H$ such that with probability at least $1-\xi$, we get

A standard form of Chernoff’s inequality is:

Letting $\delta=2$, and using the assumption $k<m$, we get $Pr(X>3E[X])<(e^{2}/27)^{E[X]}<1/n^{c},$ for any constant $c$. Hence, the probability that IncrementalSparsify succeeds, with respect to both the number of non-tree edges and the condition number, is at least $1-\xi$.

Solving using Incremental Sparsifiers

The solver of Spielman and Teng [ST06] consists of two phases. The preconditioning phase builds a chain of progressively smaller graphs ${\cal C}=\{A_{1},B_{1},A_{2},\ldots,A_{d}\}$ starting with $A_{1}=A$. The process for building ${\cal C}$ alternates between calls to a sparsification routine UltraSparsify which constructs $B_{i}$ from $A_{i}$ and a routine GreedyElimination (following below) which constructs $A_{i+1}$ from $B_{i}$. The preconditioning phase is independent from the $b$-side of the system $L_{A}x=b$.

The solve phase passes $\cal C$, $b$ and a number of iterations $t$ (depending on a desired error $\epsilon$) to the recursive preconditioning algorithm R-P-Chebyshev, described in Section 9. The time complexity of the solve phase depends on $\epsilon$, but more crucially on the quality of $\cal C$, which is a function of the sparsifier quality.

Let $\kappa(n)$ be a monotonically non-decreasing function of $n$. Let ${\cal C}=\{A=A_{1},B_{1},A_{2},\ldots,A_{d}\}$ be a chain of graphs, and denote by $n_{i}$ and $m_{i}$ the numbers of nodes and edges of $A_{i}$ respectively. We say that $\cal C$ is $\kappa(n)$-good for $A$, if:

$A_{i}\preceq B_{i}\preceq\kappa(n_{i})A_{i}$.

$A_{i+1}=\textsc{GreedyElimination}(B_{i})$.

$m_{i}/m_{i+1}\geq c_{r}\sqrt{\kappa(n_{i})}$, for some constant $c_{r}$.

Spielman and Teng analyzed a recursive preconditioned Chebyshev iteration and showed that a $\kappa(n)$-good chain for $A$ can be used to solve a system on $L_{A}$. This is captured by the following Lemma, adapted from Theorem 5.5 in [ST06].

Given a $\kappa(n)$-good chain for $A$, a vector ${x}$ such that $||{x}-L_{A}^{+}b||_{A}<\epsilon||L_{A}^{+}b||_{A}$ can be computed in ${O}(m_{d}^{3}m_{1}\sqrt{\kappa(n_{1})}\log(1/\epsilon))$ expected time.

For our solver, we follow the approach of Spielman and Teng. The main difference is that we replace their routine UltraSparsify with our routine IncrementalSparsify, which is not only faster but also constructs a better chain which translates into a faster solve phase. We are now ready to state our algorithm for building the chain. In what follows we write $v:=O(g(n_{i}))$ to mean ‘$v:=f(n_{i})$ for some explicitly known function $f(n)\in O(g(n))$’.

Proof Assume that $B_{i}$ has $n_{i}-1+m_{i}/k^{\prime}$ edges. A key property of GreedyElimination is that if $G$ is a graph with $n-1+j$ edges, GreedyElimination$(G)$ has at most $2j-2$ vertices and $3j-3$ edges [ST06]. Hence $\textsc{GreedyElimination}(B_{i})$ has at most $3m_{i}/k^{\prime}$ edges. It follows that $m_{i}/m_{i+1}\geq k^{\prime}/3$. Then, in order to satisfy the second requirement, we must have $A_{i}\preceq B_{i}\preceq c^{\prime}k^{\prime 2}A_{i}$, for some sufficiently small constant $c^{\prime}$.

The probability that $B_{i}$ has the above properties is by construction at least $1-p/(2\log n)$ if $n_{i}>\log n$ and $1-p/(2\log\log n)$ otherwise. The probability that the requirements hold for all $i$ is then at least

Combining Lemmas 7.2 and 7.3 proves our main Theorem.

Comments / Extensions

Unraveling the analysis of our bound for the condition number of the incremental sparsifier, it can been that one $\log n$ factor is due to the number of samples required by the Rudelson and Vershynin theorem. The second $\log n$ factor is due to the average stretch of the low-stretch tree.

It is quite possible that the low-stretch construction and perhaps the associated lower bound can be bypassed -at least for some graphs- by a simpler approach similar to that of [KM07]. Consider for example the case of unweighted graphs. With a simple ball-growing procedure we can concede in our incremental sparsifier a $1/\log n$ fraction of the edges, while keeping within clusters of diameters $O(\log^{2}n)$ the rest of the edges. The design of low-stretch trees may be simplified within the small diameter clusters. This diameter-restricted local sparsification is a natural idea to pursue, at least in an actual implementation of the algorithm.

References

Appendix: The Complete Solver

The purpose of this section is to provide a few more algebraic details about the chain of preconditioners, and the recursive preconditioned Chebyshev method which consists the solve phase of the solver. The material is not new and we include it only for completeness. We focus on pseudocode. We refer the reader to [ST06] for a more detailed exposition along with proofs.

Direct methods - Cholesky factorization. If $A$ is a symmetric and positive definite (SPD) matrix, it can be written in the form $A=LL^{T}$, a product known as the Cholesky factorization of $A$. This extends to Laplacians, with some care for the null space. The Cholesky factorization can be computed via a symmetric version of Gaussian elimination. Given the decomposition, solving the systems $Ly=b$ and $L^{T}x=y$ yields the solution to the system $Ax=b$; the key here is that solving with $L$ and $L^{T}$ can be done easily via forward and back substitution. A partial Cholesky factorization with respect to the first $k$ variables of $A$, puts it into the form

where $I_{k}$ denotes the $k\times k$ identity matrix, and $A_{k}$ is known as the Schur complement of $A$ with respect to the elimination of the $k$ first variables. The matrix $A_{k+1}$ is the Schur complement of $A_{k}$ with respect the the elimination of its first variable.

Given a matrix $A$, the graph $G_{A}$ of $A$ is defined by identifying the vertices of $G_{A}$ with the rows and columns of $A$ and letting the edges of $G_{A}$ encode the non-zero structure of $A$ in the obvious way.

It is instructive to take a graph-theoretic look at the partial Cholesky factorization when $k=1$. In this case, the graph $G_{A_{1}}$ contains a clique on the neighbors of the first node in $G_{A}$. In addition, the first column of $L$ is non-zero on the corresponding coordinates. This problem is known as fill. It then becomes obvious that the complexity of computing the Cholesky factorization depends crucially on the ordering of $A$. Roughly speaking, a good ordering has the property that the degrees of the top nodes of $A,A_{1},A_{2},\ldots,A_{k}$ are as small as possible. The best known algorithm for positive definite systems of planar structure runs in time $O(n^{1.5})$ and it is based on the computation of good orderings via nested dissection [Geo73, LRT79, AY].

There are two fairly simple but important facts considering the partial Cholesky factorization of equality 9.2 [ST06]. First, if the top nodes of $A,A_{1},\ldots,A_{k-1}$ have degrees $1$ or $2$, then back-substitution with $L$ requires only $O(n)$ time. Second, if $A$ is a Laplacian, then $A_{k}$ is a Laplacian. Such an ordering and the corresponding Laplacian $A_{k}$ can be found in linear time via GreedyElimination, described in Section 7. The corresponding factor $L$ can also be computed easily.

Iterative methods. Unless the system matrix is very special, direct methods do not yield nearly-linear time algorithms. For example, the nested dissection algorithm is known to be asymptotically optimal for the class of planar SPD systems, within the envelope of direct methods. Iterative methods work around the fill problem by producing a sequence of approximate solutions using only matrix-vector multiplications and simple vector-vector operations. For example Richardson’s iteration generates an approximate solution $x_{i+1}$ from $x_{i}$, by letting

The solver in this paper, as well as the Spielman and Teng solver [ST06], are based on the very well studied Chebyshev iteration [Axe94]. The preconditioned Chebyshev iteration (P-Chebyshev) is the Chebyshev iteration applied to the system $B^{+}Ax=B^{+}b$, where $A,B$ are SPD matrices, and $B$ is known as the preconditioner. The preconditioner $B$ needs not be explicitly known. The iteration requires matrix-vector products with $A$ and $B^{+}$. A product of the form $B^{+1}z$ is equivalent to solving the system $By=c$. Therefore (P-Chebyshev) requires access to only a function $f_{B}(c)$ returning $B^{+1}c$. In addition it requires a lower bound $\lambda_{\min}$ on the minimum eigenvalue of $(A,B)$ and an upper bound $\lambda_{\max}$ on the maximum generalized eigenvalue of $(A,B)$.

A well known fact about the Chebyshev method is that after $O(\sqrt{\lambda_{\max}/\lambda_{\min}}\log 1/\epsilon)$ iterations the return vector $x$ satisfies $\left\|{x}-A^{+}b\right\|_{A}\leq\epsilon\left\|A^{+}b\right\|_{A}$ [Axe94].

Hybrid methods. One of the key ideas in Vaidya’s approach was to combine direct and iterative methods into a hybrid method by exploiting properties of Laplacians. [Vai91]. For the rest of this section we will identify graphs and their Laplacians, using their natural 1-1 correspondence.

Let $A_{1}$ be a Laplacian. The incremental sparsifier $B_{1}$ of $A_{1}$ is a natural choice as preconditioner. With proper input parameters, IncrementalSparsify returns a $B_{1}$ that contains enough degree $1$ and $2$ nodes, so that GreedyElimination can make enough progress reducing $B_{1}$ to a matrix of the form

where $A_{2}$ is the output of algorithm GreedyElimination. Let $I_{j}$ denote the identity of dimension $j$ and

Recall that P-Chebyshev requires the solution of $By=c$, which is given by

The two matrix-vector products with $L_{1}^{-1},L_{1}^{-T}$ can be computed in time $O(n)$ via forward and back substitution. Therefore, we can solve a system in $B$ by solving a linear system in $A_{2}$ and performing $O(n)$ additional work. Naturally, in order to solve systems on $A_{2}$ we can recursively apply preconditioned Chebyshev iterations on it, with a new preconditioner $B_{2}$. This defines a preconditioning chain $\cal C$ that consists of progressively smaller graphs $A=A_{1},B_{1},A_{2},B_{2},\ldots,A_{d}$, along with the corresponding matrices $L_{i},\Pi_{i},Q_{i}$ for $1\leq i\leq d-1$. So, to be more precise than in Section 7, routine BuildChain has the following specifications.

We are now ready to give pseudocode for the recursive preconditioned Chebyshev iteration.

The complete solver. Finally, the pseudocode for the complete solver is as follows.