Preconditioners for symmetrized Toeplitz and multilevel Toeplitz matrices

When solving linear systems with nonsymmetric Toeplitz or multilevel Toeplitz matrices using Krylov subspace methods, the coefficient matrix may be symmetrized. The preconditioned MINRES method can then be applied to this symmetrized system, which allows rigorous upper bounds on the number of MINRES iterations to be obtained. However, effective preconditioners for symmetrized (multilevel) Toeplitz matrices are lacking. Here, we propose novel ideal preconditioners, and investigate the spectra of the preconditioned matrices. We show how these preconditioners can be approximated and demonstrate their effectiveness via numerical experiments.


Introduction. Linear systems
where A n ∈ R n×n is a Toeplitz or multilevel Toeplitz matrix, and b ∈ R n arise in a range of applications. These include the discretization of partial differential and integral equations, time series analysis, and signal and image processing [7,28]. Additionally, demand for fast numerical methods for fractional diffusion problemswhich have recently received significant attention-has renewed interest in the solution of Toeplitz and Toeplitz-like systems [10,27,32,33,46]. Preconditioned iterative methods are often used to solve systems of the form (1.1). When A n is Hermitian, CG [18] and MINRES [30] can be applied, and their descriptive convergence rate bounds guide the construction of effective preconditioners [7,28]. On the other hand, convergence rates of preconditioned iterative methods for nonsymmetric Toeplitz matrices are difficult to describe. Consequently, preconditioners for nonsymmetric problems are typically motivated by heuristics.
As described in [36] for Toeplitz matrices, and discussed in subsection 2.2 for the multilevel case, A n is symmetrized by the exchange matrix so that (1.1) can be replaced by with the symmetric coefficient matrix Y n A n . Although we can view Y n as a preconditioner, its role is not to accelerate convergence, and there is no guarantee that (1.3) is easier to solve than (1.1) [12,24]. Instead, the presence of Y n allows us to use preconditioned MINRES, with its nice properties and convergence rate bounds, to solve (1.3). We can then apply a secondary preconditioner P n ∈ R n×n to improve the spectral properties of Y n A n , and therefore accelerate convergence. An additional benefit is that MINRES may be faster than GMRES [38] even when iteration numbers are comparable, since it requires only short-term recurrences. Preconditioned MINRES requires a symmetric positive definite preconditioner P n , but it is not immediately clear how to choose this matrix when A n is nonsymmetric. In [36] it was shown that absolute value circulant preconditioners, which we describe in the next section, give fast convergence for many Toeplitz problems. However, for some problems there may be more effective alternatives based on Toeplitz matrices (see, e.g. [4,17]). Moreover, multilevel circulant preconditioners are not generally effective for multilevel Toeplitz problems [40,41,42]. Thus, alternative preconditioners for (1.3) are needed.
In this paper, we describe ideal preconditioners for symmetrized (multilevel) Toeplitz matrices and show how these can be effectively approximated. To set the scene, we present background material in section 2. Sections 3 and 4 describe the ideal preconditioners for Toeplitz and multilevel Toeplitz problems, respectively. Numerical experiments in section 5 verify our results and show how the ideal preconditioners can be efficiently approximated by circulant matrices or multilevel methods. Our conclusions can be found in section 6.
2. Background. In this section we collect pertinent results on Toeplitz and multilevel Toeplitz matrices.
In many applications, the matrix A n is associated with a generating function f ∈ L 1 ([−π, π]) via its Fourier coefficients We use the notation A n (f ) when we wish to stress that a Toeplitz matrix A n is associated with the generating function f . An important class of generating functions is the Wiener class, which is the set of functions satisfying Many properties of A n (f ) can be determined from f . For example, if f is real than A n (f ) is Hermitian and its eigenvalues are characterized by f [16, pp. 64-65].
On the other hand, if f is complex-valued then A n (f ) is non-Hermitian for at least some n and its singular values are characterized by |f | [2,34,47].
The absolute value circulant [5,29,36] derived from a circulant C n is the matrix Closely related to Toeplitz matrices are Hankel matrices H n ∈ R n×n , which have constant anti-diagonals. It is well known that a Toeplitz matrix can be converted to a Hankel matrix, or vice versa, by flipping the rows (or columns), i.e., via Y n in (1.2). Since Hankel matrices are necessarily symmetric, this means that any nonsymmetric Toeplitz matrix A n can be symmetrized by applying Y n , so that Alternatively, we may think of A n being self-adjoint with respect to the bilinear form induced by Y n [14,35].
and is skew-centrosymmetric if G n Y n = −Y n G n . Thus, (2.5) shows that symmetric Toeplitz matrices are centrosymmetric. It is clear from (2.6) that the inverse of a nonsingular centrosymmetric matrix is again centrosymmetric. Furthermore, nonsingular centrosymmetric matrices have a centrosymmetric square root [22, Corollary 1] 1 .

Multilevel Toeplitz and Hankel matrices.
Multilevel Toeplitz matrices are generalizations of Toeplitz matrices. To define a generating function for a multilevel Toeplitz matrix, let j = (j 1 , . . . , j p ) ∈ Z p be a multi-index and consider a p-variate function f ∈ L 1 ([−π, π] p ), f : [−π, π] p → C. The Fourier coefficients of f are defined as where θ, j = p k=1 θ k j k and dθ = dθ 1 . . . dθ p is the volume element with respect to the p-dimensional Lebesgue measure.
If n = (n 1 , . . . , n p ) ∈ N p , with n i > 1, i = 1, . . . , p, and π(n) = n 1 . . . n p , then f is the generating function of the multilevel Toeplitz matrix A n (f ) ∈ R π(n)×π(n) , where Here, J (k) r ∈ R r×r is the matrix whose (i, j)th entry is one if i − j = k and is zero otherwise.
Similarly, we can define a multilevel Hankel matrix as ∈ R r×r is the matrix whose (i, j)th entry is one if i + j = k + 1 and is zero otherwise. Although a multilevel Hankel matrix does not necessarily have constant anti-diagonals, it is symmetric.
Multilevel Toeplitz matrices can also be symmetrized by the exchange matrix Y n ∈ R π(n)×π(n) , Y n = Y n1 ⊗ · · · ⊗ Y np . To see this we use an approach analogous to that in the proof of [11,Lemma 5]. The key point is that Y r J where b (j1,...,jp) = a (n1−j1,...,np−jp) . Thus, Y n A n (f ) is a multilevel Hankel matrix, and hence is symmetric.
2.3. Assumptions and notation. Throughout, we assume that all Toeplitz or multilevel Toeplitz matrices A n are real, and are associated with generating functions in L 1 ([−π, π] p ). We denote the real and imaginary parts of f by f R and f I , respectively, so that f = f R + if I . We assume that the symmetric part of A n , given by A R = (A n + A T n )/2, is positive definite, which is equivalent to requiring that f R is essentially positive. Similarly, we assume that A n (|f |) is positive definite, and that |f | ≥ δ > 0 for some constant δ.

Preconditioning Toeplitz matrices.
In this section we introduce our ideal preconditioners for (1.3) when A n is a Toeplitz matrix, and analyse the spectrum of the preconditioned matrices. Although these preconditioners may be too expensive to apply exactly, they can be approximated by, e.g., a circulant matrix or multigrid solver.
3.1. The preconditioner A R . The first preconditioner we consider is the symmetric part of A n , namely A R = (A n + A T n )/2, which was previously used to precondition the nonsymmetric problem (1.1) (see [21]). When applied to the symmetrized system (1.3), spectral information can be used to bound the convergence rate of preconditioned MINRES. Accordingly, in this section we characterize the eigenvalues of A We begin by stating a powerful result that characterizes the spectra of preconditioned Hermitian Toeplitz matrices in terms of generating functions. ). Let f, g ∈ L 1 ([−π, π]) be real-valued functions with g essentially positive. Let A n (f ) and A n (g) be the Hermitian Toeplitz matrices with generating functions f and g, respectively. Then, A n (g) is positive definite and the eigenvalues of .
Lemma 3.1 shows that in the Hermitian case we can bound the extreme eigenvalues of preconditioned Toeplitz matrices using scalar quantities. If bounds on the eigenvalues nearest the origin are also available it is possible to estimate the convergence rate of preconditioned MINRES applied to the Toeplitz system. Unfortunately, this result does not carry over to nonsymmetric matrices, nor is an eigenvalue inclusion region alone sufficient to bound the convergence rate of a Krylov subspace method for nonsymmetric problems [1,15]. However, in the following we show that by symmetrizing the Toeplitz matrix A n we can obtain analogous results to Lemma 3.1, even for nonsymmetric A n . As a first step, we quantify the perturbation of the (nonsymmetric) , and let f = f R + if I , where f R and f I are realvalued functions with f R essentially positive. Additionally, let A n := A n (f ) ∈ R n×n be the Toeplitz matrix associated with f . Then Proof. It is easily seen from (2.2) that A n (f ) = A n (f R ) + iA n (f I ). Moreover, from Lemma 3.1 we also know that A R = A n (f R ) is symmetric positive definite and A n (f I ) is Hermitian. It follows that To bound := E n 2 = E n 2 , note that since E n is Hermitian, E n 2 is equal to the spectral radius of E n . Applying Lemma 3.1 thus gives that which completes the proof.
The above result tells us that the nonsymmetric preconditioned matrix will be close to the identity when the skew-Hermitian part of A n is small, as expected. Although this enables us to bound the singular values of A R , these cannot be directly related to the convergence of e.g., GMRES. In contrast, the following result will enable us to characterize the convergence rate of MINRES applied to (1.3).
, and let f = f R + if I , where f R and f I are realvalued functions with f R essentially positive. Additionally, let A n := A n (f ) ∈ R n×n be the Toeplitz matrix associated with f . Then the symmetric positive definite matrix where Y n is the exchange matrix in (1.2) and R is centrosymmetric (see (2.6) and [22] Since Y n is orthogonal, Y n E n 2 = E n 2 , and the result follows from Lemma 3.2.
Applying Weyl's theorem [20,Theorem 4.3.1] to (3.1) shows that the eigenvalues ]. However, as grows, eigenvalues could move close to the origin, and hamper MINRES convergence. We will show that this cannot happen in the following result.
Proof. We know from Lemma 3.3 that where Y n E n 2 < and Y n has eigenvalues ±1. Thus, as discussed above, the eigen- Hence, all that remains is to improve the bounds on the eigenvalues nearest the origin. Our strategy for doing so will be to apply successive similarity transformations to Y n +Y n E n ; as a byproduct, we will characterize the eigenvalues of A Before applying our first similarity transform, we recall from the proofs of Lemmas 3.2 and 3.3 that iA n (f I ) is skew-symmetric and Toeplitz, while A − 1 2 R is symmetric and centrosymmetric. It follows that Y n E n is symmetric but skew-centrosymmetric. Skew-centrosymmetry implies that whenever (λ, x), λ = 0 is an eigenpair of Y n E n , then so is (−λ, Y n x) [19,37,44]. Additionally, any eigenvectors of Y n E n corresponding to a zero eigenvalue can be expressed as a linear combination of vectors of the form x ± Y n x, x ∈ R n [44, Theorem 17]. Therefore, Y n E n has eigendecomposition where n = 2m 1 + m 2 and m 2 = m 3 + m 4 . Since Y n E n is symmetric, we may assume that U n is orthogonal. We can now apply the first similarity transform, namely Using the orthogonality of the columns of U n , it is straightforward to show that Here, Γ 2m1 = diag(1, −1, . . . , 1, −1) and the kth column of Q ∈ R n×2m1 is given by with e j ∈ R n the jth unit vector. Consequently, our second similarity transform gives Hence, letting Since the eigenvalues of Z + Σ k are ± 1 + λ 2 k we see from (3.6) that the eigenvalues of A . . , m 1 , and possibly one or both of 1 and −1. Hence, the eigenvalues are at least 1 in magnitude. This completes the proof.

Theorem 3.4 characterizes the eigenvalues of
R , and hence the convergence rate of preconditioned MINRES, in terms of the scalar quantity in (3.2). Thus, we expect that the preconditioner A R will perform best when A n is nearly symmetric, and we investigate this in section 5. However, irrespective of the degree of nonsymmetry of A n , Theorem 3.4 shows that the eigenvalues of the preconditioned matrix are at least bounded away from the origin.
3.2. The preconditioner A M . We saw in the previous section that A R is an effective preconditioner when the degree of nonsymmetry of A n is not too large. For problems that are highly nonsymmetric, however, a different preconditioner may be more effective. Here, motivated by the success of absolute value preconditioning, we consider the preconditioner A M = A n (|f |) instead. The following result describes the asymptotic eigenvalue distribution of where f = f /|f | and E n 2 = o(n) as n → ∞. Moreover, the eigenvalues of Y n A n ( f ) lie in [−1, 1].
Proof. The conditions on |f | guarantee that A n (|f |) is invertible and that its eigenvalues (singular values) are bounded away from 0. Thus, by Proposition 5 in [9], where E n 2 = o(n) as n → ∞. Since A M is Hermitian and Toeplitz, both A M and its inverse are centrosymmetric. It follows that Since Y n is unitary and Y n A n ( f ) is symmetric, the absolute values of the eigenvalues of Y n A n ( f ) coincide with the singular values of A n ( f ), which in turn are bounded above by one [47]. This proves the result.
A consequence of Theorem 3.5 is that the eigenvalues of A n (|f |) −1 Y n A n lie in [−1− , 1+ ], where for large enough n the parameter is small. Although eigenvalues may be close to the origin, most cluster at −1 and 1, in line with Theorem 3.4 in [24].
To conclude this section, we show how A M can be approximated by circulant preconditioners. First, recall from subsection 2.1 that C n (f ) is the preconditioner with eigenvalues λ j = f (2πj/n), j = 0, . . . , n − 1. For large enough dimension n, we have that A M = C n (|f |) + E n + R n , where E n has small norm and R n has small rank [13, pp. 108-110], so that C n (|f |) is a good approximation to A M for large n.
The matrix C(|f |) can in turn be approximated by the Strang absolute value circulant preconditioner |C (S) n | [5,29,36], where if C (S) n is the Strang circulant preconditioner [43] for A n , with eigenvalues λ j , j = 1, . . . , n, then the corresponding absolute value circulant preconditioner |C (S) n | has eigenvalues |λ j |, j = 1, . . . , n. For this preconditioner, we obtain the following result. n | → C n (|f |) as n → ∞. Proof. Assume that n, the dimension of A n , is n = 2m + 1. (The idea can be extended to the case of even n, as in [6, p. 37].) Then, C Since f is in the Wiener class, (D m f )(θ) converges absolutely, hence uniformly to f (θ). Thus, Since the eigenvalues of |C n (D m f )| approach those of C n (|f |) as n → ∞ we obtain the result.

Multilevel Toeplitz problems.
We now extend the results of section 3 to multilevel Toeplitz matrices.

4.1.
The preconditioner A R . The results of subsection 3.1 carry over straightforwardly to the multilevel case. They depend on the following generalization of Lemma 3.1, which can be found, e.g., in Proposition 2 in [23], and is an extension of Theorem 3.1 in [39].

Theorem 4.2 characterizes the eigenvalues of
R , which are bounded away from the origin. In turn, this allows us to bound the convergence rate of preconditioned MINRES in terms of easily-computed quantity in (4.1).

The preconditioner A M .
We can also extend the results in subsection 3.2 to the multilevel case. However, for multilevel problems this preconditioner is more challenging to approximate. Matrix algebra, e.g., block circulant, preconditioners will generally result in iteration counts that increase as the dimension increases, as previously discussed. On the other hand, constructing effective banded Toeplitz, or efficient multilevel, algorithms is challenging since it is generally necessary to compute elements of A M . Nonetheless, we present the following result for completeness. It directly generalizes the result for Toeplitz matrices, so is presented without proof. Theorem 4.3. Let f : L ∞ ([−π, π] p ), with 0 < δ < |f (θ)| for all θ ∈ [−π, π] p . Then, if A M = A n (|f |) is the multilevel Toeplitz matrix generated by |f |,

Numerical experiments.
In this section we investigate the effectiveness of the preconditioners described above, and approximations to them, for the symmetrized system (1.3). We also compare the proposed approach to using nonsymmetric preconditioners for (1.1) within preconditioned GMRES and LSQR [31]. All code is written in Matlab (version 9.4.0) and is run on a quad-core, 62 GB RAM, Intel i7-6700 CPU with 3.20GHz 2 . We apply Matlab versions of LSQR and MINRES, and a version of GMRES that performs right preconditioning. (Note that LSQR requires two matrix-vector products with the coefficient matrix and two preconditioner solves per iteration.) We take as our initial guess x 0 = (1, 1, . . . , 1) T / √ n, and we stop all methods when r k 2 / r 0 2 < 10 −8 . When more than 200 iterations are required, we denote this by '-' in the tables.
When A R or A M are too expensive to apply directly we use either a circulant or multigrid approximation. The multigrid preconditioner consists of a single V-cycle with damped Jacobi smoothing and Galerkin projections, namely linear or bilinear interpolation and restriction by full-weighting. The coarse matrices are also built by projection. The number of smoothing steps and the damping factor ω are stated below for each problem. The damping parameter is chosen by trial-and-error to minimize the number of iterations needed for small problems. When applying circulant preconditioners to (1.3) we use the absolute value preconditioner in (2.4) based on the Strang [43], optimal [8] or superoptimal [45] circulant preconditioner.
Example 5.1. Our first example is from [21,Example 2], where numerical experiments indicated that A R is an effective preconditioner for the nonsymmetric system (1.1) when GMRES is applied. The Toeplitz coefficient matrix A n = A n (f ) is formed from the generating function f (θ) = (2−2 cos(θ))(1+iθ). Since computing the Fourier coefficients for larger problems is time-consuming, smaller problems examined here. The right-hand side is a random vector (computed using the Matlab function randn).
The preconditioner A R := A n (f R ) is positive definite, since f R (θ) = 2 − 2 cos(θ) is essentially positive. Indeed, A R is the second-order finite difference matrix, namely the tridiagonal matrix with 2 on the diagonal and −1 on the sub-and super-diagonals. Accordingly, A R can be applied directly with O(n) cost. For comparison we also apply the optimal circulant preconditioner C n and its absolute value counterpart |C n |. (The optimal circulant outperformed the Strang and superoptimal circulant preconditioners for this problem.) The absolute value circulant approximates A M . Table 5.1 shows that A R requires fewer iterations than C n for MINRES and LSQR, and that MINRES with A R is the fastest method overall. The good performance of A R with MINRES can be explained by the clustered eigenvalues of A −1 R A n . Theorem 3.4 tells us that these eigenvalues lie in [−1 − π, −1] ∪ [1, 1 + π], and Figure 5.1 (b) shows that these bounds are tight. As discussed in [21], the eigenvalues of A −1 R A n are also nicely clustered (see Figure 5.1 (a)), with real part 1 and imaginary part in [−π, π]. Although we cannot rigorously link this eigenvalue characterization to the rate of Table 5.1: Iteration numbers and CPU times (in parentheses) for the optimal circulant preconditioner C n and tridiagonal preconditioner A R for Example 5.1.  GMRES convergence, Table 5.1 indicates that in this case A R is also a reasonable preconditioner for GMRES. We now consider A M , which is dense since |f (θ)| = (2 − 2 cos(θ)) √ 1 + θ 2 . Accordingly, as well as applying A M exactly-to confirm our theoretical results-we approximate A M via our V-cycle multigrid method with 2 pre-and 2 post-smoothing steps, the coarsest grid of dimension 15, and ω = 0.1 for GMRES, ω = 0.4 for LSQR and ω = 0.5 for MINRES. For LSQR, multigrid with A M gave lower timings and iteration counts than multigrid with A n , and so was used instead.
Iteration counts and CPU times (excluding the time to construct A M but including the time to set up the multigrid preconditioner) are given in Table 5.2. Both A M and its multigrid approximation give lower iteration counts than A R , with the multigrid method especially effective for MINRES applied to the symmetrized system. However, timings are higher than for A R since the O(n log(n)) multigrid method is more expensive than the O(n) solve with A R . The eigenvalues of A −1 M Y n A n , when n = 2047, are as expected from Theorem 3.5 (see Figure 5.2), since all eigenvalues lie in [−1, 1]. Indeed, most cluster at the endpoints of this interval. The eigenvalues of A −1 M A n are also localized, but not as clustered, indicating that the spectrum of A −1 M A n may differ significantly from that of A −1 M Y n A n . Example 5.2. We now examine the linear system obtained by discretising a fractional diffusion problem from [3], which we alter so as to make it nonsymmetric. The Table 5  problem is to find u(x, t) that satisfies where α ∈ (1, 2), and d + and d − are nonnegative constants. We impose the absorbing boundary conditions u(x ≤ 0, t) = u(x ≥ 1, t) = 0, t ∈ [0, 1], while u(x, 0) = 80 sin(20x) cos(10x), x ∈ [0, 1]. The Riemann-Liouville derivatives in (5.1) are where n is the integer for which n − 1 < α ≤ n.
Discretising by the shifted Grünwald-Letnikov method in space, and the backward Euler method in time [25,26], gives the linear system where g α,k = (−1) k α k , ν = τ h α and h = 1 n+1 . We set τ = 1/ n α , which makes ν constant, so that all the theory of section 3 can be directly applied, but comparable results are obtained τ = 1/n. Stated CPU times and iteration counts in this example are for the first time step. (Iteration counts and timings decrease at later time steps.) CPU times include the preconditioner setup time and solve time.
Entries of A in (5.2) are generated by [10] The real part of ϕ is essentially positive, so A R = (A n + A T n )/2 is positive definite. However, since A R is dense we approximate it by our V-cycle multigrid method (analysed in [33]) with the coarsest grid of dimension 127, 2 pre-and 2 post-smoothing steps, and ω = 0.7 for all Krylov solvers. The matrix A M is also dense and positive definite, and we approximate it using two different approaches. The first is the absolute value Strang preconditioner discussed at the end of subsection 3.2. The second is multigrid (with the same parameters as for A R , except that we use 1 pre-and postsmoothing step) applied to a banded Toeplitz approximation of A M . Specifically, if r and c are the first row and column of A M , when α = 1.25 we compute the first 50 elements in r and c, and when α > 1.25 we take the first β(1.1) log 2 (n+1) elements in r and c, where β = 40 when α = 1.5 and β = 100 when α = 1.75. This balances the time to compute these coefficients, and the resulting MINRES iteration count.
We see from Table 5.3 that our approximations to A R and A M are robust with respect to n, but both require slightly more iterations for larger α. The multigrid preconditioner for A R requires fewer iterations than the circulant, but the latter results in a lower CPU time because the preconditioner application is cheap, and indeed the absolute value preconditioner with MINRES is the fastest method overall. Of the multigrid methods, the approximation to A R with MINRES is fastest for α ≤ 1.5, while the multigrid approximation of A n with GMRES is slightly faster for large α.
In Table 5.4 we investigate the effect of d + and d − , i.e., of nonsymmetry, on the preconditioners. The results are unchanged when d + and d − are swapped, so we tabulate results for d + ≤ d − only. As expected, our approximation to A R is best suited to problems for which d + and d − do not differ too much. The hardest problem for A R is when d + = 0, since in this case A n is a Hessenberg matrix, and hence highly nonsymmetric. However, even here the iteration numbers are fairly low, since the eigenvalues are bounded away from the origin independently of n. The circulant and multigrid preconditioners based on A M are not greatly affected by altering d + and d − .
The low iteration numbers and mesh-size independent results for A R in Table 5.4 are explained by Theorem 3.4 and the relatively small upper bound (3.2), which describes how far eigenvalues of A −1 R Y n A n can deviate from 1 in magnitude. This bound is 0 when d + = d − , or when α = 2, since then A n is symmetric. However, Table 5.5 shows that even when A n is nonsymmetric the bound is quite small. Additionally, it does not change when the values of d + and d − are swapped.    Example 5.3. We now solve a two-level Toeplitz problem that also arises from fractional diffusion and is based on the symmetric problem in [3]. We seek u(x, y, t) in the domain Ω = (0, 1) 2 × (0, 1] that satisfies where α, β ∈ (1, 2), and d + , d − , e + and e − are nonnegative constants. We impose absorbing boundary conditions, and the initial condition is u(x, 0) = 100 sin(10x) cos(y)+ sin(10t)xy. We again discretize by the shifted Grünwald-Letnikov method in space, and the backward Euler method in time [25,26], which leads to the following linear system: Here n x and n y are the number of spatial degrees of freedom in the x and y directions, respectively; we choose n x = n y = n. Also, where L α is given by (5.3), and h x = 1/(n x + 1) and h y = 1/(n y + 1) are the mesh widths in the x and y directions. Unless α = β, τ /h α x and τ /h β y cannot both be independent of n; we choose τ = 1/ n α x . Note that the theory for A R still applies in this case. Stated CPU times and iteration counts are again for the first time step.
It is too costly to approximate A M by a banded Toeplitz matrix or a multigrid method, simply because it is expensive to obtain the Fourier coefficients of |f |, and so we present results for a multigrid approximation to A R only below. We also apply the nonsymmetric block circulant C n = I nxny − I ny ⊗ C x − C y ⊗ I nx preconditioner, and symmetric positive definite block circulant |C n | = I nxny + I ny ⊗ |C x | + |C y | ⊗ I nx , where C x and C y are Strang circulant approximations to L x and L y , respectively. Our multigrid method comprises 4 pre-and 4 post-smoothing steps, and a damping parameter of 0.9. The coarsest grid has n x = n y = 7.
The results in Table 5.6 show that the multigrid approximation of A R gives meshsize independent iteration counts, and that MINRES with this preconditioner is the fastest method for larger problems. For the block circulant preconditioners we see different behaviour depending on whether α > β. Specifically, when α > β, τ /h β y → 0 as n → ∞, which makes this problem easier to solve in some sense. On the other hand, when α < β the problems become harder to solve as n increases, and the block circulant with LSQR and MINRES suffer from growing iteration counts.
6. Conclusions. In this paper we presented two novel ideal preconditioners for (multilevel) Toeplitz matrices by considering the generating function f . The first, A R is formed using the real part of f . While it works best when the (multilevel) Toeplitz matrix is close to symmetric, it is reasonably robust with respect to the degree of nonsymmetry. This performance is likely attributable to the eigenvalue distribution, which remains bounded away from the origin. Our second preconditioner, A M , is based on |f |. Its performance is less affected by nonsymmetry, but it is more challenging to construct efficient approximations to A M in the multilevel case.
Our numerical results not only illustrate the effectiveness of the preconditioners, they highlight the value of symmetrization, which enables us to compute bounds on convergence rates that depend only on the scalar function f . Additionally, the combi- Table 5.6: Iteration numbers and CPU times (in parentheses) for the circulant preconditioners C n and |C n |, and multigrid preconditioners when d + = 2, d − = 0.5, e + = 0.3 and e − = 1 for Example 5.3. nation of symmetrization and preconditioned MINRES can be more computationally efficient than applying GMRES or LSQR to these problems.