The Best Circulant Preconditioners for Hermitian Toeplitz Systems.
In this paper, we propose a new family of circulant preconditioners for ill-conditioned Hermitian Toeplitz systems A x=  b. The preconditioners are constructed by convolving the generating function f of A with the generalized Jackson kernels. For an n-by-n Toeplitz matrix A, the construction of the preconditioners requires only the entries of A and does not require the explicit knowledge of f. When f is a nonnegative continuous function with a zero of order 2p, the condition number of A is known to grow as O(n2p). We show, however, that our preconditioner is positive definite and the spectrum of the preconditioned matrix is uniformly bounded except for at most 2p+1 outliers. Moreover, the smallest eigenvalue is uniformly bounded away from zero. Hence the conjugate gradient method, when applied to solving the preconditioned system, converges linearly. The total complexity of solving the system is therefore of O(n log n) operations. In the case when f is positive, we show that the convergence is superlinear. Numerical results are included to illustrate the effectiveness of our new circulant preconditioners.
Introduction
An n-by-n matrix A n with entries a ij is said to be Toeplitz if a a i\Gammaj . Toeplitz systems
of the form A n occur in a variety of applications in mathematics and engineering
E-mail: rchan@math.cuhk.edu.hk. Department of Mathematics, The Chinese University of Hong
Kong, Shatin, Hong Kong. Research supported in part by Hong Kong Research Grants Council Grant No.
CUHK 4207/97P and CUHK DAG Grant No. 2060143.
y
E-mail: mhyipa@hkusua.hku.hk. Department of Mathematics, The University of Hong Kong, Pokfu-
lam Road, Hong Kong.
z E-mail: mng@maths.hku.hk. Department of Mathematics, The University of Hong Kong, Pokfulam
Road, Hong Kong. Research supported in part by HKU CRCG grant no. 10201939.
[7]. In this paper, we consider the solution of Hermitian positive definite Toeplitz systems.
There are a number of specialized fast direct methods for solving such systems in O(n 2 )
operations, see for instance [22]. Faster methods requiring O(n log 2 n) operations have
also been developed, see [1].
Strang in [21] proposed using the preconditioned conjugate gradient method with circulant
matrices as preconditioners for solving Toeplitz systems. The number of operations
per iteration is of order O(n log n) as circulant systems can be solved efficiently by fast
Fourier transforms. Several successful circulant preconditioners have been introduced and
analyzed; see for instance [11, 5]. In these papers, the given Toeplitz matrix A n is assumed
to be generated by a generating function f , i.e., the diagonals a j of A n are given by the
Fourier coefficients of f . It was shown that if f is a positive function in the Wiener class
(i.e., the Fourier coefficients of f are absolutely summable), then these circulant preconditioned
systems converge superlinearly [5]. However, if f has zeros, the corresponding
Toeplitz systems will be ill-conditioned. In fact, for the Toeplitz matrices generated by a
function with a zero of order 2p, their condition numbers grow like O(n 2p ), see [19]. Hence
the number of iterations required for convergence will increase like O(n p ), see [2, p.24].
Tyrtyshnikov [23] has proved that the Strang [21] and the T. Chan [11] preconditioners
both fail in this case.
To tackle this problem, non-circulant type preconditioners have been proposed, see
[6, 4, 18, 16]. The basic idea behind these preconditioners is to find a function g that
matches the zeros of f . Then the preconditioners are constructed based on the function
g. These approaches work when the generating function f is given explicitly, i.e., all
Fourier coefficients
of f are available. However, when we are given only a finite
n-by-n Toeplitz system, i.e., only fa j g jjj!n are given and the underlying f is unknown,
then these preconditioners cannot be constructed. In contrast, most well-known circulant
preconditioners, such as the Strang and the T. Chan preconditioners, are defined using only
the entries of the given Toeplitz matrix. Di Benedetto in [3] has proved that the condition
numbers of the preconditioned matrices by sine transform preconditioners are uniformly
bounded. However, the preconditioners themselves may be singular or indefinite in general.
Our aim in this paper is to develop a family of positive definite circulant preconditioners
that work for ill-conditioned Toeplitz systems and do not require the explicit knowledge
of f , i.e., they require only fa j g jjj!n for an n-by-n Toeplitz system.
Our idea is based on the unified approach proposed in Chan and Yeung [9], where
they showed that circulant preconditioners can be derived in general by convolving the
generating function f with some kernels. For instance, convolving f with the Dirichlet
kernel D bn=2c gives the Strang preconditioner. They proved that for any positive 2-
periodic continuous function f , if C n is a kernel such that the convolution product C n   f
tends to f uniformly on [\Gamma-], then the corresponding circulant preconditioned matrix
n A n will have clustered spectrum. In particular, the conjugate gradient method will
converge superlinearly when solving the preconditioned system. This result turns the
problem of finding a good preconditioner to the problem of approximating f with C n   f .
Notice that D bn=2c   f , being the partial sum of f , depends solely on the first bn=2c
Fourier coefficients fa j g jjj!bn=2c of f . Thus the Strang preconditioner, and similarly for
other circulant preconditioners constructed through kernels, does not require the explicitly
knowledge of f .
In this paper, we construct our preconditioners by approximating f with the convolution
product K m;2r   f that matches the zeros of f and depends only on fa j g jjj!n . Here
K m;2r is chosen to be the generalized Jackson kernels, see [15]. Since K m;2r are positive
kernels, our preconditioners are positive definite for all n. In comparison, the Dirichlet
kernel D n is not positive and hence the Strang preconditioner is indefinite in general. We
will prove that if f has a zero of order 2p, then K m;2r   f matches the zero of f when
r ? p. Using this result, we can show that the spectra of the circulant preconditioned
matrices are uniformly bounded except for at most 2p outliers, and that their smallest
eigenvalues are bounded uniformly away from zero. It follows that the conjugate gradient
method, when applied to solving these circulant preconditioned systems, will converge
linearly. Since the cost per iteration is O(n log n) operations, see [7], the total complexity
of solving these ill-conditioned Toeplitz systems is of O(n log n) operations. In the case
when f is positive, we show that the spectra of the preconditioned matrices are clustered
around 1 and thus the method converges superlinearly. The case where f has multiple
zeros is more involved and will be considered in a future paper.
This paper is an expanded version of the proceedings paper [10] where some of the
preliminary results were reported. Recently Potts and Steidl [17] have proposed skew-
circulant preconditioners for ill-conditioned Toeplitz systems. Their idea is also to use
convolution products that match the zeros of f to construct the preconditioners. In par-
ticular, they have used the generalized Jackson kernels and the B-spline kernels proposed
in [8] in their construction. However, in order to guarantee that the preconditioners are
positive definite, the position of the zeros of f is required which in general may not be
readily available. In contrast, our circulant preconditioners can be constructed without
any explicit knowledge of the zeros of f .
The outline of the paper is as follows. In x2, we give an efficient method for computing
the eigenvalues of the preconditioners. In x3 we show that K m;2r   f matches the zeros of f .
We then analyze the spectrum of the preconditioned matrices in x4. Numerical results are
given in x5 to illustrate the effectiveness of our preconditioners in solving ill-conditioned
Toeplitz systems. Concluding remarks are given in x6.
2 Construction of Circulant Preconditioners
Let C 2- be the space of all 2-periodic continuous real-valued functions. The Fourier
coefficients of a function f in C 2- are given by
a
\Gamma-
f(')e \Gammaik' d';
a \Gammak for all k. Let A n [f ] be the n-by-n Hermitian Toeplitz matrix with the
j)th entry given by a i\Gammaj , We will use C
2- to denote the space of all
nonnegative functions in C 2- which are not identically zero. We remark that the Toeplitz
matrices A n [f ] generated by f
are positive definite for all n, see [6, Lemma 1].
Conversely, if f 2 C 2- takes both positive and negative values, then A n [f ] will be non-
definite. In this paper, we only consider f
being positive definite
Hermitian Toeplitz matrices.
We say that ' 0 is a zero of f of order p if f(' 0 is the smallest positive
integer such that f (p) (' 0 ) 6= 0 and f (p+1) (') is continuous in a neighborhood of ' 0 . By
Taylor's theorem,
p!
for all ' in that neighborhood. Since f is nonnegative, c ? 0 and p must be even. We
remark that the condition number of A n [f ] generated by such an f grows like O(n p ), see
[19]. In this paper, we will consider f having a single zero. The general case where f has
multiple zeros is more complicated and will be considered in a future paper.
The systems A n [f will be solved by the preconditioned conjugate gradient
method with circulant preconditioners. It is well known that all n-by-n circulant matrices
can be diagonalized by the n-by-n Fourier matrix F n , see [7]. Therefore, a circulant matrix
is uniquely determined by its set of eigenvalues. For a given function f , we define the
circulant preconditioner C n [f ] to be the n-by-n circulant matrix with its j-th eigenvalue
given by
We note that C n [f
Hence the matrix-vector multiplication
[f ]v, which is required in each iteration of the preconditioned conjugate
gradient method, can be done in O(n log n) operations by fast Fourier transforms. Clearly
if f is a positive function, then C n [f ] is positive definite.
In the following, we will use the generalized Jackson kernel functions
sin( m'
to construct our circulant preconditioners. Here k m;2r is a normalization constant such
that
\Gamma- K m;2r 1. It is known that k m;2r is bounded above and below by constants
independent of m, see [15, p.57] or (11) below. We note that K m;2 (') is the Fej'er kernel
and K m;4 (') is the Jackson kernel [15, p.57].
For any m, the Fej'er kernel K m;2 (') can be expressed as
where
see for instance [9]. Note that
\Gamma- K m;2
1. By (2), we see that K m;2r (')
is the r-th power of K m;2 (') up to a scaling. Hence we have
where the coefficients b (m;2r)
k can be obtained by convolving the vector (b (m;2)
times and this can be done by fast Fourier transforms, see
[20, pp.294-296]. Thus the cost of computing the coefficients fb (m;2r)
k g for all jkj -
is of order O(rm log m) operations. In order to guarantee that
\Gamma- K m;2r
can normalize b (m;2r)
0 to 1=(2-) by dividing all coefficients b (m;2r)
k by 2-b (m;2r)
The convolution product of two arbitrary functions
defined as
(g   h)(') j
\Gamma-
When we are given an n-by-n Toeplitz matrix A n [f ], our proposed circulant preconditioner
is
By (3) and (4), since
a k e ik' , the convolution product of K m;2r   f is given by
(K m;2r
a k b (m;2r)
where
0; otherwise:
Clearly, K m;2r   f depends only on a k for jkj ! n, i.e., only on the entries of the
given n-by-n Toeplitz matrix A n [f ]. Notice that by (1), to construct our preconditioner
we only need the values of K m;2r   f at 2-j=n for n. By (6), these
values can be obtained by taking one fast Fourier transform of length n. Thus the cost of
constructing C n [K m;2r   f ] is of O(n log n) operations.
We remark that the Strang [21] and the T. Chan circulant preconditioners [11] for
are just equal to C n [D bn=2c   f ] and C n [K n;2   f ] respectively where D bn=2c is the
Dirichlet kernel and K n;2 (') is the Fej'er kernel, see [9].
3 Properties of the Kernel K m;2r
In this section, we study some properties of K m;2r in order to see how good the approximation
of f by K m;2r   f will be. These properties are useful in the analysis of our circulant
preconditioners in x4. First we claim that our preconditioners are positive definite.
Lemma 3.1
2- . The preconditioner C n [K m;2r   f ] is positive definite for all
positive integers m, n and r.
Proof: By (2), K m;2r (') is positive except at
2- is nonnegative and not identically zero, the function
(K m;2r   f)(') j
\Gamma-
is clearly positive for all ' 2 [\Gamma-]. Hence by (1), the preconditioners C n [K m;2r   f ] are
positive definite.
In the following, we will use ' to denote the function ' defined on the whole real line R.
For clarity, we will use ' 2- to denote the periodic extension of ' on [\Gamma-], i.e. ' 2-
1 below). It is clear that ' 2- 2 C
2- . We
first show that K m;2r   ' 2p
2- matches the order of the zero of ' 2p
2- at
Lemma 3.2 Let p and r be positive integers with r ? p. Then
\Gamma-
where
Proof: The first two equalities in (7) are trivial by the definition of ' 2- . For the last
equality, since '=- sin('=2) - '=2 on [0; -], we have by (2)
\Gamma-
Z -sin 2r
Z m-
0sin 2r u
aeZ 1sin 2r u
Z 1sin 2r u
oe
aeZ 1u 2p du
Z 11
oe
On the other hand, we also have
\Gamma-
Z -sin 2r
Z 1sin 2r u
By setting
\Gamma-
Putting (11) back into (9) and (10), we then have (8).
We remark that using the same arguments as in (10), we can show that
i.e., the T. Chan preconditioner does not match the order of the zeros of ' 2p at
when p - 1. We will see in x5 that the T. Chan preconditioner does not work for Toeplitz
matrices generated by functions with zeros of order greater than or equal to 2.
Next we estimate (K m;2r   ' 2p
2- )(OE) for OE 6= 0. In order to do so, we first have to replace
the function ' 2p
2- in the convolution product by ' 2p defined on R.
Lemma 3.3 Let p be a positive integer. Then
\Theta

and
Proof: To prove (13), we first claim that
By the definition of
2- , we have (see Figure 1)
For \Gamma-=2], we have
For \Gamma-=2], we have
Thus we have (16).
By (16), we see that
\Gamma-
\Gamma-
\Theta K m;2r   ' 2p ('
Similarly, we also have
5-

Thus, we have (13).
To prove (14), we just note that
As for (15), we have
With Lemmas 3.2 and 3.3, we show that K m;2r   ' 2p
2- and ' 2p
are essentially the same
away from the zero of ' 2p
2- .

Figure

1: The functions
Theorem 3.4 Let p and r be positive integers with r ? p and dn=re. Then there
exist positive numbers ff and fi independent of n such that for all sufficiently large n,
Proof: We see from Lemma 3.3 that for different values of OE, (K m;2r   ' 2p
2- )(OE) can be
replaced by different functions. Hence, we proceed the proof for different ranges of values
of OE.
We first consider OE 2 [-=n; -=2]. By the binomial expansion,
\Gamma-
\Gamma-
For odd k,
\Gamma- K m;2r (t)t k
OE \Gamma2k
\Gamma-
K m;2r (t)t 2k dt:
By (7), we then have \Gamma K m;2r   ' 2p
where by (8), c k;2r are bounded above and below by positive constants independent of m
by (5), -=r -m=n - OEm, we have
Thus by (18),
Hence by (14), (17) follows for OE 2 [-=n; -=2].
The case with OE 2 [\Gamma-=2; \Gamma-=n] is similar to the case where OE 2 [-=n; -=2].
Next we consider the case OE 2 [-=2; -]. Note that
\Theta K m;2r   ' 2p
\Gamma-
\Gamma-
dt
where
is a degree 4p polynomial without the constant term. By (7), we have
\Gamma-
c 2j;2r
Thus
\Theta

c 2j;2r
2- for OE 2 [-=2; -], we have
\Theta

which is clearly bounded independent of n. For the lower bound, we use the fact that
2- for OE 2 [-=2; -] in (19), then we have
\Theta K m;2r   ' 2p
c 2j;2r
c 2j;2r
for sufficiently large n (and hence large m), the last expression is bounded uniformly
from below say by - 2p =2. Combining (20), (21) and (15), we see that (17) holds
for OE 2 [-=2; -] and for n sufficiently large.
The case where OE 2 [\Gamma-; \Gamma-=2] can be proved in a similar way as above.
Using the fact that
\Gamma-
we obtain the following corollary which deals with functions having a zero at fl 6= 0.
Corollary 3.5 Let fl 2 [\Gamma-], p and r be positive integers with r ? p and
Then there exist positive numbers ff and fi, independent of n, such that for all sufficiently
large n,
Now we can extend the results in Theorem 3.4 to any functions in C
2- with a single
zero of order 2p.
Theorem 3.6 Let f 2 C
2- and have a zero of order 2p at fl 2 [\Gamma-]. Let r ? p be any
integer and dn=re. Then there exist positive numbers ff and fi, independent of n,
such that for all sufficiently large n,
Proof: By the definition of zeros (see x2),
2- g(') for some positive continuous
function g(') on [\Gamma-]. Write
(K m;2r   f) (OE)
\Deltag(OE)
Clearly the last factor is uniformly bounded above and below by positive constants. By
Corollary 3.5, the same holds for the second factor when -=n As for the
first factor, by the Mean Value Theorem for integrals, there exists a i 2 [\Gamma-] such that
Hence
where g min and g max are the minimum and maximum of g respectively. Thus the theorem
follows.
So far we have considered only the interval -=n -
now show that the convolution product K m;2r   f matches the order of the zero of f at
the zero of f .
Theorem 3.7 Let f 2 C
2- and have a zero of order 2p at fl 2 [\Gamma-]. Let r ? p be any
integer and dn=re. Then for any jOE \Gamma flj -=n, we have
(K m;2r   f)
Proof: We first prove the theorem for the function
. By the binomial theorem,
\Gamma-
\Gamma-
Since
\Gamma- K m;2r (t)t j we have for jOEj -=n,
\Gamma-
\Gamma-
By (7), (8) and (5), we then have
Hence by (14),
On the other hand, from (22), (8) and (5),
we have
(K m;2r   ' 2p )(OE) -
\Gamma-
O
Hence by (14) again,
Thus the theorem holds for
2- .
In the general case where
2- g(') for some positive function g 2 C 2- , by
the Mean Value Theorem for integrals, there exists a i 2 [\Gamma-] such that
(K m;2r
Hence
min \Delta (K m;2r   ' 2p
for all OE 2 [\Gamma-]. Here g min and g max are the minimum and maximum of g respectively.
From the first part of the proof, we already see that (K m;2r   ' 2p
is of O
1=n 2p
for all jOE \Gamma flj -=n, hence the theorem follows.
4 Spectral Properties of the Preconditioned Matrices
4.1 Functions with a Zero
In this subsection, we analyze the spectra of the preconditioned matrices when the generating
function has a zero. We will need the following lemma.
Lemma 4.1 [4, 16] Let f 2 C
2- . Then A n [f ] is positive definite for all n. Moreover if
2- is such that 0 ! ff - f=g - fi for some constants ff and fi, then for all n,
x   A n [g]x
Next, we have our first main theorem which states that the spectra of the preconditioned
matrices are essentially bounded.
Theorem 4.2 Let f 2 C
2- and have a zero of order 2p at fl. Let r ? p and
Then there exist positive numbers ff ! fi, independent of n, such that for all sufficiently
large n, at most 2p
are outside the interval [ff; fi].
Proof: For any function g 2 C 2- , we let ~
[g] to be the n-by-n circulant matrix with the
j-th eigenvalue given by
there is at most one j such that j2-j=n \Gamma flj ! -=n, by (1),
~
is a matrix of rank at most 1.
By assumption, positive function g in C 2- . We
use the following decomposition of the Rayleigh quotient to prove the theorem:
x   A n [f ]x
x   A n
sin 2p
x
x   A n
sin 2p
x
x   ~
sin 2p
x
x   ~
sin 2p
x
x   ~
x   ~
x   ~
x   ~
We remark that by Lemma 4.1 and the definitions (1) and (23), all matrices in the factors
in the right hand side of (24) are positive definite.
As g is a positive function in C 2- , by Lemma 4.1, the first factor in the right hand
side of (24) is uniformly bounded above and below. Similarly, by (23), the third factor is
also uniformly bounded. The eigenvalues of the two circulant matrices in the fourth factor
differ only when j2-j=n \Gamma flj -=n. But by Theorem 3.6, the ratios of these eigenvalues
are all uniformly bounded when n is large. The eigenvalues of the two circulant matrices
in the last factor differ only when j2-j=n -=n. But by Theorem 3.7, their ratios
are also uniformly bounded.
It remains to handle the second factor. Define s 2p (') j sin 2p ( '\Gammafl
i.e., s 2p (') is a p-th degree trigonometric polynomial in '. Recall that for any function
the convolution product of the Dirichlet kernel D n with h is just
equal to the nth partial sum of h, i.e., (D n
j=\Gamman b j e ij' . Hence for n - 2p,
(D bn=2c   s 2p
Since C n [D bn=2c   s 2p (')] is the Strang preconditioner for A n [s 2p (')], see [9], C n [s 2p (')]
will be the Strang preconditioner for A n [s 2p (')] when n - 2p. As s 2p (') is a p-th degree
trigonometric polynomial, A n [s 2p (')] is a band Toeplitz matrix with half bandwidth p+ 1.
Therefore when n - 2p, by the definition of the Strang preconditioner,
R
where R p is a p-by-p matrix, see [21]. Thus A n [s 2p
where the n-by-n
matrix R n is of rank at most 2p + 1.
Putting this back into the numerator of the second factor in (24), we have
x   A n [f ]x
x   ~
x   ~
x   ~
x   ~
x   ~
x   A n [f ]x
x   R n x
Notice that for all sufficiently large n, except for the last factor, all factors above are
uniformly bounded below and above by positive constants. We thus have
x   A n [f ]x
when n large, where
Hence for large n,
x
If R n has q positive eigenvalues, then by Weyl's theorem [13, p.184], at most q eigenvalues
of C n [K m;2r   f are larger than ff max . By using a similar argument, we can prove
that at most 2p are less than ff min . Hence the
theorem follows.
Finally we prove that all the eigenvalues of the preconditioned matrices are bounded
from below by a constant independent of n. Hence the computational cost for solving this
class of n-by-n Toeplitz systems will be of O(n log n) operations.
Theorem 4.3 Let f 2 C
2- and have a zero of order 2p at fl. Let r ? p and
Then there exists a constant c independent of n, such that for all n sufficiently large, all
eigenvalues of the preconditioned matrix C
are larger than c.
Proof: In view of the proof of Theorem 4.2, it suffices to get a lower bound of the
second Rayleigh quotient in the right hand side of (24). Equivalently, we have to get
an upper bound of ae(A \Gamma1
denotes the spectral radius and
We note that by the definition (23), ~
the zero matrix or is given by
F
diag
for some j such that j2-j=n \Gamma flj ! -=n. Thus
By Lemma 4.1, A \Gamma1
positive definite. Thus the matrix
A
is similar to the symmetric matrix
A \Gamma1=2
Hence we have
ae
A
A \Gamma1=2
A \Gamma1=2
A \Gamma1=2
A
By [6, Theorem 1], we have
Hence the last term in (26) is of
O(1).
It remains to estimate the first term in (26). According to (25), we partition A \Gamma1
as
A
are p-by-p matrices. Then by (25),
ae
A
'-
where the last equality follows because the 3-by-3 block matrix in the equation has vanishing
central column blocks. In [3, Theorem 4.3], it has been shown that R p , B 11 , B 13 and
all have bounded ' 1 -norms and ' 1 -norms. Hence using the fact that ae(\Delta) - k
we see that (27) is bounded and the theorem follows.
By combining Theorems 4.2 and 4.3, the number of preconditioned conjugate gradient
iterations required for convergence is of O(1), see [3]. Since each PCG iteration
requires O(n log n) operations (see [7]) and so is the construction of the preconditioner
(see x2), the total complexity of the PCG method for solving Toeplitz systems generated
by
2- is of O(n log n) operations.
4.2 Positive Functions
In this subsection, we consider the case where the generating function is strictly positive.
We note that the spectrum of A n [f ] is contained in [f min
are the minimum and maximum values of f , see [6, Lemma 1]. Since f min ? 0, A n [f
is well-conditioned. In [9], it was shown that for such f , the spectrum of C
f ]A n [f ] is clustered around 1 and the PCG method converges superlinearly. Recall that
is just the T. Chan circulant preconditioner. In the following, we generalize
this result to other generalized Jackson kernels. First, it is easy to show that
(K m;2r   f) (OE) - f max . Thus the whole spectrum of C
is contained in
[f min =f i.e. the preconditioned system is also well-conditioned. We now
show that its spectrum is clustered around 1.
Theorem 4.4 Let f 2 C 2- be positive. Then the spectrum of C
clustered around 1 for sufficiently large n. Here
Proof: We first prove that K m;2r   f converges to f uniformly on [\Gamma-]. For - ? 0,
be the modulus of continuity of f . It has the
property that
see [15, p.43].
By the uniform continuity of f , for each " ? 0, there exists a ffi ? 0 such that !(f; ffi) ! ".
\Gamma-
\Gamma-
\Gamma-
\Gamma-
\Gamma-
bounded by a constant independent of n (cf. the proof of
Lemma 3.2 for Therefore, K m;2r   f converges uniformly to f . By [9, Theorem
1], the spectrum of C
clustered around 1 for sufficiently large n.
As an immediate consequence, we can conclude that when f is positive and C n [K m;2r \Lambdaf
is used as the preconditioner, the PCG method converges superlinearly, see for instance
[5].
5 Numerical Experiments
In this section, we illustrate by numerical examples the effectiveness of the preconditioner
solving Toeplitz systems. For comparisons, we also test the Strang [21]
and the T. Chan [11] circulant preconditioners. In the following, m is set to dn=re.
Example 1: The first set of examples is on mildly ill-conditioned Toeplitz systems where
the condition numbers of the systems grow like O(n ' ) for some ' ? 0. They correspond
to Toeplitz matrices generated by functions having zeros of order ', see [19]. Because of
the ill-conditioning, the conjugate gradient method will converge slowly and the number
of iterations required for convergence grows like O(n '=2 ) [2, p.24]. However, we will see
that using our preconditioner C n [K m;2r   f ] with 2r ? ', the preconditioned system will
converge linearly, i.e., the number of iterations required for convergence is independent of
n.
We solve Toeplitz systems A n [f by the preconditioned conjugate gradient
method for twelve nonnegative test functions. Since the functions are nonnegative, the
so generated are all positive definite. We remark that if f takes negative values, then
A n [f ] will be non-definite for large n. As mentioned in x2, the construction of our preconditioners
for an n-by-n Toeplitz matrix requires only the n diagonal entries fa j g jjj!n
of the given Toeplitz matrix. knowledge of f is required. In the tests, the
right-hand side vectors b are formed by multiplying random vectors to A n [f ]. The initial
guess is the zero vector and the stopping criteria is jjr q jj 2 =jjr 0 jj is the
residual vector after q iterations.

Tables

1-4 show the numbers of iterations required for convergence for different choices
of preconditioners. In the table, I denotes no preconditioner, S is the Strang preconditioner
[21], K m;2r are the preconditioners from the generalized Jackson kernel K m;2r defined in
(2) and is the T. Chan preconditioner [11]. Iteration numbers more than 3,000
are denoted by "y". We note that S in general is not positive definite as the Dirichlet
kernel D n is not positive, see [9]. When some of its eigenvalues are negative, we denote
the iteration number by "-" as the PCG method does not apply to non-definite systems
and the solution thus obtained may be inaccurate.
The first two test functions in Table 1 are positive functions and therefore correspond
to well-conditioned systems. Notice that the iteration number for the non-preconditioned
systems tends to a constant when n is large, indicating that the convergence is linear.
In this case, we see that all preconditioners work well and the convergence is fast, see
Theorem 4.4 and [9].
I

Table

1: Numbers of iterations for well-conditioned systems.
The four test functions in Table 2 are nonnegative functions with single or multiple
zeros of order 2 on [\Gamma-]. Thus the condition numbers of the Toeplitz matrices are
growing like O(n 2 ) and hence the numbers of iterations required for convergence without
using any preconditioners is increasing like O(n). We see that for these functions, the T.
Chan preconditioner does not work. This is to be expected from the fact that the order of
does not match that of ' 2 at see (12). However, we see that K m;4 , K m;6
and K m;8 all work very well as predicted from our convergence analysis in x4.
When the order of the zero is 4, like the two test functions in Table 3, the condition
number of the Toeplitz matrices will increase like O(n 4 ) and the matrices will be very
ill-conditioned even for moderate n. We see from the table that both the Strang and the
T. Chan preconditioners fail. For the T. Chan preconditioner, the failure is also to be
expected from the fact that the order of K m;2   ' 4 does not match that of ' 4 at
(12). As predicted by our theory, K m;6 and K m;8 still work very well. The numbers of
iterations required for convergence are roughly constant independent of n.
In

Table

4, we test functions that our theory does not cover. The first two functions are
not differentiable at their zeros. The last two functions are functions with slowly decaying
Fourier coefficients. We found numerically that the minimum values of
I 36 79 170 362 753 1544 53 141 293 547 1113 2213
I

Table

2: Numbers of iterations for functions with order 2 zeros.
I
26 42 71 161 167 247 24 35 58 106 144 196
KN;4 15 17 20 24 26 26 15
KN;8

Table

3: Numbers of iterations for functions with order 4 zeros.
I
I

Table

4: Numbers of iterations for other functions.
and
jkj!1024jkj 0:5 +1 e ik' are approximately equal to 0.3862 and 0.4325 respectively. Hence
the last two test functions are approximately zero at some points in [\Gamma-]. Table 4
shows that the K m;2r preconditioners still perform better than the Strang and the T.
Chan preconditioners.
To further illustrate Theorems 4.2 and 4.3, we give in Figures 2 and 3 the spectra of
the preconditioned matrices for all five preconditioners for
We see that the spectra of the preconditioned matrices for K m;6 and K m;8 are in a small
interval around 1 except for one to two large outliers and that all the eigenvalues are well
separated away from 0. We note that the Strang preconditioned matrices in both cases
have negative eigenvalues and they are not depicted in the figures.
Example 2: In image restoration, because the blurring is an averaging processing, the
resulting matrix is usually strongly ill-conditioned in the sense that its condition number
grows exponentially with respect to its size n. In contrast, the condition numbers of the
mildly ill-conditioned matrices considered in Example 1 are increasing like polynomials
of n only. Regularization techniques have been used for some time in mathematics and
engineering to treat these strongly ill-conditioned systems. The idea is to restrict the
solution in some smooth function spaces [14]. This approach has been adopted in the
circulant preconditioned conjugate gradient method and is very successful when applied
to ground-based astronomy [7].
To illustrate the idea, we use a "prototype" image restoration problem given in [12].
Strang Preconditioner (has negative eigenvalues)
T. Chan Preconditioner
K m,4 Jackson Preconditioner
K m,6 Jackson Preconditioner
K m,8 Jackson Preconditioner

Figure

2: Spectra of preconditioned matrices for
Strang Preconditioner (has negative eigenvalues)
T. Chan Preconditioner
K m,4 Jackson Preconditioner
K m,6 Jackson Preconditioner
K m,8 Jackson Preconditioner

Figure

3: Spectra of preconditioned matrices for
Consider a 100-by-100 Toeplitz matrix A with (i; entries given by
ae 0; if
where
-oe
Blurring matrices of this form (called the truncated Gaussian blur) occur in many image
restoration contexts and are used to model certain degradations in the recorded image.
The condition number of A is approximately 2.3\Theta10 6 . Thus if no regularization is used,
the result obtained will be very inaccurate.
In our experiment, we solve the regularized least squares problem min x
as suggested in [12]. The problem is equivalent to the normal equations (ffI
which we solve by the preconditioned conjugate gradient method. We choose
the solution vector x with its entries given by
see [12], and then we compute noise vector is added to b where each component
of the noise vector is taken from a normal distribution with mean zero and standard
deviation . The stopping criteria is jjr q jj 2 =jjr 0 jj is the residual
vector after q iterations.
We choose the optimal regularization parameter ff   such that it minimizes the relative
error between the computed solution x(ff) of the normal equations and the original solution
x given in (28), i.e. ff   minimizes . By trial and error, it is found to be
\Gamma6 up to one digit of accuracy. The preconditioner we used for the normal equations
is of the form ff   I is chosen to be S, T , K m;4 , K m;6 , and K m;8 . The
corresponding numbers of iterations required for convergence are equal to 21; 33; 22; 22,
and 23 respectively. The number of iterations without preconditioning is 171. The relative
error of the regularized solution is about 3:1 \Theta 10 \Gamma1 . In contrast, it is about 6:9 \Theta 10 +2
if no regularization is used. Thus we see that our preconditioners also work for strongly
ill-conditioned systems after it is regularized.
6 Concluding Remarks
We remark that even for mildly ill-conditioned matrices with condition number of order
then the matrix A n will be very ill-conditioned already for moderate
regularization is also needed in this case. Once the system is
regularized, our preconditioner C n [K m;8   f ] will work even if p ? 6, cf. Example 2 in x5
for instance. Hence in general, the circulant preconditioner C n [K m;8   f ] should be able
to handle all cases, whether the matrix A n is well-conditioned, mildly ill-conditioned, or
very ill-conditioned but regularized.



--R

Superfast solution of real positive definite Toeplitz systems
Finite Element Solution of Boundary Value Problems
Analysis of preconditioning techniques for ill-conditioned Toeplitz matrices

Circulant Preconditioners for Hermitian Toeplitz Systems
Toeplitz Preconditioners for Toeplitz Systems with Nonnegative Generating Functions
Conjugate Gradient Methods for Toeplitz Systems
Circulant Preconditioners from B-Splines
Circulant Preconditioners Constructed from Kernels
Circulant Preconditioners for Ill-Conditioned Hermitian Toeplitz Matrices
An Optimal Circulant Preconditioner for Toeplitz Systems
An algorithm for the regularization of ill-conditioned
Matrix Analysis
The Theory of Tikhonov Regularization for Fredholm Equations of the First Kind
Approximation of Functions
Preconditioners for Ill-Conditioned Toeplitz Matrices
Preconditioners for Ill-Conditioned Toeplitz Systems Constructed from Positive Kernels
Preconditioning Strategies for Hermitian Toeplitz Systems with Nondefinite Generating Functions
On the extreme eigenvalues of Hermitian (block) Toeplitz matrices
Introduction to Applied Mathematics
A Proposal for Toeplitz Matrix Calculations
An Algorithm for the Inversion of Finite Toeplitz Matrices
Circulant Preconditioners with Unbounded Inverses
--TR

--CTR
Weiming Cao , Ronald D. Haynes , Manfred R. Trummer, Preconditioning for a Class of Spectral Differentiation Matrices, Journal of Scientific Computing, v.24 n.3, p.343-371, September 2005
preconditioned conjugate gradient method;kernel functions;toeplitz systems;circulant preconditioner