A Sparse Approximate Inverse Preconditioner for Nonsymmetric Linear Systems.
This paper is concerned with a new approach to preconditioning for large, sparse linear systems. A procedure for computing an incomplete factorization of the inverse of a nonsymmetric matrix is developed, and the resulting factorized sparse approximate inverse is used as an explicit preconditioner for conjugate gradient--type methods. Some theoretical properties of the preconditioner are discussed, and numerical experiments on test matrices from the Harwell--Boeing collection and from Tim Davis's collection are presented. Our results indicate that the new preconditioner is cheaper to construct than other approximate inverse preconditioners. Furthermore, the new technique insures convergence rates of the preconditioned iteration which are comparable with those obtained with standard implicit preconditioners.
Introduction
. In this paper we consider the solution of nonsingular linear systems
of the form
(1)
where the coefficient matrix A 2 IR n\Thetan is large and sparse. In particular, we are concerned
with the development of preconditioners for conjugate gradient-type methods. It is well-known
that the rate of convergence of such methods for solving (1) is strongly influenced by
the spectral properties of A. It is therefore natural to try to transform the original system
into one having the same solution but more favorable spectral properties. A preconditioner
is a matrix that can be used to accomplish such a transformation. If G is a nonsingular
Dipartimento di Matematica, Universit'a di Bologna, Italy and CERFACS, 42 Ave. G. Coriolis, 31057
Toulouse Cedex, France (benzi@cerfacs.fr). This work was supported in part by a grant under the scientific
cooperation agreement between the CNR and the Czech Academy of Sciences.
y Institute of Computer Science, Academy of Sciences of the Czech Republic, Pod vod'arenskou v-e-z'i 2,
(tuma@uivt.cas.cz). The work of this author was supported in
part by grants GA CR No. 201/93/0067 and GA AS CR No. 230401 and by NSF under grant number
INT-9218024.
Michele Benzi and Miroslav T-uma
matrix which approximates A linear system
(2)
will have the same solution as system (1) but the convergence rate of iterative methods
applied to (2) may be much higher. Problem (2) is preconditioned from the left, but right
preconditioning is also possible. Preconditioning on the right leads to the transformed linear
system
Once the solution y of (3) has been obtained, the solution of (1) is given by
The choice between left or right preconditioning is often dictated by the choice of the
iterative method. It is also possible to use both forms of preconditioning at once (split
preconditioning), see [3] for further details.
Note that in practice it is not required to compute the matrix product GA (or AG)
explicitly, because conjugate gradient-type methods only necessitate the coefficient matrix
in the form of matrix-vector multiplies. Therefore, applying the preconditioner within a
step of a gradient-type method reduces to computing the action of G on a vector.
Loosely speaking, the closer G is to the exact inverse of A, the higher the rate of
convergence of iterative methods will be. Choosing yields convergence in one
step, but of course constructing such a preconditioner is equivalent to solving the original
problem. In practice, the preconditioner G should be easily computed and applied, so that
the total time for the preconditioned iteration is less than the time for the unpreconditioned
one. Typically, the cost of applying the preconditioner at each iteration of a conjugate
gradient-type method should be of the same order as the cost of a matrix-vector multiply
involving A. For a sparse A, this implies that the preconditioner should also be sparse with
a density of nonzeros roughly of the same order as that of A.
Clearly, the effectiveness of a preconditioning strategy is strongly problem and architecture
dependent. For instance, a preconditioner which is expensive to compute may become
viable if it is to be reused many times, since in this case the initial cost of forming the
preconditioner can be amortized over several linear systems. This situation occurs, for in-
stance, when dealing with time-dependent or nonlinear problems, whose numerical solution
gives rise to long sequences of linear systems having the same coefficient matrix (or a slowly
varying one) and different right-hand sides. Furthermore, preconditioners that are very efficient
in a scalar computing environment may show poor performance on vector and parallel
machines, and conversely.
Approximate Inverse Preconditioning 3
A number of preconditioning techniques have been proposed in the literature (see, e.g.,
[2],[3] and the references therein). While it is generally agreed that the construction of efficient
general-purpose preconditioners is not possible, there is still considerable interest in
developing methods which will perform well on a wide range of problems and are well-suited
for state-of-the-art computer architectures. Here we introduce a new algebraic preconditioner
based on an incomplete triangular factorization of A \Gamma1 . This paper is the natural
continuation of [8], where the focus was restricted to symmetric positive definite systems
and to the preconditioned conjugate gradient method (see also [5],[7]).
The paper is organized as follows. In x2 we give a quick overview of implicit and explicit
preconditioning techniques, considering the relative advantages as well as the limitations of
the two approaches. In x3 we summarize some recent work on the most popular approach
to approximate inverse preconditioning, based on Frobenius norm minimization. In x4 we
introduce the new incomplete inverse triangular decomposition technique and describe some
of its theoretical properties. A graph-theoretical characterization of fill-in in the inverse
triangular factorization is presented in x5. In x6 we consider the use of preconditioning
on matrices which have been reduced to block triangular form. Implementation details
and the results of numerical experiments are discussed in xx7 and 8, and some concluding
remarks and indications for future work are given in x9. Our experiments suggest that the
new preconditioner is cheaper to construct than preconditioners based on the optimization
approach. Moreover, good rates of convergence can be achieved by our preconditioner,
comparable with those insured by standard ILU-type techniques.
2. Explicit vs. implicit preconditioning. Most existing preconditioners can be
broadly classified as being either of the implicit or of the explicit kind. A preconditioner
is implicit if its application, within each step of the chosen iterative method, requires the
solution of a linear system. A nonsingular matrix M - A implicitly defines an approximate
applying G requires solving a linear system with coefficient
matrix M . Of course, M should be chosen so that solving a system with matrix M is easier
than solving the original problem (1). Perhaps the most important example is provided by
preconditioners based on an Incomplete LU (ILU) decomposition. Here
U where
L and -
U are sparse triangular matrices which approximate the exact L and U factors of
A. Applying the preconditioner requires the solution of two sparse triangular systems (the
forward and backward solves). Other notable examples of implicit preconditioners include
the ILQ, SSOR and ADI preconditioners, see [3].
4 Michele Benzi and Miroslav T-uma
In contrast, with explicit preconditioning a matrix G - A \Gamma1 is known (possibly as the
product of sparse matrices) and the preconditioning operation reduces to forming one (or
more) matrix-vector product. For instance, many polynomial preconditioners belong to this
class [37]. Other explicit preconditioners will be described in the subsequent sections.
Implicit preconditioners have been intensively studied, and they have been successfully
employed in a number of applications. In spite of this, in the last few years an increasing
amount of attention has been devoted to alternative forms of preconditioning, especially of
the explicit kind. There have been so far two main reasons for this recent trend. In the first
place, shortly after the usage of modern high-performance architectures became widespread,
it was realized that straightforward implementation of implicit preconditioning in conjugate
gradient-like methods could lead to severe degradation of the performance on the new
machines. In particular, the sparse triangular solves involved in ILU-type preconditioning
were found to be a serial bottleneck (due to the recursive nature of the computation), thus
limiting the effectiveness of this approach on vector and parallel computers. It should be
mentioned that considerable effort has been devoted to overcoming this difficulty. As a
result, for some architectures and types of problems it is possible to introduce nontrivial
parallelism and to achieve reasonably good performance in the triangular solves by means
of suitable reordering strategies (see, e.g., [1],[38],[54]). However, the triangular solves
remain the most problematic aspect of the computation, both on shared memory [33] and
distributed memory [10] computers, and for many problems the efficient application of an
implicit preconditioner in a parallel environment still represents a serious challenge.
Another drawback of implicit preconditioners of the ILU-type is the possibility of break-downs
during the incomplete factorization process, due to the occurrence of zero or exceedingly
small pivots. This situation typically arises when dealing with matrices which are
strongly unsymmetric and/or indefinite, even if pivoting is applied (see [11],[49]), and in
general it may even occur for definite problems unless A exhibits some degree of diagonal
dominance. Of course, it is always possible to safeguard the incomplete factorization process
so that it always runs to completion, producing a nonsingular preconditioner, but there
is also no guarantee that the resulting preconditioner will be of acceptable quality. Fur-
thermore, as shown in [23], there are problems for which standard ILU techniques produce
unstable incomplete factors, resulting in useless preconditioners.
Explicit preconditioning techniques, based on directly approximating A \Gamma1 , have been
developed in an attempt to avoid or mitigate such difficulties. Applying an explicit preconditioner
only requires sparse matrix-vector products, which should be easier to parallelize
Approximate Inverse Preconditioning 5
than the sparse triangular solves, and in some cases the construction of the preconditioner
itself is well-suited for parallel implementation. In addition, the construction of an approximate
inverse is sometimes possible even if the matrix does not have a stable incomplete LU
decomposition. Moreover, we mention that sparse incomplete inverses are often used when
constructing approximate Schur complements (pivot blocks) for use in incomplete block
factorization and other two-level preconditioners, see [2],[3],[12],[15].
Of course, explicit preconditioners are far from being completely trouble-free. Even if
a sparse approximate inverse G is computed, care must be exercised to ensure that G is
nonsingular. For nonsymmetric problems, the same matrix G could be a good approximate
inverse if used for left preconditioning and a poor one if used for right preconditioning,
see [36, p. 96],[45, p. 66],[48]. Furthermore, explicit preconditioners are sometimes not
as effective as implicit ones at reducing the number of iterations, in the sense that there
are problems for which they require a higher number of nonzeros in order to achieve the
same rate of convergence insured by implicit preconditioners. One of the reasons for this
limitation is that an explicit preconditioner attempts to approximate A \Gamma1 , which is usually
dense, with a sparse matrix. Thus, an explicit preconditioner is more likely to work well if
A \Gamma1 contains many entries which are small (in magnitude). A favorable situation is when A
exhibits some form of diagonal dominance, but for such problems implicit preconditioning
is also likely to be very effective. Hence, for problems of this type, explicit preconditioners
can be competitive with implicit ones only if explicitness is fully exploited. Finally, explicit
preconditioners are usually more expensive to compute than implicit ones, although this
difference may become negligible in the common situation where several linear systems with
the same coefficient matrix and different right-hand sides have to be solved. In this case
the time for computing the preconditioner is often only a fraction of the time required for
the overall computation. It is also worth repeating that the construction of certain sparse
approximate inverses can be done, at least in principle, in a highly parallel manner, whereas
the scope for parallelism in the construction of ILU-type preconditioners is more limited.
3. Methods based on Frobenius norm minimization. A good deal of work has
been devoted to explicit preconditioning based on the following approach: the sparse approximate
inverse is computed as the matrix G which minimizes kI \Gamma GAk (or kI \Gamma AGk
for right preconditioning) subject to some sparsity constraint (see [4], Ch. 8 of [2],[16],[43],
[44],[32],[31],[11],[30]). Here the matrix norm is usually the Frobenius norm or a weighted
variant of it, for computational reasons. With this choice, the constrained minimization
problem decouples into n independent linear least squares problems (one for each row, or
6 Michele Benzi and Miroslav T-uma
column of G), the number of unknowns for each problem being equal to the number of
nonzeros allowed in each row (or column) of G. This immediately follows from the identity
is the ith unit vector and g i is the ith column of G. Clearly, there is considerable
scope for parallelism in this approach. The resulting sparse least squares problems can
be solved, in principle, independently of each other, either by direct methods (as in [44],
[31],[30]) or iteratively ([11],[42]).
In early papers (e.g. [4],[32],[43]) the sparsity constraint was imposed a priori, and the
minimizer was found relative to a class of matrices with a predetermined sparsity pattern.
For instance, when A is a band matrix with a good degree of diagonal dominance, a banded
approximation to A \Gamma1 is justified, see [18]. However, for general sparse matrices it is very
difficult to guess a good sparsity pattern for an approximate inverse, and several recent
papers have addressed the problem of adaptively defining the nonzero pattern of G in order
to capture "large" entries of the inverse [31],[30]. Indeed, by monitoring the size of each
residual it is possible to decide whether new entries of g i are to be retained or
discarded, in a dynamic fashion. Moreover, the information on the residuals can be utilized
to derive rigorous bounds on the clustering of the singular values of the preconditioned
matrix and therefore to estimate its condition number [31]. It is also possible to formulate
conditions on the norm of the residuals which insure that the approximate inverse will
be nonsingular. Unfortunately, such conditions appear to be of dubious practical value,
because trying to fulfill them could lead to a very dense approximate inverse [16],[11].
A disadvantage of this approach is that symmetry in the coefficient matrix cannot be
exploited. If A is symmetric positive definite (SPD), the sparse approximate inverse will
not be symmetric in general. Even if a preset, symmetric sparsity pattern is enforced, there
is no guarantee that the approximate inverse will be positive definite. This could lead to a
breakdown in the conjugate gradient acceleration. For this reason, Kolotilina and Yeremin
[43],[44] propose to compute an explicit preconditioner of the form
is lower triangular. The preconditioned matrix is then GLAG T
, which is SPD, and the
conjugate gradient method can be applied. The matrix GL is the solution of a constrained
minimization problem for the Frobenius norm of I \Gamma LGL where L is the Cholesky factor
of A. In [43] it is shown how this problem can be solved without explicit knowledge of
any of the entries of L, using only entries of the coefficient matrix A. The same technique
can also be used to compute a factorized approximate inverse of a nonsymmetric matrix by
Approximate Inverse Preconditioning 7
separately approximating the inverses of the L and U factors. As it stands, however, this
technique requires that the sparsity pattern of the approximate inverse triangular factors
be specified in advance, and therefore is not suitable for matrices with a general sparsity
pattern.
There are additional reasons for considering factorized approximate inverses. Clearly,
with the approximate inverse G expressed as the product of two triangular factors it is
trivial to insure that G is nonsingular. Another argument in favor of this approach is given
in [11], where it is observed that factorized forms of general sparse matrices contain more
information for the same storage than if a single product was stored.
The serial cost for the construction of this type of preconditioner is usually very high,
although the theoretical parallel complexity can be quite moderate [44],[30]. The results of
numerical experiments reported in [44] demonstrate that factorized sparse approximate inverse
preconditioners can insure rapid convergence of the preconditioned conjugate gradient
iteration when applied to certain finite element discretizations of 3D PDE problems
arising in elasticity theory. However, in these experiments the preconditioning strategy is
not applied to the coefficient matrix directly, but rather to a reduced system (Schur comple-
ment) which is better conditioned and considerably less sparse than the original problem.
When the approximate inverse preconditioner is applied directly to the original stiffness
matrix A, the rate of convergence of the PCG iteration is rather disappointing.
A comparison between a Frobenius norm-based sparse approximate inverse preconditioner
and the ILU(0) preconditioner on a number of general sparse matrices has been made
in [30]. The reported results show that the explicit preconditioner can insure rates of convergence
comparable with those achieved with the implicit ILU-type approach. Again, it is
observed that the construction of the approximate inverse is often very costly, but amenable
to parallelization.
Factorized sparse approximate inverses have also been considered by other authors, for
instance by Kaporin [39],[40],[41], whose approach is also based on minimizing a certain matrix
functional and is closely related to that of Kolotilina and Yeremin. In the next sections
we present an alternative approach to factorized sparse approximate inverse preconditioning
which is not grounded in optimization, but is based instead on a direct method of matrix
inversion. As we shall see, the serial cost of forming a sparse approximate inverse with this
technique is usually much less than with the optimization approach, while the convergence
rates are still comparable, on average, with those obtained with ILU-type preconditioning.
8 Michele Benzi and Miroslav T-uma
4. A method based on inverse triangular factorization. The optimization approach
to constructing approximate inverses is not the only possible one. In this section we
consider an alternative procedure based on a direct method of matrix inversion, performed
incompletely in order to preserve sparsity. This results in a factorized sparse G - A \Gamma1 .
Being an incomplete matrix factorization method, our procedure resembles classical ILU-
type implicit techniques, and indeed we can draw from the experience accumulated in years
of use of ILU-type preconditioning both at the implementation stage and when deriving
theoretical properties of the preconditioner G. This paper continues the work in [8], where
the symmetric positive definite case was studied (see also [5],[7]).
The construction of our preconditioner is based on an algorithm which computes two
sets of vectors fz i g n
, which are A-biconjugate, i.e. such that w T
only if i 6= j. Given a nonsingular matrix A 2 IR n\Thetan , there is a close relationship between
the problem of inverting A and that of computing two sets of A-biconjugate vectors fz i g n
and fw i g n
. If
is the matrix whose ith column is z i and
is the matrix whose ith column is w i , then
It follows that W and Z are necessarily nonsingular and
Hence, the inverse of A is known if two complete sets of A-biconjugate vectors are known.
Note that there are infinitely many such sets. Matrices W and Z whose columns are A-
biconjugate can be explicitly computed by means of a biconjugation process applied to the
columns of any two nonsingular matrices W (0) , Z (0) 2 IR n\Thetan . A computationally convenient
choice is to let W the biconjugation process is applied to the unit basis
vectors. In order to describe the procedure, let a T
and c T
denote the ith row of A and A T ,
Approximate Inverse Preconditioning 9
respectively (i.e., c i is the ith column of A). The basic A-biconjugation procedure can be
written as follows.
THE BICONJUGATION ALGORITHM
:= a T
z (i)
z (i\Gamma1)
and
Return
This algorithm is essentially due to L. Fox, see Ch. 6 of [25]. Closely related methods
have also been considered by Hestenes and Stiefel [35, pp. 426-427],[34] and by Stewart [52].
A more general treatment is given in the recent paper [14]. Geometrically, the procedure
can be regarded as a generalized Gram-Schmidt orthogonalization with oblique projections
and nonstandard inner products, see [6],[14].
Several observations regarding this algorithm are in order. In the first place we note
that the above formulation contains some redundancy, since in exact arithmetic
Therefore, at step i the computation of the dot product q (i\Gamma1)
i may be replaced
by the assignment q (i\Gamma1)
. Another observation is the fact that the procedure, as
it stands, is vulnerable to breakdown (division by zero), which occurs whenever any of the
Michele Benzi and Miroslav T-uma
quantities p (i\Gamma1)
happens to be zero. It can be shown that in exact arithmetic, the
biconjugation algorithm will not break down if and only if all the leading principal minors of
A are nonzero (see below). For any nonsingular matrix A there exists a permutation matrix
(or Q) such that the procedure applied to PA (or to AQ) will not break down. As in
the LU decomposition with pivoting, such permutation matrices represent row (or column)
interchanges on A which can be performed, if needed, in the course of the computation.
If the biconjugation process can be carried to completion without interchanges, the
resulting Z and W matrices are upper triangular, 1 they both have all diagonal entries equal
to one, and satisfy the identity
We recognize in (5) the familiar LDU decomposition A = LDU , where L is unit lower
triangular, U is unit upper triangular and D is the diagonal matrix with the pivots down
the main diagonal. Because this factorization is unique, we have that the biconjugation
algorithm explicitly computes
and the matrix D, which is exactly the same in (5) and in A = LDU . Hence, the process
produces an inverse triangular decomposition of A or, equivalently, a triangular decomposition
(of the UDL type) of A \Gamma1 . The p i 's returned by the algorithm are the pivots in the
LDU factorization of A. If we denote by \Delta i the ith leading principal minor of A (1 - i - n)
and let the identity (5) implies that
showing that the biconjugation algorithm can be performed without breakdowns if and only
if all leading principal minors of A are non-vanishing. In finite precision arithmetic, pivoting
may be required to promote numerical stability.
Once Z, W and D are available, the solution of a linear system of the form (1) can be
computed, by (4), as
1 Note that this is not necessarily true when a matrix other than the identity is used at the outset, i.e. if
Approximate Inverse Preconditioning 11
In practice, this direct method for solving linear systems is not used on account of its cost:
for a dense n \Theta n matrix, the biconjugation process requires about twice the work as the
LU factorization of A. Notice that the cost of the solve phase using (6) is the same as for
the forward and backward solves with the L and U factors.
If A is symmetric, the number of operations in the biconjugation algorithm can be
halved by observing that W must equal Z. Hence, the process can be carried out using
only the rows of A, the z-vectors and the p (i\Gamma1)
. The columns of the resulting Z form a set
of conjugate directions for A. If A is SPD, no breakdown can occur (in exact arithmetic),
so that pivoting is not required and the algorithm computes the L T DL factorization of
A \Gamma1 . This method was first described in [26]. Geometrically, it amounts to Gram-Schmidt
orthogonalization with inner product hx; yi := x T Ay applied to the unit vectors e
It is sometimes referred to as the conjugate Gram-Schmidt process . The method is still
impractical as a direct solver because it requires about twice the work of Cholesky for dense
matrices. However, as explained in [5] and [6], the same algorithm can also be applied to
nonsymmetric systems, resulting in an implicit LDU factorization where only
are computed. Indeed, it is possible to compute a solution to (1) for any right-hand
side using just Z, D and part of the entries of A. This method has the same arithmetic
complexity as Gaussian elimination when applied to dense problems. When combined with
suitable sparsity-preserving strategies the method can be useful as a sparse direct solver, at
least for some types of problems (see [5],[6]).
For a sparse symmetric and positive definite A, the Z matrix produced by the algorithm
tends to be dense (see the next section), but it can be observed experimentally that very
often, most of the entries in Z have small magnitude. If fill-in in the Z matrix is reduced by
removing suitably small entries in the computation of the z-vectors, the algorithm computes
a sparse matrix -
Z and a diagonal matrix -
D such that
(i.e., a factorized sparse approximate inverse of A). Hence, G can be used as an explicit
preconditioner for the conjugate gradient method. A detailed study of this preconditioning
strategy for SPD problems can be found in [8], where it is proven that the incomplete
inverse factorization exists if A is an H-matrix (analogously to ILU-type factorizations).
The numerical experiments in [8] show that this approach can insure fast convergence of the
PCG iteration, almost as good as with implicit preconditioning of the incomplete Cholesky
type. The construction of the preconditioner itself, while somewhat more expensive than the
computation of the incomplete Cholesky factorization, is still quite cheap. This is in contrast
Michele Benzi and Miroslav T-uma
with the least squares approach described in the previous section, where the construction of
the approximate inverse is usually very time consuming, at least in a sequential environment.
In the remainder of this paper we consider an explicit preconditioning strategy based
on the biconjugation process described above. Sparsity in the Z and W factors of A \Gamma1 is
preserved by removing "small" fill in the z- and w-vectors. A possibility would be to drop
all newly added fill-in elements outside of a preset sparsity pattern above the main diagonal
in Z and W ; however, for general sparse matrices it is very difficult to guess a reasonable
sparsity pattern, and a drop tolerance is used instead. A trivial extension of the results
in [8] shows that the incomplete biconjugation process (incomplete inverse factorization)
cannot break down, in exact arithmetic, if A is an H-matrix. For more general matrices
it is necessary to safeguard the computation in order to avoid breakdowns. This requires
pivot modifications and perhaps some form of pivoting -we postpone the details to x7.
The incomplete biconjugation algorithm computes sparse unit upper triangular matrices
W - W and a nonsingular diagonal matrix -
D - D such that
is a factorized sparse approximate inverse of A which can be used as an explicit preconditioner
for conjugate gradient-type methods for the solution of (1).
We conclude this section with a few remarks on properties of the approximate inverse
preconditioner G just described. If A is not an H-matrix, as already mentioned, the construction
of the preconditioner could break down due to the occurrence of zero or extremely
small pivots. However, following [46], we note that there always exists ff ? 0 such that
ffI is diagonally dominant, and hence an H-matrix. Therefore, if the incomplete bicon-
jugation algorithm breaks down, one could try to select ff ? 0 and re-attempt the process
on the shifted matrix A should be large enough to insure the existence
of the incomplete inverse factorization, but also small enough so that A 0 is close to A. This
approach has several drawbacks: for ill-conditioned matrices, the quality of the resulting
preconditioner is typically poor; furthermore, the breakdown that prompts the shift may
occur near the end of the biconjugation process, and the preconditioner may have to be
recomputed several times before a satisfactory value of ff is found. A better strategy is to
perform diagonal modifications only as the need arises, shifting pivots away from zero if
their magnitude is less than a specified threshold (see x7 for details).
If A is an M-matrix, it follows from the results in [8] that G is a nonnegative matrix.
Approximate Inverse Preconditioning 13
Moreover, it is easy to see that componentwise the following inequalities hold:
where DA is the diagonal part of A. Furthermore, if G 1 and G 2 are two approximate
inverses of the M-matrix A produced by the incomplete biconjugation process and the drop
tolerance used for G 1 is greater than or equal to the drop tolerance used for G 2 , then
The same is true if sparsity patterns are used to determine the nonzero structure in -
Z and
W and the patterns for G 2 include the patterns for G 1 . This monotonicity property is
shared by other sparse approximate inverses, see for example Ch. 8 in [2]. We note that
property (7) is important if the approximate inverse is to be used within an incomplete
block factorization of an M-matrix A, because it insures that all the intermediate matrices
produced in the course of the incomplete factorization preserve the M-matrix property (see
[2, pp. 263-264]).
Finally, after discussing the similarities, we point to a difference between our incomplete
inverse factorization and the ILU-type factorization of a matrix. The incomplete factorization
of an M-matrix A induces a splitting
which is a regular splitting, and
therefore convergent: ae(I \Gamma -
denotes the spectral radius of a
[47],[55]). The same is not true, in general, for our incomplete factorization.
If one considers the induced splitting splitting need
not be convergent. An example is given by the symmetric M-matrix
A =B @
For this matrix, the incomplete inverse factorization with a drop tolerance
intermediate fill-in is dropped if smaller than T in absolute value) produces an approximate
inverse G such that ae(I \Gamma GA) - 1:215 ? 1. This shows that the approximate decomposition
cannot be obtained, in general, from an incomplete factorization of A. In this sense, the
incomplete inverse factorization is not algebraically equivalent to an incomplete LDU factorization
performed on A.
14 Michele Benzi and Miroslav T-uma
5. Fill-in in the biconjugation algorithm. In this section we give a characterization
of the fill-in occurring in the factorized inverse obtained by the biconjugation algorithm.
These results may serve as a guideline to predict the structure of the factorized approximate
inverse, and have an impact on certain aspects of the implementation.
It is well-known that structural nonzeros in the inverse matrix A \Gamma1 can be characterized
by the paths in the graph of the original matrix A (see [24],[29]). The following lemma states
necessary and sufficient conditions for a new entry (fill-in) to be added in one of the z-vectors
at the ith step of the biconjugation algorithm. A similar result holds for the w-vectors. We
make use of the standard no-cancellation assumption.
Lemma 5.1. Let 1
z (i\Gamma1)
if and only if l - i, z (i\Gamma1)
li and, at the same time, at least one of the two following
conditions holds:
Proof. Suppose that z (i\Gamma1)
0: Directly from the update formula for the
z-vectors we see that z (i\Gamma1)
li 6= 0 and l - i, since z (i\Gamma1)
lj becomes
nonzero in the ith step then clearly p (i\Gamma1)
j must be nonzero. But
z (i\Gamma1)
kj a ik
and we get the result. The opposite implication is trivial. 2

Figures

5.1 through 5.6 provide an illustration of the previous lemma. Figure 5.1 shows
the nonzero structure of the matrix FS760 1 of order 760 from the Harwell-Boeing
collection [21]. Figures 5.2-6 show the structure of the factor Z at different stages of the
biconjugation algorithm. These pictures show that in the initial steps, when most of the
entries of Z are still zero, the nonzeros in Z are induced by nonzeros in the corresponding
positions of A. A similar situation occurs, of course, for the process which computes W .
In

Figure

5.7 we show the entries of Z which are larger (in absolute
in

Figure

5.8 we show the incomplete factor -
Z obtained with drop tolerance It
can be seen how well the incomplete process is able to capture the "large" entries in the
complete factor Z. The figures were generated using the routines for plotting sparse matrix
patterns from SPARSKIT [50].
Approximate Inverse Preconditioning 15

Figure

5.1-2: Structure of the matrix FS760 1 (left) and of the factor Z
(right) after 20 steps of the biconjugation process.

Figure

5.3-4: Structure of Z after 70 steps (left) and 200 steps (right)
of the biconjugation process.
Michele Benzi and Miroslav T-uma

Figure

5.5-6: Structure of Z after 400 steps (left) and 760 steps (right)
of the biconjugation process.

Figure

5.7-8: Structure of entries in Z larger than 10 \Gamma10 (left) and structure
of incomplete factor -
Z with drop tolerance
Approximate Inverse Preconditioning 17
A sufficient condition to have a fill-in in the matrix Z after some steps of the biconju-
gation algorithm is given by the following Lemma.
Lemma 5.2. Let E) be a bipartite graph with and such that
If for some indices i l there is a path
in B, then z (i p )
Proof. We use induction on p. Let
0: Of course,
z (i 1 \Gamma1)
and from Lemma 5.1 we get z (i 1 )
Suppose now that Lemma 5.2 is true for all l ! p. Then, z (i
a
using the no-cancellation assumption
we also have z (i p )
The following theorem gives a necessary and sufficient condition for a nonzero entry to
appear in position (l; j), l ! j, in the inverse triangular factor.
Theorem 5.3. Let 1 only if for some p - 1 there are
Proof. We first show that the stated conditions are sufficient. By Lemma 5.1, the nonzeros
a
imply that z (i 1 )
l 1 j is also nonzero. If are done. Otherwise, z (i 2 \Gamma1)
and a
0: Taking into account that z l 2 we get that z (i 2 )
is
nonzero. Repeating these arguments inductively we finally get z (i p )
l
Consequently,
z (i)
Assume now that z lj 6= 0. Lemma 5.1 implies that at least one of the following two
conditions holds: either there exists li 0 6= 0, or there
exist indices i such that a i 00 k 6= 0, z (i 00 \Gamma1)
li 00 6= 0: In the
former case we have the necessary conditions. In the latter case we can apply Lemma 5.1
inductively to z (i 00 \Gamma1)
After at most j inductive steps we obtain the conditions. 2
Clearly, the characterization of fill-in in the inverse triangular factorization is less transparent
than the necessary and sufficient condition which characterize nonzeros in the non-
factorized inverse.
Michele Benzi and Miroslav T-uma
6. Preconditioning for block triangular matrices. Many sparse matrices arising
in real-world applications may be reduced to block triangular form (see Ch. 6 in [20]). In
this section we discuss the application of preconditioning techniques to linear systems with
a block (lower) triangular coefficient matrix, closely following [30].
The reduction to block triangular form is usually obtained with a two-step procedure,
as outlined in [20]. In the first step, the rows of A are permuted to bring nonzero entries on
the main diagonal, producing a matrix PA. In the second step, symmetric permutations
are used to find the block triangular form [53]. The resulting matrix can be represented as
A k1 A k2 \Delta \Delta \Delta A kkC C C C A
where the diagonal blocks A ii are assumed to be irreducible. Because A is nonsingular, the
diagonal blocks A ii must also be nonsingular.
Suppose that we compute approximate inverses of the diagonal blocks A
the incomplete biconjugation algorithm, so that A \Gamma1
Z ii
ii
ii
the inverse of A is approximated as follows (cf. [30]):
A
22
A k1 A k2 \Delta \Delta
QP:
The preconditioning step in a conjugate gradient-type method requires the evaluation
of the action of G on a vector, i.e. the computation of z = Gd for a given vector d, at each
step of the preconditioned iterative method. This can be done by a back-substitution of the
z
where
d =B @
z =B @
Approximate Inverse Preconditioning 19
with the partitioning of - z and -
d induced by the block structure of Q(PA)Q T : The computation
of which is required by certain preconditioned iterative methods, is
accomplished in a similar way.
With this approach, fill-in is confined to the approximate inverses of the diagonal blocks,
often resulting in a more sparse preconditioner. Notice also that the approximate inverses
G ii can be computed in parallel. The price to pay is the loss of part of the explicitness
when the approximate inverse preconditioner is applied, as noted in [30].
For comparison purposes, we apply the same scheme with ILU preconditioning. Specif-
ically, we approximate A as
A 21
A k1 A k2
where each diagonal block A ii is approximated by an ILU decomposition -
U ii . Applying
the preconditioner requires the solution of a linear system d at each step of the
preconditioned iteration. This can be done with the back-substitution
where
with the same partitioning of - z and -
d as above. The use of transposed ILU preconditioning
is similar.
this type of ILU block preconditioning we introduce some explicitness in the
application of the preconditioner. Again, note that the ILU factorizations of the diagonal
blocks can be performed in parallel.
We will see in the section on numerical experiments that reduction to the block triangular
form influences the behavior of the preconditioned iterations in different ways depending
on whether approximate inverse techniques or ILU-type preconditioning are used.
Michele Benzi and Miroslav T-uma
7. Implementation aspects. It is possible to implement the incomplete inverse factorization
algorithm in x4 in at least two distinct ways. The first implementation is similar in
spirit to the classical submatrix formulation of sparse Gaussian elimination as represented,
for instance, in [19],[57]. This approach relies on sparse incomplete rank-one updates of
the matrices -
Z and -
applied in the form of outer vector products. These updates are
the most time-consuming part of the computation. In the course of the updates, new fill-in
elements whose magnitude is less than a prescribed drop tolerance T are dropped. In this
approach, dynamic data structures have to be used for the -
Z and -
matrices. Note that
at step i of the incomplete inverse factorization, only the ith row a T
and the ith column c T
are required. The matrix A is stored in static data structures both by rows and by columns
(of course, a single array is needed for the numerical values of the entries of A).
For this implementation to be efficient, some additional elbow room is necessary. For
instance, in the computation of the incomplete -
Z factor the elbow room was twice the
space anticipated for storing the nonzeros in the factor itself. As we are looking for a
preconditioner with about the same number of nonzeros as the original matrix, the estimated
number of nonzeros in -
Z is half the number of nonzeros in the original matrix A. For each
column of -
Z we give an initial prediction of fill-in based on the results of x5. Thus, the
initial structure of -
Z is given by the structure of the upper triangular part of A. Of course,
W is handled similarly. If the space initially allocated for a given column is not enough, the
situation is solved in a way which is standard when working with dynamic data structures,
by looking for a block of free space at the end of the active part of the dynamic data
structure large enough to contain the current column, or by a garbage collection (see [57]).
Because most of the fill-in in -
Z and -
W appears in the late steps of the biconjugation process,
we were able to keep the amount of dynamic data structure manipulations at relatively low
levels. In the following, this implementation will be referred to as the DDS implementation.
Despite our efforts to minimize the amount of symbolic manipulations in the DDS im-
plementation, some of its disadvantages such as the nonlocal character of the computations
and a high proportion of non-floating-point operations still remain. This is an important
drawback of submatrix (right-looking, undelayed) algorithms using dynamic data structures
when no useful structural prediction is known and no efficient block strategy is used. Even
when all the operations are performed in-core, the work with both the row and column
lists in each step of the outer cycle is rather irregular. Therefore, for larger problems, most
operations are still scattered around the memory and are out-of-cache. As a consequence,
it is difficult to achieve high efficiency with the code, and any attempt to parallelize the
Approximate Inverse Preconditioning 21
computation of the preconditioner in this form will face serious problems (see [57] for a
discussion of the difficulties in parallelizing sparse rank-one updates).
For these reasons we considered an alternative implementation (hereafter referred to
as SDS) which only makes use of static data structures, based on a left-looking, delayed
update version of the biconjugation algorithm. This amounts to a rearrangement of the
computations, as shown below. For simplicity we only consider the Z factor, and assume
no breakdown occurs:
(1) Let z (0)= e 1
z (0)
(j \Gamma1)
:= a T
z (j \Gamma1)
z (j)
:= z (j \Gamma1)
(j \Gamma1)
(j \Gamma1)
z (j \Gamma1)
:= a T
z (i\Gamma1)
This procedure can be implemented with only static data structures, at the cost of increasing
the number of floating-point operations. Indeed, in our implementation we found
it necessary to recompute the dot products p (j \Gamma1)
z (j \Gamma1)
if they are used more than
once for updating subsequent columns. This increase in arithmetic complexity is more or
less pronounced, depending on the problem and on the density of the preconditioner. On
the other hand, this formulation greatly decreases the amount of irregular data structure
manipulations. It also appears better suited to parallel implementation, because the dot
products and the vector updates in the innermost loop can be done in parallel. Notice that
with SDS, it is no longer true that a single row and column of A are used at each step of the
outer loop. It is worth mentioning that numerically, the DDS and SDS implementations of
the incomplete biconjugation process are completely equivalent.
The SDS implementation is straightforward. Suppose the first steps have been
completed. In order to determine which columns of the already determined part of -
Z play
22 Michele Benzi and Miroslav T-uma
a role in the rank-one updates used to form the jth column of -
Z we only need a linked
list scanning the structure of the columns of A. This linked list is coded similarly to the
mechanism which determines the structure of the jth row of the Cholesky factor L in the
numerical factorization in SPARSPAK (see [27],[13]).
In addition to the approximate inverse preconditioner, we also coded the standard row
implementation of the classical ILU(0) preconditioner (see, e.g., [50]). We chose a no-fill
implicit preconditioner because we are mostly interested in comparing preconditioners with
a nonzero density close to that of the original matrix A.
On input, all our codes for the computation of the preconditioners check whether the
coefficient matrix has a zero-free diagonal. If not, row reordering of the matrix is used
to permute nonzeros on the diagonal. For both the ILU(0) and the approximate inverse
factorization, we introduced a simple pivot modification to avoid breakdown. Whenever
some diagonal element in any of our algorithms to compute a preconditioner was found to
be small, in our case less in absolute value than the IEEE machine precision ffl - 2:2
we increased it to 10 \Gamma3 . We have no special reasons for this choice, other than it worked well
in practice. It should be mentioned that in the numerical experiments, this safeguarding
measure was required more often for ILU(0) than for the approximate inverse factorization.
For the experiments on matrices which can be nontrivially reduced to block triangular
form, we used the routine MC13D from MA28 [19] to get the block triangular form.
8. Numerical experiments. In this section we present the results of numerical experiments
on a range of problems from the Harwell-Boeing collection [21] and from Tim
Davis' collection [17]. All matrices used were rescaled by dividing their elements by the
absolute value of their largest nonzero entry. No other scaling was used. The right-hand
side of each linear system was computed from the solution vector x   of all ones, the choice
used, e.g., in [57].
We experimented with several iterative solvers of the conjugate gradient type. Here we
present results for three selected methods, which we found to be sufficiently representative:
van der Vorst's Bi-CGSTAB method (denoted BST in the tables), the QMR method of
Freund and Nachtigal, and Saad and Schultz's GMRES (restarted every 20 steps, denoted
G(20) in the tables) with Householder orthogonalization [56]. See [3] for a description of
these methods, and the report [9] for experiments with other solvers.
Approximate Inverse Preconditioning 23
The matrices used in the experiments come from reservoir simulation (ORS*, PORES2,
SAYLR* and SHERMAN*), chemical kinetics (FS5414), network flow (HOR131), circuit
simulation (JPWH991, MEMPLUS and ADD*), petroleum engineering (WATT* matrices)
and incompressible flow computations (RAEFSKY*, SWANG1). The order N and number
NNZ of nonzeros for each test problem are given in Table 1, together with the number
of iterations and computing times for the unpreconditioned iterative methods. A y means
that convergence was not attained in 1000 iterations for Bi-CGSTAB and QMR, and 500
iterations for GMRES(20).
Its Time

Table

1: Test problems (N= order of matrix, NNZ= nonzeros in matrix) and convergence
results for the iterative methods without preconditioning.
Michele Benzi and Miroslav T-uma
All tests were performed on a SGI Crimson workstation with RISC processor R4000
using double precision arithmetic. Codes were written in standard Fortran 77 and compiled
with the optimization option -O4. CPU time is given in seconds and it was measured using
the standard function dtime.
The initial guess for the iterative solvers was always x The stopping criterion
used was jjr k jj is the (unpreconditioned) updated residual. Note that
because r we have that 1 - jjr 0 jj 1 - nzr where nzr denotes the maximum
number of nonzeros in a row of A.
The following tables present the results of experiments with the ILU(0) preconditioner
and with the approximate inverse preconditioner based on the biconjugation process (here-
after referred to as AIBC). Observe that the number of nonzeros in the ILU(0) preconditioner
is equal to the number NNZ of nonzeros in the original matrix, whereas for the AIBC
preconditioner fill-in is given by the total number of nonzeros in the factors -
W and -
D.
In the tables, the number of nonzeros in AIBC is denoted by F ill. Right preconditioning
was used for all the experiments.
The comparison between the implicit and the explicit preconditioner is based on the
amount of fill and on the rate of convergence as measured by the number of iterations.
These two parameters can realistically describe the scalar behavior of the preconditioned
iterative methods. Of course, an important advantage of the inverse preconditioner, its
explicitness, is not captured by this description.
The accuracy of the AIBC preconditioner is controlled by the value of the drop tolerance
T . Smaller drop tolerances result in a more dense preconditioner and very often (but not
always) in a higher convergence rate for the preconditioned iteration. For our experiments
we consider relatively sparse preconditioners. In most cases we were able to adjust the
value of T so as to obtain an inverse preconditioner with a nonzero density close to that
of A (and hence of the ILU(0) preconditioner). Due to the scaling of the matrix entries,
the choice very often the right one. We also give results for the approximate
inverse obtained with a somewhat smaller value of the drop tolerance, in order to show how
the number of iterations can be reduced by allowing more fill-in in the preconditioner. For
some problems we could not find a value of T for which the number of nonzeros in AIBC
is close to NNZ. In these cases the approximate inverse preconditioner tended to be either
very dense or very sparse.
Approximate Inverse Preconditioning 25
In

Table

2 we give the timings for the preconditioner computation, iteration counts and
timings for the three iterative solvers preconditioned with ILU(0). The same information is
given in Table 3 for the approximate inverse preconditioner AIBC. For AIBC we give two
timings for the construction of the preconditioner, the first for the DDS implementation
using dynamic data structures and the second for the SDS implementation using only static
data structures.
ILU - Its ILU - Time
MATRIX P-time BST QMR G(20) BST QMR G(20)
RAEFSKY1 2.457

Table

2: Time to form the ILU(0) preconditioner (P-time), number of iterations and time
for Bi-CGSTAB, QMR and GMRES(20) with ILU(0) preconditioning.
26 Michele Benzi and Miroslav T-uma
P-time AIBC - Its AIBC - Time
MATRIX Fill DDS SDS BST QMR G(20) BST QMR G(20)
5204
JPWH991 7063 0.31 0.26 15 27 28 0.24 0.67 0.78
48362 0.68 2.63 33 43 64 2.61 5.39 8.63
26654 0.89 2.45

Table

3: Time to form the AIBC preconditioner (P-time) using DDS and SDS implemen-
tations, number of iterations and time for Bi-CGSTAB, QMR and GMRES(20) with AIBC
Approximate Inverse Preconditioning 27
It appears from these results that the ILU(0) and AIBC preconditioners are roughly
equivalent from the point of view of the rate of convergence, with ILU(0) having a slight
edge. On many problems the two preconditioners give similar results. There are a few cases,
like PORES2, for which ILU(0) is much better than AIBC, and others (like MEMPLUS)
where the situation is reversed. For some problems it is necessary to allow a relatively
high fill in the approximate inverse preconditioner in order to have a convergence rate
comparable with that insured by ILU(0) (cf. SAYLR4), but there are cases where a very
sparse AIBC gives excellent results (see the ADD or the RAEFSKY matrices). It follows
that the timings for the iterative part of the solution process are pretty close, on average,
for the two preconditioners.
We also notice that using a more dense approximate inverse preconditioner (obtained
with a smaller value of T ) nearly always reduces the number of iterations, although this
does not necessarily mean a reduced computing time since it takes longer to compute the
preconditioner and the cost of each iteration is increased.
Concerning the matrix PORES2, for which our method gives poor results, we observed
that fill-in in the -
W factor was very high. We tried to use different drop tolerances for the
two factors (the one for -
being larger than the one used for -
Z) but this did not help. It
was observed in [31] that finding a sparse right approximate inverse for PORES2 is very
hard and a left approximate inverse should be approximated instead. Unfortunately, our
method produces exactly the same approximate inverse (up to transposition) for A and
A T , therefore we were not able to cope with this problem effectively. We experienced a
similar difficulty with the -
W factor for the matrix SHERMAN2. On the other hand, for
SHERMAN3 we did not face any of the problems reported in [30] and convergence with the
AIBC preconditioner was smooth.
As for the time required to compute the preconditioners, it is obvious that ILU(0) can be
computed more quickly. On the other hand, the computation of the AIBC preconditioner is
not prohibitive. There are problems for which computing AIBC is only two to three times
more expensive than computing ILU(0). More important, our experiments with AIBC
show that the overall solution time is almost always dominated by the iterative part, unless
convergence is extremely rapid, in which case the iteration part takes slightly less time than
the computation of the preconditioner.
This observation suggests that our approximate inverse preconditioner is much cheaper
to construct, in a sequential environment, than approximate inverse preconditioners based
28 Michele Benzi and Miroslav T-uma
on the Frobenius norm approach described in x3. Indeed, if we look at the results presented
in [30] we see that the sequential time required to construct the preconditioner accounts for a
huge portion, often in excess of 90%, of the overall computing time. It is worth emphasizing
that the approach based on Frobenius norm minimization and the one we propose seem to
produce preconditioners of similar quality, in the sense that they are both comparable with
ILU(0) from the point of view of fill-in and rates of convergence, at least on average.
As for the different implementations of AIBC, we see from the results in Table 3 that for
larger problems, the effect of additional floating-point operations in the SDS implementation
is such that the DDS implementation is actually faster. Nevertheless, as already observed
the implementation using static data structures may better suited for parallel architec-
tures. Because in this paper we only consider a scalar implementation, in the remaining
experiments we limit ourselves to the timings for the DDS implementation of AIBC.
In all the experiments (excluding the ones performed to measure the timings presented
in the tables) we monitored also the "true" residual jjb \Gamma Ax k jj 2 . In general, we found that
the discrepancy between this and the norm of the updated residual was small. However,
we found that for some very ill-conditioned matrices in the Harwell-Boeing collection (not
included in the tables) this difference may be very large. For instance, for some of the LNS*
and WEST* matrices, we found that jjr k jj for the final value of r k . This
happened both with the ILU(0) and with the approximate inverse preconditioner, and we
regarded this as a failure of the preconditioned iterative method.
We present in Tables 4 and 5 the results of some experiments on matrices which have
been reduced to block lower triangular form. We compared the number of iterations of
the preconditioned iterative methods and their timings for the block approximate inverse
preconditioner and for the block ILU(0) preconditioner as described in x6. Since some of
the matrices have only trivial block lower triangular form (one block, or two blocks with one
of the blocks of dimension one for some matrices) we excluded them from our experiments.
In

Table

4 we give for each matrix the number NBL of blocks and the results of experiments
with ILU(0). In Table 5 we give analogous results for the AIBC preconditioner. The
amount of fill-in (denoted by F ill) for AIBC is computed as the fill-in in the approximate
inverses of the diagonal blocks plus the number of nonzero entries in the off-diagonal blocks.
Approximate Inverse Preconditioning 29
Block ILU - Its Block ILU - Time

Table

4: Time to compute the block ILU preconditioner (P-time), number of iterations and
time for Bi-CGSTAB, QMR and GMRES(20) with block ILU(0) preconditioning.
Block AIBC - Its Block AIBC - Time
MATRIX Fill P-time BST QMR G(20) BST QMR G(20)

Table

5: Time to compute the block AIBC preconditioner (P-time) , number of iterations and
time in seconds for Bi-CGSTAB, QMR and GMRES(20) with block AIBC preconditioning.
It is clear that in general the reduction to block triangular form does not lead to a
noticeable improvement in the timings, at least in a sequential implementation. We observe
that when the block form is used, the results for ILU(0) are sometimes worse. This can
Michele Benzi and Miroslav T-uma
probably be attributed to the permutations, which are known to cause in some cases a
degradation of the rate of convergence of the preconditioned iterative method [22]. A
notable exception is the matrix WATT2, for which the number of iterations is greatly
reduced. On the other hand, the results for the block approximate inverse preconditioner
are mostly unchanged or somewhat better. Again, matrix WATT2 represents an exception:
this problem greatly benefits from the reduction to block triangular form. In any case,
permutations did not adversely affect the rate of convergence of the preconditioned iterative
method. This fact suggests that perhaps the approximate inverse preconditioner is more
robust than ILU(0) with respect to reorderings.
To gain more insight on how permutations of the original matrix can influence the
quality of both types of preconditioners, we did some experiments where the matrix A was
permuted using the minimum degree algorithm on the structure of A + A T (see [28]). We
applied the resulting permutation to A symmetrically to get PAP T , in order to preserve
the nonzero diagonal. Tables 6 and 7 present the results for the test matrices having trivial
block triangular form. The corresponding preconditioners are denoted by ILU(0)-MD and
AIBC-MD, respectively.
ILU-MD - Its ILU-MD - Time
MATRIX P-time BST QMR G(20) BST QMR G(20)
26 43 47 2.18 4.57 6.45

Table

Time to compute the ILU(0) preconditioner (P-time) for A permuted according
to minimum degree algorithm on A number of iterations and time for Bi-CGSTAB,
QMR and GMRES(20) with ILU(0)-MD preconditioning.
Approximate Inverse Preconditioning 31
AIBC-MD - Its AIBC-MD - Time
MATRIX Fill P-time BST QMR G(20) BST QMR G(20)
7152 0.31
43 48 0.38 1.23 1.38
19409 0.58 104 95 y 2.91 4.14 y

Table

7: Time to compute the AIBC preconditioner (P-time) for A permuted by the minimum
degree algorithm on A number of iterations and time for Bi-CGSTAB, QMR
and GMRES(20) with AIBC-MD preconditioning.
The results in Table 6 show that for some problems, especially those coming from PDEs,
minimum degree reordering has a detrimental effect on the convergence of the iterative
solvers preconditioned with ILU(0). In some cases we see a dramatic increase in the number
of iterations. This is in analogy with the observed fact (see, e.g., [22]) that when the
minimum degree ordering is used, the no-fill incomplete Cholesky decomposition of an SPD
Michele Benzi and Miroslav T-uma
matrix is a poor approximation of the coefficient matrix, at least for problems arising from
the discretization of 2D PDEs. The convergence of the conjugate gradient method with such
a preconditioner (ICCG(0)) is much slower than if the natural ordering of the unknowns
was used. Here we observe a similar phenomenon for nonsymmetric linear systems. Note
the rather striking behavior of matrix ADD20, which benefits greatly from the minimum
degree reordering (this matrix arises from a circuit model and not from the discretization
of a PDE).
It was also observed in [22] that the negative impact of minimum degree on the rate
of convergence of PCG all but disappears when the incomplete Cholesky factorization of
A is computed by means of a drop tolerance rather than by position. It is natural to ask
whether the same holds true for the approximate inverse preconditioner AIBC, which is
computed using a drop tolerance. The results in Table 7 show that this is indeed the case.
For most of the test problems the number of iterations was nearly unaffected (or better)
and in addition we note that the minimum degree ordering helps in preserving sparsity in
the incomplete inverse factors. While this is usually not enough to decrease the computing
times, the fact that it is possible to reduce storage demands for the approximate inverse
preconditioner without negatively affecting the convergence rates might become important
for very large problems.
We conclude this section with some observations concerning the choice of the drop
tolerance T . In all our experiments we used a fixed value of T throughout the incomplete
biconjugation process. However, relative drop tolerances, whose value is adapted from step
to step, could also be considered (see [57] for a thorough discussion of the issues related to
the choice of drop tolerances in the context of ILU). We have observed that the amount of
fill-in is distributed rather unevenly in the course of the approximate inverse factorization. A
large proportion of nonzeros is usually concentrated in the last several columns of -
Z and -
W .
For some problems with large fill, it may be preferable to switch to a larger drop tolerance
when the columns of the incomplete factors start filling-in strongly. Conversely, suppose
we have computed an approximate inverse preconditioner for a certain value of T , and we
find that the preconditioned iteration is converging slowly. Provided that enough storage is
available, one could then try to recompute at least some of the columns of -
Z and -
using
a smaller value of T . Unfortunately, for general sparse matrices there is no guarantee that
this will result in a preconditioner of improved quality. Indeed, allowing more nonzeros in
the preconditioner does not always result in a reduced number of iterations.
Approximate Inverse Preconditioning 33
Finally, it is worthwhile to observe that a dual threshold variant of the incomplete
inverse factorization could be adopted, see [51]. In this approach, a drop tolerance is
applied but a maximum number of nonzeros per column is specified and enforced during
the computation of the preconditioner. In this way, it is possible to control the maximum
storage needed by the preconditioner, which is important for an automated implementation.
This approach has not been tried yet, but we hope to do so in the near future.
9. Conclusions and future work. In this paper we have developed a sparse approximate
inverse preconditioning technique for nonsymmetric linear systems. Our approach is
based on a procedure to compute two sets of biconjugate vectors, performed incompletely
to preserve sparsity. This algorithm produces an approximate triangular factorization of
A \Gamma1 , which is guaranteed to exist if A is an H-matrix (similar to the ILU factorization).
The factorized sparse approximate inverse is used as an explicit preconditioner for
conjugate gradient-type methods. Applying the preconditioner only requires sparse matrix-vector
products, which is of considerable interest for use on parallel computers.
The new preconditioner was used to enhance the convergence of different iterative
solvers. Based on extensive numerical experiments, we found that our preconditioner can
insure convergence rates which are comparable, on average, with those from the standard
preconditioner. While the approximate inverse factorization is more time-consuming
to compute than ILU(0), its cost is not prohibitive, and is typically dominated by
the time required by the iterative part. This is in contrast with other approximate inverse
preconditioners, based on Frobenius norm minimization, which produce similar convergence
rates but are very expensive to compute.
It is possible that in a parallel environment the situation will be reversed, since the
preconditioner construction with the Frobenius norm approach is inherently parallel. How-
ever, there is some scope for parallelization also in the inverse factorization on which our
method is based: for instance, the approximate inverse factors -
Z and -
W can be computed
largely independent of each other. Clearly, this is a point which requires further research,
and no conclusion can be drawn until parallel versions of this and other approximate inverse
preconditioners have been implemented and tested.
Our results point to the fact that the quality of the approximate inverse preconditioner
is not greatly affected by reorderings of the coefficient matrix. This is important in practice
because it suggests that we may use permutations to increase the potential for parallelism or
to reduce the amount of fill in the preconditioner, without spoiling the rate of convergence.
34 Michele Benzi and Miroslav T-uma
The theoretical results on fill-in in x5 provide guidelines for the use of pivoting strategies for
enhancing the sparsity of the approximate inverse factors, and this is a topic that deserves
further research.
Based on the results of our experiments, we conclude that the technique introduced
in this paper has the potential to become a useful tool for the solution of large sparse
nonsymmetric linear systems on modern high-performance architectures. Work on a parallel
implementation of the new preconditioner is currently under way. Future work will also
include a dual threshold implementation of the preconditioner computation.

Acknowledgments

. We would like to thank one of the referees for helpful comments and
suggestions, and Professor Miroslav Fiedler for providing reference [24]. The first author
gratefully acknowledges the hospitality and excellent research environment provided by the
Institute of Computer Science of the Czech Academy of Sciences.



--R

Parallel Implementation of Preconditioned Conjugate Gradient Methods for Solving Sparse Systems of Linear Equations.
Iterative Solution Methods.
Templates for the Solution of Linear Systems.
Parallel algorithms for the solution of certain large sparse linear systems.
A Direct Row-Projection Method for Sparse Linear Systems
A direct projection method for sparse linear systems.
An explicit preconditioner for the conjugate gradient method.
A sparse approximate inverse preconditioner for the conjugate gradient method.
A sparse approximate inverse preconditioner for nonsymmetric linear systems.
Krylov methods preconditioned with incompletely factored matrices on the CM-2
Approximate inverse preconditioners for general sparse matrices.
Approximate inverse techniques for block-partitioned matrices
User's guide for SPARSPAK-A: Waterloo sparse linear equations package

Block preconditioning for the conjugate gradient method.
Approximate inverse preconditionings for sparse linear systems.
Sparse matrix collection.
Decay rates for inverses of band matrices.

Direct Methods for Sparse Matrices.
Users' guide for the Harwell-Boeing sparse matrix collection
The effect of ordering on preconditioned conjugate gradients.
A stability analysis of incomplete LU factorizations.
Inversion of bigraphs and connection with the Gauss elimination.
An Introduction to Numerical Linear Algebra.
Notes on the solution of algebraic linear simultaneous equations.
Computer Solution of Large Sparse Positive Definite Systems.
The evolution of the minimum degree algorithm.
Predicting structure in sparse matrix computations.
On approximate-inverse preconditioners
Parallel preconditioning with sparse approximate inverses.
Parallel preconditioning and approximate inverses on the Connection Machine.
A parallel preconditioned conjugate gradient package for solving sparse linear systems on a Cray Y-MP
Inversion of matrices by biorthogonalization and related results.
Method of conjugate gradients for solving linear systems.
The Theory of Matrices in Numerical Analysis.
Polynomial preconditioning for conjugate gradient calculations.
The efficient parallel iterative solution of large sparse linear sys- tems
Explicitly preconditioned conjugate gradient method for the solution of unsymmetric linear systems.

New convergence results and preconditioning strategies for the conjugate gradient method.
Factorized sparse approximate inverse (FSAI) preconditionings for solving 3D FE systems on massively parallel computers II: Iterative construction of FSAI preconditioners.
Factorized sparse approximate inverse preconditioning I: Theory.
Factorized sparse approximate inverse preconditioning II: Solution of 3D FE systems on massively parallel computers.
Krylov Methods for the Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations
An incomplete factorization technique for positive definite linear systems.
An iterative solution method for linear systems of which the coefficient matrix is a symmetric M-matrix
Some properties of approximate inverses of matrices.
Preconditioning techniques for nonsymmetric and indefinite linear systems.
SPARSKIT: A basic tool kit for sparse matrix computations.
ILUT: A dual threshold incomplete LU factorization.
Conjugate direction methods for solving systems of linear equations.

High performance preconditioning.
Matrix Iterative Analysis.
Implementation of the GMRES method using Householder transformations.
Computational Methods for General Sparse Matrices.
--TR

--CTR
Kai Wang , Jun Zhang, Multigrid treatment and robustness enhancement for factored sparse approximate inverse preconditioning, Applied Numerical Mathematics, v.43 n.4, p.483-500, December 2002
Claus Koschinski, New methods for adapting and for approximating inverses as preconditioners, Applied Numerical Mathematics, v.41 n.1, p.179-218, April 2002
Stephen T. Barnard , Luis M. Bernardo , Horst D. Simon, An MPI Implementation of the SPAI Preconditioner on the T3E, International Journal of High Performance Computing Applications, v.13 n.2, p.107-123, May       1999
N. Guessous , O. Souhar, Multilevel block ILU preconditioner for sparse nonsymmetric M-matrices, Journal of Computational and Applied Mathematics, v.162 n.1, p.231-246, 1 January 2004
Matthias Bollhfer , Volker Mehrmann, Some convergence estimates for algebraic multilevel preconditioners, Contemporary mathematics: theory and applications, American Mathematical Society, Boston, MA, 2001
Michele Benzi , Miroslav Tma, A parallel solver for large-scale Markov chains, Applied Numerical Mathematics, v.41 n.1, p.135-153, April 2002
Mansoor Rezghi , S. Mohammad Hosseini, An ILU preconditioner for nonsymmetric positive definite matrices by using the conjugate Gram-Schmidt process, Journal of Computational and Applied Mathematics, v.188 n.1, p.150-164, 1 April 2006
M. H. Koulaei , F. Toutounian, On computing of block ILU preconditioner for block tridiagonal systems, Journal of Computational and Applied Mathematics, v.202 n.2, p.248-257, May, 2007
Michele Benzi, Preconditioning techniques for large linear systems: a survey, Journal of Computational Physics, v.182 n.2, p.418-477, November 2002
approximate inverses;incomplete factorizations;sparse linear systems;conjugate gradient-type methods;sparse matrices;preconditioning