COMBINATORIAL PRECONDITIONERS SIVAN TOLEDO AND HAIM AVRON TEL-AVIV UNIVERSITY 1. Introduction The Conjugate Gradient (CG) method is an iterative algorithm for solving linear systems of Ax = b, where A is symmetric and positive denite. The convergence of the method A; when its eigenvalues are clustered,√the method converges rapidly. In particular, CG converges to within a xed tolerance in O( κ), where κ = κ(A) = λmax (A)/λmin (A) is the spectral condition number of A. When the spectrum of A is not clustered, a preconditioner can accelerate convergence. The equations depends on the spectrum of Preconditioned Conjugate Gradients (PCG) method applies the CG iteration to the linear −1/2 system (B AB −1/2 )(B 1/2 x) = B −1/2 b using a clever transformation that only requires −1 applications of A and B in every iteration; B also needs to be symmetric positive denite. −1/2 The convergence of PCG is determined by the spectrum of (B AB −1/2 ), which is the same −1 −1 as the spectrum of B A. If a representation of B can be constructed quickly and applied −1 quickly, and if B A has a clustered spectrum, the method is very eective. There are also variants of PCG that require only one of A and B to be positive denite, and variants that allow them to be singular, under some technical conditions on their null spaces. Combinatorial preconditioning is a technique that relies on graph algorithms to construct eective preconditioners. The simplest applications of combinatorial preconditioning target a class of matrices that are isomorphic to weighted undirected graph. The coecient matrix is viewed as its isomorphic graph GB GA . such that the isomorphic matrix A specialized graph algorithm constructs another graph B is a good preconditioner for aims to achieve two goals: the inverse of B It turns out that the spectrum of of properties of the graphs GA GB ; A. The graph algorithm should be easy to apply, and the spectrum of B −1 A can be bounded in terms B −1 A should be clustered. and A in particular, the quality of embeddings of GA in GB (and sometimes vice versa) plays a fundamental role in these spectral bounds. This chapter focuses on explaining the relationship between the spectrum of B −1 A and quantitative properties of embeddings of the two graphs. The last section surveys algorithms that construct combinatorial preconditioners. We omit most proofs from this chapter; some are trivial, and the others appear in the paper cited in the statement of the theorem or lemma. : January, 2010. This research was supported in part by an IBM Faculty Partnership Award and by grant 1045/09 from the Israel Science Foundation (founded by the Israel Academy of Sciences and Humanities). Date 1 COMBINATORIAL PRECONDITIONERS 2 2. Symmetric Diagonally-Dominant Matrices and Graphs We begin by exploring the structure of diagonally-dominant matrices and their relation to graphs. 2.1. Incidence Factorizations of Diagonally-Dominant Matrices. Denition 2.1. 1, 2, . . . n we have A ∈ Rn×n A square matrix Aii ≥ P i6=j is called diagonally-dominant U i= |Aij |. Symmetric diagonally dominant matrices have symmetric factorizations that each column of if for every A = UUT such has at most two nonzeros, and all nonzeros in each column have the same absolute values. We now establish a notation for such columns. Denition 2.2. and the Let 1 ≤ i, j ≤ n, i 6= j . The length-n negative edge vector hi, ji are dened by +1 k = i −1 k = j hi, −jik = 0 otherwise. and positive edge vector denoted +1 k = i +1 k = j . hi, jik = 0 otherwise. The reason for the assignment of signs to edge vectors will become apparent later. A vector hii is the unit vector hiik = hi, −ji vertex +1 k = i 0 otherwise. A symmetric diagonally dominant matrix can always be expressed as a sum of outer products of edge and vertex vectors, and therefore, as a symmetric product of a matrix whose columns are edge and vertex vectors. Lemma 2.3. ([6]) Let A ∈ Rn×n be a diagonally dominant symmetric matrix. We can decompose A as follows A = X |Aij | hi, ji hi, jiT i<j Aij >0 + X |Aij | hi, −ji hi, −jiT i<j Aij <0 n n X X T + |Aij | Aii − hii hii i=1 j=1 j6=i Matrix decompositions of this form play a prominent role in support theory, so we give them a name: COMBINATORIAL PRECONDITIONERS Denition 2.4. of the forms U U T where 3 A matrix whose columns are scaled edge and vertex vectors (that is, vectors c hi, −ji, c hi, ji, and c hii) is called U is an incidence matrix is called an an incidence matrix. A factorization A = incidence factorization. An incidence factorization with no zero columns, with at most one vertex vector for each index most one edge vector for each index pair form c hmin(i, j), − max(i, j)i is called a i, j , i, with at and whose positive edge vectors are all of the canonical incidence factorization. Lemma 2.5. Let A ∈ Rn×n be a diagonally dominant symmetric matrix. Then A has an incidence factorization A = U U T , and a unique canonical incidence factorization. 2.2. Graphs and Their Laplacians Matrices. We now dene the connection between undirected graphs and diagonally-dominant symmetric matrices. Denition 2.6. G = ({1, 2, . . . n}, E, c, d) be a weighted undirected graph on the vertex set {1, 2, . . . , n} with no self loops and no parallel edges, and with weight functions c : E → R \ {0} and d : {1, . . . , n} → R+ ∪ {0}. That is, the edge set consists of unordered pairs of n×n unequal integers (i, j) such that 1 ≤ i, j ≤ n. The Laplacian of G is the matrix A ∈ R such that P d(i) + (i,k)∈E |c(i, k)| i = j Aij = −c(i, j) (i, j) ∈ E 0 otherwise. A vertex i such that d(i) > 0 is called a strictly dominant vertex. If c = 1, the graph is not considered weighted. If c > 0 and is not always 1, the graph is weighted. If some weights Let are negative, the graph is signed. Lemma 2.7. The Laplacians of the graphs dened in Denition 2.6 are symmetric and diagonally dominant. Furthermore, these graphs are isomorphic to symmetric diagonallydominant matrices under this Laplacian mapping. We prefer to work with vertex weights rather than allowing self loops because edge and vertex vectors are algebraically dierent. In algorithms, given an explicit representation of a diagonally-dominant matrix easily compute an explicit representation of an incidence factor U A, we can (including the canonical incidence factor if desired). Sparse matrices are often represented by a data structure that stores a compressed array of nonzero entries for each row or each column of the matrix. Each entry in a row (column) array stores the column index (row index) of the nonzero, and the value of the nonzero. From such a representation of representation of U by columns. We traverse each row of nonzero in the upper (or lower) part of the d(i)'s. A. A we can easily construct A, creating a column of U a sparse for each During the traversal, we can also compute all The conversion works even if only the upper or lower part of A is represented explicitly. A as an implicit representation of U , with each odiagonal nonzero of A representing an edge-vector column of U . If A has no strictly-dominant rows, that is all. If A has strictly dominant rows, we need to compute their weights using a one-pass traversal of A. We can use the explicit representation of COMBINATORIAL PRECONDITIONERS 4 3. Support Theory Support theory is a set of tools that aim to bound the generalized eigenvalues satisfy λ that Ax = λBx from above and below. If B is nonsingular, these eigenvalues are also the B −1 A, but the generalized representation allows us to derive bounds that also eigenvalues of apply to singular matrices. Denition 3.1. λ is a v 6= 0 such that Av = λBv and Bv 6= 0. We say that ∞ is a innite generalized eigenvalue of (A, B) if there exists a vector v 6= 0 such that Bv = 0 but Av 6= 0. Note that ∞ is an eigenvalue of (A, B) if and only if 0 is an eigenvalue of (B, A). The nite and innite eigenvalues of a Let A and nite generalized eigenvalue pencil are B n-by-n be complex matrices. of the matrix pencil (pair) determined eigenvalues (A, B) We say that a scalar if there is a vector (the eigenvector uniquely determines the eigenvalue). If Av = Bv = 0 for a vector v 6= 0, we say that v Av = λBv for any scalar λ. both is an indeterminate eigenvector, because The tools of support theory rely on symmetric factorizations A = UUT and B = V V T; this is why incidence factorizations are useful. In fact, the algebraic tools of support theory are particularly easy to apply when 3.1. 2 U and V are incidence matrices. From Generalized Eigenvalues to Singular Values. T If A = UUT then Λ(A) = T Σ (U ), where Λ(A) is the set of eigenvalues of A and Σ(U ) is the set of singular values of U T , and Σ2 is the set of the squares of the singular values. The following lemma extends this trivial result to generalized eigenvalues. Lemma 3.2. ([2])Let A = U U T and B = V V T with null(B) = null(A) = S. We have and Λ (A, B) = Σ2 V + U Λ (A, B) = Σ−2 U + V . In these expressions, Σ(·) is the set of nonzero singular values of the matrix within the parentheses, Σ` denotes the same singular values to the `th power, and V + denotes the MoorePenrose pseudoinverse of V . The lemma characterizes all the generalized eigenvalues of the pair (A, B), but for large matrices, it is not particularly useful. Even if U and V are highly structured (e.g., they are + + incidence matrices), U and V are usually not structured and are expensive to compute. The next section shows that if we lower our expectations a bit and only try to bound Λfrom above and below, then we do not need the pseudo-inverses. 3.2. The Symmetric Support Lemma. V U provide complete information on the generalized eigenvalues of (A, B). Wopt = V + U , we have singular values of If we denote In the previous section we have seen that the + V Wopt = V V + U = U. COMBINATORIAL PRECONDITIONERS It turns out that any eigenvalues of W VW = U such that 5 provides some information on the generalized (A, B). Lemma 3.3. ([9]) Let A = U U T and let B = V V T , and assume that null(B) ⊆ null(A). Then max {λ | Ax = λBx, Bx 6= 0} = min kW k22 | U = V W . This lemma is fundamental to support theory and preconditioning, because it is often W possible to prove that 3.3. Norm bounds. U =VW such that exists and to give a-priori bounds on its norm. The Symmetric Product Support Lemma bounds generalized eigen- values in terms of the 2-norm of some matrix simple way to construct such a W, W such that U = V W. Even if we have a we still cannot easily derive a corresponding bound on the spectrum from the Symmetric-Support Lemma. The diculty is that there is no simple closed form expression for the W 2-norm of a matrix, since it is not related to the entries of in a simple way. It is equivalent to the largest singular value, but this must usually be computed numerically. Fortunately, there are simple (and also some not-so-simple) functions of the elements of the matrix that yield useful bounds on its 2-norm. The following bounds are standard and are well known and widely used (see [5, Fact 9.8.10.ix] and [5, Fact 9.8.15]). Lemma 3.4. The two norm of W ∈ Ck×m is bounded by (3.1) kW k22 ≤ kW k2F = k X m X Wij2 , i=1 i=1 (3.2) kW k22 ≤ kW k1 kW k∞ = m max j=1 k X ! |Wij | k max i=1 i=1 m X ! |Wij | . j=1 The next two bounds are standard and well known; they follow directly from kW k22 and from the fact that kSk1 = kSk∞ for a symmetric S . kW W T k2 = Lemma 3.5. The two norm of W ∈ Ck×m is bounded by (3.3) (3.4) kW k22 ≤ kW W T k1 = kW W T k∞ , kW k22 ≤ kW T W k1 = kW T W k∞ . The following bounds are more specialized. They all exploit the sparsity of bounds that are usually tighter than the bounds given so far. Lemma 3.6. ([11]) The two norm of W ∈ Ck×m is bounded by W to obtains COMBINATORIAL PRECONDITIONERS kW k22 (3.5) ≤ max j kW k22 (3.6) ≤ max i X kWi, : k22 X kW : ,j k22 m X X = max j i:Wi,j 6=0 2 , Wi,c i:Wi,j 6=0 c=1 k X X = max i j:Wi,j 6=0 6 2 Wr,j . j:Wi,j 6=0 r=1 The bounds in this lemma are a renement of the bound kW k22 ≤ kW k2F . nious norm, which bounds the two norm, sums the squares of all The Frobe- the elements of W. The bounds (3.5) and (3.6) sum only the squares in some of the rows or some of the columns, unless the matrix has a row or a columns with no zeros. 2 There are similar renements of the bound kW k2 ≤ kW k1 kW k∞ . Lemma 3.7. ([11]) The two norm of W ∈ Ck×m is bounded by kW k22 ≤ max (3.7) j kW k22 (3.8) 3.4. Support numbers. ≤ max i X |Wi,j | · m X i:Wi,j 6=0 c=1 X k X |Wi,j | · j:Wi,j 6=0 ! |Wi,c | , ! |Wr,j | . r=1 Support numbers generalize the notion of the maximal eigenvalue of a matrix pencil. Denition 3.8. A matrix B dominates a matrix A if for any vector x we have xT (B −A)x ≥ 0. We denote domination by Denition 3.9. The B A. support number for a matrix pencil σ(A, B) = min {t|τ B A, (A, B) for all is τ ≥ t} . is symmetric positive denite and A is symmetric, then the support number is always T T T T nite, because x Bx/x x is bounded from below by min Λ(B) > 0 and x Ax/x x is bounded If B from above by max Λ(A), which is nite. In other cases, there may not be any the formula; in such cases, we say that Example 3.10. τ >0 we have satisfying σ(A, B) = ∞. x ∈ null(B) and that A is positive denite. x (τ B − A)x = −xT Ax < 0. Therefore, σ(A, B) = ∞. Suppose that t Then for any T Example 3.11. T is not positive semidenite, then there is some x for which x Bx < 0. T This implies that for any A and for any large enough τ , x (τ B − A)x < 0. Therefore, If B σ(A, B) = ∞. The next result, like the Symmetric Support Lemma, bounds generalized eigenvalues. COMBINATORIAL PRECONDITIONERS 7 Theorem 3.12. ([9]) Let A and B be symmetric, let B also be positive semidenite. If null(B) ⊆ null(A), then σ(A, B) = max {λ|Ax = λBx, Bx 6= 0} . A primary motivation for support numbers is to bound (spectral) condition numbers. For κ(A) = λmax (A)/λmin (A) holds. Let κ(A, B) denote the −1 pencil (A, B), that is, κ(B A) when B is non-singular. symmetric matrices, the relation condition number of the matrix Theorem 3.13. When A and B are symmetric positive denite, then κ(A, B) = σ(A, B)σ(B, A). A common strategy in support theory is to bound condition numbers by bounding both σ(A, B) and σ(B, A). Typically, one direction is easy and the other is hard. Applications usually do not solve singular systems. Nonetheless, it is often convenient to analyze preconditioners in the singular context. For example, nite element systems are often singular until boundary conditions are imposed, so we can build A B from the singular part of and then impose the same boundary-constraints on both matrices. 3.5. Splitting. Support numbers are convenient for algebraic manipulation. One of their most powerful properties is that they allow us to split complicated matrices into simpler pieces (matrices) and analyze these separately. B = B1 + B2 + · · · + Bq . A = A1 + A2 + · · · Aq , and similarly, pairs (Ai , Bi ) and consider the support Let We can then match up number for each such pencil separately. Lemma 3.14. ([17, Lemma 4.7]) Let A = A1 + A2 + · · · Aq ,and similarly, B = B1 + B2 + · · · + Bq , where all Ai and Bi are symmetric and positive semidenite. Then σ(A, B) ≤ max σ(Ai , Bi ) i Proof. Let σ = σ(A, B), let σi = σ(Ai , Bi ), and let σmax = maxi σi . Then for any x ! xT (σmax B − A)x = xT σmax X Bi − X i = X ≥ X Ai x i xT (σmax Bi − Ai ) x i xT (σi Bi − Ai ) x i ≥ 0. Therefore, σ ≤ σmax . The splitting lemma is quite general, and can be used in many ways. In practice we want to break both A and B into simpler matrices that we know how to analyze. The term simpler can mean sparser, or lower rank, and so on. In order to get a good upper bound on the support number, the splitting must be chosen carefully. Poor splittings give poor bounds. Here is an example. COMBINATORIAL PRECONDITIONERS 8 3 −2 2 −2 1 0 A = = + = A1 + A2 , and B = −2 2 −2 2 0 0 2 −1 1 −1 1 0 = + = B1 +B2 . Then the Splitting Lemma says σ(A, B) ≤ −1 2 −1 1 0 1 max {σ(A1 , B1 ), σ(A2 , B2 )}. It is easy to verify that σ(A1 , B1 ) = 2 and that σ(A2 , B2 ) = 1; hence σ(A, B) ≤ 2. Note that B1 can not support A2 , so correct pairing of the terms in A and B is essential. The exact support number is σ(A, B) = λmax (A, B) = 1.557. Example 3.15. Let 4. Embeddings and Combinatorial Support Bounds A To bound σ(A, B) using the Symmetric Support Lemma, we need to factor A and B into = U U T and B = V V T , and we need to nd a W such that U = V W . We have seen that if A and B are diagonally dominant, then there is an almost trivial way to factor A B such that U and V are about as sparse as A and B . But how do we nd a W such U = V W ? In this chapter, we show that when A and B are weighted (but not signed) Laplacians, we can construct such W using an embedding of the edges of GA into paths in GB . Furthermore, when W is constructed from an embedding, the bounds on kW k2 can be and that interpreted as combinatorial bounds on the quality of the embedding. 4.1. W Dening W using Path Embeddings. such that We start with the construction of a matrix U = V W. Lemma 4.1. Let (i1 , i2 , . . . , i` ) be a sequence of integers between 1 and n, such that ij 6= ij+1 for j = 1, . . . ` − 1.Then hi1 , −i` i = `−1 X hij , −ij+1 i , j=1 where all the edge vectors are length n. W . Suppose c of U by To see why this lemma is important, consider the role of a column of the columns of U and V are all positive edge vectors. Denote column that U : ,c = hmin(i1 , i` ), − max(i1 , i` )i = (−1)i1 >i` hi1 , −i` i , (−1)i1 >i` evaluates to −1 if i1 > i` and to 1 otherwise. This column corresponds to the edge (i1 , i` ) in GU U T . Now let (i1 , i2 , . . . , i` ) be a simple path in GV V T (a simple path is a sequence of vertices (i1 , i2 , . . . , i` ) such that (ij , ij+1 ) is an edge in the graph for 1 ≤ j < ` and such that any vertex appears at most once on the path). If U = V W , then where the U : ,c = V W : ,c = k X V : ,r Wr,c . r=1 r1 , r2 , . . . , r`−1 be the columns of V that corresponds to the edges of the path (i1 , i2 , . . . , i` ), order. That is, V : ,r1 = hmin(i1 , i2 ), − max(i1 , i2 )i, V : ,r2 = hmin(i2 , i3 ), − max(i2 , i3 )i, and Let in COMBINATORIAL PRECONDITIONERS 9 so on. By the lemma, U : ,c = (−1)i1 >i` hi1 , −i` i i1 >i` = (−1) `−1 X hij , −ij+1 i j=1 i1 >i` = (−1) `−1 X (−1)ij >ij+1 V : ,rj . j=1 It follows that if we dene Wr,c = W : ,c (−1)i1 >i` (−1)ij >ij+1 r = rj for 0 otherwise, U : ,c = V W : ,c . U = V W. then we have satises to be some 1≤j<` We can construct all the columns of W in this way, so that W A path of edge vectors that ends in a vertex vector supports the vertex vector associated with the rst vertex of the path. Lemma 4.2. Let (i1 , i2 , . . . , i` ) be a sequence of integers between 1 and n, such that ij 6= ij+1 for j = 1, . . . ` − 1.Then hi1 i = hi` i + `−1 X hij , −ij+1 i , j=1 where all the edge and vertex vectors are length n. The following theorem generalizes these ideas to scaled positive edge vectors and to scaled vertex vectors. The theorem also states how to construct all the columns of W. The theorem summarizes results in [9, 4]. Theorem 4.3. Let A and B be weighted (unsigned) Laplacians and let U and V be their canonical incidence factors. Let π be a path embedding of the edges and strictly-dominant vertices of GA into GB , such that for an edge (i1 , i` ) in GA , i1 < i` , we have π(i1 , i` ) = (i1 , i2 , . . . , i` ) for some simple path (i1 , i2 , . . . , i` ) in GB , and such that for a strictly-dominant i1 in GA , π(i1 ) = (i1 , i2 , . . . , i` ) for some simple path (i1 , i2 , . . . , i` ) in GB that ends in a strictly-dominant vertex i` in GB . Denote by cV (ij , ij+1 ) the index of the column of V that is a scaling of hij , −ij+1 i. That is, V : ,cV (ij ,ij+1 ) q = −Bij ,ij+1 hmin(ij , ij+1 ), − max(ij , ij+1 )i . COMBINATORIAL PRECONDITIONERS 10 Similarly, denote by cV (ij ) the index of the column of V that is a scaling of hij i, V : ,cV (ij ) v u n X u Bi ,i hij i , u = tBij ,ij − k j ik =1 ik 6=ij and similarly for U . We dene a matrix W as follows. For a column index cU (i1 , i` ) with i1 < i` we dene Wr,cU (i1 ,i` ) p (−1)ij >ij+1 Ai1 ,i` /Bij ,ij+1 if r = cV (ij , ij+1 ) for = some edge (ij , ij+1 ) in π(i1 , i` ) 0 otherwise. For a column index cU (i1 ), we dene Wr,cU (i1 ) = r P Ai1 ,i1 − j6=i |Ai1 ,j | P 1 Bi` ,i` − j6=i |Bi` ,j | r` P (−1)ij >ij+1 if r = cV (i` ) |Ai1 ,k | if r = cV (ij , ij+1 ) for |Bij ,ij+1 | some edge (ij , ij+1 ) in π(i1 ) Ai1 ,i1 − k6=i1 otherwise. 0 Then U = V W . Proof. For scaled edge-vector columns in V W : ,cU (i1 ,i` ) = X U we have V : ,r Wr,cU (i1 ,i` ) r X = V : ,r Wr,cU (i1 ,i` ) r=cV (ij ,ij+1 ) for some edge (ij ,ij+1 ) in π(i1 ,i` ) `−1 q X Bi ,i hmin(ij , ij+1 ), − max(ij , ij+1 )i (−1)ij >ij+1 = j j+1 j=1 = `−1 q X |Ai1 ,i` | hij , −ij+1 i j=1 = U : ,cU (i1 ,i` ) . s Ai1 ,i` Bij ,ij+1 COMBINATORIAL PRECONDITIONERS For scaled vertex-vector columns in V W : ,cU (i1 ) = X U 11 we have V : ,r Wr,cU (i1 ) r = V : ,cV (i` ) WcV (i` ),cU (i1 ,i` ) + X V : ,r Wr,cU (i1 ) r=cV (ij ,ij+1 ) for some edge (ij ,ij+1 ) in π(i1 ) s P s X Ai1 ,i1 − j6=i1 |Ai1 ,j | P = Bi` ,i` − |Bik ,i` | hi` i B i` ,i` − j6=i` |Bi` ,j | j6=i ` + `−1 q X Bi j ,ij+1 hmin(ij , ij+1 ), − max(ij , ij+1 )i j=1 s P A − i ,i 1 1 j6=i1 |Ai1 ,j | ·(−1)ij >ij+1 Bi ,i j j+1 s s `−1 X X X = Ai1 ,i1 − |Ai1 ,j | hi` i + Ai1 ,i1 − |Ai1 ,j | hij , −ij+1 i j6=i1 j6=i1 j=1 s X Ai1 ,i1 − = |Ai1 ,j | hi1 i j6=i1 = U : ,cU (i1 ) . The generalization of this theorem to signed Laplacians is more complex, because a path from i1 to i` supports an edge (i1 , i` ) only if the parity of positive edges in the path and in the edge is the same. In addition, a cycle with an odd number of positive edges spans all the vertex vectors of the path. For details, see [6]. Theorem 4.3 plays a fundamental role in many applications of support theory. embedding π that can be used to construct W exists if and only if the graphs of A A path and B have related in a specic way, which the next lemma species. Lemma 4.4. Let A = U U T and B = V V T be weighted (but not signed) Laplacians with arbitrary symmetric-product factorizations. The following conditions are necessary for the equation U = V W to hold for some matrix W (by Theorem 4.3, these conditions are also sucient). (1) For each edge (i, j) in GA , either i and j are in the same connected component in GB , or the two components of GB that contain i and j both include a strictly-dominant vertex. (2) For each strictly-dominant vertex i in GA , the component of GB that contains i includes a strictly-dominant vertex. COMBINATORIAL PRECONDITIONERS Proof. is a 12 Suppose for contradiction that one of the conditions is not satised, but that there W that satises U = V W. Without loss of generality, we assume that the vertices are ordered such that vertices that belong to a connected component in GB are consecutive. Under that assumption, V1 V2 V = .. , . Vk and B1 B2 B= .. V1 V1T = . V2 V2T .. . Vk VkT Bk The blocks of V . are possibly rectangular, whereas the nonzero blocks of B are all diagonal and square. We now prove the necessity of the rst condition. i and j belong to dierent connected components of Suppose for some edge GB (i, j) in GA , (without loss of generality, to the rst two components), and that one of the components (w.l.o.g. the rst) does not have a strictly-dominant vertex. Because this component does not have a strictly-dominant vertex, T T 1 = ~0, so V1 must be rank decient. the row sums in V1 V1 are exactly zero. Therefore, V1 V1 ~ Since (i, j) is in GA , the vector hi, −ji is in the column space of the canonical incidence T factor of A, and therefore in the column space of any U such that A = U U . If U = V W , then the vector hi, −ji V , so for some x x1 V1 x1 x2 V2 x2 . = . . .. .. must also be in the column space of V1 V2 hi, −ji = V x = .. . Vk Therefore, V1 x1 is a vertex vector. xk By Theorem 4.3, if V1 Vk xk spans a vertex vector, it spans all the vertex vectors associated with the vertices of the connected component. This implies that V1 is full rank, a contradiction. The necessity of the second condition follows from a similar argument. Suppose that vertex i is strictly dominant in GA and that it belongs to a connected component in the rst) that does not have a vertex that is strictly dominant in some GB . GB (w.l.o.g. This implies that for y V1 y1 V1 y1 y2 V2 y2 . = . . .. .. V2 hii = V y = .. . Vk Again V1 y1 is a vertex vector, so V1 V1T has zero row sums. V1 yk Vk yk must be full rank, but it cannot be full rank because COMBINATORIAL PRECONDITIONERS Not every W 13 U = V W corresponds to a path embedding, even if U and V are A and B .In particular, a column of W can correspond to a multiple paths. Also, even if W does correspond to a path embedding, such that the canonical incidence factors of linear combination of the paths are not necessarily simple. A linear combination of scaled positive edge vectors that correspond to a simple cycle can be identically zero, so the coecients of such linear combinations can be added to W V W . However, it seems that of W , so cycles are unlikely to without aecting the product adding cycles to a path embedding cannot reduce the 2-norm improve support bounds. 4.2. Combinatorial Support Bounds. T A into A = U U 2-norm of W from Lemma, we factor and bound the , B To bound σ(A, B) using the Symmetric Support T into B = V V , nd a matrix W such that U = V W , above. We have seen how to factor weighted Laplacians) and how to construct an appropriate GA . W A and B (if they are from an embedding of We now show how to use combinatorial metrics of the path embeddings to bound GB in kW k2 . 2-norm directly is hard, because the 2-norm is not related in a simple way to W . But the 2-norm can be bounded using simpler norms, such as the Frobenius norm, the innity norm, and the 1-norm (see section 3.3). These simpler norms have natural and useful combinatorial interpretations when W represents a path embedding. To keep the notation and the denition simple, we now assume that A and B are weighted Bounding the the entries of Laplacians with zero row sums. We will show later how to deal with positive row sums. We also assume that W corresponds to a path embedding π. The following denitions provide a combinatorial interpretation of these bounds. Denition 4.5. GB The weighted dilation of an edge of GA in an path embedding is s dilationπ (i1 , i2 ) X = (j1, ,j2 ) (j1, ,j2 )∈π(i1 ,i2 ) The weighted congestion of an edge of GB Ai1 ,i2 . Bj1 ,j2 is s congestionπ (j1 , j2 ) X = (i1, ,i2 ) (j1, ,j2 )∈π(i1 ,i2 ) The weighted stretch of an edge of GA is stretchπ (i1 , i2 ) X = (j1, ,j2 ) (j1, ,j2 )∈π(i1 ,i2 ) The weighted crowding of an edge in Ai1 ,i2 . Bj1 ,j2 GB Ai1 ,i2 . Bj1 ,j2 is crowdingπ (j1 , j2 ) = X (i1, ,i2 ) (j1, ,j2 )∈π(i1 ,i2 ) Ai1 ,i2 . Bj1 ,j2 π of GA into COMBINATORIAL PRECONDITIONERS 14 Note that stretch is a summation of the squares of the quantities that constitute dilation, and similarly for crowding and congestion. Papers in the support-preconditioning literature are not consistent about these terms, so check the denitions carefully when you consult a paper that deals with congestion, dilation, and similar terms. Lemma 4.6. ([9, 28]) Let A and B be weighted Laplacians with zero row sums, and let π be a path embedding of GA into GB . Then stretchπ (i1 , i2 ) X σ(A, B) ≤ (i1 ,i2 )∈GA crowdingπ (j1 , j2 ) X σ(A, B) ≤ (j1 ,j2 )∈GB max dilationπ (i1 , i2 ) σ(A, B) ≤ (i1 ,i2 )∈GA · max (j1 ,j2 )∈GB congestionπ (j1 , j2 ) . We now describe one way to deal with matrices with some positive row sums. The matrix is a preconditioner. In many applications, the matrix B B is not given, but rather constructed. A is to deneπ(i1 ) = (i1 ). That is, vertex A are mapped into the same vertex vectors in the incidence factor of B . In other words, we construct B to have exactly the same row sums as A. With such a construction, the rows of W that correspond to vertex vectors in U are One simple way to deal with positive row sums in vectors in the canonical incidence factor of columns of the identity. Lemma 4.7. Let A and B be weighted Laplacians with the same row sums, let π be a path embedding of GA into GB , and let ` be the number of rows with positive row sums in A and B . Then stretchπ (i1 , i2 ) X σ(A, B) ≤ ` + (i1 ,i2 )∈GA σ(A, B) ≤ max 1, max dilationπ (i1 , i2 ) (i1 ,i2 )∈GA · max 1, max congestionπ (j1 , j2 ) . (j1 ,j2 )∈GB Proof. Under the hypothesis of the lemma, the rows and columns of a block matrix W = where WZ WZ 0 0 I`×` The sparse bounds on the W can be permuted into , represents the path embedding of the edges of follow from the structure of W GA into paths in GB . The bounds and from the proof of the previous lemma. 2-norm of a matrix lead to tighter combinatorial bounds. COMBINATORIAL PRECONDITIONERS 15 Lemma 4.8. Let A and B be weighted Laplacians with zero row sums, and let π be a path embedding of GA into GB . Then σ(A, B) ≤ (j1 ,j2 )∈GB σ(A, B) ≤ stretchπ (i1 , i2 ) , X max (i1 ,i2 )∈GA (j1 ,j2 )∈π(i1 ,i2 ) crowdingπ (j1 , j2 ) . X max (i1 ,i2 )∈GA (j1 ,j2 )∈GA (j1 ,j2 )∈π(i1 ,i2 ) We can derive similar bounds for the other sparse 4.3. Subset Preconditioners. σ(A, B) and σ(B, A). not in V, κ(A, B), we need a bound on both σ(B, A) is trivial. Many support of GA , with the same weights. That is, V is U . If we denote by V̄ the set of columns of To obtain a bound on GB to be a subgraph constructed to have a subset of the columns in that are bounds. But in one common case, bounding preconditioners construct U 2-norm we have B = VVT A = UUT = V V T + V̄ V̄ T = B + V̄ V̄ T . xT Ax ≥ xT Bx This immediately implies 4.4. Combinatorial p Trace Bounds. for any x, so λmin (A, B) ≥ 1. The Preconditioned Conjugate Gradients (PCG) al- Θ( κ(A, B) iterations only when the generalized eigenvalues are distributed poorly between λmin (A, B) and λmax (A, B). For example, if there are only two distinct eigenvalues, PCG converges in two iterations. It turns out that a bound on trace(A, B) = P λi (A, B) also yields p a bound on the number of iterations, and in some cases this bound is sharper than the O( κ(A, B) bound. gorithm requires Lemma 4.9. ([29]) The Preconditioned Conjugate Algorithm converges to within a xed tolerance in s O 3 trace(A, B) λmin (A, B) ! iterations. A variation of the Symmetric Support Lemma bounds the trace, and this leads to a combinatorial support bound. Lemma 4.10. Let A = U U T ∈ Rn×n and let B = V V T ∈ Rn×n , and assume that null(B) = null(A). Then trace(A, B) = min kW k2F | U = V W . COMBINATORIAL PRECONDITIONERS Proof. Let W be a matrix such that U = V W. For every 16 x∈ / null(B) we can write xT Ax xT V W W T V T x = xT Bx xT V V T x yT W W T y = , yT y T n y =T V x. Let S ⊆ R be a subspace orthogonal to null(B) of dimesion k , and dene T to null(B) we have dim(TS ) = k . The sets S T= V T x | x ∈ S .Because TS is orthogonal T x Ax/x Bx | x ∈ S and y W W y | y ∈ TS are identical so their minimums are equal as well. The group of subspaces {TS | dim(S) = k, S ⊥ null(B)} is a subset of all subspaces of dimension k , thererfore where xT U U T x = max min T T dim(S) = k x ∈ S x V V x S ⊥ null(B) x = 6 0 yT W W T y yT y max dim(S) = k S ⊥ null(B) min y ∈ TS y 6= 0 max yT W W T y . yT y ≤ dim(T )=k min y∈T y 6= 0 According to the Courant-Fischer Minimax Theorem, λn−k+1 (W W T ) = max dim(T )=k min y∈T y 6= 0 yT W W T y yT y and by the generalization of Courant-Fischer in [3], λn−k+1 (U U T , V V T ) = xT U U T x . max min T T dim(S) = k x∈S x V V x V ⊥ null(B) k = 1, . . . , rank(B) we have λn−k+1 (A, B) ≤ λn−k+1 (W W T ) so trace(A, B) ≤ 2 2 trace(W W ) = kW kF . This shows that trace(A, B) ≤ min kW kF | U = V W . + According to Lemma 3.2 the minimum is attainable at W = V U . Therefore, for T To the best of our knowledge, Lemma 4.10 is new. It was inspired by a specialized bound on the trace from [29]. The next theorem generalizes the result form [29]. The proof is trivial given Denition 4.5 and Lemma 4.10 COMBINATORIAL PRECONDITIONERS 17 Theorem 4.11. Let A and B be weighted Laplacians with the same row sums, let π be a path embedding of GA into GB . Then X trace(A, B) ≤ stretchπ (i1 , i2 ) , (i1 ,i2 )∈GA X trace(A, B) ≤ crowdingπ (j1 , j2 ) . (j1 ,j2 )∈GB 5. Combinatorial Preconditioners The earliest graph algorithm to construct a preconditioner was proposed by Vaidya [31] (see also [4]). He proposed to use a so-called augmented maximum spanning tree of a weighted Laplacian as a preconditioner. This is a subset preconditioner that drops some of the edges −1 in GA while maintaining the weights of the remaining edges. When B is applied using 7/4 a sparse Cholesky factorization, this construction leads to a total solution time of O(n ) 6/5 for Laplacians with a bounded degree and O(n ) when the graph is planar. For regular unweighted meshes in 2 and 3 dimensions, special constructions are even more eective [20]. −1 Vaidya also proposed to use recursion to apply B , but without a rigorous analysis; in this scheme, sparse Gaussian elimination steps are performed on B as long as the reduced matrix has a row with only two nonzeros. At that point, the reduced system is solved by constructing a graph preconditioner, and so on. Vaidya's preconditioners are quite eective in practice [12]. A generalization to complex matrices proved eective in handling a certain kind of ill conditioning [19]. Vaidya's bounds used a congestion-dilation product. The research on subset graph preconditioners continued with an observation that the sum of the stretch can also yield a spectral bound, and that so-called low-stretch trees would give better worst-case bounds than the maximum-spanning trees that Vaidya used [8]. Constructing low-stretch trees is more complicated than constructing maximum spanning trees. When the utility of low-stretch trees was discovered, one algorithm for constructing them was known [1]; better algorithms were discovered later [15]. Low-stretch trees can also be augmented, and this leads to precondition4/3 ers that can solve any Laplacian system with m nonzeros in O(m ), up to some additional polylogarithmic factors [26, 29]. By employing recursion, the theoretical running time can be reduced to close to linear in m [27]. An experimental comparison of simplied versions of these sophisticated algorithms to Vaidya's algorithm did not yield a conclusive result [30]. Heuristic subset graph preconditioners have also been proposed [16]. Gremban and Miller proposed a class of combinatorial preconditioners called preconditioners [18, 17]. support-tree Their algorithms construct a weighted tree whose leaves are the Therefore, the graph GT of the preconditioner T has additional vertices. −1 They show that applying T to extensions of residuals is equivalent to using a preconditioner vertices of B GA . that is the Schur complement of condition number κ(A, B) T with respect to the original vertices. Bounding the is more dicult than bounding the condition number of subset preconditioners, because the Schur complement is dense. strategy have been proposed later [23, 21, 22]. More eective versions of this COMBINATORIAL PRECONDITIONERS 18 Eorts to generalize these constructions to matrices that are not weighted Laplacians followed several paths. Gremban showed how to transform a linear system whose coecient matrix is a signed Laplacian to a linear system of twice the size whose matrix is a weighted 2-by-2 block matrix with diagonal blocks with the same matrix A and with identity o-diagonal blocks. A dierent Laplacian. The coecient matrix is a sparsity pattern as the original approach is to extend Vaidya's construction to signed graphs [6]. The class of symmetric T matrices with a symmetric factorization A = U U where columns of U have at most 2 nonzeros contains not only signed graphs, but also gain graphs, which are not diagonally dominant [7]; it turns out that these matrices can be scaled to diagonal dominance, which allows graph preconditioners to be applied to them [14]. The matrices that arise in nite-element discretization of elliptic partial dierential equations (PDEs) are positive semi-denite, but in general they are not diagonally dominant. However, when the PDE is scalar (e.g., describes a problem in electrostatics), the matrices can sometimes be approximated by diagonally dominant matrices. In this scheme, the coecient matrix GD A is rst approximated by a diagonally-dominant matrix is used to construct the graph GB of the preconditioner B. D, and then For large matrices of this class, the rst step is expensive, but because nite-element matrices have a natural representation as a sum of very sparse matrices, the diagonally-dominant approximation can be constructed for each term in the sum separately. There are at least three ways to construct these approximations: during the nite-element discretization process [10], algebraically [2], and geometrically [32]. A slightly modied construction that can accommodate terms that do not have a close diagonally-dominant approximation works well in practice [2]. Another approach for constructing combinatorial preconditioners to nite element problems is to rely on a graph that describes the relations between neighboring elements. This graph is the dual of the nite-element mesh; elements in the mesh are the vertices of the graph. Once the graph is constructed, it can be sparsied much like subset preconditioners. This approach, which is applicable to vector problems like linear elasticity, was proposed in [24]; this paper also showed how to construct the dual graph algebraically and how to construct the nite-element problem that corresponds to the sparsied dual graph. The rst eective preconditioner of this class was proposed in [13]. It is not yet known how to weigh the edges of the dual graph eectively, which limits the applicability of this method. However, in applications where there is no need to weigh the edges, the method is eective [25]. References [1] Noga Alon, Richard M. Karp, David Peleg, and Douglas West. A graph-theoretic game and its application to the k -server problem. SIAM Journal on Computing, 24:78100, 1995. [2] Haim Avron, Doron Chen, Gil Shklarski, and Sivan Toledo. Combinatorial preconditioners for scalar elliptic nite-element problems. SIAM Journal on Matrix Analysis and Applications, 31(2):694720, 2009. [3] Haim Avron, Esmond Ng, and Sivan Toledo. Using perturbed QR factorizations to solve linear leastsquares problems. SIAM Journal on Matrix Analysis and Applications, 31(2):674693, 2009. [4] Marshall Bern, John R. Gilbert, Bruce Hendrickson, Nhat Nguyen, and Sivan Toledo. Support-graph preconditioners. SIAM Journal on Matrix Analysis and Applications, 27:930951, 2006. COMBINATORIAL PRECONDITIONERS 19 [5] Denis S. Bernstein. Matrix Mathematics: Theory, Facts, and Formulas with Applications to Linear Systems Theory. Princeton University Press, 2005. [6] Erik G. Boman, Doron Chen, Bruce Hendrickson, and Sivan Toledo. Maximum-weight-basis preconditioners. Numerical Linear Algebra with Applications, 11:695721, 2004. [7] Erik G. Boman, Doron Chen, Ojas Parekh, and Sivan Toledo. On the factor-width and symmetric H-matrices. Numerical Linear Algebra with Applications, 405:239248, 2005. [8] Erik G. Boman and Bruce Hendrickson. On spanning tree preconditioners. Unpublished manuscript, Sandia National Laboratories, 2001. [9] Erik G. Boman and Bruce Hendrickson. Support theory for preconditioning. SIAM J. Matrix Anal. Appl., 25(3):694717, 2003. [10] Erik G. Boman, Bruce Hendrickson, and Stephen Vavasis. Solving elliptic nite element systems in near-linear time with support preconditioners. SIAM Journal on Numerical Analysis, 46(6):32643284, 2008. [11] Doron Chen, John R. Gilbert, and Sivan Toledo. Obtaining bounds on the two norm of a matrix from the splitting lemma. Electronic Transactions on Numerical Analysis, 21:2846, 2005. [12] Doron Chen and Sivan Toledo. Vaidya's preconditioners: Implementation and experimental study. Electronic Transactions on Numerical Analysis, 16:3049, 2003. [13] Samuel I. Daitch and Daniel A. Spielman. Support-graph preconditioners for 2-dimensional trusses. Mar 2007. [14] Samuel I. Daitch and Daniel A. Spielman. Faster approximate lossy generalized ow via interior point algorithms. In STOC '08: Proceedings of the 40th annual ACM Symposium on Theory of Computing, pages 451460, New York, NY, USA, 2008. ACM. [15] Michael Elkin, Yuval Emek, Daniel A. Spielman, and Shang-Hua Teng. Lower-stretch spanning trees. In Proceedings of the 37th annual ACM symposium on Theory of computing (STOC), pages 494503, Baltimore, MD, 2005. ACM Press. [16] A. Frangioni and C. Gentile. New preconditioners for KKT systems of network ow problems. SIAM Journal on Optimization, 14:894913, 2004. [17] Keith D. Gremban. Combinatorial Preconditioners for Sparse, Symmetric, Diagonally Dominant Linear Systems. PhD thesis, School of Computer Science, Carnegie Mellon University, October 1996. Available as Technical Report CMU-CS-96-123. [18] Keith D. Gremban, Gary L. Miller, and Marco Zagha. Performance evaluation of a new parallel preconditioner. In Proceedings of the 9th International Parallel Processing Symposium, pages 6569. IEEE Computer Society, 1995. A longer version is available as Technical Report CMU-CS-94-205, CarnegieMellon University. [19] Victoria E. Howle and Stephen A. Vavasis. An iterative method for solving complex-symmetric systems arising in electrical power modeling. SIAM Journal on Matrix Analysis and Applications, 26(4):1150 1178, 2005. [20] Anil Joshi. Topics in Optimization and Sparse Linear Systems. PhD thesis, Department of Computer Science, University of Illinois at Urbana-Champaign, 1997. [21] Ioannis Koutis and Gary L. Miller. A linear work, O(n1/6 ) time, parallel algorithm for solving planar Laplacians. In Nikhil Bansal, Kirk Pruhs, and Cliord Stein, editors, Proceedings of the Eighteenth Annual ACM-SIAM Symposium on Discrete Algorithms, SODA 2007, New Orleans, Louisiana, USA, January 7-9, 2007, pages 10021011. SIAM, 2007. [22] Ioannis Koutis and Gary L. Miller. Graph partitioning into isolated, high conductance clusters: theory, computation and applications to preconditioning. In SPAA '08: Proceedings of the twentieth annual Symposium on Parallelism in Algorithms and Architectures, pages 137145, New York, NY, USA, 2008. ACM. COMBINATORIAL PRECONDITIONERS 20 [23] Bruce M. Maggs, Gary L. Miller, Ojas Parekh, R. Ravi, and Shan Leung Maverick Woo. Finding eective support-tree preconditioners. In SPAA '05: Proceedings of the seventeenth annual ACM symposium on Parallelism in algorithms and architectures, pages 176185. ACM Pres, 2005. [24] Gil Shklarski and Sivan Toledo. Rigidity in nite-element matrices: Sucient conditions for the rigidity of structures and substructures. SIAM Journal on Matrix Analysis and Applications, 30(1):740, 2008. [25] Gil Shklarski and Sivan Toledo. Computing the null space of nite element problems. Computer Methods in Applied Mechanics and Engineering, 198(37-40):30843095, August 2009. [26] Daniel A. Spielman and Shang-Hua Teng. Solving sparse, symmetric, diagonally-dominant linear systems in time 0(m1.31 ). In Proceedings of the 44th Annual IEEE Symposium on Foundations of Computer Science, pages 416427, October 2003. [27] Daniel A. Spielman and Shang-Hua Teng. Nearly-linear time algorithms for graph partitioning, graph sparsication, and solving linear systems. In STOC '04: Proceedings of the thirty-sixth annual ACM symposium on Theory of computing, pages 8190, New York, NY, USA, 2004. ACM Press. [28] Daniel A. Spielman and Shang-Hua Teng. Nearly-linear time algorithms for preconditioning and solving symmetric, diagonally dominant linear systems. Unpublished manuscript available online at http://arxiv.org/abs/cs/0607105, 2009. [29] Daniel A. Spielman and Jaeoh Woo. A note on preconditioning by low-stretch spanning trees. Mar 2009. [30] Uri Unger. An experimental evaluation of combinatorial preconditioners. Master's thesis, Tel-Aviv University, July 2007. [31] Pravin M. Vaidya. Solving linear equations with symmetric diagonally dominant matrices by constructing good preconditioners. Unpublished manuscript. A talk based on this manuscript was presented at the IMA Workshop on Graph Theory and Sparse Matrix Computations, Minneapolis, October 1991. [32] Meiqiu Wang and Vivek Sarin. Parallel support graph preconditioners. In Yves Robert, Manish Parashar, Ramamurthy Badrinath, and Viktor K. Prasanna, editors, High Performance Computing - HiPC 2006, volume 4297, chapter 39, pages 387398. Springer Berlin Heidelberg, Berlin, Heidelberg, 2006.

- Similar pages
- Italian User's Manual