On Two Variants of an Algebraic Wavelet Preconditioner.
A recursive method of constructing preconditioning matrices for the nonsymmetric stiffness matrix in a wavelet basis is proposed for solving a class of integral and differential equations. It is based on a level-by-level application of the wavelet scales decoupling the different wavelet levels in a matrix form just as in the well-known nonstandard form. The result is a powerful iterative method with built-in preconditioning leading to two specific algebraic multilevel iteration algorithms: one with an exact Schur preconditioning and the other with an approximate Schur preconditioning. Numerical examples are presented to illustrate the efficiency of the new algorithms.
Introduction
The discovery of wavelets is usually described as one of the most important advances in mathematics
in the twentieth century as a result of joint eorts of pure and applied mathematicians. Through
the powerful compression property, wavelets have satisfactorily solved many important problems
in applied mathematics e.g. signal and image processing; see [23, 20, 34, 38] for a summary.
There remain many mathematical problems to be tackled before wavelets can be used for
solution of dierential and integral equations in a general setting. The traditional wavelets were
designed mainly for regular domains and uniform meshes. This was one of the reasons why wavelets
may not be immediately applicable to arbitrary problems. The introduction of the lifting idea,
interpolatory wavelets [35, 25, 1] and adaptivity [17] provides a useful way of constructing wavelets
functions in non-regular domains and in high dimensions.
However, the algebraic (sparse) structure of the matrix generated by a wavelet method is usually
a nger-like one that is a di-cult sparse pattern to deal with; refer to [10, 13, 15, 18]. Firstly
direct solution of a linear system with such a matrix is either not feasible or ine-cient. Secondly
iterative solution requires a suitable preconditioner and this choice of preconditioner is usually dependent
of the smoothness of the underlying operator (in addition to the assumptions for wavelets
compression). Often a diagonal preconditioner is not su-cient. For some particular problems, several
preconditioning techniques have been suggested. For instance, the nger matrix from wavelets
representation of a Calderon-Zygmund operator plus a non-constant diagonal matrix cannot be
preconditioned eectively by a diagonal matrix. In this case, one can use the idea of two-stage
preconditioning as proposed in [13] or to use other modied wavelets methods such as a centring
algorithm [15]; see also [38] for another modied algorithm and [19] for using approximate inverses.
there exists a gap in realizing the full e-ciency oered by wavelet bases for model prob-
lems. That is to say, a generally applicable iterative algorithm is still lacking. For a recent and
general survey of iterative methods, refer to [31].
This paper proposes two related and e-cient iterative algorithms based on the wavelet formulation
for solving an operator equation with conventional arithmetic. Both algorithms use the
Schur complements recursively but dier in how to use coarse levels to solve Schur complements
equations. In the rst algorithm, we precondition a Schur complement by using coarse levels while
in the second we use approximate Schur complements to construct a preconditioner. We believe
that our algorithms can be adapted to higher dimensional problems more easily than previous work
in the subject.
The motivation of this work follows from the observation that any 1-scale compressed results
(matrices) can be conveniently processed before applying the next scale. In this way, regular patterns
created by past wavelet scales are not destroyed by the new scales like in the non-standard (NS)
form [10] and unlike in the standard wavelet bases; we dene the notation and give further details
in Section 2. This radical but simple idea will be combined in Section 3 with the Schur complement
method and Richardson iterations in a multi-level iterative algorithm. Moreover the Richardson iterations
can be replaced by a recursive generalized minimal residuals (GMRES) method [30]. The
essential assumption for this new algorithm to work is the invertibility of an approximate band
matrix; in the Appendix we show that for a class of Calderon-Zygmund and pseudo-dierential
operators such an invertibility is ensured. In practice we found that our method works equally well
for certain operators outside the type for which we can provide proofs. In Section 4, we present
an alternative way of constructing the preconditioner by using approximate Schur complements.
Section 5 discusses the complexity issues while Section 6 presents several numerical experiments to
illustrate the eectiveness of the two new algorithms.
We remark that our rst algorithm is similar to the framework of a NS form reformulation of
the standard wavelets bases (based on the pyramid algorithm) but does not make use of the NS
form itself, although our algorithm avoids a nger matrix (just like a NS form method) that could
arise from overlapping dierent wavelet scales. The NS form work was by Beylkin, Coifman and
Rokhlin [10] often known as the BCR paper. As a by-product, the NS form reduces the
ops from
O(n log n) to O(n). However, the NS form does not work with conventional arithmetic although
operations with the underlying matrix (that has a regular sparse pattern) can be specially designed;
in fact the NS form matrix itself is simply singular in conventional arithmetic. The recent work in
[22] has attempted to develop a direct solution method based on the NS form that requires a careful
choice of a threshold; here our method is iterative. In the context of designing recursive sparse
preconditioners, it is similar to the ILUM type preconditioner to a certain extent [33]. Our second
algorithm is similar to the algebraic multi-level iteration methods (AMLI) that were developed
for nite elements [4, 2, 3, 36]; here our method uses wavelets and does not require estimating
eigenvalues.
Wavelets splitting of an operator
This section will set up the notation to be used later and motivate the methods in the next
sections. We rst introduce the standard wavelet method. For simplicity, we shall concentrate on
the Daubechies' order m orthogonal wavelets with low pass lters c
lters (such that d In fact, the ideas and expositions in this
paper apply immediately to the more general bi-orthogonal wavelets [20].
Following the usual setting of [10, 20, 22, 11, 34], the lter coe-cients c j 's and d j 's dene the
scaling function (x) and the wavelet function (x). Further, dilations and translations of (x)
and (x) dene a multi-resolution analysis for L 2 in d-dimensions, in particular,
(R d
where the subspaces satisfy the relations
In numerical realisations, we select a nite dimension space V 0 (in the nest scale) as our approximation
space to the innite decomposition of L 2 in (1) i.e. eectively use
to approximate L 2 (R d ). Consequently for a given operator its innite and exact
operator representation in wavelet bases
is approximated in space V 0 by
where are both projection operators.
For brevity, dene operators
Then one can observe that
A further observation based on T
is that the wavelet
coe-cients of T j 1 will be equivalently generated by the block operator
Now we change the notation and consider the discretization of all continuous operators. Dene
A as the representation of operator T 0 in space V 0 . Assume that A on the nest
level is of dimension  the dimension of matrices on a coarse level j is
'. The operator splitting in (3) for the case of dimensions
can be discussed similarly [10, 22]) corresponds to the two-dimensional wavelet transform
e
where the one level transform from j 1 to j (for any
with rectangular matrices P j and Q j (corresponding to operators P j and Q j ) dened respectively
as
For a class of useful and strongly elliptic operators i.e. Calderon-Zygmund and pseudo-dierential
operators, it was shown in BCR [10] that matrices A
k;i are indeed
'sparse' satisfying the decaying property
c m;j
m;j is a generic constant depending on m and j only.
To observe a relationship between the above level-by-level form and the standard wavelet rep-
resentation, dene a square matrix of size  0
I  j
. Then the standard wavelet transform can be
written as
that transforms matrix A into e

Figure

1: The level-by-level form (left) versus the non-standard wavelet form [10] (right)
22464160208Thus the diagonal blocks of e
A are the same as A j 's of a level-by-level form. However the o-
diagonal blocks of the former are dierent from B j and C j of the latter. To gain some insight
into the structure of the o-diagonal blocks of matrix e
A with the standard wavelet transform, we
consider the following case of wavelets. Firstly after level 1
transform, we obtain
e
nn
Secondly after level 2 transform, we get
e
e
nn
Finally after level 3 transform, we arrive at
e
e
I 3n=4
nn
Clearly the o-diagonal blocks of e
A 3 are perturbations of that of the level-by-level form o-diagonal
blocks in fact the one-sided transforms for the o-diagonal blocks are responsible for
the resulting (complicated) sparsity structure. This can be observed more clearly for a typical
example with Fig.1 where the left plot shows the level-by-level
representation set-up that will be used in this paper and in Fig.2 where the left plot shows the
standard wavelet representation as in (8).

Figure

2: The standard wavelet form representation (left) versus an alternative centering form [15]
(right) for the example in Figure 1
12864112Motivated by the exposition in (8) of the standard form, we shall propose a preconditioning and
iterative scheme that operates on recursive 1-level transforms. Thus it will have the advantage of
making full use of the NS form idea and its theory while avoiding the problem of a non-operational
NS form matrix.
Remark 1 Starting from the level-by-level set-up, taking T ' and the collection of all triplets
1j' as a sparse approximation for T 0 is the idea of the NS form [10, 22]. By way of
comparison, in Fig.1, the NS form representation versus the level-by-level form are shown. It turns
out that this work uses the identical set-up to the NS form without using the NS form formulation
itself because we shall not use the proposed sparse approximation. Note that the centering algorithm
[15] (see the right plot in Fig.2) is designed as a permutation of the standard wavelet form (see
the left plot of Fig.2) and is only applicable to a special class of problems where its performance is
better.
3 An exact Schur preconditioner with level-by-level wavelets
We now present our rst and new recursive method for solving the linear system dened
on the nest scale V 0 i.e.
A is of size  0   discussed in the previous section, and x
Instead of considering a representation of T 0 in the decomposition space (2) and then the resulting
linear system, we propose to follow the space decomposition and the intermediate linear system in a
level-by-level manner. A sketch of this method is given in Fig.3 (the left plot) where we try to show
a relationship between the multi-resolution (MR for wavelet representation) and the multi-level
(ML for preconditioning via Schur) ideas from the nest level (top) to the coarsest level (bottom).

Figure

3: Illustration of Algorithms 1 (left) and 2 (right). Here we take
the nest level and 3 the coarsest level), use '2' to indicate a DWT step (one level of wavelets) and
' to denote the direct solution process on the coarsest level. The arrows denote the sequence
of operations (written on the arrowed lines) with each algorithm interacting the two states (two
columns on the plots) of multi-resolution wavelets and multi-level Schur decomposition. The left
plot shows that for Algorithm 1, a Richardson step (or GMRES) takes the results of a DWT step to
the next level via the Schur decomposition while the right plot shows that the Schur decomposition
takes the results of a DWT step to the next level via a Schur approximation.2 Richardson
Schur LU
Schur LU
Schur LU
Direct Solution
MR/Wavelets ML/Schur2
Approximation
Schur LUApproximation
Schur LUSchur LU
Direct Solution
MR/Wavelets ML/Schur
Firstly at level 0, we consider V
and the wavelet transform (4) yields
e
where e x
e
nn
following the general result in (5), it is appropriate to consider the approximation of A
band matrices. To be more precise, let B  (D) denote a banded matrix of D with semi-bandwidth
where integer   0. Dene A suitable  (to be
specied later). Then matrix e
can be approximated by
nn
Or equivalently matrix
nn
is expected to be small in some norm (refer to the Appendix). Write equation (10) as
Consequently we propose to use M our preconditioner to equation (14). This preconditioner
can be used to accelerate iterative solution; we shall consider two such methods: the
Richardson method and the GMRES method [30].
The most important step in an iterative method is to solve the preconditioning equation:
or in a decomposed form  A 1
y (1)y (2)!

r (1)r (2)!
Using the Schur complement method we obtain
y (2)= z 2
y (1)
Here the third equation of (17), unless its dimension is small (i.e. V 1 is the coarsest scale), has to
be solved by an iterative method with the preconditioner T 1 ; we shall denote the preconditioning
step by
where T 1 is of size  1   1 . This sets up the sequence of a multilevel method where the main
characteristic is that each Schur complements equation in its exact form is solved iteratively with
a preconditioner that involves coarse level solutions.
At any level j (1  j < '), the solution of the following linear system
with T j of size  j   j and through solving
e
can be similarly reduced to that of
with T j+1 of size  j+1   j+1 . The solution procedure from the nest level to the coarsest level can
be illustrated by the following diagram (for
Transform
Schur
Precondition
j+1   j+1
A j+1 z
y (2)
y (1)
where as with (12) and (13)
A j+1 B j+1
nn
nn
The coarsest level is set up in the previous section, where a system like (20) is solved
by a direct elimination method. As with conventional multi-level methods, each ne level iteration
leads to many coarse level iteration cycles. This can be illustrated in Fig. 4 where
(bottom plot) are assumed and at the coarsest level (
direct solution is
used. In practice, a variable may be used to achieve certain accuracy for the
preconditioning step i.e. convergence up to a tolerance is pursued whilst by way of comparison
smoothing rather convergence is desired in an usual multilevel method. Our experiments have
shown that are often su-cient to ensure the overall convergence.
We now summarise the formulation as an algorithm. The iterative solver for (19) at level i can
be the Richardson method
or the GMRES method [30] for solving e
(or actually a combination of the two). For
simplicity and generality, we shall use the word \SOLVE" to denote such an iterative solver (either
Richardson or GMRES).
Algorithm 1 (Recursive I)
1. and start on the nest level.
2. Apply one level DWT to T j x to obtain ~
3. Use  j steps of SOLVE for ~
4. In each step, implement the preconditioner
Restrict to the coarse level:
A j+1 z
y (2)
y (1)

Figure

4: Iteration cycling patterns of Algorithm 1 with levels: top for
3. In each case, one solves a ne level equation (starting from the nest level 0) by iteratively
solving coarser level equations  times; on the coarsest level
(here level 3) a direct solution is
carried out.2020
5. Use SOLVE for the above third equation with the preconditioner T j+1 i.e. solve T j+1 x
b j+1 .
7. If (on the coarsest level), apply a direct solver to T j x
and proceed with Step 8; otherwise return to Step 2.
8.
9. Interpolate the coarse level solution to the ne level j:
x (2)
x (1)
x (2)
x (1)
10. Apply one level inverse DWT to e
y j to obtain y j .
11. If (on the nest level), check the residual error | if small enough accept the solution
x 0 and stop the algorithm. If j > 0, check if  j steps (cycles) have been carried out; if not,
return to Step 2 otherwise continue with Step 8 on level j.
The rate of convergence of this algorithm depends on how well the matrix T j approximates
j and this approximation is known to be accurate for a suitable  and for a class of Calderon-
Zygmund and pseudo-dierential operators [10]. For this class of problems, it remains to discuss
the invertibility of matrix A j which is done in the Appendix; a detailed analysis on T j  T j may
be done along the same lines as Lemma 1. For other problem classes, the algorithm may not work
at all for the simple reason that A j may be singular e.g. the diagonal of matrix may
have zero entries. Some extensions based on the idea of [13] may be applied as discussed in Section
6.
Remark 2 We remark that for a class of general sparse linear systems, Saad, Zhang, Botta,
Wubs et al [28, 12, 32, 33] have proposed a recursive multi-level preconditioner (named as ILUM)
similar to this Algorithm 1. The rst dierence is that we need to apply one level of wavelets to
achieve a nearly sparse matrix while these works start from a sparse matrix and permute it to obtain
a desirable pattern suitable for Schur decomposition. The second dierence is that we propose an
iterative step before calling for the Schur decomposition while these works try to compute the exact
Schur decomposition approximately. Therefore it is feasible to rene our Algorithm 1 to adopt the
ILUM idea (using independent sets) for other problem types. However one needs to be careful in
selecting the dimensions of the leading Schur block if a DWT is required for compression purpose.
4 An approximate Schur preconditioner with level-by-level wavelets
In the previous algorithm, we use coarse level equations to precondition the ne level Schur complement
equation. We now propose an alternative way of constructing a preconditioner for a ne level
equation. Namely we approximate and compute the ne level Schur complement before employing
coarse levels to solve the approximated Schur complement equation. A sketch of this method is
shown in Fig.3 showing the natural coupling of wavelet representation (level-by-level form) and
Schur complement. symbol '2'. To dierentiate from Algorithm 1, we change the notation for all
matrices.
At any level k for consider the solution (compare to (19))
A
Applying 1-level of DWT, we obtain
A
A (k)
Note that we have the block LU decomposition
A (k)
A (k)
I
#"
I A (k)1
A (k)0 S (k)
22 A (k)
A (k)
12 is the true Schur complement. To approximate this Schur
complement, we must consider approximating the second term in the above S (k) . We propose to
form band matrix approximations
For level these approximations are possible for a small bandwidth ; see Appendix. Seeking
a band approximation to the inverse of A (k)
makes sense because A (k)
is expected to have a
decaying property (refer to (5)). Let S denote the set of all matrices that have the sparsity pattern
of a band  matrix B
). The formation of a sparse approximate inverse (SPAI) is to nd a
band matrix B 11 2 S such that
min
Refer to [8, 24, 16]. Brie
y as with most SPAI methods, the use of F-norm decouples the minimisation
into least squares (LS) problems for individual columns c j of B 11 . More precisely, owning
to
the j-th LS problem is to solve A (k)
which is not expensive since c j is sparse. Once B 11 is
found, dene an approximation to the true Schur complement S (k) as
22 A (k)
and set A This generates a sequence of matrices A (k) .
Now comes the most important step about the new preconditioner. Setting M
the ne level preconditioner M (0) is dened recursively by
I
#"
I B (k)
is an approximation to the true Schur complement S (k) of A (k) . Here
Observe that this preconditioner is dened through the V-cycling pattern
recursively using the coarse levels.
To go beyond the V-cycling, we propose a simple residual correction idea. We view the solution
y [j]
k of the preconditioning equation (compare to (15) and (9))
as an approximate solution to the equation
Then the residual vector is
k . This calls for a repeated solution M
r k and
gives the correction and a new approximate solution to (25):
In practice we take  see the top
plot in Fig.4) as our experiments suggest that   2 is su-cient to ensure the overall convergence.
Thus an essential feature of this method dierent from Algorithm 1 is that every approximated
Schur complement matrix needs to be transformed to the next wavelet level in order to admit the
matrix splitting (23) while inverse transforms are needed to pass coarse level information back to
a ne level as illustrated in Fig.3. Ideally we may wish to use A
22 A (k)
the approximate Schur complement but taking A (k)
21 and A (k)
12 as full matrices would jeopardize the
e-ciency of the overall iterative method. We proposed to use band matrices (or thresholding) to
approximate these two quantities just as in (13) with Algorithm 1.
To summarise, the solution of the preconditioning equation from the nest level to the coarsest
level and back up is illustrated in Fig.5 for means the same entry
and exit point. The general algorithm for solving M (0) y can be stated as follows:
Algorithm 2 (Recursive II)
1. Apply one level DWT to T k to obtain A (k)
(k)2. Find the approximate inverse B (k)
3. Generate the matrix T
22 A (k)
12 .
Solution Stage
1. and start on the nest level.
2. Apply one level DWT to r k and consider ~
3. Solve the preconditioning equation M (k) ~ y by
Restricting to the coarse level:
r (2)
4. Solve for the above third equation at the next level T k+1 y
5.
(on the coarsest level), apply a direct solver to T k y
and proceed with Step 8; otherwise return to Step 2.
7. Set k := k 1.
8. Interpolate the coarse level k solution to the ne level k:
e
y (2)
y (1)
i:e:
e
y (1)
e
y (2)
e
y (1)
9. Apply one level inverse DWT to e
y k to obtain y k .
10. When (on the nest level), check the residual error | if small enough accept the
solution y 0 and stop the algorithm.
When k > 1, check if  cycles have been carried out; if not, nd the residual vector and return
to Step 2 otherwise continue with Step 7 on level k.
Remark 3 It turns out that this algorithm is similar to the algebraic multi-level iteration methods
(AMLI) that was developed for a class of symmetric positive denite nite element equations in a
hierarchical basis [4, 2, 3, 36]. In fact, let  and then S (k) is implicitly dened by the following
I P
A (k+1)
where P  is a degree  polynomial satisfying
As with AMLI, for the valid choice P 1 gives rise to
the V-cycling pattern and For  > 1, the polynomial P  (t) is chosen to improve the preconditioner;
ideally
A (k)
O(1) asymptotically. Refer to these original papers about how to work
out the coe-cients of P  (t) based on eigenvalue estimates. However we do not use any eigenvalue
estimates to construct P  .
We also remark that an alternative denition of a recursive preconditioner dierent from AMLI
is the ILUM method as mentioned in Remark 2, where in a purely algebraic way (using independent
sets of the underlying matrix graph) A (k)
11 is dened as a block diagonal form after a suitable
permutation. This would give rise to another way of approximating the true Schur complement
A (k). However the sparsity structures of blocks A (k)and A (k)will aect
the density of nonzeros in matrix S (k) and an incomplete LU decomposition has to be pursued as
in [28, 12, 33].
5 Complexity analysis
Here we mainly compare the complexity of Algorithms 1-2 in a full cycle. Note that both algorithms
can be used by a main driving iterative solver where each step of iteration will require n 2
ops
(one
op refers to 1 multiplication and 1 addition) unless some fast matrix vector multiplication
methods are used. One way to reduce this
op count is to use a small threshold so that only a
sparse form of e
A is stored. Also common to both algorithms are the DWT steps which are not
counted here.
To work out
op counts for diering steps of the two algorithms, we list the main steps as follows
Algorithm 1 Algorithm 2
Main 3 band solves: 3n i  2 2 band-band multiplications: 8n i  2
step 4 band-vector multiplications: 8n i  4 band-vector multiplications: 8n i
3 o-band-vector multiplications (R i
Therefore an -cycling (iteration) across all levels would require these
ops (assuming
and
after ignoring the low order terms. Therefore we obtain F II =F I
9
for a typical situation
with I  10. Thus we expect Algorithm I to be cheaper than II if the same
number of iteration steps are recorded. Here by Algorithm I we meant the use of a Richardson
iteration (in SOLVE of Algorithm 1); however if a GMRES iteration is used for preconditioning
then the
op count will increase. Of course as is well known,
op count is not always a reliable
indicator for execution speed; especially if parallel computing is desired a lot of other factors have
to be considered.
Remark 4 For sparse matrices, all
op counts will be much less as DWT matrices are also
sparse. The above complexity analysis is done for a dense matrix case. Even in this case, setting
up preconditioners only adds a few equivalent iteration steps to a conventional iteration solver. One
can usually observe overall speed-up. For integral operators, the proper implementation is to use the
biorthorgonal wavelets as trial functions to yield sparse matrices A directly (for a suitable
threshold); in this case a dierent complexity analysis is needed as all matrices are sparse and
experiments have shown that although this approach is optimal a much larger complexity constant
(independent of n) is involved. For structured dense matrices where FFT is eective, wavelets may
have to be applied implicitly to preserve FFT representation. Further work is in progress.
6 Numerical experiments
Here we shall compare the new algorithms with two previous methods: the WSPAI method by
Chan-Tang-Wan [14] (CTW) and the two stage method by Chan-Chen [13] (CC). Further comparisons
with other methods such as SPAI and ILU type can be found in [14] and [39]. Note that
WSPAI, where applicable, is faster that SPAI and ILU.
We now present numerical experiments from two sets of problems solved by these 4 methods:
M1 | Algorithm 1 with the SOLVE step replaced by a Richardson iteration method for
steps on each level.
M2 | Algorithm 1 with the SOLVE step replaced by a GMRES iteration method on each
level; in particular GMRES(25) for the main nest level iteration and GMRES() for coarser
levels.
M3 | Algorithm 1 with the SOLVE step replaced by a GMRES(25) iteration method on the
nest level (outer iteration method) and a Richardson iteration method for  steps on coarser
levels.
M4 | Algorithm 2 with the main iteration method being a GMRES(25)iteration method on
the nest level and  steps of residual correction on all levels for preconditioning.
The test problems are the following:
{ Example 1. Symmetric case [10]: A
{ Example 2. Symmetric case
log ji Ljlog jj Lj
6 otherwise
Set 2:
{ Example 3. Unsymmetric case: A
{ Example 4. An anisotropic PDE problem in both x and y directions:
where the coe-cients are dened as ([39, 13, Ch.5])
100 (x; y) 2 [0; 0:5]  [0; 0:5] or [0:5; 1]  [0:5; 1]
100 (x; y) 2 [0; 0:5]  [0:5; 1] or [0:5; 1]  [0; 0:5]

Table

1: Number of iteration steps to reach 1). Here  is for
cycling pattern, n is the problem size and ' denotes the number of levels used.
Method Size Levels Case of Case of
used n ' Steps Steps
{ Example 5. A discontinuous coe-cient PDE problem as tested in [14, 39, 13]:
where the coe-cients are dened as ([39, Ch.3])
Here set 1 problems fall into the class for which our algorithms are expected to work; we mainly
test the dependence of the variants of Algorithms 1-2 on the parameter choices. This is because for
symmetric positive denite (SPD) matrices, their principal submatrices are also SPD and invertible
so the assumptions of Lemma 1 (Appendix) are satised. Set 2 problems, although outside the
class of problems that we have provided an analysis, are solvable by previously studied wavelets
preconditioners so we include these for comparison purpose. For all cases, we stop an iterative
method when the residual in the 2-norm is reduced by . The coarsest level is xed at

Tables

1-2 display the numerical results for dierent problem sizes from solving set 1
problems, where 'Steps' means the number of iteration steps required. Clearly one can observe
that 'Steps' is approximately independent of the problem size which is usually expected of a good
preconditioner.
For set 2 problems, we do not expect M1 { M4 to work. To extend these methods, we propose
to combine one step of Stage-1 preconditioning as proposed in [13] to smooth out the given matrix
(w.r.t. level 0 to 1). Specically we select a diagonal matrix such that the sum of
diagonal entries of B 1 and C 1 in
is minimised. Then our algorithms can be applied to the new linear system D 1
points to one way of possible further work. Other approaches similar to stage 1 preconditioning

Table

2: Number of iteration steps to reach 1). Here  is for
cycling pattern, n is the problem size and ' denotes the number of levels used.
Method Size Levels Case of Case of
used n ' Steps Steps
512
43
may be considered e.g. for some indenite problems we have found that the permutation idea by
Du [21] can be combined with our algorithms here; the idea was recently used in [9] to extend the
applicability of sparse approximate inverse preconditioners.
In

Table

3, we use 'Stage-1' to indicate if such a step has been carried out; whenever such
an one-o stage 1 diagonal preconditioning is used, we put \Yes" in this column. The symbol '*'
denotes a case where no convergence has been achieved after 100 steps. Clearly the simple idea
of Stage-1 preconditioning does improve the performance of M1{M4 except for M1 (although M2,
M4 appear to be more robust than M3). Therefore we shall only compare M2{M4 with the work
of CTW [14] and CC [13] next.
Finally in Table 4, we show results from solving Examples 4 and 5 where the CPU seconds
are obtained from a Sun Ultra-2 workstation (using Matlab 5) and the other notation is the same
as in

Table

3. Here 'diag' refers to the diagonal preconditioner which does not work for these
two examples and data with a little 'f' in column 5 are used in constructing Figs.6 and 7. Table
4 demonstrates that when combined with stage-1 preconditioning, the new algorithms can out
perform the previous methods. In particular, it appears that M3 is the fastest method in the
table. To compare the residuals and CPU time of M3 (the best case), CTW [14] and CC [13] in
solving Examples 4 and 5, we plot respectively in Fig.6 and Fig.7 the convergence history and CPU
time (all data are taken from Table 4), where one can observe that M3 (New I) out performs the
others. This again conrms that the new algorithm M3 converges the fastest. As remarked already,
comparisons with other non-wavelets preconditioners (e.g. SPAI and ILU) can be found in [14] and
[39] where it was concluded that WSPAI is faster.
As also remarked earlier, the proposed multi-level algorithms can be potentially developed much
further by incorporating the ideas from [6, 13, 9, 21]. More importantly as they are closely based on
wavelets compression theory, generalizations to multi-dimensions appear to be more straightforward
than similar and known wavelet preconditioners; these aspects are currently being investigated.

Table

3: Number of iteration steps to reach Note that this
example is outside the scope of M1 M4. Here  is for cycling pattern, n is the problem size and
' denotes the number of levels used.
Method With Size Levels Case of Case of
used Stage-1 n ' Steps Steps

Table

4: Comparison of the new algorithms with previous work for Examples 4 5 (set 2). Data
indicated by 'f' in column 5 are used in Figs.6-7.
Problem Stage-1 Method  Convergence Steps CPU
Diag
Diag
CC [13] 54 f 302
Conclusions
This paper has presented two related algorithms implementing an algebraic wavelet preconditioner.
The rst one is similar to the set up of a NS form representation of a wavelet basis while the second
resembles the AMLI preconditioner designed for nite elements. Both algorithms are observed to
give excellent performance for a class of symmetric positive denite Calderon-Zygmund and pseudo-
dierential operators. Combined with a minor stage 1 preconditioning step, they are immediately
applicable to problems outside this class of problems. We note that there are several methods that
are designed to deal with anisotropic and highly indenite elliptic problems [6, 9, 21] and alternative
methods of constructing a sparse preconditioner [33, 12]; these should be investigated in the near
future to further extend the multi-level preconditioner to an even wider class of problems.

Acknowledgements

This work is partially supported by grants NSF ACR 97-20257, NASA Ames NAG2-1238, Sandia
Lab LG-4440 and UK EPSRC GR/R22315. The second author wishes to thank the Department
of Mathematics, UCLA for its hospitality during his visits conducting part of this work.



--R

A wavelet-based approach for the compression of kernel data in large scale simulations of 3D integral problems
Algebraic multilevel preconditioning methods I
Algebraic multilevel preconditioning methods II
Algebraic multilevel iteration method for Stieljes matrices
The algebraic multilevel iteration methods
On the additive version of the algebraic multilevel iteration method for anisotropic elliptic problems
Hierarchical bases and the

Preconditioning highly inde
Fast wavelet transforms and numerical algorithms I
Fast wavelet transforms for matrices arising from boundary element methods
Matrix renumbering ILU: An e

Wavelet sparse approximate inverse preconditioners
Discrete wavelet transforms accelerated sparse preconditioners for dense boundary element systems
An analysis of sparse approximate inverse preconditioners for boundary integral equa- tions
Adaptive wavelet methods for elliptic operator equations - convergence rates
Wavelet methods for second-order elliptic problems
Wavelet adaptive method for second order elliptic problems: boundary conditions and domain decomposition
Wavelet and multiscale methods for operator equations

LU factorization of non-standard forms and direct multiresolution solvers
Multiresolution representation and numerical algorithms: a brief review
On a family of two-level preconditionings of the incomplete block factorization type

Algebraic multilevel iteration preconditioning technique
Preconditioning of inde
ILUM: a multi-elimination ILU preconditioner for general sparse matrices
Iterative Methods for Sparse Linear Systems
GMRES: a generalized minimal residual algorithm for solving unsymmetric linear systems
Iterative solution of linear systems in the 20th century
BILUM: Block versions of multielimination and multilevel ILU preconditioner for general sparse linear systems
Enhanced multi-level block ILU preconditioning strategies for general sparse linear systems
Wavelets and Filter Banks
The lifting scheme: a construction of second generation of wavelets
On two ways of stabilizing the hierarchical basis multilevel methods
Nearly optimal iterative methods for solving
Wavelet Transforms and PDE Techniques in Image Compression
Scalable and multilevel iterative methods
On the multi-level splitting of nite element spaces
--TR
schur complements;multiresolution;sparse approximate inverse;multilevel preconditioner;wavelets;level-by-level transforms