A note on parallel preconditioning for the all-at-once solution of Riesz fractional diffusion equations

The governing Riesz fractional diffusion equations ¹¹1Here it is worth noting that if the Riesz fractional derivative is replaced with the fractional Laplacian, we can just follow the work in [50] for the spatial discretization of one- or two-dimensional model problem and our proposed numerical techniques in the next sections are easy and simple to be adapted for such an extension. of the anomalous diffusion process can be written as

where $u(x,t)$ may represent, for example, a solute concentration, and constant $\kappa_{\gamma}>0$ the diffusion coefficient. This equations is a superdiffusive model largely used in fluid flow analysis, financial modelling and other applications. The Riesz fractional derivative $\partial^{\gamma}u(x,t)/\partial|x|^{\gamma}$ is defined by [51]

in which ${}_{a}D^{\gamma}_{x}$ and ${}_{x}D^{\gamma}_{b}$ are the left and right Riemann-Liouville fractional derivatives of order $\gamma\in(1,2)$ given, respectively, by [5]

and

As $\gamma\rightarrow 2$, it notes that Eq. (2.1) degenerates into the classical diffusion equation [51].

Next, we focus on the numerical solution of Eq. (2.1). We consider a rectangle $\bar{Q}_{T}=\{(x,t):~{}a\leq x\leq b,~{}0\leq t\leq T\}$ discretized on the mesh $\varpi_{h\tau}=\varpi_{h}\times\varpi$, where $\varpi_{h}=\{x_{i}=a+ih,~{}i=0,1,\cdots,N;~{}h=(b-a)/N\}$, and $\varpi_{\tau}=\{t_{k}=k\tau,~{}j=0,1,\cdots,M;~{}\tau=T/M\}$. We denote by $v=\{v_{i}~{}|~{}0\leq i\leq N\}$ any grid function.

A considerable amount of work has been devoted in the past years to the development of fast methods for the approximation of the Riesz and Riemann-Liouville (R-L) fractional derivatives, such as the first-order accurate shifted Grünwald approximation [3, 51] and the second-order accurate weighted-shifted Grünwald difference (WSGD) approximation [9, 52]. Due to the relationship between the two kinds of fractional derivatives, all the numerical schemes proposed for the R-L fractional derivatives can be easily adapted to approximate the Riesz fractional derivative. Although the solution approach described in this work can accomodate any spatial discretized method, we choose the so-called fractional centred difference formula [8] of the Riesz fractional derivatives for clarity.

For any function $u(x)\in L^{1}(\mathbb{R})$, we denote

where the $\gamma$-dependent weight coefficient is defined as

As noted in [8], $\omega^{(\gamma)}_{\ell}=\mathcal{O}(\ell^{-1-\gamma})$ and the fractional centred difference formula $\Delta^{\gamma}_{h}u(x)$ exists for any $u(x)\in L^{1}(\mathbb{R})$. Some properties of the coefficient $\omega^{(\gamma)}_{\ell}$ and the operator $\Delta^{\gamma}_{h}u(x)$ are presented in the following lemmas.

At this stage, let $u(x,t)\in\mathcal{C}^{4,3}_{x,t}([a,b]\times[0,T])$ be a solution to the problem (2.1) and consider Eq. (2.1) at the set of grid points $(x,t)=(x_{i},t_{k})\in\bar{Q}_{T}$ with $i=1,2,\cdots,N-1,~{}k=1,\cdots,M$. The first-order time derivative at the point $t=t_{k}$ is approximated by the second-order backward difference formula, i.e.,

then we define $U^{k}_{i}=u(x_{i},t_{k})$ and $f^{k}_{i}=f(x_{i},t_{k})$ to obtain

where $\{R^{k}_{i}\}$ are small and satisfy the inequality

Note that we cannot omit the small term in the derivation of a two-level difference scheme that is not selfstarting from Eq. (2.7) due to the unknown information of $u(x_{i},t_{1})$. One of the most popular strategies to compute the first time step solution $u^{1}_{i}$ is to use a first-order backward Euler scheme (which is only used once). This yields the following implicit difference scheme [37, 10] for Eq. (2.1):

where

In order to implement the proposed scheme, here we define ${\bm{u}}^{k}=[u(x_{1},t_{k}),u(x_{2},t_{k}),\cdots,u(x_{N-1},t_{k})]^{\top}$ and ${\bm{f}}^{k}=[f(x_{1},t_{k}),f(x_{2},t_{k}),\cdots,f(x_{N-1},t_{k})]^{\top}$; then $\kappa_{\gamma}\Delta^{\gamma}_{h}u^{k}_{i}$ can be written into the matrix-vector product form $A{\bm{u}}^{k}$, where

It is easy to prove that $T_{x}$ is a symmetric positive definite (SPD) Toeplitz matrix (see [8]). Therefore, it can be stored with only $N-1$ entries and the fast Fourier transforms (FFTs) can be applied to carry out the matrix-vector product in $\mathcal{O}(N\log N)$ operations. The matrix-vector form of the 2-step BDF method with start-up backward Euler method for solving the model problem (2.1) can be formulated as follows:

We refer to matrix $A$ as the Jacobian matrix. In particular, we note that the numerical scheme (2.10) has been proved to be unconditionally stable and second-order convergent in the discrete $L^{2}$ norm [10]. Note that it is computationally more efficient than the C-N scheme presented in [8] as it requires one less matrix-vector multiplication per time step.

Instead of computing the solution of (2.10) step-by-step, we try to get an all-at-once approximation by solving the following linear system:

where $I_{s}$ and $I_{t}$ are two identity matrices of order $N_{s}~{}(=N-1)$ and $N_{t}~{}(=M)$, respectively. We denote the above linear system as follows:

where $\mathcal{A}=C\otimes I_{s}-\tau I_{t}\otimes A$ with

and it is clear that $\mathcal{A}$ is invertible, because its all diagonal blocks (i.e, either $I_{s}-\tau A$ or $\frac{3}{2}I_{s}-\tau A$) are invertible [33, 34].

In this subsection, we investigate the properties of the discrete all-at-once formulation (2.11). This will guide us to discuss the design of an efficient solver for such a large linear system.

Proof. Consider the following matrix splitting and then perform the matrix factorization,

where the vector ${\bm{e}}_{1}=[1,0,\cdots,0]^{T}\in\mathbb{R}^{N_{t}}$. According to the Sherman-Morrison formula [53], we can write

On the other hand, we can compute the inverse of $\hat{C}$ as follows,

Therefore, we know that

and

The explicit expression of $C^{-1}$ has the following form:

Since $\hat{C}^{-1}$ is a lower triangular Toeplitz matrix, $\|C^{-1}\|_{1}$ is the absolute sum of elements of its last row, i.e.,

Analogously, $\|C^{-1}\|_{\infty}$ is the absolute sum of elements of its first column, i.e.,

Therefore, the above estimates are proved. $\Box$

Proof. Since $Re(\lambda(A))\leq 0$, we may claim that $\|\mathcal{A}{\bm{z}}\|\geq\|(C\otimes I_{s}){\bm{z}}\|$ for any vector ${\bm{z}}$. In particular, using the properties of the Kronecker product, we may estimate $\sigma_{min}(\mathcal{A})\geq\sigma_{min}(C)$. The latter is computed using the spectral norm of $C^{-1}$, which can be bounded by the following inequality:

According to Lemma 2.3, we have

which proves the second bound.

The first estimate is derived similarly. By applying the triangle inequality, we can write

whereas

so the spectral norms are not greater than 4. Finally, the proof is completed by recalling that $\tau=T/N_{t}$. $\Box$

The condition number can be expected to vary linearly with the main properties of the system, such as number of time steps, length of time interval, and norm of the Jacobian matrix [39, 40, 54]. More refined estimates could be provided by taking into account particular properties of $A$. In Theorem 2.1, the time interval $[0,T]$ may not necessarily be equal to the whole observation range of an application. The global time scheme (2.11) could be applied by splitting the desired interval $[0,\hat{T}]$ into a sequence of subintervals $[0,T],[T,2T],\cdots,[\hat{T}-T,\hat{T}]$ for “large” time steps of size $T$ each, solving by (2.11) for $[(q-1)T,qT],~{}q=1,\cdots,\hat{T}/T$, extracting the last snapshot $x_{Nt}$ and restarting the method using $x_{Nt}$ as the initial state for the next interval. The optimal value of $T$ should provide the fastest computation [39]. Another way to accelerate the solution is to reduce the condition number of the linear system by preconditioning. This standard computational aspect will be considered in the next section.

Since the matrix $C$ defined in (2.13) can be viewed as the sum of a Toeplitz matrix and a rank-1 matrix, it is natural to define our first structuring-circulant preconditioner as

where

is a $\alpha$-circulant matrix of Strang type that can be diagonalized as

with

The matrix $\mathbb{F}_{N_{t}}$ in (3.4b) is the discrete Fourier matrix, and is unitary. In addition, it is worth noting that if we choose $\alpha$ very small, then the condition number of $V_{\alpha}$ increases quickly and therefore large roundoff errors may be introduced in our proposed method; see [36, 19] for a discussion of these issues. Using the property of the Kronecker product, we can factorize $\mathcal{P}_{\alpha}$ as

and this implies that we can compute ${\bm{z}}=\mathcal{P}^{-1}_{\alpha}{\bm{v}}$ via the following three steps:

where ${\bm{z}}_{j}=({\bm{z}}_{j,1}^{\top},{\bm{z}}_{j,2}^{\top},\dots,{\bm{z}}_{j,N_{t}}^{\top})^{\top}$ (with $j=1,2$) and $\lambda^{(\alpha)}_{n}$ is the $n$-th eigenvalue of $C_{\alpha}$. The first and third steps only involve matrix-vector multiplications with $N_{t}$ subvectors, and these can be computed simultaneously by the FFTs²²2It notes that there are also parallel FFTs available at http://www.fftw.org/parallel/parallel-fftw.html. in $N_{t}$ CPUs [54]. The major computational cost is spent in the second step, however this step is naturally parallel for all the $N_{t}$ time steps.

Next, we study the invertibility of the matrix $\mathcal{P}_{\alpha}$. According to Eq. (2.13) and because matrix $A$ is negative definite, the following result assures the invertibility of $\mathcal{P}_{\alpha}$³³3It is unnecessary to keep that $r_{0}\geq|r_{1}|+|r_{2}|$, which is introduced in [35, Remark 4]..

In practice, we often choose $\alpha\in(0,1)$ [35]; then it holds that $\mathcal{R}e(\lambda^{(\alpha)}_{n})>0$ according to Proposition 3.1. In Fig. 1, we plot the complex quantities $\{\lambda^{(\alpha)}_{n}\}^{N_{t}}_{n=1}$ on the complex plane. We see that for $\alpha\in(0,1)$, it holds $\mathcal{R}e(\lambda^{(\alpha)}_{n})>0$ and consequently all the linear systems involved at Step-(b) are positive definite, thus not difficult to solve [31]. As shown in [36], a multigrid method using Richardson iterations as a smoother with an optimized choice of the damping parameter is very effective. By the way, a detailed convergence analysis with multigrid precondioning in the context of the $\tau$ algebra considered in this paper is presented in [38].

In addition, it is not hard to derive the following result.

Theorem 3.1 implies that the generalized minimal residual (GMRES) method by Saad and Schultz [56] can compute the exact solution of the preconditioned system $\mathcal{P}^{-1}_{\alpha}\mathcal{A}$ (if it is diagonalizable) in at most $2N_{t}+1$ iterations [56, 23, 33]. Although it should be noted that such a convergence behavior of the preconditioned GMRES is not sharp when $N_{t}$ is not small, the result can give useful clues on the effectiveness of the preconditioner $\mathcal{P}_{\alpha}$ to approximate the all-at-once matrix $\mathcal{A}$. Moreover, in Section 4 we will show numerically that most of the eigenvalues of the preconditioned matrix $\mathcal{P}^{-1}_{\alpha}\mathcal{A}$ actually cluster at point 1 of the spectrum, supporting the theoretical findings of this section on the good potential of the proposed preconditioner $\mathcal{P}_{\alpha}$ to accelerate the iterative solution of $\mathcal{P}^{-1}_{\alpha}\mathcal{A}$.