A parallel solver for large-scale Markov chains.
We consider the parallel computation of the stationary probability distribution vector of ergodic Markov chains with large state spaces by preconditioned Krylov subspace methods. The parallel preconditioner is obtained as an explicit approximation, in factorized form, of a particular generalized inverse of the generator matrix of the Markov process. Graph partitioning is used to parallelize the whole algorithm, resulting in a two-level method.Conditions that guarantee the existence of the preconditioner are given, and the results of a parallel implementation are presented. Our results indicate that this method is well suited for problems in which the generator matrix can be explicitly formed and stored.
Introduction
Discrete Markov chains with large state spaces arise in many applications,
including for instance reliability modeling, queueing network analysis, large
scale economic modeling and computer system performance evaluation. The
stationary probability distribution vector of an ergodic Markov process with
n \Theta n transition probability matrix P is the unique 1 \Theta n vector - which
satisfies
Letting the computation of the stationary vector
reduces to finding a nontrivial solution to the homogeneous linear system
0. The ergodicity assumption means that P (and therefore A) is irreducible.
Perron-Frobenius theory [9] guarantees that A has rank n \Gamma 1, and that the
(one-dimensional) null space N (A) of A is spanned by a vector x with positive
entries. Upon normalization in the ' 1 -norm, this is the stationary distribution
vector of the Markov process.
The coefficient matrix A is a singular M-matrix, called the generator of the
Markov process. 3 The matrix A is nonsymmetric, although it is sometimes
structurally symmetric. See [26] for a good introduction to Markov chains and
their numerical solution.
Due to the very large number n of states typical of many real-world applica-
tions, there has been increasing interest in recent years in developing parallel
algorithms for Markov chain computations; see [2], [5], [10], [17], [19], [24].
Most of the attention so far has focused on (linear) stationary iterative meth-
ods, including block versions of Jacobi and Gauss-Seidel [10], [19], [24], and on
iterative aggregation/disaggregation schemes specifically tailored
to stochastic matrices [10], [17]. In contrast, little work has been done with
parallel preconditioned Krylov subspace methods. Partial exceptions are [5],
where a symmetrizable stationary iteration (Cimmino's method) was accelerated
using conjugate gradients on a Cray T3D, and [19], where an out-of-core,
parallel implementation of Conjugate Gradient Squared (with no precondi-
tioning) was used to solve very large Markov models with up to 50 million
states. The suitability of preconditioned Krylov subspace methods for solving
Markov models has been demonstrated, e.g., in [25], although no discussion
of parallelization issues is given there.
In this paper we investigate the use of a parallel preconditioned iterative
method for large, sparse linear systems in the context of Markov chain com-
Strictly speaking, the generator matrix is We work with A
instead of Q to conform to the familiar notation of numerical linear algebra.
putations. The preconditioning strategy is a two-level method based on sparse
approximate inverses, first introduced in [3]. However, due to the singularity
of the generator matrix A, the applicability of approximate inverse techniques
in this context is not obvious. That this is indeed possible is a consequence of
the fact that A is a (singular) M-matrix.
The paper is organized as follows. In section 2 we discuss the problem of pre-conditioning
singular equations in general, and we establish a link between
some standard preconditioners and generalized inverses. Sections 3-5 are devoted
to AINV preconditioning for Markov chain problems, including a discussion
of the parallel implementation and a theoretical analysis of the existence
of the preconditioner. Numerical tests are reported in section 6, and some
conclusions are presented in section 7.
Preconditioning Markov chain problems
In the Markov chain context, preconditioning typically amounts to finding
an easily invertible nonsingular matrix M (the preconditioner) which is a
good approximation to A; a Krylov subspace method is then used to solve
preconditioning). Notice that even if A itself is singular, the preconditioner
must be nonsingular so as not to change the solution set, i.e., the null space
(A) of A. Preconditioners can be generated by means of splittings
N , such as those used in stationary iterative methods including Jacobi, Gauss-
Seidel, SOR and block versions of these schemes; see [26]. Also in this class
are the popular incomplete LU (ILU) factorization preconditioners. ILU-type
methods have been successfully applied to Markov chain problems by Saad
[25] in a sequential environment. The existence of incomplete factorizations
for nonsingular M-matrices was already proved in [20]; an investigation of the
existence of ILU factorizations for singular M-matrices can be found in [11].
Incomplete factorization methods work quite well on a wide range of problems,
but they are not easily implemented on parallel computers. For this and other
reasons, much effort has been put in recent years into developing alternative
preconditioning strategies that have natural parallelism while being comparable
to ILU methods in terms of robustness and convergence rates. This work
has resulted in several new techniques known as sparse approximate inverse
preconditioners; see [7] for a recent survey and extensive references. Sparse
approximate inverse preconditioners are based on directly approximating the
inverse of the coefficient matrix A with a sparse matrix G - A \Gamma1 . The application
of the preconditioner only requires matrix-vector products, which
are easily parallelized. Until now, these techniques have been applied almost
exclusively to nonsingular systems of equations b. The only exception
seems to be [13], where the SPAI preconditioner [16] was used in connection
with Fast Wavelet Transform techniques on singular systems stemming from
discretizations of the Neumann problem for Poisson's equation.
The application of approximate inverse techniques in the singular case raises
several interesting theoretical and practical questions. Because the inverse of
A does not exist, it is not clear what matrix G is an approximation of. It
should presumably be some generalized inverse of A, but which one? Note
that this question can be asked of M \Gamma1 for any preconditioner M - A. In
[26], page 143, it is stated that M \Gamma1 should be an approximation of the group
generalized inverse A ] , and that an ILU factorization A -
U implicitly yields
such an approximation: ( -
. As we will see, this interpretation is not
entirely correct and is somewhat misleading. The group inverse (see [12]) is
only one of many possible generalized inverses. It is well known [21] that the
group inverse plays an important role in the modern theory of finite Markov
chains. However, it is seldom used as a computational tool, in part because
its computation requires knowledge of the stationary distribution vector -.
As it turns out, different preconditioners result (implicitly or explicitly) in
approximations M \Gamma1 to different generalized inverses of A, which are typically
not the group inverse A ] . Let us consider ILU preconditioning first. If A is a n\Theta
irreducible, singular M-matrix, then A has the LDU factorization
where L and U are unit lower and upper triangular matrices (respectively) and
D is a diagonal matrix of rank
(see [26]). Notice that L and U are nonsingular M-matrices; in particular, L \Gamma1
and U \Gamma1 have nonnegative entries. Define the matrix
It can be easily verified that A \Gamma satisfies the first two of Penrose's four conditions
[12]:
The first identity states that A \Gamma is an inner inverse of A and the second
that A \Gamma is an outer inverse of A. A generalized inverse satisfying these two
conditions is called a (1; 2)-inverse of A or an inner-outer inverse. Another
term that is found in the literature is reflexive inverse; see [12]. Because A \Gamma
does not necessarily satisfy the third and fourth Penrose conditions, it is not
the Moore-Penrose pseudoinverse A y of A in general. Because A y is obviously
a (1; 2)-inverse, this kind of generalized inverse is non-unique. Indeed, there
are infinitely many such (1; 2)-inverses in general. Each pair R, N of subspaces
of IR n that are complements of the null space and range of A (respectively)
uniquely determines a (1; 2)-inverse GN;R of A with null space N (GN;R
and range R(GN;R see [12]. In the case of A \Gamma it is readily verified that
denotes the i-th unit
basis vector in IR n . It is easy to see that R is complementary to N (A) and N
is complementary to R(A). The pseudoinverse A y corresponds to
The (1; 2)-inverse A \Gamma is also different from the group inverse A ] , in general.
This can be seen from the fact that in general AA \Gamma 6= A \Gamma A, whereas the group
inverse always satisfies AA A. Also notice that for a singular irreducible
M-matrix A the (1; 2)-inverse A nonnegative ma-
trix, which is not true in general for either the group or the Moore-Penrose
Let now A -
U be an incomplete LDU factorization of A, with -
unit lower triangular, -
U - U unit upper triangular and -
diagonal matrix with positive entries on the main diagonal. Then clearly
Hence, an ILU factorization of A yields an implicit approximation to A \Gamma rather
than to A ] . This can also be seen from the fact that ( -
always nonnegative, which is not true of A ] .
It is straightfoward to check that A \Gamma A is the oblique projector onto
along N (A) and that AA \Gamma is the oblique projector onto
along g. Therefore A \Gamma A has eigenvalues 0 with multiplicity
1, and 1 with multiplicity likewise for AA \Gamma . Hence, it makes good sense
to construct preconditioners based on approximating A \Gamma (either implicitly or
explicitly), since in this case most eigenvalues of the preconditioned matrix
will be clustered around 1.
Next we consider the approximate inverse preconditioner AINV; see [4], [6].
This method is based on the observation that if Z and W are matrices whose
columns are A-biorthogonal, then W T diagonal matrix. When all
the leading principal minors of A (except possibly the last one) are nonzeros,
Z and W can be obtained by applying a generalized Gram-Schmidt process
to the unit basis vectors e 1 . In this case Z and W are unit upper
triangular. It follows from the uniqueness of the LDU factorization that
U \Gamma1 and LDU is the LDU factorization of A. The
diagonal matrix D is the same in both factorizations. An approximate inverse
in factorized form
W T can be obtained by dropping small entries
in the course of the generalized Gram-Schmidt process. Similar to ILU, this
incomplete inverse factorization is guaranteed to exist for nonsingular M -
matrices [4]; see the next section for the singular M-matrix case. In either
case,
W T is a nonnegative matrix. Hence, this preconditioner can
be interpreted as a direct (explicit) approximation to the (1; 2)-inverse A \Gamma of
Lastly, we take a look at sparse approximate inverse techniques based on
Frobenius norm minimization; see, e.g., [16] and [14]. With this class of meth-
ods, an approximate inverse G is computed by minimizing the functional
subject to some sparsity constraints. Here jj \Delta jj F denotes
the Frobenius matrix norm. The sparsity constraints could be imposed
a priori, or dynamically in the course of the algorithm. In either case, it is natural
to ask what kind of generalized inverse is being approximated by G when
A is a singular matrix. It can be shown that the Moore-Penrose pseudoinverse
A y is the matrix of smallest Frobenius norm that minimizes jjI \Gamma AXjj F ;
see for instance [23], page 428. Hence, in the singular case the SPAI preconditioner
can be seen as a sparse approximate Moore-Penrose inverse of A.
This is generally very different from the approximate (1; 2)-inverses obtained
by either ILU or AINV. For instance, SPAI will not produce a nonnegative
preconditioner in general.
In the next section we restrict our attention to the AINV preconditioner and
its application to Markov chain problems.
3 The AINV method for singular matrices
The AINV preconditioner [4], [6] is based on A-biorthogonalization. This is
a generalized Gram-Schmidt process applied to the unit basis vectors e i ,
n. In this generalization the standard inner product is replaced by
the bilinear form h(x; Ay. This process is well defined, in exact arith-
metic, if the leading principal minors of A are nonzero, otherwise some form
of pivoting (row and/or column interchanges) may be needed. If A is a nonsingular
M-matrix, all the leading principal minors are positive and the process
is well defined with no need for pivoting. This is perfectly analogous to the LU
factorization of A, and indeed in exact arithmetic the A-biorthogonalization
process computes the inverses of the triangular factors of A. When A is a
singular irreducible M-matrix, all the leading principal minors of A except
the n-th one (the determinant of are positive, and the process can still be
completed.
In order to obtain a sparse preconditioner, entries (fill-ins) in the inverse factors
Z and W less than a given drop tolerance in magnitude are dropped in the
course of the computation, resulting in an incomplete process. The stability
of the incomplete process for M-matrices was analyzed in [4]. In particular, if
d i denotes the i-th pivot, i.e., the i-th diagonal entry of -
D in the incomplete
process, then (Proposition 3.1 in [4]) -
Because
for an M-matrix, no pivot breakdown can occur.
Exactly the same argument applies to the case where A is an irreducible,
singular M-matrix. In this case there can be no breakdown in the first
steps of the incomplete A-biorthogonalization process, since the first
leading principal minors are positive, and the pivots in the incomplete process
cannot be smaller than the exact ones. And even if the n-th pivot -
happened
to be zero, it could simply be replaced by a positive number in order to have
a nonsingular preconditioner. The argument in [4] shows that -
d n must be a
nonnegative number, and it is extremely unlikely that it will be exactly zero
in the incomplete process.
Another way to guarantee the nonsingularity of the preconditioner is to perturb
the matrix A by adding a small positive quantity to the last diagonal
entry. This makes the matrix a nonsingular M-matrix, and the incomplete
A-biorthogonalization process can then be applied to this slightly perturbed
matrix to yield a well defined, nonsingular preconditioner. In practice, how-
ever, this perturbation is not necessary, since dropping in the factors typically
has an equivalent effect (see Theorem 3 below).
The AINV preconditioner has been extensively tested on a variety of symmetric
and nonsymmetric problems in conjunction with standard Krylov subspace
methods like conjugate gradients (for symmetric positive definite matrices)
and GMRES, Bi-CGSTAB and TFQMR (for unsymmetric problems). The
preconditioner has been found to be comparable to ILU methods in terms
of robustness and rates of convergence, with ILU methods being somewhat
faster on average on sequential computers. The main advantage of AINV over
the ILU-type methods is that its application within an iterative process only
requires matrix-vector multiplies, which are much easier to vectorize and to
parallelize than triangular solves.
Unfortunately, the computation of the preconditioner using the incomplete A-
biorthogonalization process is inherently sequential. One possible solution to
this problem, adopted in [8], is to compute the preconditioner sequentially on
one processor and then to distribute the approximate inverse factors among
processors in a way that minimizes communication costs while achieving good
load balancing. This approach is justified in applications, like those considered
in [8], in which the matrices are small enough to fit on the local memory of
one processor, and where the preconditioner can be reused a number of times.
In this case the time for computing the preconditioner is negligible relative
to the overall costs. In the Markov chain setting, however, the preconditioner
cannot be reused in general and it is imperative that set-up costs be minimized.
Furthermore, Markov chain problems can be very large, and it is desirable to
be able to compute the preconditioner in parallel.
4 The parallel preconditioner
In the present section we describe how to achieve a fully parallel precon-
ditioner. The strategy used to parallelize the preconditioner construction is
based on the use of graph partitioning. This approach was first proposed in
[3] in the context of solving nonsingular linear systems arising from the discretization
of partial differential equations.
The idea can be illustrated as follows. If p processors are available, graph
partitioning can be used to decompose the adjacency graph associated with
the sparse matrix A (or with A A is not structurally symmetric) in
p subgraphs of roughly equal size in such a way that the number of edge
cuts is approximately minimized. Nodes which are connected by cut edges
are removed from the subgraphs and put in the separator set. By numbering
the nodes in the separator set last, a symmetric permutation Q T AQ of A is
obtained. The permuted matrix has the following structure:
The diagonal blocks A correspond to the interior nodes in the graph
decomposition, and should have approximately the same order. The off-diagonal
blocks connections between the subgraphs, and the diagonal
block A S the connections between nodes in the separator set. The order of
A S is equal to the cardinality of the separator set and should be kept as small
as possible. Note that because of the irreducibility assumption, each block A i
is a nonsingular M-matrix, and each of the LDU factorizations A
LDU be the LDU factorization of Q T AQ. Then it is easy to see
that
S is the inverse of the unit lower
triangular factor of the Schur complement matrix
In the next section we show that S is a singular, irreducible M-matrix, hence
it has a well defined LDU factorization
Likewise,
U
S is the inverse of the unit upper
triangular factor of S. It is important to observe that L \Gamma1 and U \Gamma1 preserve
a good deal of sparsity, since fill-in can occur only within the nonzero blocks.
The matrix D is simply defined as
note that all diagonal entries of D are positive except for the last one, which
is zero. The (1; 2)-inverse D \Gamma of D is defined in the obvious way.
Hence, we can write the (generally dense) generalized inverse (Q T AQ) \Gamma of
as a product of sparse matrices L In practice, however,
the inverse factors L \Gamma1 and U \Gamma1 contain too many nonzeros. Since we are only
interested in computing a preconditioner, we just need to compute sparse
approximations to L \Gamma1 and U \Gamma1 .
This is accomplished as follows. With graph partitioning, the matrix is distributed
so that processor P i holds A One of the
processors, marked as P S , should also hold A S . Each processor then computes
sparse approximate inverse factors -
W i such that -
using the AINV algorithm. Once this is done, each processor computes the
product
. Until this point the computation
proceeds in parallel with no communication. The next step is the accumulation
of the approximate Schur complement -
This accumulation is done in steps with a fan-in across the
processors. In the next section we show that although the exact Schur complement
S is singular, the approximate Schur complement -
S is a nonsingular
M-matrix under rather mild conditions.
As soon as -
S is computed, processor P S computes a factorized sparse approximate
using the AINV algorithm. This is a sequential
bottleneck, and explains why the size of the separator set must be kept small.
Once the approximate inverse factors of -
are computed, they are broadcast
to all remaining processors. (Actually, the preconditioner application can be
implemented in such a way that only the -
needs to be broadcast.)
Notice that because only matrix-vector products are required in the application
of the preconditioner, there is no need to
explicitly. In this way, a factorized sparse approximate
(1; 2)-inverse of Q T AQ is obtained.
This is a two-level preconditioner, in the sense that the computation of the
preconditioner involves two phases. In the first phase, sparse approximate inverses
of the diagonal blocks A i are computed. In the second phase, a sparse
approximate inverse of the approximate Schur complement -
S is computed.
this second step the preconditioner would reduce to a block Jacobi
method with inexact block solves (in the terminology of domain decomposition
methods, this is additive Schwarz with inexact solves and no overlap).
It is well known that for a fixed problem size, the rate of convergence of
this preconditioner tends to deteriorate as the number of blocks (subdomains)
grows. Hence, assuming that each block is assigned to a processor in a parallel
computer, this method would not be scalable. However, the approximate
Schur complement phase provides a global exchange of information across the
processors, acting as a "coarse grid" correction in which the "coarse grid"
nodes are interface nodes (i.e., they correspond to vertices in the separator
set). As we will see, this prevents the number of iterations from growing as
the number of processors grows. As long as the cardinality of the separator
set is small compared to the cardinality of the subdomains (subgraphs), the
algorithm is scalable in terms of parallel efficiency. Indeed, in this case the
application of the preconditioner at each step of a Krylov subspace method
like GMRES or Bi-CGSTAB is easily implemented in parallel with relativeley
little communication needed.
5 The approximate Schur complement
In this section we investigate the existence of the approximate (1; 2)-inverse
of the generator matrix A. The key role is played by the (approximate) Schur
complement.
First we briefly review the situation for the case where A is a nonsingular
M-matrix. Assume A is partitioned as
A 21 A 22
Then it is well-known that the Schur complement
is also a nonsingular M-matrix; see, e.g., [1]. Moreover, the same is true of
any approximate Schur complement
provided that O - X 11 - A \Gamma1
, where the inequalities hold componentwise;
see [1], page 264.
In the singular case, the situation is slightly more complicated. In the following
we will examine some basic properties of the exact Schur complement of a
singular, irreducible M-matrix A corresponding to an ergodic Markov chain.
Recall that is the irreducible row-stochastic transition
probability matrix.
be an irreducible row-stochastic matrix partitioned as
Assume that
A 21 A 22
Then the Schur complement S of A 11 in A is a singular, irreducible M-matrix
with a one-dimensional null space.
Proof: Consider the stochastic complement [22] \Sigma of P 22 in
Note that I \Gamma P 11 is invertible since P is irreducible. From the theory developed
in [22], we know that \Sigma is row stochastic and irreducible since P is. Consider
now the Schur complement S of A 11 in A:
22
Clearly, S is an irreducible singular M-matrix. It follows (Perron-Frobenius
theorem) that S has a one-dimensional null space. 2
The previous lemma is especially useful in cases where the exact Schur complement
is used. In the context of preconditioning it is often important to
know properties of approximate Schur complements. As shown in the previous
section, graph partitioning induces a reordering and block partitioning of
the matrix A in the form (1) where
and
We are interested in properties of the approximate Schur complement -
obtained
by approximating the inverses of the diagonal blocks A i with AINV:
In particular, we are interested in conditions that guarantee that -
S is a non-singular
M-matrix, in which case the AINV algorithm can be safely applied to
S, resulting in a well defined preconditioner. We begin with a lemma. Recall
that a Z-matrix is a matrix with nonpositive off-diagonal entries [9].
Lemma 2 Let S be a singular, irreducible M-matrix. Let C - O; C 6= O be
such that -
S is a nonsingular M-matrix.
Proof: Since S is a singular M-matrix, we have
where ae(B) denotes the spectral radius of B. Let
is a Z-matrix, we have that B - C \Gamma D. Therefore we can write the modified
S as -
We distinguish the two following simple cases: Assume first that
S is a nonsingular M-matrix since it can be written as ae(B)I \Gamma (B \Gamma C) with
C). The last inequality follows from the irreducibility of B
and properties of nonnegative matrices; see [9], page 27, Cor. 1.5 (b). Assume,
on the other hand, that
B. Note that by
assumption, at least one of the diagonal entries d ii of D must be positive. Let
denote the largest such positive diagonal entry. Since B + ffiI is irreducible,
similar to the previous case. It follows
that -
D) is a nonsingular M-matrix.
Finally, if both D 6= O and C 6= D the result follows by combining the two
previous arguments. 2
In the context of our parallel preconditioner, this lemma says that if the inexactness
in the approximate inverses of the diagonal blocks A i results in an
approximate Schur complement -
S that is still a Z-matrix and furthermore if
nonnegative and nonzero, then -
S is a nonsingular M-matrix.
The following theorem states sufficient conditions for the nonsingularity of the
approximate Schur complement. B i  and B \Lambdaj denote the i-th row and the j-th
column of matrix B, respectively.
Theorem 3 Let
A 21 A 22
be a singular irreducible M-matrix, with A 11 2 IR m\Thetam . Assume that ~
A
11 is an
approximation to A \Gamma1
11 such that O - ~
A
A
11 . Furthermore,
assume that there exist indices such that the following three
conditions are satisfied:
A
the approximate Schur complement -
~
A
A 12 is a nonsingular M-matrix

Proof: From Lemma 1 we know that the exact Schur complement S of A 11
in A is a singular, irreducible M-matrix. Note that the approximate Schur
complement -
S induced by an approximation ~
A
11 of the block A \Gamma1
11 and the
exact Schur complement S are related as follows:
~
A
A \Gamma1)A
with C - O since A 21 - O; A 12 - O; and A \Gamma1
A
S is a
Z-matrix by its definition. From the assumption that ~
A \Gamma16= A \Gamma1and that
there exists at least one nonzero entry ff ij of ~
A
11 not equal to (in fact, strictly
less than) the corresponding entry of A \Gamma1
11 , we see that if the corresponding
row i of A 21 and column j of A 12 are both nonzero, there is a nonzero entry
in C: The result now follows from Lemma 2. 2
Let us now apply these results to the preconditioner described in the previous
section. In this case, A 11 is block diagonal and the AINV algorithm is used
to approximate the inverse of each diagonal block separately, in parallel. The
approximate Schur complement (2) is the result of subtracting p terms of the
We refer to these terms as (Schur
complement) updates. Each one of these updates is nonnegative and approximates
the exact update C
from below in the (entrywise) nonnegative
ordering, since O -
see [6]. Theorem 3 says that as long
as at least one of these updates has an entry that is strictly less than the
corresponding entry in C
, the approximate Schur complement -
S is a
nonsingular M-matrix.
In practice, these conditions are satisfied as a result of dropping in the approximate
inversion of the diagonal blocks A i . It is nevertheless desirable to have
rigorous conditions that ensure nonsingularity. The following theorem gives a
sufficient condition for having a nonsingular approximate Schur complement
as a consequence of dropping in AINV. Namely, it specifies conditions under
which any dropping forces -
S to be nonsingular. Note that the conditions of
this theorem do not apply to the global matrix (1), since A 11 is block diagonal
and therefore reducible. However, they can be applied to any individual
Schur complement update for which the corresponding diagonal block A i is
irreducible, making the result fairly realistic.
Theorem 4 Let A 11 2 IR m\Thetam and the singular M-matrix
A 21 A 22
be both irreducible. Let A be the LU factorization of A 11 . Assume
that in each column of L 11 (except the last one) and in each row of U 11 (except
the last one) there is at least one nonzero entry in addition to the diagonal
one:
and
Denote by -
Z 11
11 the factorized sparse approximate inverse of A 11
obtained with the AINV algorithm. Then the approximate Schur complement
S is a nonsingular M-matrix provided that -
Z 11 6= Z 11 and -
Proof: First note that the two conditions for nonzero entries in L 11 and U 11 are
implied by similar conditions for the entries in the lower and upper triangular
parts of A 11 . Namely, it is easy to see that (barring fortuitous cancellation)
and for 0:
Here tril(B) and triu(B) denote the lower and upper triangular part of matrix
B, respectively. These conditions are easier to check than the weaker ones on
the triangular factors of A 11 . Conditions (3) and (4) imply that there is a path
in the graph of L T
11 and there is a path in the graph of U 11 for
Because the matrix is irreducible, then A 12 6= O and A 21 6= O. This means
that there exist indices s such that rs 6= 0: The
existence of the previously mentioned paths implies that
where W 11 and Z 11 are the exact inverse factors of A 11 . The approximate
inverse factors from the AINV algorithm satisfy [6]
O -
Z 11
Hence, the conditions of Theorem 3 are satisfied and the result is proved. 2
It is instructive to consider two extreme cases. If A 11 is diagonal, then the
approximate Schur complement is necessarily equal to the exact one, and is
therefore singular. In this case, of course, the conditions of the last theorem
are violated. On the other hand if A 11 is irreducible and tridiagonal, its inverse
factors are completely dense and by the last theorem it is enough to drop a
single entry in each inverse factor to obtain a nonsingular approximate Schur
complement.
The purpose of the theory developed here is to shed light on the observed
robustness of the proposed preconditioner rather than to serve as a practical
tool. In other words, it does not seem to be necessary to check these conditions
in advance. Indeed, thanks to dropping the approximate Schur complement
was always found to be a nonsingular M-matrix in actual computations.
6 Numerical experiments
In this section we report on results obtained with a parallel implementation of
the preconditioner on several Markov chain problems. The underlying Krylov
subspace method was Bi-CGSTAB [27], which was found to perform well for
Markov chains in [15]. Our FORTRAN implementation uses MPI and dynamic
memory allocation. The package METIS [18] was used for the graph
partitioning, working with the graph of A + A T whenever A was not structurally
symmetric.
The test problems arise from real Markov chain applications and were provided
by T. Dayar. These matrices have been used in [15] to compare different
methods in a sequential environment. A description of the test problems is
provided in Table 1 below. Here n is the problem size and nnz the number of
nonzeros in the matrix. All the test problems are structurally nonsymmetric
except ncd and mutex. Most matrices are unstructured.

Tables

2-11 contain the test results. All runs were performed on an SGI Origin
2000 at Los Alamos National Laboratory (using up to 64 processors), except
for those with matrices leaky, ncd and 2d which were performed on an Origin
2000 at the Helsinki University of Technology (using up to 8 processors). In
all cases, the initial guess was a constant nonzero vector; similar results were
obtained with a randomly generated initial guess. In the tables, "P-time"
denotes the time to compute the preconditioner, "P-density" the ratio of the
number of nonzeros in the preconditioner to the number of nonzeros in the
matrix A, "Its" denotes the number of iterations needed to reduce the ' 2 -
norm of the initial residual by eight orders of magnitude, "It-time" the time
to perform the iterations, and "Tot-time" the sum of "P-time" and "It-time."
All timings are in seconds. Furthermore, "Sep-size" is the cardinality of the
Information on test problems.
Matrix n nnz Application
hard 20301 140504 Complete buffer sharing in ATM networks
Multiplexing model of a leaky bucket
2d 16641 66049 A two-dimensional Markov chain model
telecom 20491 101041 A telecommunication model
ncd 23426 156026 NCD queueing network
mutex 39203 563491 Resource sharing model
qn 104625 593115 A queueing network
separator set (i.e., the order of the Schur complement matrix) and "Avg-
dom" the average number of vertices in a subdomain (subgraph) in the graph
partitioning of the problem. The drop tolerance - in the AINV algorithm was
the same at both levels of the preconditioner (approximate inversion of A i for
approximate inversion of the approximate Schur complement
S), except for the mutex problem (see below).

Tables

present results for the matrix hard, using three different values of
the drop tolerance in the AINV algorithm. It can be seen that changing the
value of - changes the density of the preconditioner and the number of itera-
tions. However, the total timings are scarcely affected, especially if at least 8
processors are being used. See [8] for a similar observation in a different con-
text. It is also clear from these runs that good speed-ups are obtained so long as
the size of the separator set is small compared to the average subdomain size.
As soon as the separator set is comparable in size to the average subdomain or
larger, the sequential bottleneck represented by the Schur complement part of
the computation begins to dominate and performance deteriorates. The number
of iterations remains roughly constant (with a slight downward trend) as
the number of processors grows. This is due to the influence of the approximate
Schur complement.
The same problem was also solved using Bi-CGSTAB with diagonal precondi-
tioning. This required approximately 700 iterations and 16.4 seconds on one
processor. If implemented in parallel, this method would probably give results
only slightly worse than those obtained with AINV. A similar observation
applies to matrices qn and mutex. On the other hand, diagonally preconditioned
Bi-CGSTAB did not converge on the telecom problem. Hence, AINV
is a more robust approach. Furthermore, the ability to reduce the number of
iterations, and therefore the total number of inner products, is an advantage
on distributed memory machines, on which inner products incur an additional
penalty due to the need for global communication.
Matrix hard,
P-time 2.35 1.19 0.65 0.40 0.34
P-density 6.21 5.90 5.67 5.14 4.52
Its
It-time 9.43 2.91 1.42 0.99 0.88
Tot-time 11.8 4.10 2.07 1.39 1.22
Sep-size 156 321 540 900 1346
Avg-dom 10073 5155 2470 1213 592

Table
Matrix hard,
P-time 1.19 0.60 0.35 0.22 0.21
P-density 3.10 3.02 2.95 2.78 2.52
Its 106 109 99 98 97
It-time 6.85 2.90 1.59 1.05 1.19
Tot-time 8.04 3.50 1.94 1.27 1.40

Table
Matrix hard,
P-time 0.73 0.38 0.22 0.15 0.14
P-density 1.39 1.36 1.34 1.28 1.19
Its 170 167 159 151 153
It-time 6.03 2.52 1.50 1.30 1.64
Tot-time 6.76 2.90 1.72 1.45 1.78
Results for matrices leaky and 2d are reported in Tables 5 and 6. These two
matrices are rather small, so only up to 8 processors were used. Note that the
speed-ups are better for 2d than for leaky. Also notice that the preconditioner
is very sparse for leaky, but rather dense for 2d.

Tables

7 and 8 refer to the telecom test problem. Here we found that very
small values of - (and, consequently, very dense preconditioners) are necessary
in order to achieve convergence in a reasonable number of iterations. This
problem is completely different from the matrices arising from the solution of
elliptic partial differential equations. Notice the fairly small size of the sep-
Matrix leaky,
P-time 0.26 0.17 0.09
P-density
No. its 134 134 132
It-time 1.39 0.84 0.62
Tot-time 1.65 1.01 0.71
Sep-size 48 144 335
Avg-dom 4105 2028 990

Table
Matrix 2d,
P-time 0.38 0.16 0.09
P-density 8.60 7.12 6.54
No. its 33 36 37
It-time 1.46 0.56 0.38
Tot-time 1.84 0.72 0.47
Sep-size 129 308 491
Avg-dom 8256 4083 2018

Table
Matrix telecom,
P-time 24.9 10.0 3.35 1.21 1.05 1.37
P-density 167 116 70 44 28
Its 11 12 14 13 12 12
It-time 36.5 14.0 6.04 0.96 0.38 0.43
Tot-time 61.4 24.0 9.39 2.17 1.43 1.80
Sep-size 34 97 220 471 989 1603
Avg-dom 10229 5099 2534 1251 609 295
arator set, which causes the density of the preconditioner to decrease very
fast as the number of processors (and corresponding subdomains) grows. As
a result, speed-ups are quite good (even superlinear) up to 32 processors. For
a sufficiently high number of processors, the density of the preconditioner be-
Matrix telecom,
P-time 7.04 3.78 1.37 0.61 0.49 0.56
P-density
No. its
It-time 59.4 27.0 12.8 2.18 1.36 1.71
Tot-time 66.4 30.8 14.2 2.79 1.86 2.27

Table
Matrix
P-time 1.42 0.69 0.31
P-density 4.13 2.65 1.90
No. its 292 288 285
It-time 17.0 8.45 6.38
Tot-time 18.4 9.14 6.69
Sep-size 3911 6521 12932
Avg-dom 9758 4226 1379

Table
Matrix
P-time 1.19 0.40 0.10
P-density 0.14 0.14 0.14
No. its
It-time 1.64 1.19 1.51
Tot-time 2.83 1.59 1.61
Sep-size 13476 17749 20654
Avg-dom 12864 5363 2319
comes acceptable, and the convergence rate is the same or comparable to that
obtained with a very dense preconditioner on a small number of processors.

Tables

results for matrices ncd and mutex, respectively. For the
first matrix we see that the separator set is larger than the average subdomain
already for nevertheless, it is possible to use effectively up
to 8 processors. Matrix mutex exhibits a behavior that is radically different
Matrix qn,
P-time 4.12 2.37 2.26 2.82
P-density 1.27 1.23 1.18 1.14
No. its
It-time 13.4 6.33 4.58 5.68
Tot-time 17.5 8.70 6.84 8.50
Sep-size 2879 6579 13316 20261
Avg-dom 50873 24511 11414 5273
from that of matrices arising from PDE's in two or three space dimensions.
The separator set is huge already for 2. This is due to the fact that the
problem has a state space (graph) of high dimensionality, leading to a very
unfavorable surface-to-volume ratio in the graph partitioning. In order to solve
this problem, we had to use two different values of - in the two levels of AINV;
at the subdomain level we used when forming the approximate
Schur complement we dropped everything outside the main diagonal, resulting
in a diagonal -
S. In spite of this, convergence was very rapid. Nevertheless, it
does not pay to use more than processors.
In

Table

11 we report results with the largest example in our data set, qn.
This model consists of a network of three queues, and is analogous to a three-dimensional
problem. Because of the fairly rapid growth of the separator
set, it does not pay to use more than processors.
The test problems considered so far, although realistic, are relatively small.
Hence, it is difficult to make efficient use of more than 16 processors, with
the partial exceptions of matrices hard and telecom. To test the scalability of
the proposed solver on larger problems, we generated some simple reliability
problems analogous to those used in [2] and [5]; see also [26], page 135. These
problems have a closed form solution. In Table 12 we show timing results for
running 100 preconditioned Bi-CGSTAB iterations on a reliability problem of
entries. This problem is sufficiently
large to show the good scalability of the algorithm up to processors.
We conclude this section on numerical experiments by noting that in virtually
all the runs, the preconditioner construction time has been quite modest and
the total solution time has been dominated by the cost of the iterative phase.
Reliability model,
P-time 9.53 4.86 2.49 1.31 0.90 0.90
P-density 4.05 4.02 3.97 3.90 3.80 3.70
It-time 138.8 70.5 37.2 16.8 9.47 7.96
Tot-time 148.3 75.4 39.7 18.1 10.4 8.86
Sep-size 542 1229 2186 3268 4919 7336
Avg-dom 124729 62193 30977 15421 7659 3792
Conclusions
We have investigated the use of a parallel preconditioner for Krylov subspace
methods in the context of Markov chain problems. The preconditioner is a
direct approximation, in factorized form, of a (1; 2)-inverse of the generator
matrix A, and is based on an A-biorthogonalization process. Parallelization
is achieved through graph partitioning, although other approaches are also
possible. The existence of the preconditioner has been justified theoretically,
and numerical experiments on a parallel computer have been carried out in
order to assess the effectiveness and scalability of the proposed technique.
The numerical tests indicate that the preconditioner construction costs are
modest, and that good scalability is possible provided that the amount of
work per processor is sufficiently large compared to the size of the separator
set.
The method appears to be well suited for problems in which the generator
can be explicitly formed and stored. Parallelization based
on graph partitioning is usually effective, with the possible exception of problems
with a state space of high dimensionality (i.e., a large state descriptor
set). For such problems, a different parallelization strategy is needed in order
to achieve scalability of the implementation.

Acknowledgements

. We would like to thank Professors M. Gutknecht and
W. Sch-onauer for their kind invitation to take part in this commemoration of
our friend and colleague R-udiger Weiss. We are indebted to Tu-grul Dayar for
providing the test matrices used in the numerical experiments and for useful
information about these problems, as well as for his comments on an early
version of the paper. Thanks also to Carl Meyer for his valuable input on
generalized inverses.



--R


The arithmetic mean method for finding the stationary vector of Markov chains

A sparse approximate inverse preconditioner for the conjugate gradient method
A parallel block projection method of the Cimmino type for finite Markov chains
A sparse approximate inverse preconditioner for nonsymmetric linear systems
A comparative study of sparse approximate inverse preconditioners
Approximate inverse preconditioning in the parallel solution of sparse eigenproblems
Nonnegative Matrices in the Mathematical Sciences (Academic Press
Distributed steady state analysis using Kronecker algebra
Incomplete factorization of singular M-matrices
Generalized Inverses of Linear Transformations (Pitman Publishing Ltd.
Fast wavelet iterative solvers applied to the Neumann problem
A priory sparsity patterns for parallel sparse approximate inverse preconditioners
Comparison of partitioning techniques for two-level iterative solvers on large
Parallel preconditioning with sparse approximate inverses
Asynchronous iterations for the solution of Markov systems
A fast and high quality multilevel scheme for partitioning irregular graphs
Distributed disk-based solution techniques for large Markov models
An iterative solution method for linear systems of which the coefficient matrix is a symmetric M-matrix
The role of the group generalized inverse in the theory of finite Markov chains
the theory of nearly reducible systems
Matrix Analysis and Applied Linear Algebra (SIAM
Experimental studies of parallel iterative solutions of Markov chains with block partitions
Preconditioned Krylov subspace methods for the numerical solution of Markov chains
Introduction to the Numerical Solution of Markov Chains (Princeton University Press
BiCGSTAB: A fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems
--TR
Incomplete factorization of singular M-matrices
Stochastic complementation, uncoupling Markov chains, and the theory of nearly reducible systems
BI-CGSTAB: a fast and smoothly converging variant of BI-CG for the solution of nonsymmetric linear systems
Iterative solution methods
A Sparse Approximate Inverse Preconditioner for the Conjugate Gradient Method
Parallel Preconditioning with Sparse Approximate Inverses
A Sparse Approximate Inverse Preconditioner for Nonsymmetric Linear Systems
A Fast and High Quality Multilevel Scheme for Partitioning Irregular Graphs
A comparative study of sparse approximate inverse preconditioners
Matrix analysis and applied linear algebra
Comparison of Partitioning Techniques for Two-Level Iterative Solvers on Large, Sparse Markov Chains
A Priori Sparsity Patterns for Parallel Sparse Approximate Inverse Preconditioners

--CTR
Aliakbar Montazer Haghighi , Dimitar P. Mishev, A parallel priority queueing system with finite buffers, Journal of Parallel and Distributed Computing, v.66 n.3, p.379-392, March 2006
Ilias G. Maglogiannis , Elias P. Zafiropoulos , Agapios N. Platis , George A. Gravvanis, Computing the success factors in consistent acquisition and recognition of objects in color digital images by explicit preconditioning, The Journal of Supercomputing, v.30 n.2, p.179-198, November 2004
Nicholas J. Dingle , Peter G. Harrison , William J. Knottenbelt, Uniformization and hypergraph partitioning for the distributed computation of response time densities in very large Markov models, Journal of Parallel and Distributed Computing, v.64 n.8, p.908-920, August 2004
Michele Benzi, Preconditioning techniques for large linear systems: a survey, Journal of Computational Physics, v.182 n.2, p.418-477, November 2002
singular matrices;generalized inverses;Bi-CGSTAB;graph partitioning;AINV;parallel preconditioning;discrete Markov chains;iterative methods