Multilevel Tau preconditioners for symmetrized multilevel Toeplitz systems with applications to solving space fractional diffusion equations

In recent years, there has been growing interest in developing effective preconditioners for symmetrized Toeplitz systems. The underlying idea stems from [32], in which the original (real) nonsymmetric Toeplitz matrix $T_{m}\in\mathbb{R}^{m\times m}$ is symmetrized by premultiplying it with the anti-identity matrix

i.e., $[Y_{m}]_{j,k}=1$ if and only if $j+k=m+1$ and $[Y_{m}]_{j,k}=0$ elsewhere. For the now symmetric system with $Y_{m}T_{m}$, the minimal residual (MINRES) method can be used. Its convergence behaviour is related only to the eigenvalues and therefore effective preconditioners can be constructed exploiting the spectral information.

While the generalized minimal residual (GMRES) method is a common preconditioned solver applied to the original nonsymmetric Toeplitz system with $T_{m}$, it has a significant drawback: its preconditioning is largely heuristic, as discussed in [39]. Besides, in general, the convergence behaviors of GMRES cannot be rigorously analyzed solely through the knowledge of the eigenvalues, according to [14].

For the symmetrized Toeplitz matrix $Y_{m}T_{m}$, various preconditioning techniques have been proposed. Absolute value circulant preconditioners were first proposed in [32], with the preconditioned matrix eigenvalues shown to cluster around $\pm 1$ under certain conditions. Additionally, Toeplitz and band-Toeplitz preconditioners have been proposed in [31] and [18], respectively.

The same symmetrization preconditioning approach has also been employed in solving evolutionary differential equations. Initially introduced for solving ordinary and partial differential equations (PDEs) in [26] and [27], respectively, this approach has since been extended to a number of PDE problems, as discussed in [17] and [16], for example.

black Notably, this symmetrization approach was further generalized to the multilevel Toeplitz case in [31]. The author proposed a multilevel Toeplitz preconditioner as an ideal preconditioner for a class of symmetrized multilevel Toeplitz systems. Although this approach was shown to be effective in some cases, as demonstrated numerically in the same work, its overall preconditioning effect and the issues related to fast implementation were not fully explored.

Motivated by the need to address these issues, in this work, we first rigorously demonstrate that this ideal preconditioner can indeed achieve optimal convergence when employing the MINRES method, with a convergence rate that is independent of the mesh size. Then, to further illustrate its applicability and develop a fast implementation strategy, we consider solving Riemann-Liouville (R.-L.) fractional diffusion equations as an application. The findings demonstrate that Tau preconditioning constitutes an optimal strategy for symmetrized multilevel Toeplitz systems that emerge from the discretization of fractional diffusion equations.

It should be noted that, although some commonly used preconditioners, such as the aforementioned circulant and Toeplitz preconditioners, have been proposed for symmetrized multilevel Toeplitz systems, the performance of Tau preconditioners remains unexplored in this context. This is noteworthy considering their demonstrated effectiveness in symmetric Toeplitz systems, as studied in several works (see, for example, [6, 7, 33]).

As a model application, we consider the following space R.-L. fractional diffusion equation

where $\Omega=\prod_{i=1}^{d}(a_{i},b_{i})$ is an open hyper-rectangle, $\partial\Omega$ denotes the boundary of $\Omega$, $\alpha_{i}\in(1,2)$ are the fractional derivative orders, $f(x{,t})$ is the source term, $u_{0}(x)$ is a given function, and the diffusion coefficients $d_{i,\pm}$ are nonnegative constants. The left- and the right-sided R.-L. fractional derivatives in (1.2) are defined by

respectively, where ${\Gamma(\cdot)}$ denotes the gamma function.

Denote the set of all positive integers and the set of all nonnegative integers by $\mathbb{N}^{+}$ and $\mathbb{N}$, respectively. For any $m,k\in\mathbb{N}$ with $m\leq k$, define the set $m\wedge k:=$ $\{m,m+1,\ldots,k-1,k\}$. For $n_{i}\in\mathbb{N}^{+},i\in 1\wedge d$, denote by $h_{i}=\left(b_{i}-a_{i}\right)/\left(n_{i}+1\right)$, the stepsize along the $i$-th direction. Denote

black The Crank–Nicolson method is applied for the temporal derivative, and the \colorblack weighted and shifted Grünwald scheme \colorblack [37, 38] is used for the space fractional derivatives. \colorblack The time interval $[0,T]$ is partitioned into $M$ sub-intervals, each with a time step size of $\tau:=T/M$. This results in a scheme that is \colorblack second-order accurate in both time and space, leading to the linear system

where $\nu=\frac{1}{\tau}$ and $I_{m}$ represents an $m\times m$ identity matrix. The vector \colorblack $\mathbf{f}^{k+\frac{1}{2}}$ is known from the numerical scheme used and denotes the sampling of $f$ at the grid points. The matrix $A_{n}$ is a nonsymmetric multilevel Toeplitz matrix. More precisely, the coefficient matrix $A_{n}$ can be expressed as \colorblack

in which ‘ $\otimes$ ’ denotes the Kronecker product. The matrix $L_{n_{i}}^{(\alpha_{i})}$ is defined as

black where the coefficients $w_{0}^{(\alpha_{i})}=\frac{\alpha_{i}}{2}g_{0}^{(\alpha_{i})},w_{k}^{(\alpha_{i})}=\frac{\alpha_{i}}{2}g_{k}^{(\alpha_{i})}+\frac{2-\alpha_{i}}{2}g_{k-1}^{(\alpha_{i})}$ for $k\geq 1$. Note that the coefficients $g_{k}^{(\alpha_{i})}$ are given by

for $k\geq 1$.

From [38] and following the notation to be introduced in Section 3, we know that $L_{n_{i}}^{(\alpha_{i})}=T_{n_{i}}[g_{\alpha_{i}}]$ is a Toeplitz matrix generated by the complex-valued function

black

Thus, the generating function of $A_{n}$ is associated by

where $\bm{\alpha}=(\alpha_{1},\dots,\alpha_{d})$.

black Alternatively, if a backward Euler method is applied for the temporal derivative, and the shifted Grünwald scheme [24, 25] is used for the space fractional derivatives. This result in a scheme that is first-order accurate in both time and space, leading to the linear system

In this case, the nonsymmetric multilevel Toeplitz matrix $\widetilde{A}_{n}$ can be expressed as

The matrix $\widetilde{L}_{n_{i}}^{(\alpha_{i})}$ is defined as

where the coefficients $\widetilde{g}_{k}^{(\alpha_{i})}$ given by $\widetilde{g}_{k}^{(\alpha_{i})}=(-1)^{k}\binom{\alpha_{i}}{k}=\frac{(-1)^{k}}{k!}\alpha_{i}(\alpha_{i}-1)\cdots(\alpha_{i}-k+1),~{}k\geq 0,$ where $\binom{\alpha_{i}}{0}=1$.

From [8], we know that $\widetilde{L}_{n_{i}}^{(\alpha_{i})}=T_{n_{i}}[\widetilde{g}_{\alpha_{i}}]$ is a Toeplitz matrix generated by the complex-valued function

Thus, the generating function of $\widetilde{A}_{n}$ is associated by

Instead of solving \colorblack (1.3) (or (1.9)) directly, we employ the MINRES method to solve the following equivalent system denoted by \colorblack

In an endeavor to expedite the convergence of the MINRES algorithm applied to the symmetrized system $Y_{n}A_{n}\mathbf{x}=Y_{n}\mathbf{b}$ \colorblack (or $Y_{n}\widetilde{A}_{n}{\mathbf{x}}=Y_{n}{\mathbf{b}})$, the absolute value circulant preconditioner proposed in [32] was considered as a potential solution. However, as demonstrated in the numerical experiments of [31], its effectiveness has been shown to be unsatisfactory. Furthermore, a comprehensive analysis was provided in [18] showing that circulant preconditioners do not generally ensure rapid convergence. This is particularly true for ill-conditioned Toeplitz systems, where the preconditioned matrix may exhibit eigenvalue outliers that are very close to zero. Moreover, circulant preconditioners are known to be suboptimal in the context of multilevel Toeplitz systems, as discussed in [34, 35]. As an alternative, [18] proposed a band-Toeplitz preconditioner for symmetrized Toeplitz systems, specifically when the generating functions exhibit zeros of even order. Yet, this approach is not applicable for the matrix $Y_{n}A_{n}$ considered in this work. This is because the function $f_{\bm{\alpha}}$ associated with $Y_{n}A_{n}$ has a zero of fractional order $\bm{\alpha}$ as discussed in [8], a scenario where the band-Toeplitz preconditioner does not apply.

Recently, the use of the symmetric part of the multilevel Toeplitz coefficient matrix as a preconditioner for the target nonsymmetric multilevel Toeplitz systems was proposed in [31]. While this ideal \colorblack symmetric part preconditioner can facilitate fast MINRES convergence, the computationally expensive nature of its implementation limits practical application. Moreover, the effectiveness of the preconditioner has been observed numerically but has not yet been theoretically proven.

In this context, the main contributions of our work are twofold:

1. 1.

   We show that the ideal preconditioner proposed in [31] can indeed lead to a MINRES convergence rate that is independent of the mesh size for \colorblack a class of symmetrized multilevel Toeplitz systems. Moreover, building on this preconditioner, we develop a practical and optimal preconditioned MINRES method for these symmetrized systems (see Corollary 2).

2. 2.
   \color

   black To illustrate our preconditioning approach, we propose a multilevel Tau preconditioner as an practical preconditioner for $Y_{n}A_{n}\mathbf{x}=Y_{n}\mathbf{b}$ \colorblack (or $Y_{n}\widetilde{A}_{n}{\mathbf{x}}=Y_{n}{\mathbf{b}})$. This preconditioner balances the effectiveness of preconditioning with computational feasibility. Our approach offers a practical strategy for preconditioning fractional diffusion equations and achieves mesh-independent MINRES convergence. To the best of our knowledge, this study is the first to show that Tau preconditioning is an optimal choice for \colorblack a range of symmetrized multilevel Toeplitz systems stemming from space R.-L. fractional diffusion equations, \colorblack where both first and second order schemes are considered.

The paper is organized as follows: Our proposed preconditioner is defined in Section 2. Section 3 reviews some preliminary results on multilevel Toeplitz matrices. Section 4 presents our main results on the effectiveness of our proposed preconditioners. Numerical examples in Section 5 demonstrate the expected performance of our proposed preconditioners, and support the theoretical results concerning the eigenvalues of the associated preconditioned matrices.

In this section, our multilevel Tau preconditioner $P_{n}$ for (1.13) is presented.

For a symmetric Toeplitz matrix $T_{m}\in\mathbb{R}^{m\times m}$ with $(t_{1},t_{2},...,t_{m})^{\top}\in\mathbb{R}^{m}$, define its $\tau$-matrix [1] approximation as

where $H_{m}$ is the Hankel matrix with $(t_{3},t_{4},...,t_{m},0,0)^{\top}$ as its first column and $(0,0,t_{m},...,t_{4},t_{3})^{\top}$ as its last column. A crucial property of the Tau matrix defined in (2.1) is that it is diagonalizable by sine transform matrix, i.e.,

where $Q_{m}=[{\rm diag}(q_{i})]_{i=1}^{m}$ is a diagonal matrix with

is a sine transform matrix. It is easy to verify that $S_{m}$ is a symmetric orthogonal matrix, i.e., $S_{m}=S_{m}^{\top}=S_{m}^{-1}$. The product of matrix $S_{m}$ and a given vector of length $m$ can be fast computed in $\mathcal{O}(m\log m)$ operations using discrete sine transform (DSTs) [2]. Let ${\bf e}_{m,i}\in\mathbb{R}^{m}$ denotes the $i$-th column of the $m\times m$ identity matrix. We also note that the $m$ numbers $\{q_{i}\}_{i=1}^{m}$ defined in (2.3) can be computed by

From the equation above, we know that the computation of $\{q_{i}\}_{i=1}^{m}$ requires only $\mathcal{O}(m\log m)$ operations.

For a real square matrix $Z$, denote the symmetric part of $Z$ as

Then, our multilevel Tau preconditioner $P_{n}$ for $Y_{n}A_{n}$ in (1.13) is defined as follows

black Similarly, in the first-order case, a multilevel Tau preconditioner $\widetilde{P}_{n}$ for $Y_{n}\widetilde{A}_{n}$ in (1.13) can be defined as follows

From (2.2) & (2.3) and properties of the one-dimensional sine transform matrix $S_{m}$, we know that $P_{n}$ is diagonalizable by a $d$-dimension sine transform matrix, namely,

where $\Lambda_{i}$ contains the eigenvalues of $\tau\left(\mathcal{H}\left(L_{n_{i}}^{\left(\alpha_{i}\right)}\right)\right)$. Thus, the product of $P_{n}$ and a given vector can be efficiently computed in $\mathcal{O}(n\log n)$ operations using DSTs.

In Section 4, we will show that the eigenvalues of the preconditioned matrix $P_{n}^{-1}Y_{n}A_{n}$ \colorblack (or $\widetilde{P}_{n}^{-1}Y_{n}\widetilde{A}_{n}$) are contained in a disjoint interval enclosing $\pm 1$, which leads to theoretically guaranteed matrix size/mesh-independent convergence when MINRES is applied.

In this section, we provide some useful background knowledge regarding multilevel Toeplitz and Hankel matrices.

Now consider the Banach space $L^{1}([-\pi,\pi]^{k})$ of all complex-valued Lebesgue integrable functions over $[-\pi,\pi]^{k}$, equipped with the norm

where ${\rm d}{\theta}={\rm d}{\theta_{1}}\cdots{\rm d}{\theta_{k}}$ denotes the volume element with respect to the $k$-dimensional Lebesgue measure.

Let $f:$ $[-\pi,\pi]^{k}\to\mathbb{C}$ be a function belonging to $L^{1}([-\pi,\pi]^{k})$ and periodically extended to $\mathbb{R}^{k}$. The multilevel Toeplitz matrix $T_{{n}}[f]$ of size $n\times n$ with $n=n_{1}n_{2}\dots n_{k}$ is defined as

where

are the Fourier coefficients of $f$ and $J^{j}_{m}$ is the $m\times m$ matrix whose $(l,h)$-th entry equals 1 if $(l-h)=j$ and $0$ otherwise. The function $f$ is called the generating function of $T_{{n}}[f]$.

It is easy to prove that (see e.g., [28, 3, 4, 13]) if $f$ is real-valued, then $T_{{n}}[f]$ is Hermitian; if $f$ is real-valued and nonnegative, but not identically zero almost everywhere, then $T_{{n}}[f]$ is Hermitian positive definite; if $f$ is real-valued and even, $T_{{n}}[f]$ is (real) symmetric.

Throughout this work, we assume that $f\in L^{1}([-\pi,\pi]^{k})$ and is periodically extended to $\mathbb{R}^{k}$.

Similar to a multilevel Toeplitz matrix, we can define a multilevel Hankel matrix as

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. Clearly, a multilevel Hankel matrix is symmetric.

Multilevel Toeplitz matrices can be symmetrized by the permutation matrix $Y_{n}\in\mathbb{R}^{n\times n}$, namely, $Y_{n}=Y_{n_{1}}\otimes\dots\otimes Y_{n_{k}}$. Knowing that $Y_{r}J_{r}^{(k)}=K_{r}^{(r-k)}$, we can easily show that

where $b_{(j_{1},\dotsc,j_{k})}=a_{(n_{1}-j_{1},\dotsc,n_{k}-j_{k})}$. Hence, $Y_{n}T_{n}[f]$ is a symmetric multilevel Hankel matrix. For more properties regarding multilevel Hankel matrices, see [10].

A crucial aspect of developing effective preconditioners for $Y_{n}T_{n}[f]$ is understanding its asymptotic spectral distribution associated with $f$. This was established for the uni-level case in [22, 12, 15] and later generalized to the multilevel case in [23, 11].

black

The main results of this work are divided into the following subsections.