Sparse inverse incidence matrices for Schilders' factorization applied to resistor network modeling

Schilders' factorization can be used as a basis for preconditioning indefinite linear systems which arise in many problems like least-squares, saddle-point and electronic circuit simulations. Here we consider its application to resistor network modeling. In that case the sparsity of the matrix blocks in Schilders' factorization depends on the sparsity of the inverse of a permuted incidence matrix. We introduce three different possible permutations and determine which permutation leads to the sparsest inverse of the incidence matrix. Permutation techniques are based on types of sub-digraphs of the network of an incidence matrix.

arise in electronic circuit simulations and many other applications where A is symmetric and positive (semi) definite, and B T is of maximal column rank m (≤ n). Preconditioning techniques to solve (1) have become very important, especially for problems which arise from Stokes equation resulting in saddle point problems [6]. Schilders' factorization [7] can be used as a basis for such preconditioners. A possibly permutedÃ where the blocks L 1 , L 2 & M depend onB −T 1 . The permutation has to be such that B T is permuted into a lower trapezoidal form where the top partB T 1 is a lower triangular matrix of size m × m andB T 2 is an (n − m) × m matrix. Generally, A is sparse and so isÃ. But (2) can have blocks which are more dense because ofB −T 1 . In other words, sparsity of the blocks in (2) depends on the sparsity ofB −T 1 . In order to illustrate the involvement ofB −T 1 while computing the blocks in (2), we consider (4) and (5). The details for deriving these formulas can be found in [7].

1Ã 11B
−1 1 can become dense whereB −T 1 is dense. Eventually this will lead to even more dense blocks L 1 and M . This entails to have a permutation such thatB −T 1 is sparse. Therefore, we develop an algorithm which permutes an incidence matrixB T into a lower trapezoidal form (3) such that it results sparseB −T 1 . In this paper we consider Schilders' factorization applied to resistor network modeling in which A is a diagonal matrix with entries of resistance values, and B T is an incidence matrix. The modeling can be done by applying Kirchhoff's current law and Ohm's law for resistors [5]. In order to explain the definition of an incidence matrix, and to understand some of its properties, a brief summary on graph theory and incidence matrices is given in Section 2. In Section 3, aspects of the sparsity ofB −T 1 are discussed. Then we introduce three different possible permutations and determine which permutation leads to the sparsest inverse of the incidence matrix. Section 4 gives the permutation algorithm which leads to a sparse inverse,B −T 1 . Numerical experiments performed on industrial resistor networks are presented in Section 5.

2.
A brief summary on graph theory and incidence matrices. A graph G consists of a finite set V = {η 0 , η 1 , . . . , η m } called the vertex set and a finite set E = {ξ 1 , ξ 2 , . . . , ξ n } called the edge set. An η ∈ V is called a vertex (or node) of G and a ξ ∈ E an edge of G. An edge ξ is represented by an unordered pair of nodes {η i , η j }, which are said to be adjacent to each other, and called the end points of ξ.
A directed graph (or digraph) G consists of a node set V and an edge set E, where each edge ξ is an ordered pair (η i , η j ) known as an arc (or directed edge). We write η i → η j to represent that the arc ξ connects the two nodes in the direction from η i to η j . Here, call η i the initial node and η j the terminal node of ξ. If the direction between the two nodes is unknown, we shall write it as η i η j to avoid ambiguity.

SPARSE INVERSE INCIDENCE MATRICES 229
A path in a graph G is an alternating sequence of distinct nodes and edges. For example the path between η 0 and η r is given by η 0 , ξ 1 , η 1 , ξ 2 , . . . , ξ r , η r . If there is a path between every pair of nodes, then G is said to be connected. A digraph is called weakly connected if its underlying graph is connected. A digraph is strongly connected if for every pair distinct nodes η i and η j , there is a directed path starting from η i to η j . A Hamiltonian path in a digraph is a path in a single direction that visits each node exactly once. A tournament is a digraph in which every pair of nodes is connected by a single directed edge [1]. A digraph is called loopless if it contains none of (η i , η i ).
To obtain a reduced incidence matrix B T of size n×m fromB T , we call one node the reference or ground node, and delete its related column. This causes all rows of B T , related to the arcs connected to the reference node, to only have one nonzero entry. This is necessary for the construction of a permutation to a lower trapezoidal matrix. Moreover, systems of the form (1) resulting from resistor network modelings are made consistent by grounding a node and removing the corresponding column [5]. Theorem 2.3 gives an important property of an incidence matrix which relates to its underlying digraph.

Theorem 2.3. [2]
Suppose G is a digraph with m + 1 nodes. Then the column rank of the incidence matrix B T is m if and only if G is weakly connected.
3. Sparsity of inverse incidence matrix. Before exploring different types of permutations, we first give an overview of interdependence between the sparsity of inverse and nilpotency degree of a lower triangular matrix which contains maximally two nonzero entries in each row. For this we consider a unit lower bidiagonal matrix H and compute its inverse using the expansion of (I + N H ) −1 , where I is the identity matrix and N H is the strictly lower triangular part of H (the nilpotent part of H). The inverse of unit lower bidiagonal matrices can also be computed by using Gaussian elimination [8]. However, the former gives an insight on how the interdependence between the sparsity of inverse and nilpotency degree can be conceived. A unit lower bidigonal matrix H of size n × n is of the form it is trivial to observe that That is every subdiagonal entry of H −1 is a product of some β i 's and not equal to zero. So H −1 is full and N H has nilpotency degree k = n. On the other hand if a lower triangular matrix with maximally two nonzero entries in each row is of the form then again it is easy to obtain which is as sparse as G and the nilpotent part N G has degree k = 2. Based on these two extreme cases, we consider an example of an incidence matrix B T which is constructed from a digraph G as shown in Figure 1. It is permuted to a lower trapezoidal formB T (3) in two different ways giving two differentB T 1 's of size 4 × 4. Sparsity ofB −T 1 depends on the type of underlying digraph. If the permutation is based on a Hamiltonian path (Figure 1 (a)), thenB −T 1 is a full lower triangular matrix. If it is based on a star digraph (Figure 1 (b)), then one can obtain the most sparseB −T 1 . The former has nilpotency degree k = 4, while the latter has k = 2. Therefore we are also induced to look for a permutation that leads to a small nilpotency degree. It is always possible to find star sub-digraphs in large networks. from a digraph G whose ground node is shaded.
In order to illustrate further the sparsity ofB −T 1 by looking at examples of larger sizes, we generated incidence matrices of full column rank. The easiest way to generate such incidence matrix is to consider a digraph containing a Hamiltonian path. This is being indicated by Theorems 2.1 and 2.3 (see also in [4]). Construct a tournament T from a randomly generated m + 1 nodes η 0 , η 1 , . . . , η m , where η 0 is the ground node. And then choose randomly a sub-digraph G ⊆ T of n arcs such that it contains the Hamiltonian path. The required incidence matrix B T can be constructed from G using Algorithm 1. Note that Algorithm 1 can be used only for generating a random test incidence matrix. It ensures to construct an incidence matrix B T of full column rank. It should also be noted that n ≤ (m + 1)m/2 since a tournament from m + 1 nodes contains only (m + 1)m/2 arcs. However if G is already known to be fulfilling the requirements of Theorem 2.3, then it is not necessary for it to contain a Hamiltonian path. For such G, we can directly construct B T by using the steps 10 − 14 of Algorithm 1.
Using the above technique, a random network with an incidence matrix B T of size 100 × 60 is generated. Then three different permutations are applied on B T to reduce it to a lower trapezoidal formB T as shown in Figure 2. Permutation I is done by choosing a Hamiltonian path in G. It is evident from Figure 2, that permutation I does not promise good results of sparse blocks in Schilders' factorization since it gives full inverse incidence matrices. Permutation II is done randomly without taking any other considerations than the lower trapezoidal form. This permutation may sometimes give sparse inverse incidence matrices. However, it is quite difficult to predict the sparsity in advance since the permutation does not emphasize any special connectivity of the underlying digraph. Permutation III is done by prioritizing on the star sub-digraphs of G that are connected to each other starting from the ground node. It leads to the sparsestB −T 1 and the smallest nilpotency degree comparing to the other two. Thus we propose Permutation III in the next section. 4. Permutation which leads to small nilpotency degree and sparseB −T 1 . Let G be the digraph of a node-arc incidence matrixB T of size n × (m + 1). Let the node set of G be V = {η j } m j=0 , where η 0 is the reference node, and the arc set be E = {ξ i } n i=1 . After deleting the column ofB T corresponding to η 0 , we obtain a reduced incidence matrix B T of size n × m. Define Without loss of generality, assume that B T has maximal column rank m and m ≤ n. Then by Theorem 2.3, G is at least weakly connected. Note that each row of B T contains maximally two non-zero entries from {−1, 1}. This leads us to represent the column numbers of B T by the elements of V and its row numbers by elements of E. Let R [is] be an n × n row permutation matrix such that the rows s and i s are permuted. Similarly denote an m × m column permutation matrix by C [js] indicating that the columns s and j s are permuted.
Suppose η j1 η j0 by the arc ξ i1 , where η j0 = η 0 , η j1 ∈ V , and ξ i1 ∈ E. Permuting the 1 st row with i th 1 row and the 1 st column with j th 1 column of B T would give , where B T (1) is an (n − 1) × (m − 1) reduced incidence matrix, u is an (n − 1) × 1 column vector, and b i1j1 ∈ {−1, 1}. Similarly, continue the process until we find the last remaining node, η jr η j0 by the arc ξ ir , where 2 ≤ r ≤ m and 1 ≤ i r ≤ n. So that permuting the r th row with i th r row and the r th column with j th r column would give B T with the following structure.
where B T (r) is now of size (n − r) × (m − r), and v 1 , . . . , v r are the column vectors of size (n − r) × 1. A diagonal entry b isjs ∈ {−1, 1}, 1 ≤ s ≤ r, is the only nonzero entry in s th row. This completes permuting the columns and rows related to the nodes and arcs of the star sub-digraph whose central node is η j0 . Now, suppose each of η jr+1 , . . . , η jp η j1 by the arcs ξ ir+1 , . . . , ξ ip , respectively, where r + 1 ≤ p ≤ m. If there exists no such connection for η j1 , then there will be at least one such connection for one of η j2 , . . . , η j r . Otherwise G is disconnected. As soon as one such connection is found, permute every corresponding r + 1 th , . . . , p th row with i th r+1 , . . . , i th p row and every corresponding r + 1 th , . . . , p th column with j th r+1 , . . . , j th p column. This is the permutation with respect to another star sub-digraph whose central node is one of the nodes η j1 , . . . , η j r , connected to the previous star sub-digraph. Repeating the same process for p + 1, p + 2, . . . , m, we obtain B T with a lower trapezoidal form: The top m × m lower triangular part isB T 1 , and the remaining (n − m) × m rectangular matrix isB T 2 = [w 1 , . . . , w r , . . . , w p . . . , w m ], where w s , 1 ≤ s ≤ m, is an (n − m) × 1 column vector. Here, ξ i N = ξ im+1 , . . . , ξ in , which gives the row numbers ofB T 2 . In fact we have shown that if B T is an n × m incidence matrix with full column rank, containing maximally two nonzero entries in each row, then there exist an n × n row permutation matrix R and an m × m column permutation matrix C such that whereB T 1 is an m × m lower triangular andB T 2 is an (n − m) × m rectangular matrix.
Consider an example of a digraph shown in Figure 3, which has 9 nodes (including the ground node η 0 ) and 12 arcs. Using Definition 2.2, we construct a node-arc ξ 12 ξ 10 Figure 3. A digraph.
incidence matrixB T of size 12 × 9 as follows: After deleting the column related to ground node, η 0 , a reduced incidence matrix B T of size 12 × 8 is obtained as follows: Applying the new permutation technique, we obtain the permuted incidence matrix B T as follows:B The row and column permutation matrices, respectively are given by R = [e 9 , e 11 , e 12 , e 4 , e 7 , e 3 , e 2 , e 6 , e 1 , e 10 , e 5 , e 8 ] T and C = [e 5 , e 7 , e 8 , e 4 , e 6 , e 2 , e 1 , e 3 ], where e i denotes the i th unit vector of size 12 × 1 for R and 8 × 1 for C. It is easy to see that the nilpotent part of B T 1 has a degree k = 3.
Algorithm 1 Generates a random full column rank incidence matrix.
Require: Number of nodes m, and number of arcs n such that m ≤ n. Ensure: Node-arc incidence matrixB T . 1: generate {η 1 , η 2 , . . . , η m } 2: for j = 1 to m − 1 do 3: for i = 1 to m − j do 4: end for 6: end for 7: for i = 1 to n do 8: G = {ξ i : ξ i ∈ T } ⊃ weakly connected path 9: end for 10: for i = 1 to n do 11:B T (i, G(i, 1)) = 1 12:B T (i, G(i, 2)) = −1 13: end for 14: delete the ground node column ofB T to obtain a reduced incidence matrix B T Algorithm 2 Permutes an incidence matrix to a lower trapezoidal form which leads to a sparse inverse.
Require: Incidence matrix B T . Ensure: Lower trapezoidal incidence matrixB T , row and column permutation matrices P R and P C . 1: P R = I n×n , P C = I m−1×m−1 , k = 1 2: for i = 1 to n do 3: find column index r such that B T (i, r) is the only non-zero entry in i th row 4: permute rows i and k of B T and P R

5:
permute columns r and k of B T and P C 6: k = k + 1 (only after each permutation) 7: end for 8: s = k (as updated above) 9: j = 1 10: while s ≤ m do  Figure 4). Component R 1 is banded and has star sub-digraphs with central nodes connected to only a few number of arcs as shown in Figure 4 (a). Due to this banded structure, the two diagonals ofB T 1 are quite close to each other creating more fill-ins in the inversẽ B −T 1 . Component R 2 in Figure 4 (b) has star sub-digraphs with many number of arcs connected to each central node, which results large distance between the two diagonals ofB T 1 giving sparseB −T 1 . Thus sparsity ofB −T 1 depends on the number of arcs connected to each central node of the underlying star sub-digraphs as explained in Section 3. Sparsity of the blocks from other components is shown in Table 1. 6. Concluding remarks. The sparsity of blocks in Schilders' factorization depends on the sparsity of the inverse incidence matrixB −T 1 . Permutation algorithm to achieve a potential optimum nilpotency degree ofB T 1 is important. We have examined three different permutations and proposed the one which leads to the sparsest inverse. The technique of permutation is based on the star sub-digraphs connected to each other starting from the ground node. Sparsity of the inverseB −T 1 depends on the number of arcs connected to each central node of the underlying star sub-digraphs. We have developed an algorithm which takes care efficiently all the above conditions.

Resistor
Size PermutationB T