Multigrid Method for Ill-Conditioned Symmetric Toeplitz Systems.
In this paper, we consider solutions of Toeplitz systems where the Toeplitz matrices An are generated by nonnegative functions with zeros. Since the matrices An are ill conditioned, the convergence factor of classical iterative methods, such as the damped Jacobi method, will approach one as the size n of the matrices becomes large. Here we propose to solve the systems by the multigrid method. The cost per iteration for the method is of O(n log n) operations. For a class of Toeplitz matrices which includes weakly diagonally dominant Toeplitz matrices, we show that the convergence factor of the two-grid method is uniformly bounded below one independent of n, and the full multigrid method has convergence factor depending only on the number of levels. Numerical results are given to illustrate the rate of convergence.
Introduction
In this paper we discuss the solutions of ill-conditioned symmetric Toeplitz systems A n by the
multigrid method. The n-by-n matrices A n are Toeplitz matrices with generating functions f that
are nonnegative even functions. More precisely, the matrices A n are constant along their diagonals
with their diagonal entries given by the Fourier coefficients of f :
\Gamma-
Since f are even functions, we have [A n are symmetric.
In [10, pp.64-65], it is shown that the eigenvalues - j lie in the range of f('), i.e.
min
Moreover, we also have
lim
f(') and lim
Department of Mathematics, Chinese University of Hong Kong, Shatin, Hong Kong. Research supported by
HKRGC grants no. CUHK 178/93E and CUHK 316/94E.
y Institute of Applied Mathematics, Chinese Academy of Science, Beijing, People's Republic of China.
z Department of Mathematics, Chinese University of Hong Kong, Shatin, Hong Kong.
Consequently, if f(') is nonnegative and vanishes at some points ' 0 2 [\Gamma-], then the condition
number is unbounded as n tends to infinity, i.e. A n is ill-conditioned. In fact, if the
zeros of f are of order -, then instance [4].
Superfast direct methods for Toeplitz matrices have been developed around 1980. They can
solve n-by-n Toeplitz systems in O(n log 2 n) operations, see for instance [1]. However, their stability
properties for ill-conditioned Toeplitz matrices are still unclear. Iterative methods based on
the preconditioned conjugate gradient method were proposed in 1985, see [11, 13]. With circulant
matrices as preconditioners, the methods require O(n log n) operations per iteration. For Toeplitz
systems generated by positive functions, these methods have shown to converge superlinearly. How-
ever, circulant preconditioners in general cannot handle Toeplitz matrices generated by functions
with zeros, see the numerical results in x6. The band-Toeplitz preconditioners proposed in [4, 5]
can handle functions with zeros, but are restricted to the cases where the order of the zeros are
even numbers. Thus they are not applicable for functions like
Classical iterative methods such as the Jacobi or Gauss-Seidel methods are also not applicable
when the generating functions have zeros. Since lim n!1 the convergence factor is
expected to approach 1 for large n. In [8, 9], Fiorentino and Serra proposed to use multigrid method
coupled with Richardson method as smoother for solving Toeplitz systems. Their numerical results
show that the multigrid method gives very good convergence rate for Toeplitz systems generated by
nonnegative functions. The cost per iteration of the multigrid method is of O(n log n) operations
which is the same as the preconditioned conjugate gradient methods.
However, in [8, 9], the convergence of the two-grid method (TGM) on first level is only proved for
the so-called band - matrices. These are band matrices that can be diagonalized by sine transform
matrices. A typical example is the 1-dimensional discrete Laplacian matrix diag[\Gamma1; 2; \Gamma1]. In
general, - matrices are not Toeplitz matrices and vice versa. The proof of convergence of the TGM
for Toeplitz matrices was not given there.
From the computational point of view, the matrix on the coarser grid in TGM is still too
expensive to invert. One therefore usually does not use TGM alone but instead applies the idea of
TGM recursively on the coarser grid to get down to the coarsest grid. The resulting method is the
full multigrid method (MGM). We remark that the convergence of MGM for Toeplitz matrices or
for - matrices was not discussed in [8, 9].
In this paper, we consider the use of MGM for solving ill-conditioned Toeplitz systems. Our
interpolation operator is constructed according to the position of the first non-zero entry on the
first row of the given Toeplitz matrix and is different from the one proposed by Fiorentino and
Serra [8, 9]. We show that for a class of ill-conditioned Toeplitz matrices which includes weakly
diagonally dominant Toeplitz matrices, the convergence factor of TGM with our interpolation
operator is uniformly bounded below 1 independent of n. We also prove that for this class of
Toeplitz matrices, the convergence factor of MGM with V -cycles will be level-dependent. One
standard way of removing the level-dependence is to use "better" cycles such as the F - or the
W -cycles, see [12]. We remark however that our numerical results show that MGM with V -cycles
already gives level-independent convergence. Since the cost per iteration is of O(n log n) operations,
the total cost of solving the system is therefore of O(n log n) operations.
We note that the class of functions that we can handle includes functions with zeros of order 2 or
less and also functions such as which cannot be handled by band-Toeplitz preconditioners
proposed in [4, 5]. We will also give examples of functions that can be handled by multigrid method
with our interpolation operator but not with the interpolation operator proposed in [8].
The paper is organized as follows. In x2, we introduce the two-grid method and the full multigrid
method. In x3, we analyze the convergence rate of two-grid method. We first establish in x3.1 the
convergence of two-grid method on the first level for the class of weakly diagonally dominant
Toeplitz matrices. The interpolation operator for these matrices can easily be identified. Then in
x3.2, we consider a larger class of Toeplitz matrices which are not necessarily diagonally dominant.
The convergence of full multigrid method is studied in x4 by establishing the convergence of the two-
grid method on the coarser levels. In x5, we give the computational cost of our method. Numerical
results are given in x6 to illustrate the effectiveness of our method. Finally, concluding remarks are
given in x7.
Given a Toeplitz system A n we define a sequence of sub-systems on different
levels:
Here q is the total number of levels with being the finest level. Thus for
and are just the size of the matrix A m . We denote the interpolation and
restriction operators by I m
respectively. We will
choose
I m+1
The coarse grid operators are defined by the Galerkin algorithm, i.e.
Thus, if A m is symmetric and positive definite, so is A m+1 . The smoothing operator is denoted
by Typical smoothing operators are the Jacobi, Gauss-Seidel and Richardson
iterations, see for instance [3]. Once the above components are fixed, a multigrid cycling procedure
can be set up. Here we concentrate on the V -cycle scheme which is given as follows, see [3, p.48].
procedure
then u q :=
begin do i := 1 to - 1
d m+1 := I m+1
e m+1
do
Here I nm is the nm -by-n m identity matrix. If we set 2, the resulting multigrid method is the
two-grid method (TGM).
3 Convergence of TGM for Toeplitz Matrices
In this section, we discuss the convergence of TGM for Toeplitz matrices. We first give an estimate
of the convergence factor for Toeplitz matrices that are weakly diagonally dominant. Then we
extend the results to a larger class of Toeplitz matrices.
Let us begin by introducing the following notations. We say A
positive (respectively semi-positive) definite matrix. In particular, A ? 0 means that A
is positive definite. The spectral radius of A is denoted by ae(A). For A ? 0, we define the following
inner products which are useful in the convergence analysis of multigrid methods, see [12, p.77-78]:
Here h\Delta; \Deltai is the Euclidean inner product. Their respective norms are denoted by k 2.
Throughout this section, we denote the fine and coarse grid levels of the TGM as the h- and
H-levels respectively. For smoothing operator, we consider the damped-Jacobi iteration, which is
given by
see [3, p.10]. The following theorem shows that kG h k 1 - 1 if ! is properly chosen.
Theorem 1 ([12, p.84]) Suppose A ? 0. Let ff be such thatff
Then
satisfies
Inequality (6) is called the smoothing condition. We see from the theorem that the damped-Jacobi
method (4) with
For a Toeplitz matrix A generated by an even function f , we see from (1) that ae(A) -
Moreover, diag(A) is just a constant multiple of the identity matrix. Thus it
is easy to find an ff that satisfies (5). In applications where f is not known a priori, we can estimate
ae(A) by the Frobenius norm or matrix 1-norm of A. The estimate can be computed in O(n)
operations.
For TGM, the correction operator is given by
with the convergence factor given by k(G h are the
numbers of pre- and post-smoothing steps in the MGM algorithm in x2. For simplicity, we will
consider only The other cases can be established similarly as we have kG h k 1 - 1.
Thus the convergence factor of our TGM is given by kG h T h k 1 . The following theorem gives a general
estimate on this quantity.
Theorem 2 ([12, p.89]) Let chosen such that G h satisfies the
smoothing condition (6), i.e.
Suppose that the interpolation operator I h
H has full rank and that there exists a fi ? 0 such that
min
and the convergence factor of the h-H two-level TGM satisfies
r
Inequality (7) is called the correcting condition. From Theorems 1 and 2, we see that if ff is
chosen according to (5) and that the damped-Jacobi method is used as the smoother, then we only
have to establish (7) in order to get the convergence results. We start with the following class of
matrices.
3.1 Weakly Diagonally Dominant Toeplitz Matrices
In the following, we write n-by-n Toeplitz matrix A generated by f as its j-th
diagonal as a j , i.e. [A] is the j-th Fourier coefficient of f . Let ID be the class of Toeplitz
matrices generated by functions f that are even, nonnegative and satisfy
a
Given a matrix A 2 ID, let l be the first non-zero index such that a l 6= 0. If a l ! 0, we define the
interpolation operator as
I h
I l2 I l 1
I l
I l .2 I l .C C C C C C C C C C A
Here I l is the l-by-l identity matrix. If a l ? 0, we define the interpolation operator as
I h
I l
I l
I l
I l .
Theorem 3 Let A 2 ID and l be the first non-zero index where a l 6= 0. Let the interpolation
operator be chosen as in (10) or (11) according to the sign of a l . Then there exists a fi ? 0
independent of n such that (7) holds. In particular, the convergence factor of TGM is bounded
uniformly below 1 independent of n.
Proof: We will prove the theorem for the case a l ! 0. The proof for the case a l ? 0 is similar
and is sketched at the end of this proof. We first assume that n
according to (10), we have kl. For any e
where
For ease of indexing, we set e
We note that with I h
H as defined in (10) and the norm k \Delta k 0 in (3), we have
l
Thus (7) is proved if we can bound the right hand side above by fihe h ; Ae h i for some fi independent
of e h . To do so, we observe that for the right hand side above, we have
a 0
l
l
\Gammae 2il+j e (2i\Gamma1)l+j
- a 0
l
\Gammae 2il+j e (2i\Gamma1)l+j
- a 0
l
is the n h -by-n h Toeplitz matrix generated by 1 \Gamma cos l'. Thus
min
Hence to establish (7), we only have to prove that
for some fi independent of e h . To this end, we note that the n h -by-n h matrix A is generated by
a j cos j':
But by (9),
a
In particular, by (1)
Thus, by (12), we then have
2a l
2a l
Hence (7) holds with
Next we consider the case where n h is not of the form (2k+1)l. In this case, we let
We then embed the vectors e h and e H into longer vectors
e ~ h and e ~
H of size n ~ h and n ~
H by zeros. Then since
~
he ~
we see that the conclusion still holds.
We remark that the case where a l ? 0 can be proved similarly. We only have to replace the
function above by (1 Since in this case, f n h (') - 2a l (1 we then have
From this, we get (15) and hence (7) with fi defined as in (16).
3.2 More General Toeplitz Matrices
The condition on ID class matrices is too strong. For example, it excludes the matrix
However, from (12) and (13), we see that (7) can be proved if we can find a positive number fi
independent of n and an integer l such that
Since by (14) and (17), we see that (18) holds for any matrices B in ID, we immediately have the
following corollary.
Corollary 1 Let A be a symmetric positive definite Toeplitz matrix. If there exists a matrix B 2 ID
such that A - B. Then (7) holds provided that the interpolation operator for A is chosen to be the
same as that for B.
More generally, we see by (1) that if the generating function f of A satisfies
min
for some l, then (18) holds. Thus we have the following theorem.
Theorem 4 Let A be generated by an even function f that satisfies (19) for some l. Let the
interpolation operator be chosen as in (10) or (11) according to the sign of a l . Then (7) holds. In
particular, the convergence factor of TGM is uniformly bounded below 1 independent of the matrix
size.
It is easy to prove that (19) holds for any even, nonnegative functions with zeros that are of
order 2 or less. As an example, consider
and T n [1 \Gamma cos '] 2 ID, it follows from Theorem 4 that if the interpolation operator for A is chosen
to be the same as that for T n [1 \Gamma cos '], the convergence factor of the resulting TGM will be
bounded uniformly below 1. We note that T n [1 \Gamma cos '] is just the 1-dimensional discrete Laplacian:
diag[\Gamma1; 2; 1]. Our interpolation operator here is the same as the usual linear interpolation operator
used for such matrices, see [3, p.38]. However, we remark that the matrix is a dense
matrix.
As another example, consider the dense matrix T n [j'j]. Since -j'j - ' 2 on [\Gamma-], we have by
Hence T n [j'j] can also be handled by TGM with the same linear interpolation operator used for
Convergence Results for Full Multigrid Method
In TGM, the matrix A H on the coarse grid is inverted exactly. From the computational point of
view, it will be too expensive. Usually, A H is not solved exactly, but is approximated using the
TGM idea recursively on each coarser grid until we get to the coarsest grid. There the operator
is inverted exactly. The resulting algorithm is the full multigrid method (MGM). In x3, we have
proved the convergence of TGM for the first level. To establish convergence of MGM, we need to
prove the convergence of TGM on coarser levels.
Recall that on the coarser grid, the operator A H is defined by the Galerkin algorithm (2), i.e.
h A h I h
H . We note that if n will be a block-Toeplitz-
Toeplitz-block matrix and the blocks are l-by-l Toeplitz matrices. In particular, if l = 1, then A H
is still a Toeplitz matrix. However, if n h is not of the form (2k will be a sum of a
block-Toeplitz-Toeplitz-block matrix and a low rank matrix (with rank less than or equal to 2l).
We will only consider the case where n j. For then on each level
. Hence the main diagonals of the coarse-grid operators A m ,
will still be constant. Recall that from the proof of Theorem 3 that (18) implies (7).
We now prove that if (18) holds on a finer level, it holds on the next coarser level when the same
interpolation operator is used.
Theorem 5 Let a h
0 and a H
0 be the main diagonal entries of A h and A H respectively. Let the
interpolation operator I h
H be defined as in (10) or (11). Suppose that
A h - a hfi h
for some fi h ? 0 independent of n. Then
with
a
Proof: We first note that if we define the (n H
I l I l
I l I lC A ; (24)
then there exists a permutation matrix P such that
I h
I nH
(cf (10) and (11)). Moreover, for the same permutation matrix P , we have
I nH \UpsilonK t
I nH+l
By (2) and (21), we have
h A h I h
a hfi h I H
But by (25) and (26), we have
a hfi h I H
I nH \UpsilonK t
I nH+l
!/
I nH
By the definition of K in (24), we have
a hfi h
Combining this with (28), we get
a hfi h I H
Hence (27) implies (22) with (23).
Recall by (5) that we can choose ff h such that
Notice that K t K - I n h and therefore
h A h I h
I h
Thus on the coarser level, we can choose ff H as
According to (8), (30) and (23), we see that
s
s
ff h a H
s
Recursively, we can extend this result from the next coarser-level to the q-th level and hence obtain
the level-dependent convergence of the MGM:
s
s
We remark that this level-dependent result is the same as that of most MGM, see for instance
[12, 2]. One standard way to overcome level-dependent convergence is to use "better" cycles such
as the F - or W -cycles, see [12]. We note however that our numerical results in x6 shows that MGM
with V -cycles already gives level-independent convergence.
We remark that we can prove the level-independent convergence of MGM in a special case.
Theorem 6 Let f(') be such that
for some integer l and positive constants c 1 and c 2 . Then for any 1 - m - q,
r
Proof: From (31), we have
Recalling the Galerkin algorithm (2) and using (29) recursively, we then have
By the right-hand inequality and (18), we see that
c 1 a m:
and hence by the left hand side of (32)
Therefore by the definition of ff in (5), we see that
c 2 a m:
According to (8), we then conclude that
s
r
As an example, we see that MGM can be applied to T n [' 2 ] with the usual linear interpolation
operator and the resulting method will be level-independent.
Computational Cost
Let us first consider the case where j. Then on each level, n
some k. From the MGM algorithm in x2, we see that if we are using the damped-Jacobi method
(4), the pre-smoothing and post-smoothing steps become
Thus the main cost on each level depends on the matrix-vector multiplication A m y for some vector
y. If we are using one pre-smoothing step and one post-smoothing step, then we require two such
matrix vector multiplications - one from the post-smoothing and one from the computation of the
residual. We do not need the multiplication in the pre-smoothing step since the initial guess u m is
the zero vector.
On the finest level, A is a Toeplitz matrix. Hence Ay can be computed in two 2n-length FFTs,
see for instance [13]. If l = 1, then on each coarser level, A m will still be a Toeplitz matrix. Hence
A m y can be computed in two 2nm -length FFTs. When l ? 1, then on the coarser levels, A m will
be a block-Toeplitz-Toeplitz-block matrix with l-by-l Toeplitz sub-blocks. Therefore A m y can also
be computed in roughly the same amount of time by using 2-dimensional FFTs. Thus the total
cost per MGM iteration is about eight 2n-length FFTs.
In comparison, the circulant-preconditioned conjugate gradient methods require two 2n-length
FFTs and two n-length FFTs per iteration for the multiplication of Ay and C \Gamma1 y respectively. Here
C is the circulant preconditioner, see [13]. The band-Toeplitz preconditioned conjugate gradient
methods require two 2n-length FFTs and one band-solver where the band-width depends on the
order of the zeros, see [4]. Thus the cost per iteration of using MGM is about 8/3 times as that
required by the circulant preconditioned conjugate gradient methods and 4 times of that required
by the band-Toeplitz preconditioned conjugate gradient methods.
Next we consider the case when n is not of the form 1)l. In that case, on the coarser level,
A m will no longer be a block-Toeplitz-Toeplitz-block matrix. Instead it will be a sum of such a
matrix and a low rank matrix (with rank less than 2l). Thus the cost of multiplying A m y will be
increased by O(ln).
6 Numerical Results
In this section, we apply the MGM algorithm in x2 to ill-conditioned real symmetric Toeplitz
systems A n We choose as solution a random vector u such that 0 - 1. The right hand
side vector b is obtained accordingly. As smoother, we use the damped-Jacobi method (4) with
chosen as a 0 =max f(') for pre-smoother and post-smoother. We use one
pre-smoothing and one post-smoothing on each level.
The zero vector is used as the initial guess and the stopping criterion is kr
where r j is the residual vector after j iterations. In the following tables, we give the number
of iterations required for convergence using our method, see column under M . For comparison,
we also give the number of iterations required by the preconditioned conjugate gradient method
with no preconditioner (I), the Strang (S) circulant preconditioner, the T. Chan (C) circulant
preconditioner and also the band (B) preconditioners, see [6, 7, 4]. The double asterisk * signifies
more than 200 iterations are needed.
For the first example, we consider functions with single zero at the point The functions
we tried are We note that T n
ID. However, we have Therefore, according
to Corollary 1, we can use the interpolation operator (10) with

Table

1: Number of Iterations for Functions with Single Zero.
Next we consider two functions with jumps and a single zero at the point
and
We note that both matrices T n [J j (')], are not in ID. However, since J j (') -
for all ', we can still use the interpolation operator defined by 1 \Gamma cos ' for both T n [J j (')]. We
remark that for J 2 ('), since the zero is not of even order, band circulant preconditioners cannot be
constructed.

Table

2: Number of Iterations for Functions with Jumps.
Finally, we consider two functions with multiple zeros. They are
and
ID. But we
note that ' 2 (- both matrices can use the interpolation operator
in (10) with l = 2. In particular, our interpolation operator will be different from that proposed
in [8], which in this case will use the interpolation operator in (10) with l = 1. Their resulting
MGM converges very slowly with convergence factor very close to 1 (about 0.98 for both functions
when 64). For comparison, we list in Table 3, the number of MGM iterations required by such
interpolation operator under column F .

Table

3: Number of Iterations for Functions with Multiple Zeros.
7 Concluding Remarks
We have shown that MGM can be used to solve a class of ill-conditioned Toeplitz matrices. The
resulting convergence rate is linear. The interpolation operator depends on the location of the first
non-zero diagonals of the matrices and its sign.
Here we have only proved the convergence of multigrid method with damped-Jacobi as smoothing
operator. However, our numerical results show that multigrid method with some other smoothing
operators, such as the red-black Jacobi, block-Jacobi and Gauss-Seidel methods, will give better
convergence rate. As an example, for the function convergence factors of
MGM with the point- and block-Jacobi methods as smoothing operator are found to be about 0:71
and 0:32 respectively for 64 - n - 1024.

Acknowledgment

We will like to thank Prof. Tony Chan and Dr. J. Zou for their valuable comments.



--R

Superfast Solution of Real Positive Definite Toeplitz Systems
Convergence Estimates for Multigrid Algorithms without Regularity Assumptions

Toeplitz Preconditioners for Toeplitz Systems with Nonnegative Generating Func- tions
Fast Toeplitz Solvers Based on Band-Toeplitz Preconditioner
Toeplitz Equations by Conjugate Gradients with Circulant Precondi- tioner
An Optimal Circulant Preconditioner for Toeplitz Systems
Multigrid Methods for Toeplitz Matrices
Multigrid Methods for Symmetric Positive Definite Block Toeplitz Matrices with Nonnegative Generating Functions

Linear and Nonlinear Deconvolution Problems
in Multigrid Methods
A Proposal for Toeplitz Matrix Calculations
--TR

--CTR
S. Serra Capizzano , E. Tyrtyshnikov, How to prove that a preconditioner cannot be superlinear, Mathematics of Computation, v.72 n.243, p.1305-1316, July
toeplitz matrices;multigrid method