Positive approximations of the inverse of fractional powers of SPD M-matrices

TL;DR

This paper introduces a new rational approximation method for efficiently computing the inverse fractional powers of SPD M-matrices, preserving positivity and providing error bounds, with applications in fractional diffusion problems.

Contribution

It develops and analyzes positive rational approximations of fractional matrix powers using BURA, ensuring preservation of key properties and offering sharp error estimates.

Findings

The method preserves positivity of the approximation.

Error bounds are established for the rational approximations.

Numerical tests confirm theoretical properties and effectiveness.

Abstract

This study is motivated by the recent development in the fractional calculus and its applications. During last few years, several different techniques are proposed to localize the nonlocal fractional diffusion operator. They are based on transformation of the original problem to a local elliptic or pseudoparabolic problem, or to an integral representation of the solution, thus increasing the dimension of the computational domain. More recently, an alternative approach aimed at reducing the computational complexity was developed. The linear algebraic system $\cal A^\alpha \bf u=\bf f$, $0< \alpha <1$ is considered, where $\cal A$ is a properly normalized (scalded) symmetric and positive definite matrix obtained from finite element or finite difference approximation of second order elliptic problems in $\Omega\subset\mathbb{R}^d$, $d=1,2,3$. The method is based on best uniform rational approximations (BURA) of the function $t^{\beta-\alpha}$ for $0 < t \le 1$ and natural $\beta$. The maximum principles are among the major qualitative properties of linear elliptic operators/PDEs. In many studies and applications, it is important that such properties are preserved by the selected numerical solution method. In this paper we present and analyze the properties of positive approximations of $\cal A^{-\alpha}$ obtained by the BURA technique. Sufficient conditions for positiveness are proven, complemented by sharp error estimates. The theoretical results are supported by representative numerical tests.

Tables

Table 1. Table 1: Errors E α ​ ( k , k ; 1 ) subscript 𝐸 𝛼 𝑘 𝑘 1 E_{\alpha}(k,k;1) of BURA r α 1 ​ ( t ) subscript superscript 𝑟 1 𝛼 𝑡 r^{1}_{\alpha}(t) of t 1 − α superscript 𝑡 1 𝛼 t^{1-\alpha} on [ 0 , 1 ] 0 1 [0,1] .

$\alpha$	$E_{\alpha }(5,5,1)$	$E_{\alpha }(6,6,1)$	$E_{\alpha }(7,7,1)$
0.25	2.8676e-5	9.2522e-6	3.2566e-6
0.50	2.6896e-4	1.0747e-4	4.6037e-5
0.75	2.7162e-3	1.4312e-3	7.8966e-4

Table 2. Table 2: ℓ 2 superscript ℓ 2 \ell^{2} relative error ‖ 𝐮 h , r − 𝐮 h ‖ 2 ‖ 𝐟 2 ‖ 2 subscript norm subscript 𝐮 ℎ 𝑟 subscript 𝐮 ℎ 2 subscript norm subscript 𝐟 2 2 \frac{\|{\bf u}_{h,r}-{\bf u}_{h}\|_{2}}{\|{\bf f}_{2}\|_{2}} .

$h$	$\alpha =0.25$			$\alpha =0.5$			$\alpha =0.75$
	$k=5$	$k=6$	$k=7$	$k=5$	$k=6$	$k=7$	$k=5$	$k=6$	$k=7$
$2^{-5}$	6.3e-5	1.2e-4	7.5e-5	1.9e-4	3.1e-4	1.0e-4	8.5e-4	3.5e-4	2.1e-4
$2^{-6}$	6.8e-4	2.1e-4	1.2e-4	9.3e-4	6.2e-4	2.6e-4	3.0e-4	7.3e-4	1.9e-4
$2^{-7}$	4.9e-3	9.9e-4	2.2e-4	3.2e-3	1.3e-3	5.5e-4	1.6e-3	1.0e-3	6.6e-5
$2^{-8}$	1.2e-2	4.9e-3	1.0e-3	5.0e-3	2.4e-3	1.0e-3	2.7e-3	6.7e-4	4.4e-4
$2^{-9}$	3.7e-2	9.6e-3	4.9e-3	5.6e-3	1.3e-3	1.2e-3	8.2e-4	1.3e-3	1.1e-3
$2^{-10}$	2.1e-2	3.6e-2	9.5e-3	2.4e-2	8.8e-3	1.4e-3	5.8e-3	3.0e-3	1.5e-3
$2^{-11}$	1.9e-1	2.3e-2	3.6e-2	4.2e-2	1.4e-2	8.9e-3	1.6e-3	9.9e-4	4.3e-4

Equations

$$a(u,v):={\int }_{\Omega }\left(\mathbf{a}(x)\nabla u(x)\cdot \nabla v(x)+q(x)\right)dx={\int }_{\Omega }f(x)v(x)dx,\forall v\in V,$$

$$V:=\{v\in H^1(\Omega ):v(x)=0\text{mbox}on{\Gamma }_D\},$$

$$c\Vert \mathbf{z}{\Vert }^2\le {\mathbf{z}}^T\mathbf{a}(x)\mathbf{z}\le C\Vert \mathbf{z}{\Vert }^2\forall \mathbf{z}\in {\mathbb{R}}^d,\forall x\in \Omega ,$$

$${\mathcal{L}}^{\alpha }u(x)=\sum _{i=1}^{\infty }{\lambda }_i^{\alpha }{c}_i{\psi }_i(x),\text{mbox}whereu(x)=\sum _{i=1}^{\infty }{c}_i{\psi }_i(x).$$

$${\mathcal{L}}^{\alpha }u=f$$

$$a(u,v):={\int }_{\Omega }\nabla u(x)\cdot \nabla v(x)dx.$$

$$-div\left(y^{1-2\alpha }\nabla v(x,y)\right)=0,(x,y)\in \Omega \times {\mathbb{R}}_{+},$$

$$\underset{y\to \infty }{\lim }v(x,y)=0,x\in \Omega ,$$

$$\underset{y\to 0^{+}}{\lim }\left(-y^{1-2\alpha }{v}_y(x,y)\right)=f(x),x\in \Omega .$$

$$a(x)\frac{\partial u}{\partial n}+\mu (x)u=0,x\in \partial \Omega ,$$

$$u(x)=w(x,1),w(x,0)={\delta }^{-\alpha }f,$$

$$(t\mathcal{D}+\delta \mathcal{I})\frac{dw}{dt}+\alpha \mathcal{D}w=0,$$

$${\mathcal{L}}^{-\alpha }=\frac{2\sin (\pi \alpha )}{\pi }{\int }_0^{\infty }{t}^{2\alpha -1}{\left(\mathcal{I}+t^2\mathcal{L}\right)}^{-1}dt$$

$${\mathcal{A}}^{\alpha }\mathbf{u}=\mathbf{f},\text{mbox}where0<\alpha <1.$$

$$\mathbf{u}={\mathcal{A}}^{-\alpha }\mathbf{f}=\mathcal{W}{\mathcal{D}}^{-\alpha }{\mathcal{W}}^T\mathbf{f}.$$

$${\mathcal{A}}_e={\mathcal{K}}_e+{\mathcal{M}}_e,$$

$$\mathcal{R}(m,k):=\{P_m/Q_k:P_m\in {\mathcal{P}}_m,Q_k\in {\mathcal{P}}_k\},$$

$$\underset{r\in \mathcal{R}(m,k)}{\min }\underset{t\in [0,1]}{\max }\left|t^{\beta -\alpha }-r(t)\right|,$$

$$E_{\alpha }(m,k;\beta ):=\underset{t\in [0,1]}{\max }\left|t^{\beta -\alpha }-r_{\alpha }^{\beta }(t)\right|.$$

$$\Vert {\mathcal{A}}^{-\beta }{r}_{\alpha }^{\beta }(\mathcal{A})\mathbf{f}-{\mathcal{A}}^{-\alpha }\mathbf{f}{\Vert }_{{\mathcal{A}}^{\gamma +\beta }}\le E_{\alpha }(m,k;\beta )\Vert \mathbf{f}{\Vert }_{{\mathcal{A}}^{\gamma -\beta }},\forall \gamma \in \mathbb{R},\forall \mathbf{f}\in {\mathbb{R}}^N.$$

$$t^{-\beta }{r}_{\alpha }^{\beta }(t)=\sum _{j=0}^{m-k-\beta }{b}_j{t}^j+\sum _{j=1}^{\beta }\frac{c_{0,j}}{t^j}+\sum _{j=1}^k\frac{c_j}{t-d_j},$$

$${\mathbf{u}}_r:={\mathcal{A}}^{-\beta }{r}_{\alpha }^{\beta }(\mathcal{A})\mathbf{f}=\sum _{j=0}^{m-k-\beta }{b}_i{\mathcal{A}}^j\mathbf{f}+\sum _{j=1}^{\beta }{c}_{0,j}{\mathcal{A}}^{-j}\mathbf{f}+\sum _{j=1}^k{c}_j(\mathcal{A}-d_{j}I{)}^{-1}\mathbf{f}$$

$${\mathcal{A}}^{-\beta }{r}_{\alpha }^{\beta }(\mathcal{A})=\sum _{j=1}^{\beta }{c}_{0,j}{\mathcal{A}}^{-j}+\sum _{j=1}^k{c}_j(\mathcal{A}-d_{j}I{)}^{-1}.$$

$$0>{\zeta }_1>d_1>{\zeta }_2>d_2>\cdots >{\zeta }_k>d_k>-\infty$$

$$0$$

$${\eta }_j^{1-\alpha }-r_{\alpha }^1({\eta }_j)$$

$$r_{\alpha }^1(t)=b_0^{\ast }+\sum _{j=1}^k\frac{c_j^{\ast }}{t-d_j},t^{-1}{r}_{\alpha }^1(t)=\frac{c_{0,1}}{t}+\sum _{j=1}^k\frac{c_j}{t-d_j}.$$

$$c_{0,1}=E_{\alpha }(k,k;1),c_j=c_j^{\ast }/d_j,j=1,\dots ,k.$$

$$\frac{1}{t(t-d_j)}=\frac{1}{d_j}\left(\frac{1}{t-d_j}-\frac{1}{t}\right),j=1,\dots ,k.$$

$$c_{0,1}=b^{\ast}_{0}-\sum_{j=1}^{k}\frac{c^{\ast}_{j}}{d^{j}}=r^{1}_{\alpha}(0)=-\big{(}\eta_{1}^{1-\alpha}-r^{1}_{\alpha}(\eta_{1})\big{)}=E_{\alpha}(k,k;1).$$

$$\underset{k\to \infty }{\lim }{e}^{2\pi \sqrt{\alpha k}}{E}_{1-\alpha }(k,k;1)=4^{1+\alpha }|\sin \pi \alpha |$$

Full text

During last few years, several different techniques are proposed to localize the nonlocal fractional diffusion operator. They are based on transformation of the original problem to a local elliptic or pseudoparabolic problem, or to an integral representation of the solution, thus increasing the dimension of the computational domain.

More recently, an alternative approach aimed at reducing the computational complexity was developed. The linear algebraic system $\cal A^{\alpha}\bf u=\bf f$, $0<\alpha<1$ is considered, where $\cal A$ is a

properly normalized (scalded) symmetric and positive definite matrix obtained from finite element or finite difference

approximation of second order elliptic problems in

$\Omega\subset\mathbb{R}^{d}$,

$d=1,2,3$. The method is based on best uniform rational approximations (BURA) of the function $t^{\beta-\alpha}$ for $0<t\leq 1$ and natural $\beta$.

The maximum principles are among the major qualitative properties of linear elliptic operators/PDEs. In many studies and applications, it is important that such properties are preserved by the selected numerical solution method. In this paper we present and analyze the properties of positive approximations of $\cal A^{-\alpha}$ obtained by

the BURA technique. Sufficient conditions for positiveness are proven, complemented by sharp error estimates. The theoretical results are supported by representative numerical tests.

1 Introduction

This work is inspired by the recent development in the fractional calculus and its various applications, i.e., to Hamiltonian chaos, [28], anomalous diffusion in complex systems, [2], long-range interaction in elastic deformations, [23], nonlocal electromagnetic fluid flows, [18], image processing, [10]. A more recent impressive examples of anomalous diffusion models in chemical engineering are provided in [19]. Such kind of applications lead to fractional order partial differential equations that involve in general non-symmetric elliptic operators see, e.g. [13]. An important subclass of this topic are the fractional powers of self-adjoint elliptic operators, which are nonlocal but

self-adjoint.

In particular, the fractional Laplacian [21] describes an unusual diffusion process associated with random excursions. In general, the parabolic equations with fractional derivatives in time are associated with sub-diffusion, while the fractional elliptic operators are related to super-diffusion.

Let us consider the elliptic boundary value problem in a weak form: find $u\in V$ such that

[TABLE]

where

[TABLE]

$\Gamma=\partial\Omega$, and $\Gamma=\bar{\Gamma}_{D}\cup\bar{\Gamma}_{N}$. We assume that $\Gamma_{D}$ has positive measure, $q({x})\geq 0$ in $\Omega$, and ${\bf a}({x})$ is an SPD $d\times d$ matrix, uniformly bounded in $\Omega$, i.e.,

[TABLE]

for some positive constants $c$ and $C$. Also, $\Omega$ is a polygonal domain in $\mathbb{R}^{d}$, $d\in\{1,2,3\}$, and $f({x})$ is a given Lebesgue integrable function on $\Omega$ that belongs to the space $L_{2}(\Omega)$. Further, the case when ${\bf a}({x})$ does not depend on $x$ is referred to as problem in homogeneous media, while the general case models processes in non homogeneous media. The bilinear form $a(\cdot,\cdot)$ defines a linear operator ${\cal L}:V\to V^{*}$ with $V^{*}$ being the dual of $V$. Namely, for all $u,v\in V$ $a(u,v):=\langle{\cal L}u,v\rangle$, where $\langle\cdot,\cdot\rangle$ is the pairing between $V$ and $V^{*}$.

One possible way to introduce ${\cal L}^{\alpha}$, $0<\alpha<1$, is through its spectral decomposition, i.e.

[TABLE]

Here $\{\psi_{i}(x)\}_{i=1}^{\infty}$ are the eigenfunctions of $\cal L$, orthonormal in $L_{2}$-inner product and $\{\lambda_{i}\}_{i=1}^{\infty}$ are the corresponding positive real eigenvalues. This definition generalizes the concept of equally weighted left and right Riemann-Liouville fractional derivative, defined in one space dimension, to the multidimensional case. There is still ongoing research about the relations of the different definitions and their applications, see, e.g. [3].

The numerical solution of nonlocal problems is rather expensive. The following three approaches (A1 - A3) are based on transformation of the original problem

[TABLE]

to a local elliptic or pseudo-parabolic problem, or on integral representation of the solution, thus increasing the dimension of the original computational domain.

The Poisson problem is considered in the related papers refereed bellow, i.e.

[TABLE]

• A1

Extension to a mixed boundary value problem in the semi-infinite cylinder $C=\Omega\times\mathbb{R}_{+}\subset\mathbb{R}^{d+1}$

A “Neumann to Dirichlet” map is used in [6]. Then, the solution of fractional Laplacian problem is obtained by $u(x)=v(x,0)$ where $v:\Omega\times\mathbb{R}_{+}\rightarrow\mathbb{R}$ is a solution of the equation

[TABLE]

where $v(\cdot,y)$ satisfies the boundary conditions of (1) $\forall y\in\mathbb{R}_{+}$,

[TABLE]

as well as

[TABLE]

It is shown that the variational formulation of this equation is well posed in the related weighted Sobolev space. The finite element approximation uses the rapid decay of the solution $v(x,y)$ in the $y$ direction, thus enabling truncation of the semi-infinite cylinder to a bounded domain of modest size. The proposed multilevel method is based on the Xu-Zikatanov identity [27]. The numerical tests for $\Omega=(0,1)$ and $\Omega=(0,1)^{2}$ confirm the theoretical estimates of almost optimal computational complexity.

• A2

Transformation to a pseudo-parabolic problem

The problem (1) is considered in [25, 26] assuming the boundary condition

[TABLE]

which ensures ${\cal L}={\cal L}^{*}\geq\delta{\cal I}$, $\delta>0$. Then the solution of fractional power diffusion problem $u$ can be found

as

[TABLE]

where $w(x,t),0<t<1$, is the solution of pseudo-parabolic equation

[TABLE]

and ${\cal D}={\cal L}-\delta{\cal I}\geq 0$. Stability conditions are obtained for the fully discrete schemes under consideration. A further development of this approach is presented in [16] where the case of fractional order boundary conditions is studied.

• A3

Integral representation of the solution

The following representation of the solution of (1) is used in [4]:

[TABLE]

Among others, the authors introduce an exponentially convergent quadrature scheme. Then, the approximate solution of $u$ only involves evaluations of $({\cal I}+t_{i}{\cal A})^{-1}f$, where $t_{i}\in(0,\infty)$ is related to the current quadrature node, and where $\cal I$ and ${\cal A}$ stand for the identity and the finite element stiffness matrix corresponding to the Laplacian. The computational complexity of the method depends on the number of quadrature nodes. For instance, the presented analysis shows that approximately 50 auxiliary linear systems have to be solved to get accuracy of the quadrature scheme of order $O(10^{-5})$ for $\alpha\in\{0.25,0.5,0.75\}$.

A further development of this approach is available in [5], where the theoretical analysis is extended to the class of regularly accretive operators.

An alternative approach is applied in [11] where a class of optimal solvers for linear systems with fractional power of symmetric and positive definite (SPD) matrices is proposed. Let ${\cal A}\in\mathbb{R}^{N\times N}$ be a normalized SPD matrix generated by a finite element or finite difference approximation of some self-adjoint elliptic problem. An efficient method for solving algebraic systems of linear equations involving fractional powers of the matrix ${\cal A}$ is considered, namely for solving the system

[TABLE]

The fractional power of SPD matrix ${\cal A}$, similarly to the infinite dimensional counterpart ${\cal L}$, is expressed through the spectral representation of ${\bf u}$ through the eigenvalues and eigenvectors $\{(\Lambda_{i},{\bf\Psi}_{i})\}_{i=1}^{N}$ of ${\cal A}$, assuming that the eigenvectors are $l_{2}$-orthonormal, i.e. ${\bf\Psi}_{i}^{T}{\bf\Psi}_{j}=\delta_{ij}$ and $\Lambda_{1}\leq\Lambda_{2}\leq\dots\Lambda_{N}{\leq 1}$. Then ${\cal A}={\cal W}{\cal D}{\cal W}^{T}$, ${\cal A}^{\alpha}={\cal W}{\cal D}^{\alpha}{\cal W}^{T}$, where the $N\times N$ matrices ${\cal W}$ and ${\cal D}$ are defined as ${\cal W}=({\bf\Psi}_{1}^{T},{\bf\Psi}_{2}^{T},...,{\bf\Psi}_{N}^{T})$ and ${\cal D}=diag(\Lambda_{1},\dots,\Lambda_{N})$, ${\cal A}^{-\alpha}={\cal W}{\cal D}^{-\alpha}{\cal W}^{T}$, and the solution of ${\cal A}^{\alpha}{\bf u}={\bf f}$ can be expressed as

[TABLE]

Instead of the system (5), one can solve the equivalent system ${\cal A}^{\alpha-\beta}{\bf u}={\cal A}^{-\beta}{\bf f}{:=\mathbf{F}}$ with $\beta\geq 1$ an integer. Then the idea is to approximately evaluate ${\cal A}^{\beta-\alpha}{\mathbf{F}}$ using a set of equations involving inversion of ${\cal A}$ and ${\cal A}-d_{j}I$, for $j={1},\dots,k$. The integer parameter $k\geq 1$ is the number of partial fractions of the best uniform rational approximation (BURA) $r^{\beta}_{\alpha}(t)$ of $t^{\beta-\alpha}$ on the interval $(0,1]$. One can observe that the algorithm of [4], see A3, can be viewed as a particular rational approximation of ${\cal A}^{-\alpha}$. It is also important,

that in certain sense the results from [11] are more general, and are applicable to a wider class of sparse SPD matrices.

Assuming that ${\cal A}$ is a large-scale matrix, the computational complexity of the discussed methods for numerical solution of fractional diffusion problems is substantially high. Then, the parallel implementation of such methods for real life problems is an unavoidable topic.

In this context, there are some serious advantages of the last two approaches, see e.g., in [8].

The maximum principles are among the major qualitative properties of the elliptic or parabolic operators/PDEs. In general, to solve PDEs we use some numerical method, and it is a natural requirement that such qualitative properties are preserved on the discrete level. Most of the studies which deal with such topics give sufficient conditions for the discretization parameters in order to guarantee the certain maximum principle. It is easily see, that under certain such assumptions, the solutions of

(A1) and (A3)

satisfy certain maximum principle. In this paper, we study positive approximations of ${\cal A}^{-\alpha}$, obtained by BURA technique introduced in [11], under rather general assumptions for

the normalized SPD matrix ${\cal A}$.

The rest of the paper is organized as follows. In Section 2 we provide a brief introduction to the topic of monotone matrices including some basic properties of the M-matrices and their relations to FEM discretization of elliptic PDEs. Sufficient conditions for positive approximations of the inverse of a given SPD matrix, based on BURA technique are presented in Section 3. The analysis in Section 4 is devoted to a class of best rational approximations of ${\cal A}^{-\alpha}$, that satisfy such sufficient conditions. Sharp error estimates for a class of BURA approximations are also included in this section. Some numerical tests and short concluding remarks are given at the end.

2 Monotone matrices and SPD M-matrices

The maximum principles are some of the most useful properties used to solve a wide range of problems in the PDEs. For instance, their use is often essential to study the uniqueness and necessary conditions of solvability, approximation

and boundedness

of the solution, as well as, for quantities of physical interest like maximum stress, torsional stiffness, electrostatic capacity, charge density etc. Under certain regularity conditions, a classical maximum principle for elliptic problems reads as follows. Suppose that ${\cal L}u\geq 0$ in $\Omega$, then a nonnegative maximum is attained at the boundary $\partial\Omega$. Let us assume additionally that $u\geq 0$ in $\partial\Omega$. Then the positivity preserving property holds, that is, $u(x)>0$ for $x\in\Omega$ or $u\equiv 0$.

To solve PDEs we use some numerical methods, and it is a natural requirement that such qualitative properties are preserved on the discrete level. Most of the papers which deal with this topic give sufficient conditions for the discretization parameters in order to guarantee the certain maximum principle. For instance, when FEM is applied, the related results are usually described in terms of properties of the related mass and stiffness matrices.

Definition 2.1

A real square matrix ${\cal A}$ is called monotone if for all real vectors $v$, ${\cal A}{\bf v}\geq 0$ implies ${\bf v}\geq 0$, where $\geq$ is in element-wise sense.

The next property is sometimes used as an alternative definition.

Proposition 2.2

Let ${\cal A}$ be a real square matrix. ${\cal A}$ is monotone if and only if ${\cal A}^{-1}\geq 0$.

Definition 2.3

The class of Z-matrices are those matrices whose off-diagonal entries are less than or equal to zero. Let ${\cal A}$ be a $N\times N$ real Z-matrix, then ${\cal A}$ is a non-singular M-matrix if every real eigenvalue of ${\cal A}$ is positive. A symmetric M-matrix is sometimes called a Stieltjes matrix.

The M-matrices are among most often used monotone matrices. They arise naturally in some discretizations of elliptic operators.

Let ${\cal A}$ be an SPD matrix obtained after FEM approximation of (1) by linear triangle elements. Let us assume also that $\Omega\subset\mathbb{R}^{2}$ is discretized by a nonobtuse triangle mesh $\tau_{h}$, and the coefficients $a(x)=a_{e}$ and $q(x)=q_{e}$ are piecewise constants on the triangles $e\in\tau_{h}$. Then ${\cal A}$ can be assembled by the element matrices

[TABLE]

where ${\cal K}_{e}$ is the element stiffness matrix and ${\cal M}_{e}$ is the element mass matrix. Subject to a scaling factor, the off diagonal elements of the symmetric and positive semidefinite matrix ${\cal K}_{e}$ are equal to some of $-\cot(\theta_{i})\leq 0$, $i=1,2,3$, where $0<\theta_{i}\leq\pi/2$ are the nonobtuse angles of the triangle $e$. Imposing the boundary conditions we get that the global stiffness matrix is SPD M-matrix. The element mass matrix is positive diagonal matrix if a proper quadrature formula is applied. A similar result can be obtained by the standard diagonalization known as lumping the mass. Then, the global mass matrix is positive diagonal matrix and ${\cal A}$ is SPD M-matrix. A more general considerations of this kind are available in [15] including the case of coefficient anisotropy as well as the nonconforming linear finite elements. Similar representations of the element mass and stiffness matrices are derived if $\Omega\subset\mathbb{R}^{3}$ is discretized by a nonobtuse tetrahedral mesh $\tau_{h}$. We get again that ${\cal A}$ is again SPD M-matrix, see e.g., [14].

Remark 2.4

Not all monotone matrices are M-matrices, and the sum of two monotone matrices is not always monotone. The next examples prove these statements.

Example 2.5

• E1.1

${\cal A}_{1}=\left(\begin{matrix}-1&3\\ 2&-4\end{matrix}\right)$* is not M-matrix, but ${\cal A}_{1}^{-1}={\frac{1}{2}}\left(\begin{matrix}4&3\\ 2&1\end{matrix}\right)$ and therefore ${\cal A}_{1}$ is monotone.*

• E1.2

${\cal A}_{2}={\cal A}_{1}+6{\cal I}$* is a sum of two monotone matrices, but ${\cal A}_{2}^{-1}={\frac{1}{4}}\left(\begin{matrix}2&-3\\ -2&0\end{matrix}\right),$ and therefore ${\cal A}_{2}$ is not monotone.*

In what follows, we study positive approximations of ${\cal A}^{-\alpha}$ for a given normalized SPD M-matrix ${\cal A}$. It follows straightforwardly that the inverse of each such approximation will approximate ${\cal A}^{\alpha}$ in the class of monotone functions.

3 Positive approximations of the inverse of SPD

M-matrices

Explicit computation and memory storage of the fractional power ${\cal A}^{\alpha}$ in (5) for large-scale problems is expensive and impractical. Even when ${\cal A}$ is sparse, ${\cal A}^{\alpha}$ is typically dense. Therefore, we study possible positive approximations of the action of ${\cal A}^{-\alpha}$ based solely on the information of ${\cal A}$. For this purpose, we consider the class of rational functions

[TABLE]

fix a positive integer $\beta$, and search for an appropriate candidate $r$ in it, that approximates well the univariate function $t^{\beta-\alpha}$ on the unit interval $[0,1]$. Note that, due to the normalization of ${\cal A}$, this interval covers the spectrum of ${\cal A}^{\alpha}$,

Definition 3.1

Let $\alpha\in(0,1)$, and $\beta,m,k\in\mathbb{N}\setminus\{0\}$. The minimizer $r^{\beta}_{\alpha}\in\mathcal{R}(m,k)$ of the problem

[TABLE]

will be called $\beta$-Best Uniform Rational Approximation ($\beta$-BURA). Its error will be denoted by

[TABLE]

Based on classical Spectral Theory arguments (see [11, Theorem 2.1]), the univariate approximation error $E_{\alpha}(m,k;\beta)$ is an upper bound for the multivariate relative error $\|{\cal A}^{-\beta}r^{\beta}_{\alpha}({\cal A}){\bf f}-{\cal A}^{-\alpha}{\bf f}\|_{{\cal A}^{\gamma+\beta}}/\|{\bf f}\|_{{\cal A}^{\gamma-\beta}}$ for the corresponding matrix-valued BURA approximation of the exact solution ${\bf u}$ in (6). Here, $\gamma\in\mathbb{R}$ can be arbitrary,

and the Krylov norms are defined via standard energy dot product, i.e. $\|{\bf f}\|^{2}_{{\cal A}^{\gamma-\beta}}=\langle{\cal A}^{\gamma-\beta}{\bf f},{\bf f}\rangle$.

Proposition 3.2

Let ${\cal A}\in\mathbb{R}^{N\times N}$ be an SPD matrix with eigenvalues $0<\Lambda_{1}\leq\Lambda_{2}\leq\dots\leq\Lambda_{N}\leq 1$. Let $r^{\beta}_{\alpha}$ be the $\beta$-BURA for given $\alpha,\beta,m,k$. Then,

[TABLE]

For the practical computation of ${\cal A}^{-\beta}r^{\beta}_{\alpha}({\cal A}){\bf f}$ we use the partial fraction decomposition of $t^{-\beta}r^{\beta}_{\alpha}(t)$, which is of the form

[TABLE]

provided all $r^{\beta}_{\alpha}$ has no complex poles and the real ones $\{d_{j}\}_{1}^{k}$ are all of multiplicity 1. Later, we will see that for $\beta=1$ and $m=k$ the above assumption on the poles of $r^{\beta}_{\alpha}$ holds true for any $\alpha\in(0,1)$. Furthermore, in all our numerical experiments with various $\beta,m,k$ and $\alpha\in\{0.25,0.5,0.75\}$ the assumption always remains valid. Hence, it does not seem to restrict the application range of the proposed method. On the other hand, under (9) the approximate solution

[TABLE]

of ${\bf u}$ can be efficiently numerically computed via solving several linear systems, that involve ${\cal A}$ and its diagonal variations ${\cal A}-d_{j}I$, for $j=1,\dots,k$.

Definition 3.3

A real symmetric matrix ${\cal A}^{-1}$ is said to be doubly nonnegative if it is both positive definite, and entrywise nonnegative.

Our first goal is to analyze under what conditions on the coefficients and the poles in (9), the matrix ${\cal A}^{-\beta}r^{\beta}_{\alpha}({\cal A})$ remains doubly nonnegative. Clearly the matrix is symmetric whenever ${\cal A}$ is, so the main investigations are on assuring ${\cal A}^{-\beta}r^{\beta}_{\alpha}({\cal A})\geq 0$. The following proposition contains sufficient conditions for positivity.

Proposition 3.4

If ${\cal A}$ is a normalized SPD M-matrix, $m<k+\beta$, $c_{0}\geq 0$, $c\geq 0$, and ${\mathbf{d}}<0$ (entrywise), then ${\cal A}^{-\beta}r^{\beta}_{\alpha}({\cal A})$ in (10) is doubly nonnegative.

Proof: Since $m<k+\beta$, equation (10) is simplified to

[TABLE]

For every $j=1,\dots,k$, the matrix ${\cal A}-d_{j}I$ is an SPD M-matrix, as $d_{j}<0$ and the diagonal elements increase their values, i.e. become stronger dominant. Hence, $({\cal A}-d_{j}I)^{-1}\geq 0$.We have ${\cal A}^{-1}\geq 0$, thus ${\cal A}^{-j}=({\cal A}^{-1})^{j}\geq 0$, $j=1,\dots,\beta$, as each entry of ${\cal A}^{-j}$ is a sum of nonnegative summands. Finally, a linear combination of positively scaled doubly nonnegative matrices is also a doubly nonnegative matrix. $\Box$

Note that, when applying pure polynomial approximation techniques for $t^{-\alpha}$ on $[\Lambda_{1},1]$ like in [12], there is practically no chance to come up with a positive approximation of ${\cal A}^{-\alpha}$. First of all, such an approximant is a linear combination of positive degrees of ${\cal A}$ and in particular ${\cal A}$ itself appears with a nonzero coefficient. This matrix has non-positive off-diagonal entries. Furthermore, it was numerically observed that the coefficient sequence in the linear combination is sign alternating. Hence, the proposed $\beta$-BURA approach seems the right and most natural tool for constructing positive approximations of ${\cal A}^{-\alpha}$, or alternatively, monotone approximations of ${\cal A}^{\alpha}$. Another disadvantage of the former approach is the restriction on $\Lambda_{1}$ to be well-separated from zero, which is also a restriction on the condition number of ${\cal A}$.

4 Analysis of a class of best rational approximations of fractional power

of SPD M-matrices

Among all various classes of best rational approximations, the diagonal sequences $r\in{\cal R}(k,k)$ of the Walsh table of $t^{\alpha}$, $\alpha\in(0,1)$ are studied in greatest detail [20, 9, 24, 22]. There is an existence and uniqueness of the BURA elements for all $k$ and $\alpha$. The distribution of poles, zeros, and extreme points of those elements plays a central role in asymptotic convergence analysis, when $k\to\infty$, thus is well known. In this section, we show that the above diagonal class perfectly fits within our positive ${\cal A}^{-\alpha}$ approximation framework.

First, we collect some preliminary results that will be later needed for the proof of the main theorem. The following characterization lemma, which we state here without proof, is vital for our further investigations.

Lemma 4.1

[22, Lemma 2.1]** Let $m=k$ and $0<\alpha<1$.

• (a)

The best rational approximant $r^{1}_{\alpha}$ is of exact numerator and denominator degree $k$.

• (b)

All $k$ zeros $\zeta_{1},\dots,\zeta_{k}$ and poles $d_{1},\dots,d_{k}$ of $r^{1}_{\alpha}$ lie on the negative half-axis $\mathbb{R}_{<0}$ and are interlacing; i.e., with an appropriate numbering we have

[TABLE]

• (c)

The error function $t^{1-\alpha}-r^{1}_{\alpha}(t)$ has exactly $2k+2$ extreme points $\eta_{1},\dots,\eta_{2k+2}$ on $[0,1]$, and with an appropriate numbering we have

[TABLE]

The next lemma builds a bridge between the fractional decompositions of $r^{1}_{\alpha}$ and $t^{-1}r^{1}_{\alpha}$.

Lemma 4.2

Let $m=k$, $0<\alpha<1$, and

[TABLE]

Then

[TABLE]

Proof: The second part of (14) follows directly from

[TABLE]

For the first part, we combine the above identity with (12) and (13)

[TABLE]

The proof of the lemma is completed. $\Box$

Our last lemma provides an asymptotic bound on $E_{\alpha}(k,k;1)$. The proof can be found in [24].

Lemma 4.3

[24, Theorem 1]** The limit

[TABLE]

holds true for each $\alpha>0$.

Now, we are ready to formulate and prove our main result.

Theorem 4.4

Let $\beta=1$ and $m=k$. For every normalized SPD M-matrix ${\cal A}$ and every $\alpha\in(0,1)$, the matrix ${\cal A}^{-1}r^{1}_{\alpha}({\cal A})$ is doubly nonnegative and for all $\gamma\in\mathbb{R}$

[TABLE]

Proof:

Based on the results in Lemma 4.1, we can quickly derive ${\cal A}^{-1}r^{1}_{\alpha}({\cal A})\geq 0$. For this purpose, we study the sign pattern of $c$ and $d$ and assure the applicability of Proposition 3.4. From (11) we know that all the poles $\{d_{j}\}$ are real, negative, and of multiplicity 1. The same holds true for the zeros $\{\zeta_{j}\}$. Since $r^{1}_{\alpha}(t)$ is continuous on $\mathbb{R}\setminus\{d_{j}\}$, the function changes its sign $2k$ times - at each zero $\zeta_{j}$ and at each pole $d_{j}$. In Lemma 4.2, we have already computed that $r^{1}_{\alpha}(0)=c_{0,1}=E_{\alpha}(k,k;1)>0$, thus, due to interlacing, at each pole $d_{j}$ we have

[TABLE]

Since $c^{\ast}<0$ and $d<0$, from Lemma 4.2 it follows that $c_{0}>0$ and $c>0$. Hence, Proposition 3.4 gives rise to ${\cal A}^{-1}r^{1}_{\alpha}({\cal A})\geq 0$.

The error estimate (15) is a direct corollary of Proposition 3.2 and Lemma 4.3. $\Box$

Some remarks are in order. For any fixed $\alpha<1$, the relative error (15) decays exponentially as $k\to+\infty$ with order $\sqrt{(1-\alpha)k}$. When $\alpha\to 1$ the relative error decays linearly independently of $k$, since $|\sin{\pi(1-\alpha)}|\to\pi(1-\alpha)$. It is straightforward to extend the coefficient correspondence (14) to $c_{j}=c^{\ast}_{j}/d^{\beta}_{j}$ for any $(m,k,\beta)$, such that $m<k+\beta$. Therefore, a necessary condition for Proposition 3.4 to be applicable is the sequence $c^{\ast}$ to have constant sign. Due to the proof of Theorem 4.4, it implies that zeros $\{\zeta_{i}\}_{1}^{m}$ and poles $\{d_{j}\}_{1}^{k}$ of $r^{\beta}_{\alpha}$ should be interlacing, thus $|m-k|\leq 1$. Furthermore (see [11, (20)]) the following identity always holds true

[TABLE]

Hence, another necessary condition for applicability of Proposition 3.4 is $m\geq k+\beta-1$. Combining all derived constraints, we observe that ${\cal A}^{-\beta}r^{\beta}_{\alpha}({\cal A})$ could be represented as a sum of doubly nonnegative matrices only if $m=k+\beta-1$ and $m\leq k+1$, meaning that we are left with the admissible triples

[TABLE]

In [22] it is remarked that for the case $(k,k,2)$ it cannot be theoretically excluded that one root $\zeta_{i}$ and one pole $d_{i}$ of the 2-BURA $r^{2}_{\alpha}(t)$ lie outside of $\mathbb{R}_{<0}$. For the case $(k+1,k,2)$ we can show that $c_{0,2}=-E_{\alpha}(k+1,k;\beta)<0$, thus the assumptions of the proposition are again violated. In conclusion, the triple $(k,k,1)$, investigated in Theorem 4.4 is the unique choice of parameters for which one can prove positiveness of the approximation ${\cal A}^{-\beta}r^{\beta}_{\alpha}({\cal A})$ using Proposition 3.4.

5 Numerical tests

The main goal of the numerical tests is to illustrate the positive properties of the proposed approximations of inverse of fractional powers of SPD M-matrices. Complementary, we provide a short discussion related to some interpretations of the results from Section 4 in the case of fractional diffusion problems.

The test problems are in 1D. In this case we are able easier to compute the exact solutions of linear systems with fractional powers of the corresponding tridiagonal matrix. Note that this does not cause restrictions to the derived conclussions. As it was shown in [11], if $\Omega\subset\mathbb{R}^{d},d>1$, some PCG solver of optimal complexity (e.g., BoomerAMG) can be utilized for efficient solution of the arising sparse linear systems, fully preserving the accuracy and efficiency of the composite algorithm.

We consider fractional powers of the Poisson’s equation on the unit interval with Dirichlet boundary conditions:

[TABLE]

On a uniform grid with mesh parameter $h=1/(N+1)$, using central finite differences, the operator $\mathcal{L}$ is approximated by the $N\times N$ matrix ${\cal A}_{h}:=tridiag(-1,2,-1)/h^{2}$, which in turn can be rewritten as

[TABLE]

The matrix ${\cal A}$ is a normalized, SPD M-matrix, which eigenvectors and eigenvalues are explicitly known:

[TABLE]

We approximate $\mathcal{L}^{\alpha}$ by ${\cal A}_{h}^{\alpha}$ and, due to Theorem 4.4, the $\ell^{2}$ relative error is bounded by an $h$-dependent constant.

Corollary 5.1

Let ${\bf u}_{h}:={\cal A}_{h}^{-\alpha}{\bf f}=\left(\frac{h}{2}\right)^{2\alpha}{\cal A}^{-\alpha}{\bf f}$ be the exact solution of the discretized fractional Poisson’s equation ${\cal A}_{h}^{\alpha}{\bf u}={\bf f}$. Let $r^{1}_{\alpha}\in\mathcal{R}(k,k)$ be 1-BURA and denote by ${\bf u}_{h,r}:=\left(\frac{h}{2}\right)^{2\alpha}{\cal A}^{-1}r^{1}_{\alpha}({\cal A}){\bf f}$. Then

[TABLE]

Proof: Indeed, let ${\bf u}={\cal A}^{-\alpha}{\bf f}$ and ${\bf u}_{r}={\cal A}^{-1}r^{1}_{\alpha}({\cal A}){\bf f}$. Applying

[TABLE]

where $\mathrm{k}({\cal A})$ is the condition number of ${\cal A}$, we derive

[TABLE]

The result follows from (15) for $\gamma=1$. $\Box$

Due to Corollary 5.1, we can compute the minimal degree $k$ that guarantees $\|{\bf u}_{h,r}-{\bf u}_{h}\|_{2}/\|{\bf f}\|_{2}<\varepsilon$ for every given pair $(\varepsilon,h)$. Such an $\ell^{2}$ error analysis is outside of the scope of this paper, so we will not further elaborate on it.

In our numerical experiments, we choose $\alpha\in\{0.25,0.5,0.75\}$ and $k\in\{5,6,7\}$. The square root of an M-matrix is again an M-matrix [1] and it is easy to check that ${\cal A}_{h}^{3}$ is also an M-matrix. Therefore, all considered ${\cal A}_{h}^{\alpha}$ are M-matrices, their inverse matrices are doubly nonnegative (but dense!), and constructing computationally cheep approximants within the same class is of great practical importance. The univariate error estimates $E_{\alpha}(k,k;1)$ for the above choice of parameters are summarized in Table 1. Note that each of them satisfies the inequality (15) even without introducing the low-order term $\mathrm{o}(1)$ in the right-hand-side.

A modified Remez algorithm is used for the derivation of $r^{1}_{\alpha}$ [17, 7].

For $f$ in (4) we take two different positive functions, supported on the interval $[1/2,3/4]$. The first one $f_{1}$ is piecewise constant and discontinuous, while the second one $f_{2}$ is a $C^{2}$ cubic spline function, corresponding to the Irwin-Hall distribution. Together with the exact discretized solutions ${\cal A}_{h}^{-\alpha}{\bf f}$, they are illustrated on Fig. 1.

We consider mesh parameters $h=2^{-m}$, $m\in\{5,6,\dots,11\}$. On Fig. 2 the corresponding approximants ${\bf u}_{h,r}$, $h=2^{-11}$, of ${\bf u}_{h}={\cal A}_{h}^{-\alpha}{\bf f}_{1}$ are plotted. As suggested by Corollary 5.1, ${\bf u}_{h,r}$ fails

to approximate well ${\bf u}_{h}$ on a fine grid for smaller $k$

and $\alpha$ (see $k=5$, $\alpha\in\{0.25,0.5\}$). For

larger $\alpha$,

the exponential growth of the $\ell^{2}$ relative error with $h\to 0$ is less significant, as it is of order $2(1-\alpha)$, thus when $\alpha=0.75$ all the three approximants, corresponding to $k=\{5,6,7\}$ follow closely the graph of ${\cal A}_{h}^{-3/4}{\bf f}_{1}$. On Fig. 3 and in Table 2 we numerically confirm the asymptotic behavior of the relative $\ell^{2}$ error from Corollary 5.1.

The sufficient conditions from Proposition 3.4 hold true for the cases under consideration. This means that positivity of all considered approximations is guaranteed. Therefore, the discrete maximum principle is always inherited. The presented numerical results are fully aligned with the theory. What is very important is the numerical robustness of positivity with respect to both accuracy parameters $h$ and $k$ which is confirmed for all $\alpha\in\{0.25,0.5,0.75\}$. Even in the case of lower accuracy, we do not observe any oscillations. The monotonicity preservation of the data is clearly expressed, capturing their geometrical shape.

6 Concluding remarks

This study is inspired by some quite recent results in the numerical methods for fractional diffusion problems. In the Introduction, we discussed three methods based on reformulation of the original nonlocal problem into local (elliptic, pseudo parabolic, and integral) problems. In all cases, the cost is in the increased dimension of computational domain from $d$ to $d+1$.

Our approach is based on best uniform rational approximations of $t^{\beta-\alpha}$, $0\leq t\leq 1$. The primal motivation is to reduce the computational complexity. A next important step is made in this paper. Here, we provide sufficient conditions to guarantee positive approximation of the inverse of fractional powers of normalized SPD M-matrices. Therefore, we get a numerical method which preserves the maximum principle. The presented numerical results clearly confirm the monotone behaviour of the solution, without any observed oscillations. Further research has to be devoted to the topic of accuracy of mass conservation.

The currently available methods and algorithms for numerical solution of boundary value problems with fractional power of elliptic operators have a quite different nature. A serious theoretical and experimental study is required to get a comparative analysis of their advantages and disadvantages for particular classes of problems. For instance, the error analysis is in different functional spaces assuming different conditions for smoothness. The comparison of the computational complexity is also an open question.

As a part of our analysis, Theorem 1 provides a sharp error estimates for the 1-BURA based approximations $E(k,k,1)$. Then, at the beginning of Section 5, we showed how this result can be used to derive relative error estimates of the numerical solution of fractional order elliptic problems in $\ell^{2}$. The numerical tests are well aligned with this theoretical estimates. The presented approach has a strong potential for further development addressing different pairs of functional spaces in the relative error estimates, varying the smoothness assumptions, for $d=1,2,3$.

In addition, a lot of new numerical tests are needed to evaluate/confirm/compare the computational efficiency for more realistic towards real-life large-scale super diffusion problems.