Preconditioners for Nondefinite Hermitian Toeplitz Systems.
This paper is concerned with the construction of circulant preconditioners for Toeplitz systems arising from a piecewise continuous generating function with sign changes.If the generating function is given, we prove that for any $\varepsilon >0$, only ${\cal O} (\log N)$ eigenvalues of our preconditioned Toeplitz systems of size N  N are not contained in $[-1-\varepsilon, -1+\varepsilon] \cup [1 -\varepsilon, 1+\varepsilon]$. The result can be modified for trigonometric preconditioners. We also suggest circulant preconditioners for the case that the generating function is not explicitly known and show that only ${\cal O} (\log N)$ absolute values of the eigenvalues of the preconditioned Toeplitz systems are not contained in a positive interval on the real axis.Using the above results, we conclude that the preconditioned minimal residual method requires only ${\cal O} (N \log ^2 N)$ arithmetical operations to achieve a solution of prescribed precision if the spectral condition numbers of the Toeplitz systems increase at most polynomial in N. We present various numerical tests.
Introduction
Let L 2- be the space of 2-periodic Lebesgue integrable real-valued functions and let C 2- be
the subspace of 2-periodic real-valued continuous functions with norm
The Fourier coefficients of f 2 L 2- are given by
a
\Gamma-
Research supported in part by the Hong Kong-German Joint Research Collaboration Grant from the
Deutscher Akademischer Austauschdienst and the Hong Kong Research Grants Council.
and the sequence fAN (f)g 1
N=1 of (N; N)-Toeplitz matrices generated by f is defined by
real-valued the matrices AN (f) are Hermitian.
We are interested in the iterative solution of Toeplitz systems
where the generating function f 2 L 2- . To be more precise, we are looking for good preconditioning
strategies so that Krylov space methods applied to the preconditioned system converge
in a few number of iteration steps. Note that by the Toeplitz structure of AN each iteration
step requires only O(N log N) arithmetical operations by using fast Fourier transforms.
Preconditioning techiques for Toeplitz systems have been well-studied in the past 10 years.
However, most of the papers in this area are concerned with the case where the generating
function f is either positive or nonnegative, see for instance [4, 3, 18, 6, 16, 9] and the references
therein. In this paper, we consider f that has sign changes. The method we propose here will
also work for generating functions that are positive or nonnegative.
Up to now iterative methods for Toeplitz systems with generating functions having different
signs were only considered in [18, 20] and in connection with non-Hermitian systems in [7, 5].
In [7], we have constructed circulant preconditioners for non-Hermitian Toeplitz matrices with
known generating function of the form
where p is an arbitrary trigonometric polynomial and h is a function from the Wiener class
with jhj ? 0. We proved that the preconditioned matrices have singular values properly
clustered at 1. Then, if the spectral condition number of AN (f) fulfills - 2 (AN
the conjugate gradient method (CG) applied to the normal equation requires only O(log N)
iteration steps to produce a solution of fixed precision. However, in general nothing can be
said about the eigenvalues of the preconditioned matrix.
In this paper, we consider real-valued functions f 2 L 2- of the form
where
Y
is a trigonometric polynomial with a finite number of zeros of even
order 2s j and where h 2 L 2- is a piecewise continuous function with simple discontinuities at
there exist h(- j \Sigma
0))=2. Further, we assume that
In particular, we are interested in the Heavyside function h.
A similar setting was also considered in [18]. S. Serra Capizzano suggested the application
of band-Toeplitz preconditioners AN (p s ) in combination with CG applied to the normal
equation. He proved, beyond a more general result which can not directly be used for precon-
ditioning, that at most o(N) eigenvalues of the preconditioned matrix AN (p s
absolute values not contained in a positive interval on the real axis.
A result with o(N) outlyers was also obtained in [19], where the application of preconditioned
GMRES was examined.
In the following, we construct circulant preconditioners for the minimal residual method (MIN-
RES). Note that preconditioned MINRES avoids the transformation of the original system
to the normal equation but requires Hermitian positive definite preconditioners. Then, the
preconditioned matrices are again Hermitian, so that the absolute values of their eigenvalues
coincide with their singular values. If the generating function is given, we prove that for
" ? 0, only O(log N) singular values of the preconditioned matrices are not contained
in [1 \Gamma We also construct circulant preconditioners for the case that the generating
function of the Toeplitz matrices is not explicitly known. For this, we use positive reproducing
kernels with special properties previously applied by the authors in [16, 9] and show
that O(log N) singular values of the preconditioned matrices are not contained in a positive
interval on the real axis. Then, if in addition - 2 (AN preconditioned MINRES
converges in at most O(log N) iteration steps. In summary, the proposed algorithm requires
only O(N log 2 N) arithmetical operations.
This paper is organized as follows: In Section 2, we introduce circulant preconditioners for
(1.1) under the assumption that the generating function of the sequence of Toeplitz matrices is
known and prove clustering results for the eigenvalues of the preconditioned matrices. Section
3 deals with the construction of preconditioners if the generating function of the Toeplitz
matrices is not explicitly known. In Section 4, we modify the results of Section 2 with respect
to trigonometric preconditioners. The convergence of MINRES applied to our preconditioned
Toeplitz systems is considered in Section 5. Finally, we present numerical results in Section
6.
Circulant preconditioners involving generating functions
First we introduce some basic notation. By RN (M) we denote arbitrary (N; N)-matrices of
rank at most M . Let M N (g) be the circulant (N; N)-matrix
diag
F
where F N denotes the N-th Fourier matrix
e \Gamma2-ijk=N
and where F   is the transposed complex conjugate matrix of F . For a trigonometric polynomial
k=\Gamman 1
the matrices AN (q) and M N (q) are related by
(see [14]). For a function g with a finite number of zeros we define the set I N (g) by
I N (g) :=
and the points xN;l (g) (l
xN;l (g) :=
where ~ l 2 1g is the next higher index to l so that ~ l 2 I N (g). For N large enough
we can simply choose ~ l By M N;g (f) we denote the circulant matrix
diag (f(x N;l
If g has m zeros, then we have by construction that
Assume now that the sequence fAN (f)g 1
N=1 of nonsingular Toeplitz matrices is generated by
a known piecewise continuous function f 2 L 2- of the form (1.2) - (1.4). Then we suggest
the Hermitian positive definite circulant matrix M N;f (jf j) as preconditioner for MINRES.
We examine the distribution of the eigenvalues of M N;f (jf
.
The following theorem is Lemma 10 of [22] written with respect to our notation.
Theorem 2.1 Let h 2 L 2- be a piecewise continuous function having only simple discontinuities
at By FN we denote the Fej'er kernel
FN (t) :=
e ikt
cos kt (2.4)
sin
and by FN   h the cyclic convolution of FN and h. Then, for any " ? 0, there exist constants
independent of N so that the number -("; AN ) of eigenvalues of AN
absolute value exceeding " can be estimated by
In other words, we have by Theorem 2.1 that
where V N is a matrix of spectral norm - " and where
Using Theorem 2.1, we can prove the following lemma.
Lemma 2.2 Let be given by (1:2) - (1.4). Then, for any " ? 0 and sufficiently
large N , the number of singular values of M N;f
2 which are not
contained in the interval [1 \Gamma
Proof. By (2.6) and since the eigenvalues of M N;f (jhj) are restricted from below by h \Gamma , it
remains to show that for any " ? 0 and sufficiently large N , except for O(log N) eigenvalues,
all eigenvalues of M N;f (jhj) \Gamma1 M N (FN   h) have absolute values in [1 \Gamma Indeed we
will prove that there are only O(1) outlyers.
For this we follow mainly the lines of proof of Gibb's ph-anomenon. Without loss of generality
we assume that h 2 L 2- has only one jump at
First we examine FN   g, where g is given by
By (2.4) and since g has Fourier series
sin
we obtain Z xFN (t)
sin
and further by (2.5)
Z x/
sin Ntsin t! 2
Z x/
sin
sin N t' 2
Z Nx0
sin t
xand by partial integration and definition of g
where si (y) :=
y
Rsin t
dt. We are interested in the behavior of
Here dxe denotes the smallest integer - x. It is well known that lim
. Thus, if
so that
" for all N -
The same holds if we approach 0 from the left, i.e. if we consider 2-l=N for
Next we have by definition of g and h that
is a continuous function. Since FN is a reproducing kernel, for any " ? 0, there exists
~
so that for all l 2
Assume that l = we obtain by (2.7)
and (2.8) that for any " ? 0 there exists
and consequently, since jh
denote the number of zeros of f which are equal to one of the points 2-l=N
1). Then the set
contains at least absolute values of eigenvalues of M N;f (jhj) \Gamma1 M N (FN   h) and we
conclude by (2.9) that except for O(1) eigenvalues and sufficiently large N , all eigenvalues of
have absolute values contained in [1 \Gamma This completes the
proof.
Remark 2.3 In a similar way as above we can prove that for any " ? 0 and N sufficiently
large, the number of eigenvalues of AN (h) with absolute values not in the interval
is O(log N ).
Note that the property that at most o(N) eigenvalues of AN (h) have absolute values not
contained in [h simply from the fact that the singular values of AN (h) are
distributed as jhj [13, 19].
Theorem 2.4 Let be given by (1:2) - (1.4). Then, for any " ? 0 and
sufficiently large N , except for O(log N) singular values, all singular values of
contained in [1 \Gamma
Proof. The polynomial p s in (1.3) can be rewritten as
where
Y
and -
p(t) is the complex conjugate of p(t). By straightforward computation it is easy to check
that
where only the first s columns (rows) of R c (r)
N are nonzero columns (rows).
p jhj the eigenvalues of M N;f (jf coincide with the eigenvalues of
BN
Now we obtain by (2.10), (2.1) and (2.3) that
BN
By Lemma 2.2, for any " ? 0 and N sufficiently large, except for O(log N) singular values,
all singular values of M N;f
are contained in [1 \Gamma "; 1+ "]. Now the
assertion follows by (2.12) and Weyl's interlacing theorem [12, p. 184].
3 Circulant preconditioners involving positive kernels
In many applications we only know the entries a k (f) of the Toeplitz matrices AN (f ), but
not the generating function itself. In this case, we use even positive reproducing kernels
. These are trigonometric polynomials of the form
c N;k cos kt; c
satisfying KN - 0,2-
\Gamma-
KN
and the reproducing property
lim
Since
(KN
\Gamma-
a k (f) c N;k e ikx ;
the cyclic convolution of KN and f is determined by the first N Fourier coefficients of f . As
preconditioner which can be constructed from the entries of AN (f) without explicit knowledge
of f we suggest the circulant matrix M N;KN \Lambdaf (jK N   f j).
In order to obtain a suitable distribution of the eigenvalues of the preconditioned matrices,
we need kernels with a special property which is related to the order
of the zeros of p s .
The generalized Jackson kernels J m;N of degree - are defined by
determined by (3.1). Here btc denotes the largest
integer - t. In particular, we have that
i.e. there exist positive constants c 1 ; c 2 so that c 1 N 1\Gamma2m - m;N - c 2 N 1\Gamma2m . See [10, pp.
possibility for the construction of the Fourier coefficients of J m;N is prescribed
in [9].
The B-spline kernels B m;N of degree - are defined by
sinc
where Mm denotes the centered cardinal B-spline of order m and
sinc t :=
sin t
See [16, 8]. Since
cos kt
the Fourier coefficients of B m;N are given by values of centered cardinal B-splines. Note that
is just the Fej'er kernel FN .
The above kernels have the following important property:
Theorem 3.1 be given by (1:2) - (1:4). Assume that for all t j (j 2
exists a neighborhood so that f is a monotone function in this
neighborhood and moreover f(t m;N be given by
(3.2) or (3.3), where
Then there exist so that for N !1, except for O(1) points, all points of the
set
Proof. 1. First we consider the upper bound. Since p s and KN are nonnegative, we obtain
Z
\Gamma-
Z
\Gamma-
In [16, 9], we proved that m - oe implies that for all x 2 I N (p s ) ' I N (f ), there exists a
(KN   p s )(x)
Thus, since jh(x)j - h \Gamma for
(KN   p s )(x)
2. Next we deal with the lower bound.
2.1. Let x 2 I N (f) be not in the neighborhood of t j (j exist
independent of N so that
since KN is a reproducing kernel and by using the same arguments as in
the proof of Lemma 2.2 if x is in the neighborhood of some we obtain that,
for any " ? 0 there exists N("), so that except for at most a constant number of points, all
considered points x 2 I N (f) satisfy
and thus
2.2. It remains to consider the points
For simplicity we assume that
i.e. p s has only a zero of order 2s at
lim
For any fixed 0
(KN
Z
\Gammab
\Gammab
Z
\Gamma-
Z
Z
\Gammab
b+x
and since f is bounded
(KN
Z
\Gammab
b+x
By definition of KN we see that for any fixed 0 ! ~ b -
Z
KN (t) dt - const N \Gamma2m+1 ; (3.5)
so that we get for small x (e.g. x ! b=2)
(KN
Z
\Gammab
2.2.1. Assume that h has no jump at
that h(t) - h \Gamma or h(t) - \Gammah \Gamma for t 2 [\Gamma"; "]. We restrict our attention to the case h - h \Gamma .
monotone increasing on (0; -), we
obtain for x(N) 2 (0; ") " I N (f) and N sufficiently large that
Z
\Gamma"
Z
Z
with a positive constant c independent of N . On the other hand, we have by definition of
s and since by assumption . Then we
obtain by (3.6) with that for N large enough
(KN   f)(x(N))
const
with a positive constant const independent of N .
The proof for x(N) 2 (\Gamma"; follows the same lines.
2.2.2. Finally, we assume that h has a jump at
Without loss of generality let h(0 by assumption on f , there exists
that
Z
We consider points of the form
with lim
in case of Jackson kernels and fl := 1 in case of
B-spline kernels. Then we have for t 2
and consequently for sufficiently small " 1 and y, since sin is odd and monotone increasing on
(0; -=2) that
Further, by definition of the B-spline kernels
m;N (t) := N
sinc
and similarly as in (3.9) we see that
By assumption h does not change the sign in (0; " 1 ). Then we obtain by (3.8), monotonicity
of p s in (0; -) and m - s
Z
Z
y
where K 0
m;N g. Set
there exist
Z
wl
l=r\Gammak
wl
Z
and further by (3.5) and since lim
Z
Straightforward computation yields
const
Z' sin u
Hence we get for N large enough that
Z
const
and by (3.10) that
Z
const (3.11)
with positive constants const independent of N .
Now we consider x(N) 2 I N (f) with y k (N) - x(N) ! y k+1 (N ).
Z
Z
Z
Z
and since f is by assumption monotone increasing on [\Gamma"
Z
Z
Z
Z
Z
Z
and by (3.5) and since f is bounded
Z
Z
By assumption
R
const
R
and since f(y k (N) - const N \Gamma2s and m - s by (3.12), (3.11) that for N large
enough
Z
const
with a nonnegative constant const independent of N . Finally, we use (3.6) with
again to finish the proof.
To show our main result we also need the following lemma.
Lemma 3.2 Let A 2 C N;N be a Hermitian positive definite matrix having
in [a be a Hermitian matrix with
singular values in [b 1. Then at least N \Gamma 4n
of A B are contained in [\Gammaa
Proof. 1. Assume first that n i.e. A has only eigenvalues in [a
the j-th eigenvalue of the matrix B. We consider the eigenvalues of
to t 2 R. By Weyl's interlacing theorem (see [12, p. 184]) we obtain for t - 0 that
and for t ! 0 that
we obtain by (3.13) and (3.14) that - j
. On the other hand, we see by (3.13) and (3.14) that - j
. Thus, since - j
is a continuous function in t 2 R, there exists
This implies that t j
is an eigenvalue of AB. Consequently, every corresponds to an eigenvalue
of AB. (Eigenvalues are called with multiplicities.)
The examination of the same lines.
In summary, N \Gamma n 2 eigenvalues of AB are contained in [\Gammaa
2. Let n 1 eigenvalues of A be outside [a since A is positive definite, the matrix
can be splitted as
A 1=2
where ~
A 1=2
is Hermitian with all eigenvalues in [a 1=2
of rank n 1 . The eigenvalues of AB coincide with the eigenvalues of A 1=2 BA 1=2 . Hence it
remains to show that at most 4n 1 singular values of A 1=2 BA 1=2 are not contained in
we have
A 1=2 BA
A 1=2
A 1=2 BA 1=2
~
A 1=2
A 1=2
By 1. all but n 2 singular values of ~
A 1=2
A 1=2
are contained in [a
and Weyl's interlacing theorem yield the assertion.
Theorem 3.3 Let be given by (1:2) - (1:4). Assume that for all t j (j 2
exists a neighborhood so that f is a monotone function in this
neighborhood and moreover f(t m;N be given by
(3.2) or (3.3), where
By ff; fi we denote the constants from Theorem 3:1.
Then, for any " ? 0 and sufficiently large N , except for O(log N) singular values, all singular
values of M N (jK N   f
are contained in [ff \Gamma
Proof. Let BN (f) be defined by (2.11). Then we obtain by (2.12) that
The distribution of the eigenvalues of M N;f
2 is known by Lemma
2.2. It remains to examine the eigenvalues of the Hermitian positive definite matrix
These eigenvalues coincide with the reciprocal eigenvalues of M N;f (jf
By definition of M N;g and since KN is a reproducing kernel, except for O(1) eigenvalues, all
eigenvalues of M N;f (jf are given by j(K N   f)(2-l=N)j=jf(2-l=N)j
(l 2 I N (f )). Thus, by Theorem 3.1, for N !1 only O(1) eigenvalues of
are not contained in [ff; fi]. Consequently, by (3.18), Lemma
2.2, Lemma 3.2 and Weyl's interlacing theorem at most O(log N) singular values of
2 are not contained in [ff \Gamma
4 Trigonometric preconditioners
In addition to Section 2, we suppose that the Toeplitz matrices AN 2 R N;N are symmetric,
i.e. the generating function f 2 L 2- is even. This suggests the application of so-called
trigonometric preconditioners. Note that in the symmetric case the multiplication of a vector
with AN can be realized using fast trigonometric transforms instead of fast Fourier transforms
(see [14]). In this way complex arithmetic can be completely avoided in the iterative solution
of (1.1). This is one of the reasons to look for preconditioners which can be diagonalized by
trigonometric matrices corresponding to fast trigonometric transforms instead of the Fourier
matrix F N .
In practice, four discrete sine transforms (DST I - IV) and four discrete cosine transforms
(DCT I - IV) were used (see [21]). Any of these eight trigonometric transforms can be realized
with O(N log N) arithmetical operations. Likewise, we can define preconditioners with respect
to any of these transforms.
In this paper, we restrict our attention to the so-called discrete cosine transform of type
II (DCT-II) and discrete sine transform of type II (DST-II), which are determined by the
following transform matrices:
where ffl N
1). We propose the preconditioners
diag (jf(~x N;l
diag (jf(~x N;l )j) N
where
~ xN;l :=
l-
l-
~ l-
and where ~ l 2 1g is the next higher index to l such that jf(~x N;l )j ? 0. See [15].
Then we can prove in a completely similar way as in Section 2 that for any " ? 0 and
sufficiently large N except for O(log N) singular values, all singular values of
are contained in [1 \Gamma
5 Convergence of preconditioned MINRES
In order to prescribe the convergence behavior of preconditioned MINRES with our preconditioners
of the previous sections, we have to estimate the smaller outlyers for increasing
N .
Lemma 5.1 Let f 2 L 2- be defined by (1:2)-(1:4). Assume that - 2 (AN
(ff ? 0). Then the smallest absolute values of the eigenvalues of M N;f (jf
behave for N !1 as O(N \Gammaff ).
Proof. Since
and both kM N;f (jf j)k 2 and kM N;KN \Lambdaf (jK N   f j)k 2 are restricted from above, it remains to
show that there exists a constant c ? 0 independent of N so that
kAN
The above inequality follows immediately from the fact that the singular values of AN (f) are
distributed as jf j (see [13, 19]).
We want to combine our knowledge of the distribution of the eigenvalues of our preconditioned
matrices with results concerning the convergence of MINRES.
Theorem 5.2 Let A 2 C N;N be a Hermitian matrix with p and q isolated large and small
singular values, respectively:
Let for the
solution of
iteration steps to achieve precision - , i.e. jjr (k) jj 2
the k-th iterate.
The theorem can be proved by using the same technique as in [1, pp. 569 - 573]. Namely,
based on the known estimate
jjr
k denotes the space of polynomials of degree - k with p k are the
eigenvalues of A, we choose p k as product of the linear polynomials passing through the
outlyers and the modified Chebyshev polynomials
The above summand p ln 2 can be further reduced if we use polynomials of higher degree for
the larger outlyers.
Note that a similar estimate can be given for the CG method applied to the normal equation
A   b. Here we need
iteration steps to archive precision
jje (0) jj A
Note that the latter
method requires two matrix-vector multiplications in each iteration step.
By Theorem 2.4, Theorem 3.3 and Lemma 5.1 our preconditioned MINRES with both preconditioners
produces a solution of (1.1) of prescribed
precision in O(log N) iteration steps and with O(N log 2 N) arithmetical operations. The
same holds for preconditioned CG applied to the normal equation.
6 Numerical results
In this section, we test our circulant and trigonometric preconditioners in connection with
different iterative methods on a SGI O2 work station. As transform length we use
right-hand side b of (1.1) the vector consisting of N entries "1" and as start vector the zero
vector.
We begin with a comparison of MINRES applied to
N g and CGNE (Craig's method) (cf. [17, p. 239]) applied to
For both algorithms we have used MATLAB implementations of B. Fischer. See also [11]. In
particular, his implementation of preconditioned MINRES avoids the splitting (6.2).
In order to make the following computations with MINRES and CGNE comparable, we have
stopped both computations if
Example 1. We begin with Hermitian Toeplitz matrices AN (f) arising from the generating
function

Table

1 presents the number of iterations for circulant preconditioners. The first row of the
table contains the exponent n of the transform length . According to Theorem 2.4
and Theorem 5.2, the preconditioners M N (jf j; F N ) lead to very good results. As expected,
the preconditioners M N;KN \Lambdaf (jK N   f j; F N ) with the Fej'er kernels are not suitable
for (1.1) (cf. also [16]), while the preconditioners with do their job.
Further, CGNE needs half the number of iterations but twice the number of matrix-vector
multiplications per iteration than MINRES. See also Section 5.
method
MINRES I N 23 71 277 *
43
CGNE I N 11 37 164 *

Table

1:
Example 2. Next, we consider the symmetric Toeplitz matrices AN (f) arising from the
generating function
with

Tables

2 presents the number of iterations for trigonometric preconditioners. The results
are similar to those of Example 1, except that CGNE requires nearly the same number of
iterations as MINRES.
method
MINRES I N 9 17
MINRES M N;f (jf j; C II
MINRES M N;f (jf j; S II
MINRES M N;FN \Lambdaf (jF N   f j; C II
CGNE I

Table

2: f 2

Acknowledgment

. The authors wish to thank B. Fischer for for the MATLAB implementations
of PMINRES and CGNE.



--R

Iterative Solution Methods.

preconditioning for Toeplitz matrices.
Toeplitz preconditioners for Toeplitz systems with nonnegative generating functions.

Conjugate gradient methods of Toeplitz systems.
Preconditioners for non-Hermitian Toeplitz systems
Circulant preconditioners from B-splines
The best circulant preconditioners for Hermitian Toeplitz matrices.
Constructive Approximation.
Polynomial Based Iteration Methods for Symmetric Linear Systems.
Matrix Analysis.
On the distribution of singular values of Toeplitz matrices.
Optimal trigonometric preconditioners for nonsymmetric Toeplitz systems.
Preconditioners for ill-conditioned Toeplitz matrices
Preconditioners for ill-conditioned Toeplitz matrices constructed from positive kernels
Iterative Methods for Sparse Linear Systems.
Preconditioning strategies for Hermitian Toeplitz systems with nondefinite generating functions.
A unifying approach to some old and new theorems on distribution and clustering.

Fast algorithms for the discrete W transform and for the discrete Fourier transform.
Circulant preconditioners for Toeplitz matrices with piecewise continuous generating functions.
--TR
circulant matrices;preconditioners;nondefinite Toeplitz matrices;krylov space methods