Preconditioners for symmetrized Toeplitz and multilevel Toeplitz matrices

TL;DR

This paper introduces novel preconditioners for symmetrized Toeplitz and multilevel Toeplitz matrices, enabling efficient Krylov subspace methods like MINRES with rigorous iteration bounds.

Contribution

It proposes new ideal preconditioners for symmetrized Toeplitz matrices and analyzes their spectra, addressing a gap in effective preconditioning techniques.

Findings

Preconditioners improve convergence of Krylov methods.

Spectral analysis supports preconditioner effectiveness.

Numerical experiments validate theoretical results.

Abstract

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.

Figures

Tables

Table 1. Table 1 : Iteration numbers and CPU times (in parentheses) for the optimal circulant preconditioner C n subscript 𝐶 𝑛 C_{n} and tridiagonal preconditioner A R subscript 𝐴 𝑅 A_{R} for Example   5.1 .

$n$	GMRES				LSQR				MINRES
	$C_n$		$A_R$		$C_n$		$A_R$		$\left|C_n\right|$		$A_R$
1023	37	(0.058)	67	(0.067)	105	(0.15)	62	(0.05)	82	(0.057)	68	(0.029)
2047	48	(0.13)	68	(0.12)	186	(0.4)	67	(0.081)	111	(0.12)	70	(0.046)
4095	62	(0.25)	69	(0.22)	—	—	73	(0.13)	170	(0.21)	71	(0.065)
8191	—	—	72	(0.51)	—	—	78	(0.2)	—	—	72	(0.13)

Table 2. Table 2 : Iteration numbers and CPU times (in parentheses) for the exact preconditioner A M subscript 𝐴 𝑀 A_{M} and its multigrid approximation M ​ G ​ ( A M ) 𝑀 𝐺 subscript 𝐴 𝑀 MG(A_{M}) for Example   5.1 . The second column shows the time needed to compute the elements of A M subscript 𝐴 𝑀 A_{M} .

$n$	$A_M$ time	GMRES				LSQR				MINRES
		$A_n$		$MG(A_n$)		$A_n$		$MG\left(A_M\right)$		$A_M$		$MG\left(A_M\right)$
1023	6.4	1	(0.06)	39	(0.15)	1	(0.068)	33	(0.11)	11	(0.11)	24	(0.044)
2047	21	1	(0.28)	41	(0.28)	1	(0.33)	37	(0.24)	11	(0.6)	24	(0.082)
4095	73	1	(1.9)	39	(0.36)	1	(1.9)	39	(0.27)	12	(3.6)	25	(0.096)
8191	$2.5\times {10}^2$	1	(8.7)	42	(0.88)	1	(11)	43	(0.67)	12	( 22)	25	(0.21)

Table 3. Table 3 : Iteration numbers and CPU times (in parentheses) for the Strang circulant C n subscript 𝐶 𝑛 C_{n} , absolute value Strang circulant | C n | subscript 𝐶 𝑛 |C_{n}| and multigrid preconditioners when d + = 0.5 subscript 𝑑 0.5 d_{+}=0.5 and d − = 1 subscript 𝑑 1 d_{-}=1 for Example   5.2 .

$\alpha$	$n$	GMRES				LSQR				MINRES
		$C_n$		$MG\left(A_n\right)$		$C_n$		$MG\left(A_n\right)$		$\left|C_n\right|$		$MG\left(A_M\right)$		$MG\left(A_R\right)$
1.25	1023	5	(0.01)	4	(0.016)	6	(0.011)	6	(0.02)	10	(0.0084)	12	(0.18)	8	(0.014)
	4095	6	(0.017)	4	(0.045)	6	(0.016)	6	(0.053)	10	(0.013)	12	(0.18)	8	(0.043)
	16383	6	(0.065)	4	(0.17)	6	(0.066)	7	(0.22)	10	(0.054)	13	(0.32)	8	(0.17)
	65535	6	(0.25)	4	(0.66)	6	(0.25)	7	(0.76)	9	(0.19)	13	(0.82)	8	(0.6)
	262143	6	(0.99)	4	(4.4)	6	( 1)	7	( 5)	9	(0.72)	13	(4.5)	8	(4.3)
1.5	1023	6	(0.0062)	4	(0.021)	6	(0.0062)	7	(0.025)	10	(0.0048)	13	(0.37)	8	(0.013)
	4095	6	(0.018)	4	(0.046)	6	(0.018)	7	(0.061)	10	(0.015)	13	(0.5)	8	(0.044)
	16383	6	(0.062)	4	(0.17)	6	(0.067)	7	(0.21)	9	(0.05)	13	(0.76)	9	(0.19)
	65535	6	(0.24)	5	(0.7)	7	(0.28)	8	(0.8)	9	(0.19)	13	(1.4)	9	(0.66)
	262143	6	(0.93)	5	(5.2)	7	(1.1)	8	(5.7)	9	(0.72)	15	(6.1)	9	(4.7)
1.75	1023	6	(0.0085)	5	(0.043)	7	(0.0088)	7	(0.021)	9	(0.0062)	13	(1.6)	9	(0.014)
	4095	6	(0.015)	5	(0.046)	7	(0.015)	8	(0.058)	9	(0.01)	15	(2.2)	9	(0.036)
	16383	6	(0.062)	5	(0.2)	7	(0.075)	8	(0.24)	9	(0.049)	15	(3.2)	10	(0.21)
	65535	6	(0.24)	5	(0.71)	7	(0.28)	8	(0.81)	9	(0.19)	15	(4.8)	11	(0.75)
	262143	6	(0.9)	5	(5.2)	7	(1.1)	9	(6.3)	9	(0.72)	16	( 11)	11	(5.7)

Table 4. Table 4 : Iteration numbers and CPU times (in parentheses) for the Strang circulant C n subscript 𝐶 𝑛 C_{n} , absolute value Strang circulant | C n | subscript 𝐶 𝑛 |C_{n}| and multigrid preconditioners when α = 1.5 𝛼 1.5 \alpha=1.5 for Example   5.2 .

$(d_{+},d_{-})$	$n$	GMRES				LSQR				MINRES
		$C_n$		$MG\left(A_n\right)$		$C_n$		$MG\left(A_n\right)$		$\left|C_n\right|$		$MG\left(A_M\right)$		$MG\left(A_R\right)$
(0,3)	4095	5	(0.031)	5	(0.053)	7	(0.02)	5	(0.059)	10	(0.016)	10	(0.47)	13	(0.049)
	16383	4	(0.044)	5	(0.21)	7	(0.078)	5	(0.23)	10	(0.058)	10	(0.83)	13	(0.27)
	65535	4	(0.18)	6	(0.84)	7	(0.29)	6	(0.93)	10	(0.22)	10	(1.6)	14	(0.97)
	262143	4	(0.75)	6	(6.2)	7	(1.2)	6	(6.7)	11	(0.92)	11	(6.8)	14	(7.1)
(1,3)	4095	7	(0.015)	5	(0.042)	7	(0.014)	5	(0.045)	10	(0.013)	10	(0.4)	9	(0.037)
	16383	7	(0.072)	5	(0.21)	7	(0.078)	5	(0.23)	11	(0.06)	10	(0.81)	10	(0.21)
	65535	7	(0.29)	5	(0.71)	8	(0.33)	5	(0.77)	11	(0.24)	11	(1.6)	10	(0.7)
	262143	7	(1.1)	6	(6.2)	8	(1.3)	6	(6.6)	11	(0.93)	11	(6.6)	10	(5.3)
(1,1)	4095	6	(0.013)	4	(0.031)	6	(0.013)	5	(0.041)	10	(0.01)	9	(0.39)	9	(0.034)
	16383	6	(0.064)	4	(0.17)	6	(0.068)	5	(0.23)	10	(0.058)	9	(0.79)	9	(0.19)
	65535	6	(0.25)	4	(0.57)	6	(0.26)	5	(0.77)	9	(0.19)	9	(1.4)	9	(0.64)
	262143	6	(0.93)	5	(5.2)	7	(1.2)	5	(5.8)	9	(0.79)	9	(5.9)	9	(4.9)

Table 5. Table 5 : Upper bound in Eq.   11 for Example   5.2 .

$\alpha$	$(d_{+},d_{-})$
	(0,3)	(1,3)	(0.5,1)	(1,1)
1	1.13	0.67	0.25	0.00
1.25	0.70	0.39	0.17	0.00
1.5	0.42	0.23	0.11	0.00
1.75	0.20	0.11	0.05	0.00

Table 6. Table 6 : Iteration numbers and CPU times (in parentheses) for the circulant preconditioners C n subscript 𝐶 𝑛 C_{n} and | C n | subscript 𝐶 𝑛 |C_{n}| , and multigrid preconditioners when d + = 2 subscript 𝑑 2 d_{+}=2 , d − = 0.5 subscript 𝑑 0.5 d_{-}=0.5 , e + = 0.3 subscript 𝑒 0.3 e_{+}=0.3 and e − = 1 subscript 𝑒 1 e_{-}=1 for Example   5.3 .

$(\alpha ,\beta )$	$n^2$	GMRES				LSQR				MINRES
		$C_n$		$MG\left(A_n\right)$		$C_n$		$MG\left(A_n\right)$		$\left|C_n\right|$		$MG\left(A_R\right)$
(1.5,1.25)	961	16	(0.032)	5	(0.011)	23	(0.033)	5	(0.014)	42	(0.028)	12	(0.013)
	16129	15	(0.12)	5	(0.058)	21	(0.11)	6	(0.07)	39	(0.12)	12	(0.07)
	261121	14	(1.5)	5	(1.1)	18	(1.4)	6	(1.3)	34	(1.5)	12	(1.0)
(1.5,1.75)	961	21	(0.029)	4	(0.0086)	28	(0.038)	4	(0.0099)	43	(0.027)	10	(0.01)
	16129	21	(0.16)	4	(0.051)	35	(0.2)	5	(0.065)	57	(0.19)	10	(0.049)
	261121	20	(2.1)	5	(1.2)	40	(3.1)	5	(1.0)	67	(2.8)	12	(0.97)

Equations

$$A_{n}x=b,$$

$$Y_n=\left[\begin{array}{ccc}& & 1\\ & \text{iddots}& \\ 1& & \end{array}\right],$$

$$Y_n{A}_{n}x=Y_{n}b,$$

$$A_n=\left[\begin{array}{ccccc}{a}_0& a_{-1}& \dots & a_{-n+2}& a_{-n+1}\\ a_1& a_0& a_{-1}& & a_{-n+2}\\ ⋮& a_1& a_0& \ddots & ⋮\\ a_{n-2}& & \ddots & \ddots & a_{-1}\\ a_{n-1}& a_{n-2}& \dots & a_1& a_0\end{array}\right].$$

$$a_k=\frac{1}{2\pi }{\int }_{-\pi }^{\pi }f(\theta ){\mathrm{e}}^{-\mathrm{i}k\theta }\mathrm{d}\theta ,k\in \mathbb{Z}.$$

$$f(\theta )=\sum _{k=-\infty }^{\infty }{a}_k{\mathrm{e}}^{-\mathrm{i}k\theta },\sum _{k=-\infty }^{\infty }|a_k|<\infty .$$

$$C_n=\left[\begin{array}{ccccc}{c}_0& c_{n-1}& \dots & c_2& c_1\\ c_1& c_0& c_{n-1}& & c_2\\ ⋮& c_1& c_0& \ddots & ⋮\\ c_{n-2}& & \ddots & \ddots & c_{n-1}\\ c_{n-1}& c_{n-2}& \dots & c_1& c_0\end{array}\right].$$

$$(F_n{)}_{j,k}=\frac{1}{\sqrt{n}}{\mathrm{e}}^{\frac{2\pi \mathrm{i}jk}{n}},j,k=0,\dots ,n-1,$$

$${\lambda }_k=\sum _{j=0}^{n-1}{c}_j{\mathrm{e}}^{\frac{2\pi \mathrm{i}jk}{n}}.$$

$$|C_n|=F_n^{\ast }|{\Lambda }_n|F_n.$$

$$H_n=\left[\begin{array}{ccccc}{a}_1& a_2& a_3& \dots & a_n\\ a_2& a_3& & \text{iddots}& a_{n+1}\\ a_3& & \text{iddots}& \text{iddots}& ⋮\\ ⋮& \text{iddots}& \text{iddots}& & a_{2n-2}\\ a_n& a_{n+1}& \dots & a_{2n-2}& a_{2n-1}\end{array}\right],$$

$$Y_n{A}_n=A_n^T{Y}_n.$$

$$G_n{Y}_n=Y_n{G}_n$$

$$a_j=a_{(j_1,\dots ,j_p)}=\frac{1}{(2\pi {)}^p}{\int }_{[-\pi ,\pi {]}^p}f(\theta ){\mathrm{e}}^{-\mathrm{i}⟨\theta ,j⟩}\mathrm{d}\theta ,j\in {\mathbb{Z}}^p,$$

$$A_n(f)=\sum _{j_1=-n_1+1}^{n_1-1}\cdots \sum _{j_p=-n_p+1}^{n_p-1}{J}_{n_1}^{(j_1)}\otimes \cdots \otimes J_{n_p}^{(j_p)}{a}_{(j_1,\dots ,j_p)}.$$

$$H_n(f)=\sum _{j_1=1}^{2n_1-1}\cdots \sum _{j_p=1}^{2n_p-1}{K}_{n_1}^{(j_1)}\otimes \cdots \otimes K_{n_p}^{(j_p)}{a}_{(j_1,\dots ,j_p)},$$

$$Y_n{A}_n(f)=$$

$$=$$

$$r=\underset{x\in [-\pi ,\pi ]}{\mathrm{essinf}}\frac{f(\theta )}{g(\theta )},R=\underset{\theta \in [-\pi ,\pi ]}{\mathrm{esssup}}\frac{f(\theta )}{g(\theta )}.$$

$$A_R^{-\frac{1}{2}}{A}_n{A}_R^{-\frac{1}{2}}=I_n+E_n,$$

$$\Vert E_n{\Vert }_2=ϵ<\underset{\theta \in [-\pi ,\pi ]}{\mathrm{esssup}}\left|\frac{f_I(\theta )}{f_R(\theta )}\right|.$$

$$A_R^{-\frac{1}{2}}{A}_n{A}_R^{-\frac{1}{2}}=A_R^{-\frac{1}{2}}(A_R+\mathrm{i}{A}_n(f_I))A_R^{-\frac{1}{2}}=I_n+E_n,$$

$$ϵ<\underset{\theta \in [-\pi ,\pi ]}{\mathrm{esssup}}\left|\frac{f_I(\theta )}{f_R(\theta )}\right|,$$

$$A_R^{-\frac{1}{2}}(Y_n{A}_n)A_R^{-\frac{1}{2}}=Y_n+Y_n{E}_n,$$

$$\Vert Y_n{E}_n{\Vert }_2=ϵ<\underset{\theta \in [-\pi ,\pi ]}{\mathrm{esssup}}\left|\frac{f_I(\theta )}{f_R(\theta )}\right|.$$

$$A_R^{-\frac{1}{2}}(Y_n{A}_n)A_R^{-\frac{1}{2}}=Y_n(A_R^{-\frac{1}{2}}{A}_n{A}_R^{-\frac{1}{2}})=Y_n(I_n+E_n)=Y_n+Y_n{E}_n.$$

$$ϵ<\underset{\theta \in [-\pi ,\pi ]}{\mathrm{esssup}}\left|\frac{f_I(\theta )}{f_R(\theta )}\right|.$$

$$A_R^{-\frac{1}{2}}(Y_n{A}_n)A_R^{-\frac{1}{2}}=Y_n+Y_n{E}_n,$$

$$\Lambda_{n}=\blockarray{cccc}&m_{1}m_{1}m_{2}\\ \block{c[ccc]}m_{1}\Lambda_{\mathrm{pos}}\\ m_{1}-\Lambda_{\mathrm{pos}}\\ m_{2}0\\$$

$$U_{n}=\blockarray{cccc}m_{1}&m_{1}m_{3}m_{4}\\ \block{[cccc]}U_{\mathrm{pos}}Y_{n}U_{\mathrm{pos}}U_{\text{sym}}+Y_{n}U_{\text{sym}}U_{\text{skew}}-Y_{n}U_{\text{skew}}\\ ,$$

Full text

{AMS}

65F08, 65F10, 15B05, 35R11

1 Introduction

Linear systems

[TABLE]

where $A_{n}\in\mathbb{R}^{n\times n}$ is a Toeplitz or multilevel Toeplitz matrix, and $b\in\mathbb{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, 27]. Additionally, demand for fast numerical methods for fractional diffusion problems—which have recently received significant attention—has renewed interest in the solution of Toeplitz and Toeplitz-like systems [10, 26, 31, 32, 46].

Preconditioned iterative methods are often used to solve systems of the form Eq. 1. When $A_{n}$ is Hermitian, CG [18] and MINRES [29] can be applied, and their descriptive convergence rate bounds guide the construction of effective preconditioners [7, 27]. 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 [35] for Toeplitz matrices, and discussed in Section 2.2 for the multilevel case, $A_{n}$ is symmetrized by the exchange matrix

[TABLE]

so that Eq. 1 can be replaced by

[TABLE]

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 Eq. 3 is easier to solve than Eq. 1 [12, 23]. Instead, the presence of $Y_{n}$ allows us to use preconditioned MINRES, with its nice properties and convergence rate bounds, to solve Eq. 3. We can then apply a secondary preconditioner $P_{n}\in\mathbb{R}^{n\times 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 [37] 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 [35] 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 Eq. 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.

2.1 Toeplitz and Hankel matrices

Let $A_{n}\in\mathbb{R}^{n\times n}$ be the nonsingular Toeplitz matrix

[TABLE]

In many applications, the matrix $A_{n}$ is associated with a generating function $f\in L^{1}([-\pi,\pi])$ via its Fourier coefficients

[TABLE]

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

[TABLE]

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, 33, 47].

Circulant matrices are Toeplitz matrices of the form

[TABLE]

They are diagonalized by the Fourier matrix, i.e., $C_{n}=F_{n}^{\ast}\Lambda_{n}F_{n}$, where

[TABLE]

and $\Lambda_{n}=\operatorname{diag}(\lambda_{0},\dotsc,\lambda_{n-1})$, with

[TABLE]

We denote by $C_{n}(f)$ the circulant with eigenvalues $\lambda_{j}=f(2\pi j/n)$, $j=0,\dotsc,n-1$. The absolute value circulant [5, 28, 35] derived from a circulant $C_{n}$ is the matrix

[TABLE]

Closely related to Toeplitz matrices are Hankel matrices $H_{n}\in\mathbb{R}^{n\times n}$,

[TABLE]

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 Eq. 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

[TABLE]

Alternatively, we may think of $A_{n}$ being self-adjoint with respect to the bilinear form induced by $Y_{n}$ [14, 34].

A matrix $G_{n}\in\mathbb{R}^{n\times n}$ is centrosymmetric if

[TABLE]

and is skew-centrosymmetric if $G_{n}Y_{n}=-Y_{n}G_{n}$. Thus, Eq. 8 shows that symmetric Toeplitz matrices are centrosymmetric. It is clear from Eq. 9 that the inverse of a nonsingular centrosymmetric matrix is again centrosymmetric. Furthermore, nonsingular centrosymmetric matrices have a centrosymmetric square root [22, Corollary 1]111In [22] the proof is given only for a centrosymmetric matrix of even dimension. However, the extension to matrices of odd dimension is straightforward..

2.2 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},\dotsc,j_{p})\in\mathbb{Z}^{p}$ be a multi-index and consider a $p$-variate function $f\in L^{1}([-\pi,\pi]^{p})$, $f:[-\pi,\pi]^{p}\rightarrow\mathbb{C}$. The Fourier coefficients of $f$ are defined as

[TABLE]

where $\langle\theta,j\rangle=\sum_{k=1}^{p}\theta_{k}j_{k}$ and $\mathrm{d}\theta=\mathrm{d}\theta_{1}\dots\mathrm{d}\theta_{p}$ is the volume element with respect to the $p$-dimensional Lebesgue measure.

If $n=(n_{1},\dotsc,n_{p})\in\mathbb{N}^{p}$, with $n_{i}>1$, $i=1,\dotsc,p$, and $\pi(n)=n_{1}\dots n_{p}$, then $f$ is the generating function of the multilevel Toeplitz matrix $A_{n}(f)\in\mathbb{R}^{\pi(n)\times\pi(n)}$, where

[TABLE]

Here, $J^{(k)}_{r}\in\mathbb{R}^{r\times 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

[TABLE]

where $K^{(k)}_{r}\in\mathbb{R}^{r\times 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}\in\mathbb{R}^{\pi(n)\times\pi(n)}$, $Y_{n}=Y_{n_{1}}\otimes\dots\otimes Y_{n_{p}}$. 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_{r}^{(k)}=K_{r}^{(r-k)}$, so that

[TABLE]

where $b_{(j_{1},\dotsc,j_{p})}=a_{(n_{1}-j_{1},\dotsc,n_{p}-j_{p})}$. 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}([-\pi,\pi]^{p})$. We denote the real and imaginary parts of $f$ by $f_{R}$ and $f_{I}$, respectively, so that $f=f_{R}+\mathrm{i}f_{I}$. We assume that the symmetric part of $A_{n}$, given by $A_{R}=(A_{n}+A_{n}^{T})/2$, is positive definite, which is equivalent to requiring that $f_{R}$ is essentially positive. Similarly, we assume that $|f|\geq\delta>0$ for some constant $\delta$, so that $A_{n}(|f|)$ is positive definite with $\lambda_{\min}(A_{n}(|f|)\geq\delta$. Moreover, $\lambda_{\min}(A_{n}(|f|)>\delta$ if $\operatorname*{esssup}|f|>\delta=\operatorname*{essinf}|f|$ (see Lemma 3.1).

3 Preconditioning Toeplitz matrices

In this section we introduce our ideal preconditioners for Eq. 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_{n}^{T})/2$, which was previously used to precondition the nonsymmetric problem Eq. 1 (see [21]). When applied to the symmetrized system Eq. 3, spectral information can be used to bound the convergence rate of preconditioned MINRES. Accordingly, in this section we characterize the eigenvalues of $A_{R}^{-\frac{1}{2}}Y_{n}A_{n}A_{R}^{-\frac{1}{2}}$.

We begin by stating a powerful result that characterizes the spectra of preconditioned Hermitian Toeplitz matrices in terms of generating functions.

Lemma 3.1** ([39, Theorem 3.1]).**

Let $f,g\in L^{1}([-\pi,\pi])$ 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 $A_{n}^{-1}(g)A_{n}(f)$ lie in $(r,R)$, where $r<R$ and

[TABLE]

*If $r=R$ then $A_{n}^{-1}(g)A_{n}(f)=I_{n}$, the identity matrix of dimension $n$. *

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) preconditioned matrix $A_{R}^{-\frac{1}{2}}A_{n}A_{R}^{-\frac{1}{2}}$ from the identity.

Lemma 3.2**.**

Let $f\in L^{1}([-\pi,\pi])$, and let $f=f_{R}+\mathrm{i}f_{I}$, where $f_{R}$ and $f_{I}$ are real-valued functions with $f_{R}$ essentially positive. Additionally, let $A_{n}:=A_{n}(f)\in\mathbb{R}^{n\times n}$ be the Toeplitz matrix associated with $f$. Then $A_{R}=A_{n}(f_{R})=(A_{n}+A_{n}^{T})/2$ is symmetric positive definite and

[TABLE]

where

[TABLE]

Proof 3.3**.**

It is easily seen from Eq. 5 that $A_{n}(f)=A_{n}(f_{R})+\mathrm{i}A_{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

[TABLE]

where $E_{n}=\mathrm{i}{\widehat{E}_{n}}=\mathrm{i}A_{R}^{-\frac{1}{2}}A_{n}(f_{I})A_{R}^{-\frac{1}{2}}$.

To bound $\epsilon:=\|E_{n}\|_{2}=\|{\widehat{E}_{n}}\|_{2}$, note that since ${\widehat{E}_{n}}$ is Hermitian, $\|{\widehat{E}_{n}}\|_{2}$ is equal to the spectral radius of ${\widehat{E}_{n}}$. Applying Lemma 3.1 thus gives that

[TABLE]

*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}^{-\frac{1}{2}}A_{n}A_{R}^{-\frac{1}{2}}$, 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 Eq. 3.

Lemma 3.4**.**

Let $f\in L^{1}([-\pi,\pi])$, and let $f=f_{R}+\mathrm{i}f_{I}$, where $f_{R}$ and $f_{I}$ are real-valued functions with $f_{R}$ essentially positive. Additionally, let $A_{n}:=A_{n}(f)\in\mathbb{R}^{n\times n}$ be the Toeplitz matrix associated with $f$. Then the symmetric positive definite matrix $A_{R}=A_{n}(f_{R})=(A_{n}+A_{n}^{T})/2$ is such that

[TABLE]

where $Y_{n}$ is the exchange matrix in Eq. 2 and

[TABLE]

Proof 3.5**.**

Since $A_{R}$ is a symmetric Toeplitz matrix, it is centrosymmetric. Hence, $A_{R}^{-\frac{1}{2}}$ is centrosymmetric (see Eq. 9 and [22]), so that $Y_{n}A_{R}^{-\frac{1}{2}}=A_{R}^{-\frac{1}{2}}Y_{n}$. Combining this with Lemma 3.2 shows that

[TABLE]

*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 Eq. 10 shows that the eigenvalues of $A_{R}^{-\frac{1}{2}}(Y_{n}A_{n})A_{R}^{-\frac{1}{2}}$ lie in $[-1-\epsilon,-1+\epsilon]\cup[1-\epsilon,1+\epsilon]$. However, as $\epsilon$ grows, eigenvalues could move close to the origin, and hamper MINRES convergence. We will show that this cannot happen in the following result.

Theorem 3.6**.**

Let $f\in L^{1}([-\pi,\pi])$, and let $f=f_{R}+\mathrm{i}f_{I}$, where $f_{R}$ and $f_{I}$ are real-valued functions with $f_{R}$ essentially positive. Additionally, let $A_{n}:=A_{n}(f)\in\mathbb{R}^{n\times n}$ be the Toeplitz matrix associated with $f$ and let $A_{R}=A_{n}(f_{R})=(A_{n}+A_{n}^{T})/2$. Then, the eigenvalues of $A_{R}^{-\frac{1}{2}}(Y_{n}A_{n})A_{R}^{-\frac{1}{2}}$ lie in $[-1-\epsilon,-1]\cup[1,1+\epsilon]$, where

[TABLE]

Proof 3.7**.**

We know from Lemma 3.4 that

[TABLE]

where $\|Y_{n}E_{n}\|_{2}<\epsilon$ and $Y_{n}$ has eigenvalues $\pm 1$. Thus, as discussed above, the eigenvalues of $A_{R}^{-\frac{1}{2}}(Y_{n}A_{n})A_{R}^{-\frac{1}{2}}$ lie in $[-1-\epsilon,-1+\epsilon]\cup[1-\epsilon,1+\epsilon]$. 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_{R}^{-\frac{1}{2}}(Y_{n}A_{n})A_{R}^{-\frac{1}{2}}$ in terms of the eigenvalues of $Y_{n}E_{n}$.

Before applying our first similarity transform, we recall from the proofs of Lemmas 3.2 and 3.4 that $\mathrm{i}A_{n}(f_{I})$ is skew-symmetric and Toeplitz, while $A_{R}^{-\frac{1}{2}}$ is symmetric and centrosymmetric. It follows that $Y_{n}E_{n}$ is symmetric but skew-centrosymmetric. Skew-centrosymmetry implies that whenever $(\lambda,x)$, $\lambda\neq 0$ is an eigenpair of $Y_{n}E_{n}$, then so is $(-\lambda,Y_{n}x)$ [19, 36, 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\pm Y_{n}x$, $x\in\mathbb{R}^{n}$ [44, Theorem 17]. Therefore, $Y_{n}E_{n}$ has eigendecomposition $Y_{n}E_{n}=U_{n}\Lambda_{n}U_{n}^{T}$, where

[TABLE]

and

[TABLE]

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

[TABLE]

Using the orthogonality of the columns of $U_{n}$, it is straightforward to show that

[TABLE]

Thus, $U_{n}^{T}Y_{n}U_{n}=Q_{n}\Gamma_{n}Q_{n}^{T}$, where

[TABLE]

Here, $\widehat{\Gamma}_{2m_{1}}=\operatorname{diag}(1,-1,\dotsc,1,-1)$ and the $k$th column of $\widehat{Q}\in\mathbb{R}^{n\times 2m_{1}}$ is given by

[TABLE]

with $e_{j}\in\mathbb{R}^{n}$ the $j$th unit vector. Consequently, our second similarity transform gives

[TABLE]

with

[TABLE]

where if $\Lambda_{\mathrm{pos}}=\operatorname{diag}(\lambda_{1},\lambda_{2},\dotsc,\lambda_{m_{1}})$ then

[TABLE]

Hence, letting

[TABLE]

we find that

[TABLE]

*Since the eigenvalues of $Z+\Sigma_{k}$ are $\pm\sqrt{1+\lambda_{k}^{2}}$ we see from Eq. 15 that the eigenvalues of $A_{R}^{-\frac{1}{2}}(Y_{n}A_{n})A_{R}^{-\frac{1}{2}}$ are $\pm\sqrt{1+\lambda_{k}^{2}}$, $k=1,\dotsc,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.6 characterizes the eigenvalues of $A_{R}^{-\frac{1}{2}}(Y_{n}A_{n})A_{R}^{-\frac{1}{2}}$, and hence the convergence rate of preconditioned MINRES, in terms of the scalar quantity in Eq. 11. 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.6 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 $A_{M}^{-1}Y_{n}A_{n}$.

Theorem 3.8**.**

Assume that $f\in L^{\infty}([-\pi,\pi])$ with $0<\delta\leq|f(\theta)|$ for all $\theta\in[-\pi,\pi]$. Then, if $A_{M}=A_{n}(|f|)$,

[TABLE]

*where $\widetilde{f}=f/|f|$ and $\|E_{n}\|_{2}=o(n)$ as $n\rightarrow\infty$. Moreover, the eigenvalues of $Y_{n}A_{n}(\widetilde{f})$ lie in $[-1,1]$. *

Proof 3.9**.**

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],

[TABLE]

where $\|\widetilde{E}_{n}\|_{2}=o(n)$ as $n\rightarrow\infty$. Since $A_{M}$ is Hermitian and Toeplitz, both $A_{M}$ and its inverse are centrosymmetric. It follows that

[TABLE]

*where $E_{n}=Y_{n}\widetilde{E}_{n}$ and $\|E_{n}\|_{2}=\|Y\widetilde{E}_{n}\|_{2}=\|\widetilde{E}_{n}\|_{2}$. Hence $\|E_{n}\|_{2}=o(n)$ as $n\rightarrow\infty$. Since $Y_{n}$ is unitary and $Y_{n}A_{n}(\widetilde{f})$ is symmetric, the absolute values of the eigenvalues of $Y_{n}A_{n}(\widetilde{f})$ coincide with the singular values of $A_{n}(\widetilde{f})$, which in turn are bounded above by one [47]. This proves the result. *

A consequence of Theorem 3.8 is that the eigenvalues of $A_{n}(|f|)^{-1}Y_{n}A_{n}$ lie in $[-1-\epsilon,1+\epsilon]$, where for large enough $n$ the parameter $\epsilon$ is small. Although eigenvalues may be close to the origin, most cluster at $-1$ and 1, in line with Theorem 3.4 in [23].

To conclude this section, we show how $A_{M}$ can be approximated by circulant preconditioners. First, recall from Section 2.1 that $C_{n}(f)$ is the preconditioner with eigenvalues $\lambda_{j}=f(2\pi j/n)$, $j=0,\dotsc,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_{n}^{(S)}|$ [5, 28, 35], where if $C_{n}^{(S)}$ is the Strang circulant preconditioner [43] for $A_{n}$, with eigenvalues $\lambda_{j}$, $j=1,\dotsc,n$, then the corresponding absolute value circulant preconditioner $|C_{n}^{(S)}|$ has eigenvalues $|\lambda_{j}|$, $j=1,\dotsc,n$. For this preconditioner, we obtain the following result.

Theorem 3.10**.**

*Let $f:[-\pi,\pi]\rightarrow\mathbb{C}$ be in the Wiener class and let $A_{n}=A_{n}(f)\in\mathbb{R}^{n\times n}$. Then the Strang preconditioner $C_{n}^{(S)}$, is such that $|C_{n}^{(S)}|\rightarrow C_{n}(|f|)$ as $n\rightarrow\infty$. *

Proof 3.11**.**

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_{n}^{(S)}=C_{n}(\mathcal{D}_{m}\star f)$, where

[TABLE]

is the convolution of $f$ with the Dirichlet kernel $\mathcal{D}$ [6], and $a_{k}$ is as in Eq. 5.

Since $|C_{n}^{(S)}|=|C_{n}(\mathcal{D}_{m}\star f)|$ and $C_{n}(|f|)$ are both diagonalized by the Fourier matrix, they will be identical if all their eigenvalues, defined by Eq. 6, match. The eigenvalues of $C_{n}(\mathcal{D}_{m}\star f)$ are $(\mathcal{D}_{m}\ast f)(2\pi j/n)$, $j=0,\dotsc,n-1$. Hence the $j$th eigenvalue of $|C_{n}(\mathcal{D}_{m}\star f)|$ is

[TABLE]

Since $f$ is in the Wiener class, $(\mathcal{D}_{m}\star f)(\theta)$ converges absolutely, hence uniformly to $f(\theta)$. Thus,

[TABLE]

*Since the eigenvalues of $|C_{n}(\mathcal{D}_{m}\star f)|$ approach those of $C_{n}(|f|)$ as $n\rightarrow\infty$ we obtain the result. *

4 Multilevel Toeplitz problems

We now extend the results of Section 3 to multilevel Toeplitz matrices.

4.1 The preconditioner $A_{R}$

The results of Section 3.1 carry over straightforwardly to the multilevel case. They depend on the following generalization of Lemma 3.1. The result essentially appeared in Theorem 2.4222Although the result is stated for $f,g$ nonnegative, the proof also holds for indefinite $f$. in [38].

Lemma 4.1** ([38]).**

Let $f,g\in L^{1}([-\pi,\pi]^{p})$ with $g$ essentially positive. Let

[TABLE]

*Then the eigenvalues of $A_{n}^{-1}(g)A_{n}(f)$ lie in $(r,R)$ if $r<R$. If $r=R$ then $A_{n}^{-1}(g)A_{n}(f)=I_{n}$, where $I_{n}$ is the identity matrix of dimension $\pi(n)=n_{1}\cdots n_{p}$. *

With this result, Lemmas 3.2, 3.4 and 3.6 carry over directly to the multilevel case. In particular, we have the following characterization of the eigenvalues of $A_{R}^{-\frac{1}{2}}(Y_{n}A_{n})A_{R}^{-\frac{1}{2}}$.

Theorem 4.2**.**

Let $f\in L^{1}([-\pi,\pi]^{p})$, and let $f=f_{R}+\mathrm{i}f_{I}$, where $f_{R}$ and $f_{I}$ are real-valued functions with $f_{R}$ essentially positive. Additionally, let $A_{n}:=A_{n}(f)\in\mathbb{R}^{\pi(n)\times\pi(n)}$ be the multilevel Toeplitz matrix associated with $f$ and let $A_{R}=A_{n}(f_{R})=(A_{n}+A_{n}^{T})/2$. Then, the eigenvalues of $A_{R}^{-\frac{1}{2}}(Y_{n}A_{n})A_{R}^{-\frac{1}{2}}$ lie in $[-1-\epsilon,-1]\cup[1,1+\epsilon]$, where

[TABLE]

Theorem 4.2 characterizes the eigenvalues of $A_{R}^{-\frac{1}{2}}(Y_{n}A_{n})A_{R}^{-\frac{1}{2}}$, 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 Eq. 16.

4.2 The preconditioner $A_{M}$

We can also extend the results in Section 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^{\infty}([-\pi,\pi]^{p})$, with $0<\delta<|f(\theta)|$ for all $\theta\in[-\pi,\pi]^{p}$. Then, if $A_{M}=A_{n}(|f|)$ is the multilevel Toeplitz matrix generated by $|f|$,

[TABLE]

*where $\tilde{f}=f/|f|$ and $\|E_{n}\|_{2}=o(n)$ as $n\rightarrow\infty$. Moreover, the eigenvalues of $Y_{n}A_{n}(\tilde{f})$ lie in $[-1,1]$. *

5 Numerical experiments

In this section we investigate the effectiveness of the preconditioners described above, and approximations to them, for the symmetrized system Eq. 3. We also compare the proposed approach to using nonsymmetric preconditioners for Eq. 1 within preconditioned GMRES and LSQR [30]. 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.20GHz333Code is available from https://github.com/jpestana/fracdiff. 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,\dotsc,1)^{T}/\sqrt{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 $\omega$ 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 (3) we use the absolute value preconditioner in Eq. 7 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 Eq. 1 when GMRES is applied. The Toeplitz coefficient matrix $A_{n}=A_{n}(f)$ is formed from the generating function $f(\theta)=(2-2\cos(\theta))(1+\mathrm{i}\theta)$. 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}(\theta)=2-2\cos(\theta)$ 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 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_{R}^{-1}A_{n}$. Theorem 3.6 tells us that these eigenvalues lie in $[-1-\pi,-1]\cup[1,1+\pi]$, and Fig. 1 (b) shows that these bounds are tight. As discussed in [21], the eigenvalues of $A_{R}^{-1}A_{n}$ are also nicely clustered (see Fig. 1 (a)), with real part $1$ and imaginary part in $[-\pi,\pi]$. Although we cannot rigorously link this eigenvalue characterization to the rate of GMRES convergence, Table 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(\theta)|=(2-2\cos(\theta))\sqrt{1+\theta^{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 $\omega=0.1$ for GMRES, $\omega=0.4$ for LSQR and $\omega=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 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_{M}^{-1}Y_{n}A_{n}$, when $n=2047$, are as expected from Theorem 3.8 (see Fig. 2), since all eigenvalues lie in $[-1,1]$. Indeed, most cluster at the endpoints of this interval. The eigenvalues of $A_{M}^{-1}A_{n}$ are also localized, but not as clustered, indicating that the spectrum of $A_{M}^{-1}A_{n}$ may differ significantly from that of $A_{M}^{-1}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 problem is to find $u(x,t)$ that satisfies

[TABLE]

where $\alpha\in(1,2)$, and $d_{+}$ and $d_{-}$ are nonnegative constants. We impose the absorbing boundary conditions $u(x\leq 0,t)=u(x\geq 1,t)=0$, $t\in[0,1]$, while $u(x,0)=80\sin(20x)\cos(10x),\,x\in[0,1].$ The Riemann-Liouville derivatives in Eq. 17 are

[TABLE]

where $n$ is the integer for which $n-1<\alpha\leq n$.

Discretising by the shifted Grünwald-Letnikov method in space, and the backward Euler method in time [24, 25], gives the linear system

[TABLE]

[TABLE]

where $g_{\alpha,k}=(-1)^{k}\binom{\alpha}{k}$, $\nu=\frac{\tau}{h^{\alpha}}$ and $h=\frac{1}{n+1}$. We set $\tau=1/\lceil n^{\alpha}\rceil$, which makes $\nu$ constant, so that all the theory of Section 3 can be directly applied, but comparable results are obtained $\tau=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 Eq. 18 are generated by [10]

[TABLE]

The real part of $\varphi$ is essentially positive, so $A_{R}=(A_{n}+A_{n}^{T})/2$ is positive definite. However, since $A_{R}$ is dense we approximate it by our V-cycle multigrid method (analysed in [32]) with the coarsest grid of dimension 127, 2 pre- and 2 post-smoothing steps, and $\omega=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 Section 3.2. The second is multigrid (with the same parameters as for $A_{R}$, except that we use 1 pre- and post-smoothing 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 $\alpha=1.25$ we compute the first 50 elements in $r$ and $c$, and when $\alpha>1.25$ we take the first $\lceil\beta(1.1)^{\log_{2}(n+1)}\rceil$ elements in $r$ and $c$, where $\beta=40$ when $\alpha=1.5$ and $\beta=100$ when $\alpha=1.75$. This balances the time to compute these coefficients, and the resulting MINRES iteration count.

We see from Table 3 that our approximations to $A_{R}$ and $A_{M}$ are robust with respect to $n$, but both require slightly more iterations for larger $\alpha$. 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 $\alpha\leq 1.5$, while the multigrid approximation of $A_{n}$ with GMRES is slightly faster for large $\alpha$.

In Table 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_{+}\leq 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 4 are explained by Theorem 3.6 and the relatively small upper bound Eq. 11, which describes how far eigenvalues of $A_{R}^{-1}Y_{n}A_{n}$ can deviate from 1 in magnitude. This bound is 0 when $d_{+}=d_{-}$, or when $\alpha=2$, since then $A_{n}$ is symmetric. However, Table 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 $\Omega=(0,1)^{2}\times(0,1]$ that satisfies

[TABLE]

where $\alpha,\beta\in(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 [24, 25], which leads to the following linear system:

[TABLE]

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,

[TABLE]

where $L_{\alpha}$ is given by Eq. 19, 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 $\alpha=\beta$, $\tau/h_{x}^{\alpha}$ and $\tau/h_{y}^{\beta}$ cannot both be independent of $n$; we choose $\tau=1/\lceil n_{x}^{\alpha}\rceil$. 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_{n_{x}n_{y}}-I_{n_{y}}\otimes C_{x}-C_{y}\otimes I_{n_{x}}$ preconditioner, and symmetric positive definite block circulant $|C_{n}|=I_{n_{x}n_{y}}+I_{n_{y}}\otimes|C_{x}|+|C_{y}|\otimes I_{n_{x}}$, 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 6 show that the multigrid approximation of $A_{R}$ gives mesh-size 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 $\alpha>\beta$. Specifically, when $\alpha>\beta$, $\tau/h_{y}^{\beta}\rightarrow 0$ as $n\rightarrow\infty$, which makes this problem easier to solve in some sense. On the other hand, when $\alpha<\beta$ 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 combination of symmetrization and preconditioned MINRES can be more computationally efficient than applying GMRES or LSQR to these problems.

Acknowledgments

The author would like to thank Mariarosa Mazza and Stefano Serra Capizzano for their careful reading of an earlier version of this manuscript and helpful discussions, and the anonymous referees for their suggestions. All data underpinning this publication are openly available from Zenodo at http://doi.org/10.5281/zenodo.1327565.