# Fine spectral estimates with applications to the optimally fast solution of
large FDE linear systems

M. Bogoya  University of Insubria, Como, Italy
Dipartimento di Scienza e Alta Tecnologia
<EMAIL_ADDRESS>S.M. Grudsky  CINVESTAV del I.P.N., Mexico
D.F., Mexico
Departamento de Matemáticas
Southern Federal University, Rostov-on-Don, Russia
Regional Mathematical Center
<EMAIL_ADDRESS>S. Serra–Capizzano  University of Insubria, Como,
Italy
Dipartimento di Scienza e Alta Tecnologia
<EMAIL_ADDRESS>C. Tablino–Possio  University of Milano-
Bicocca, Milano, Italy
Dipartamento di Matematica e Applicazioni
<EMAIL_ADDRESS>
###### Abstract

In the present note we consider a type of matrices stemming in the context of
the numerical approximation of distributed order fractional differential
equations (FDEs): from one side they could look standard, since they are,
real, symmetric and positive definite. On the other hand they present specific
difficulties which prevent the successful use of classical tools. In
particular the associated matrix-sequence, with respect to the matrix-size, is
ill-conditioned and it is such that a generating function does not exists, but
we face the problem of dealing with a sequence of generating functions with an
intricate expression. Nevertheless, we obtain a real interval where the
smallest eigenvalue belongs, showing also its asymptotic behavior. We observe
that the new bounds improve those already present in the literature and give a
more accurate spectral information, which are in fact used in the design of
fast numerical algorithms for the associated large linear systems,
approximating the given distributed order FDEs. Very satisfactory numerical
results are presented and critically discussed, while a section with
conclusions and open problems ends the current note.
Keywords: Toeplitz sequences, algebra of matrix-sequences, generating
function, fractional operators
Subclass: MSC 15B05, 15A18, 26A33

##  1 Introduction

When considering the numerical approximation of fractional differential
equations (FDEs), due to the nonlocal nature of the involved operators, the
matrices are intrinsically dense, even in the case where the approximation
technique is of local nature (Finite Differences, Finite Elements, Finite
Volumes, Isogeometric Analysis etc.); see [14, 18, 19, 20, 22, 25, 34, 35, 36,
37] and references therein. The same situation is present also in the context
of distributed order FDEs (see [21] and references therein). When employing
equispaced gridding, the resulting structures have the advantage of being of
Toeplitz nature, so that the cost of a matrix-vector multiplication is almost
linear i.e. of $O(n\log n)$ arithmetic operations, where $n$ is the matrix
order and the constant hidden in the big $O$ is very moderate (refer to [8,
23] and to the references there reported). These computational features, joint
with the usually large dimensions of the considered linear systems, lead to
the search for suitable iterative solvers: typically the most successful
iterative procedures belong to the class of preconditioned Krylov methods, to
the algorithms of multigrid type, or to clever combinations of them [8, 10,
11, 12, 23, 31, 29].

A crucial information for an appropriate design of an efficient solver of such
a kind is the precise knowledge of the asymptotical behavior of the minimal
eigenvalue, which in the positive definite case is also related to the
asymptotic spectral conditioning: in our setting we emphasize that the
considered matrices are real, symmetric, and positive definite.

In this work, starting from [5, 21], we improve the bounds present in the
literature. We take into account the techniques already developed in [4, 6]:
the additional difficulty, in the present setting, relies on the fact that the
generating function is not fixed, since it depends on the matrix order $n$,
which is somehow a challenging novelty with respect to the classical work on
the subject [4, 6, 27].

In reality, in technical terms, we consider the real-valued symbol
$f_{n}(\theta)\coloneqq\sum_{j=0}^{n-1}h^{jh}|\theta|^{2-jh}$ with
$h\coloneqq\frac{1}{n}$. In this note we obtain a real interval where the
smallest eigenvalue of the Toeplitz matrix $A_{n}\coloneqq T_{n}(f_{n})$
belongs to, with $f_{n}$ being the generating function of $A_{n}$, for any
fixed $n$.

Based on the spectral information, few algorithmic proposals are also
discussed starting from those presented in [12, 13, 21] and are of the type
mentioned before: preconditioned conjugate gradient (PCG) algorithms,
multigrid solvers (typically the V-cycle), and combinations of them that is a
V-cycle where the smoothing iterations (one step) incorporate the proposed PCG
choices.

The work is organized as follows. Section 2 contains the setting of the
problem regarding the minimal eigenvalue of $A_{n}$ and its study and
solution. Section 3 is devoted to numerical experiments concerning the
solution of large linear systems with coefficient matrix $A_{n}$, including
few evidences of the spectral clustering at one of the proposed preconditioned
matrix-sequences.

##  2 The problem and its solution

Let $h\coloneqq\frac{1}{n}$ and consider the function
$f_{n}\colon[-\pi,\pi]\to\mathbb{R}$ given by

$f_{n}(\theta)\coloneqq\theta^{2}\sum_{j=0}^{n-1}h^{jh}|\theta|^{-jh}=\theta^{2}\frac{1-\frac{1}{n|\theta|}}{1-\big{(}\frac{1}{n|\theta|}\big{)}^{\frac{1}{n}}}.$

Note that $f_{n}$ coincides with the function $\hat{F}_{n}$ used in [5]. We
employ the notation $A_{n}\coloneqq T_{n}(f_{n})$, then $f_{n}$ is the so-
called generating function or symbol of $A_{n}$. Using that,

$\displaystyle 1-\Big{(}\frac{1}{n|\theta|}\Big{)}^{\frac{1}{n}}$
$\displaystyle=$ $\displaystyle
1-\exp\Big{\\{}-\frac{\log(n|\theta|)}{n}\Big{\\}}$ $\displaystyle=$
$\displaystyle\frac{\log(n|\theta|)}{n}\Big{\\{}1+O\Big{(}\frac{|\log(|n|\theta|)|}{n}\Big{)}\Big{\\}}\quad\mbox{as}\quad
n\to\infty,$

we obtain

$\displaystyle f_{n}(\theta)$ $\displaystyle=$
$\displaystyle\frac{n\theta^{2}}{\log(n|\theta|)}\Big{\\{}1-\frac{1}{n|\theta|}\Big{\\}}\Big{\\{}1+O\Big{(}\frac{|\log(n|\theta|)|}{n}\Big{)}\Big{\\}}$
$\displaystyle=$
$\displaystyle\frac{|\theta|(n|\theta|-1)}{|\log(n|\theta|)|}\Big{\\{}1+O\Big{(}\frac{|\log(n|\theta|)|}{n}\Big{)}\Big{\\}}$
$\displaystyle=$
$\displaystyle\frac{|\theta|(n|\theta|-1)}{\log(n|\theta|)}+O\Big{(}\theta^{2}+\frac{|\theta|}{n}\Big{)}.$

To simplify the previous expression, we consider now the function $g$ given by

$g(\sigma)\coloneqq\frac{\sigma^{2}-\sigma}{\log(\sigma)}.$ (2.1)

Then we infer the relation $nf_{n}(\theta)=g(n|\theta|)+r_{n}(\theta)$ where
$r_{n}(\theta)=O(n\theta^{2}+|\theta|)$ as $n\to\infty$.

We now establish a connection with the simple-loop method [3, 4]. Consider the
Toeplitz matrix $T_{n}(a)$ with symbol $a\colon\mathbb{T}\to\mathbb{R}$ given
by

$a(t)=-\frac{1}{t}+2-t=4\sin^{2}\Big{(}\frac{\theta}{2}\Big{)},$

where $t=\mathrm{e}^{\mathrm{i}\theta}$. For a symbol $b$, let
$\lambda_{j}(T_{n}(b))$ and $\psi_{j}(T_{n}(b))$ be its $j$th eigenvalue and
eigenfunction, respectively. Since $a$ is real-valued, its eigenvalues must be
real and we arrange them in increasing order, that is

$\lambda_{1}(T_{n}(a))\leqslant\lambda_{2}(T_{n}(a))\leqslant\cdots\leqslant\lambda_{n}(T_{n}(a)).$

If $\psi_{j}(T_{n}(a))$ is the eigenfunction corresponding to
$\lambda_{j}(T_{n}(a))$, then it is well-known that (see for example [4, 7])
it leads to the equations

$\displaystyle\lambda_{1}(T_{n}(a))$ $\displaystyle=$ $\displaystyle
4\sin^{2}\Big{(}\frac{s}{2}\Big{)},$ $\displaystyle\psi_{1}(T_{n}(a))(\theta)$
$\displaystyle=$
$\displaystyle\frac{c_{n}\,\mathrm{e}^{\frac{\mathrm{i}}{2}(n+1)\theta}}{(n+1)^{\frac{3}{2}}}\cdot\frac{\cos\big{(}\frac{(n+1)\theta}{2}\big{)}}{\sin\big{(}\frac{\theta-s}{2}\big{)}\sin\big{(}\frac{\theta+s}{2}\big{)}},$
(2.2)

where $s\coloneqq\frac{\pi}{n+1}$ and $c_{n}$ is bounded when $n\to\infty$.
The constant $c_{n}$ can be calculated from the relation
$\|\psi_{1}(T_{n}(a))\|_{2}=1$, see Figure 1.

Figure 1: The constant $c_{n}$ in Equation (2) for different values of $n$.

###  2.1 Upper bound for $\mathbf{\lambda_{1}(A_{n})}$

######  Theorem 2.1.

We have

$n\lambda_{1}(A_{n})\leqslant
k_{1}+O\Big{(}\frac{1}{n}\Big{)},\quad\mbox{as}\quad n\to\infty,$

where the constant $k_{1}$ is given by

$k_{1}=\Big{\\{}\int_{0}^{\infty}\frac{u^{2}-u}{(u^{2}-\pi^{2})^{2}}\cdot\frac{\cos^{2}(\frac{u}{2})}{\log(u)}\operatorname{d\\!}u\Big{\\}}\Big{\\{}\int_{0}^{\pi}\frac{\cos^{2}(\frac{u}{2})}{(u^{2}-\pi^{2})^{2}}\operatorname{d\\!}u\Big{\\}}^{-1},$

and can be numerically approximated as $12.9301$.

###### Proof.

Let $\langle\cdot,\cdot\rangle$ be the inner product of the Hilbert space
$L^{2}(\mathbb{T})$. Note that the essential ranges of the symbols $n^{2}a$
and $nf_{n}$ have a similar behavior in an interval of the type
$\big{[}0,O\big{(}\frac{1}{n}\big{)}\big{]}$. Hence using the well-known
formula

$\lambda_{1}(A_{n})=\inf\\{\langle
A_{n}\psi,\psi\rangle\colon\|\psi\|_{2}=1\\},$

we obtain

$\displaystyle\lambda_{1}(A_{n})$ $\displaystyle\leqslant$
$\displaystyle\langle A_{n}\psi_{1}(A_{n}),\psi_{1}(A_{n})\rangle$ (2.3)
$\displaystyle=$
$\displaystyle\frac{1}{2\pi}\int_{-\pi}^{\pi}f_{n}(\theta)|\psi_{1}(A_{n})(\theta)|^{2}\operatorname{d\\!}\theta$
$\displaystyle=$ $\displaystyle\mathcal{I}_{1}(n)+\mathcal{I}_{2}(n),$

where

$\displaystyle\mathcal{I}_{1}(n)$ $\displaystyle\coloneqq$
$\displaystyle\frac{1}{2\pi
n}\int_{-\pi}^{\pi}g(n|\theta|)|\psi_{1}(A_{n})(\theta)|^{2}\operatorname{d\\!}\theta,$
$\displaystyle\mathcal{I}_{2}(n)$ $\displaystyle\coloneqq$
$\displaystyle\frac{1}{2\pi
n}\int_{-\pi}^{\pi}r_{n}(\theta)|\psi_{1}(A_{n})(\theta)|^{2}\operatorname{d\\!}\theta.$

It is easy to see that $\lim_{n\to\infty}n\,\mathcal{I}_{1}(n)<\infty$, indeed

$\displaystyle n\,\mathcal{I}_{1}(n)$ $\displaystyle=$
$\displaystyle\frac{c_{n}^{2}}{2\pi(n+1)^{3}}\int_{-\pi}^{\pi}g(n|\theta|)\cdot\frac{\cos^{2}\big{(}\frac{(n+1)\theta}{2}\big{)}}{\sin^{2}\big{(}\frac{\theta-s}{2}\big{)}\sin^{2}\big{(}\frac{\theta+s}{2}\big{)}}\operatorname{d\\!}\theta$
(2.4) $\displaystyle\sim$
$\displaystyle\frac{16nc_{n}^{2}}{\pi}\int_{0}^{\pi}g(n\theta)\cdot\frac{\cos^{2}\big{(}\frac{n\theta}{2}\big{)}}{((n\theta)^{2}-\pi^{2})^{2}}\operatorname{d\\!}\theta$
$\displaystyle\sim$
$\displaystyle\frac{16c_{n}^{2}}{\pi}\int_{0}^{\infty}\frac{u^{2}-u}{(u^{2}-\pi^{2})^{2}}\cdot\frac{\cos^{2}\big{(}\frac{u}{2}\big{)}}{\log(u)}\operatorname{d\\!}u.$

Here the notation $f\sim g$ means $\lim_{n\to\infty}\frac{f(n)}{g(n)}=1$.
Using (2) a similar calculation produces

$\displaystyle c_{n}^{-2}$ $\displaystyle\sim$
$\displaystyle\frac{16n}{\pi}\int_{0}^{\pi}\frac{\cos^{2}\big{(}\frac{n\theta}{2}\big{)}}{((n\theta)^{2}-\pi^{2})^{2}}\operatorname{d\\!}\theta$
(2.5) $\displaystyle=$
$\displaystyle\frac{16}{\pi}\int_{0}^{\infty}\frac{\cos^{2}\big{(}\frac{u}{2}\big{)}}{(u^{2}-\pi^{2})^{2}}\operatorname{d\\!}u.$

Note that the last integral does not depend on $n$, its value can be
numerically found producing $c_{n}\approx 2.2214$, which agrees with the
Figure 1. For the second integral in (2.3) we recall that
$r_{n}(\theta)=O(n\theta^{2}+|\theta|)$ and write
$\mathcal{I}_{2}(n)=O(\Phi(n))$ where

$\displaystyle\Phi(n)$ $\displaystyle=$
$\displaystyle\frac{1}{n^{2}}\int_{-\pi}^{\pi}\\{(n|\theta|)^{2}-n|\theta|\\}|\psi_{1}(A_{n})(\theta)|^{2}\operatorname{d\\!}\theta$
$\displaystyle=$
$\displaystyle\frac{c_{n}^{2}}{n^{5}}\int_{-\pi}^{\pi}\\{(n|\theta|)^{2}-n|\theta|\\}\frac{\cos^{2}\big{(}\frac{n|\theta|}{2}\big{)}}{((n|\theta|)^{2}-\pi^{2})^{2}}\operatorname{d\\!}\theta$
$\displaystyle\sim$
$\displaystyle\frac{16c_{n}^{2}}{n^{2}}\int_{0}^{\infty}\frac{u^{2}-u}{(u^{2}-\pi^{2})^{2}}\cos\Big{(}\frac{u}{2}\Big{)}\operatorname{d\\!}u,$

which combined with (2.5) produces

$\mathcal{I}_{2}(n)=O\Big{(}\frac{1}{n^{2}}\Big{)}.$ (2.6)

Finally, combining (2.3), (2.4), (2.5), and (2.6) we obtain the thesis. ∎

###  2.2 Lower bound for $\mathbf{\lambda_{1}(A_{n})}$

In this part we will implement the trick used in [6] which works as follows.
Let $b\in C(\mathbb{T})$ and $q_{n}\colon\mathbb{T}\to\mathbb{C}$ given by

$q_{n}(t)\coloneqq\sum_{j=n}^{\infty}q_{j,n}t^{j}+\sum_{j=n}^{\infty}q_{-j,n}t^{-j}.$
(2.7)

Since $T_{n}(q_{n})$ is the zero matrix we clearly have
$T_{n}(a)=T_{n}(a+q_{n})$, thus instead of working with $T_{n}(a)$ we can use
$T_{n}(a+q_{n})$ which under a swiftly selected symbol $q_{n}$, can be
advantageous.

######  Theorem 2.2.

We have

$n\lambda_{1}(A_{n})\geqslant
k_{2}+O\Big{(}\frac{1}{n}\Big{)},\quad\mbox{as}\quad n\to\infty,$

where the constant $k_{2}$ is given by

$k_{2}=\frac{1}{\pi}\int_{0}^{\pi}\frac{\sigma^{2}-\sigma}{\log(\sigma)}\operatorname{d\\!}\sigma,$

which can be approximated numerically as $2.2945$.

###### Proof.

In our case instead of working with the symbol $nf_{n}$ we will use

$\displaystyle g_{n}(\theta)$ $\displaystyle\coloneqq$ $\displaystyle
nf_{n}(\theta)+p(n|\theta|),$ $\displaystyle=$ $\displaystyle
g(n|\theta|)+p(n|\theta|)+r_{n}(\theta),$

where $p(\sigma)\coloneqq\sum_{j=1}^{\infty}p_{j}\cos(\sigma j)$, $p_{j}$ are
real constants, and $r_{n}$ is given in (2.1). Let
$t=\mathrm{e}^{\mathrm{i}\theta}$, then we can write

$p(n|\theta|)=p(n\theta)=\sum_{j=1}^{\infty}\frac{p_{j}}{2}(t^{nj}+t^{-nj}),$

which clearly shows the form (2.7). Additionally, the function $p$ is
$2\pi$-periodic, even, and satisfies
$\int_{0}^{\pi}p(\sigma)\operatorname{d\\!}\sigma=0$. We thus deduce
$T_{n}(p(n|\theta|))=0$, and hence

$\displaystyle T_{n}(nf_{n})$ $\displaystyle=$ $\displaystyle T_{n}(g_{n})$
$\displaystyle=$ $\displaystyle
T_{n}(g(n|\theta|))+T_{n}(p(n|\theta|))+T_{n}(r_{n})$ $\displaystyle=$
$\displaystyle T_{n}(g(n|\theta|))+T_{n}(r_{n}).$

Keeping in mind that, for any real-valued symbol $\beta$, the smallest
eigenvalue of $T_{n}(\beta)$ must be greater or equal to the infimum of
$\beta$, we obtain

$\displaystyle n\lambda_{1}(A_{n})$ $\displaystyle\geqslant$
$\displaystyle\inf\\{g_{n}(\theta)\colon\theta\in[-\pi,\pi]\\}$
$\displaystyle=$
$\displaystyle\inf\\{g(n|\theta|)+p(n|\theta|)+r_{n}(\theta)\colon\theta\in[-\pi,\pi]\\}$
$\displaystyle=$
$\displaystyle\inf\Big{\\{}g(\sigma)+p(\sigma)+O\Big{(}\frac{\sigma^{2}+\sigma}{n}\Big{)}\colon\sigma\in[0,n\pi]\Big{\\}},$

where in the last line we used the variable change $\sigma=n|\theta|$. Then we
obtain

$n\lambda_{1}(A_{n})\geqslant\inf\Big{\\{}g(\sigma)+p(\sigma)+O\Big{(}\frac{\sigma^{2}+\sigma}{n}\Big{)}\colon\sigma\in[0,M]\Big{\\}},$

for a sufficiently large constant $M$. Thus we have

$n\lambda_{1}(A_{n})\geqslant\inf\\{g(\sigma)+p(\sigma)\colon\sigma\in[0,M]\\}+O\Big{(}\frac{1}{n}\Big{)}.$
(2.8)

In order to obtain a neat lower bound for $\lambda_{1}(A_{n})$, we need to
choose the coefficients $p_{j}$ in such a way that

$m\coloneqq\inf_{0\leqslant\sigma\leqslant M}\\{g(\sigma)+p(\sigma)\\}$

be maximal. Since $p$ is $2\pi$-periodic and $g$ is strictly increasing with
$g(0)=0$ and $g(\infty)=\infty$, the infimum of $g+p$ must be in the interval
$[0,\pi]$, and consequently we infer that

$m=\inf_{0\leqslant\sigma\leqslant\pi}\\{g(\sigma)+p(\sigma)\\}.$

Consider the integral

$k_{2}\coloneqq\frac{1}{\pi}\int_{0}^{\pi}\\{g(\sigma)+p(\sigma)\\}\operatorname{d\\!}\sigma=\frac{1}{\pi}\int_{0}^{\pi}g(\sigma)\operatorname{d\\!}\sigma,$

which satisfies $k_{2}\geqslant m$, and take $p$ as

$p(\sigma)=\begin{cases}k_{2}-g(\sigma),&\sigma\in[-\pi,\pi];\\\
k_{2}-g(\sigma-2\pi j),&\sigma\in[\pi j,\pi(j+2)],\quad
j\in\mathbb{Z}.\end{cases}$

It is immediate that the function $p$ is $2\pi$-periodic,
$[g(\sigma)+p(\sigma)]_{\sigma\in[-\pi,\pi]}\equiv k_{2}$,
$g(\sigma)+p(\sigma)\geqslant k_{2}$ for $\sigma\notin[-\pi,\pi]$, which
implies $m=k_{2}$. which combined with (2.8) proves the theorem. ∎

Finally, combining the Theorems 2.1 and 2.2 we obtain

$2.2945\approx k_{2}\leqslant n\lambda_{1}(A_{n})\leqslant k_{1}\approx
12.9301,$

for every sufficiently large $n$.

######  Remark 2.3.

In this remark we give a specific account on a comparison of the obtained
results with respect to those available in the literature and which amount so
far to only two works [5, 21]: our bounds are more precise as emphasized in
the next section.

##  3 Few selected numerical experiments

The current section is divided into two parts. In the first we discuss our
theoretical results, regarding the bounds on the minimal eigenvalue of
$A_{n}$, as a function of the matrix order $n$, also in comparison with the
bounds present in the literature. Regarding [21], based on [26], the only
relevant observation is that the minimal eigenvalue of
$T_{n}^{-1}(|\theta|^{2})h^{\alpha}T_{n}(|\theta|^{2-\alpha})$ is well
separated from zero, for any choice of $\alpha\in(0,2)$, and this provides a
qualitative indication that the minimal eigenvalue of ${A_{n}\over n}$
converges to zero with a speed of ${1\over n^{2}}$. Concerning [5], the latter
claim is indeed proved formally, but the constants are not computed, while in
the present note we improve the findings, by determining quite precise
lowerbounds and upperbounds (see Figure 2).

Figure 2: The normalized minimum eigenvalue $n\lambda_{1}(A_{n})$ of $A_{n}$
for different values of $n$. The left image includes the lower and upper
bounds given by Theorems 2.2 and 2.1, respectively, while the right image
shows a standard data range.

The second part is devoted to exploit the spectral information for the use and
the design of specialized iterative algorithms, when solving large linear
systems having ${A_{n}\over n}$ as coefficient matrix. In Table 1 we employ
the standard conjugate gradient (CG), since the coefficient matrix is positive
definite, but we do not expect optimality due to the ill-conditioning of the
related matrix-sequence. We then use the preconditioned conjugate gradient
(PCG) with Strang and Frobenius optimal preconditioners (see [8, 9, 23, 28]
taking the preconditioner in the algebra of the circulant matrices and in the
algebra of sine transforms, containing the Finite Difference discrete
Laplacian $T_{n}(2-2\cos\theta)$). Also in this case it is nontrivial to
obtain optimality, due to the ill-conditioning of the involved matrices, which
grows to infinity as the matrix size tends to infinity.

Few observations are in order: the considered matrix $A_{n}$ is real,
symmetric, and Toeplitz and the $\tau$ algebra is again inherently real,
symmetric, and with Toeplitz generator [9]. The latter statement represents a
qualitative explanation of the fact that the $\tau$ preconditioners perform
substantially better than the analogs in the circulant algebra (see Table 1):
the theoretical ground to the preceding qualitative remark relies in the
notion of good algebras developed in [30].

Notice that the tridiagonal $\tau$ discrete Finite Difference Laplacian
$\Delta_{n}$ is optimal: the related linear system solution is extremely cheap
both in a sequential model (via the Thomas algorithm) and in a parallel model
of computation (via e.g. classical Multigrid methods [17]). We also notice the
remarkable difference between the performance of the optimal $\tau$
preconditioner and of the optimal circulant preconditioner. The reason is
spectral (plus the good algebra argument [30]). The minimal eigenvalue of the
optimal circulant preconditioner is averaged, due to a Cesàro sum effect [9,
28], and hence it behaves as ${1\over n}$ (instead of $1\over n^{2}$) and this
explains the reason why the number of iterations grows as $\sqrt{n}$, in the
light of the classical results on conjugate gradient methods [1]. On the other
hand, the optimal $\tau$ preconditioner matches very well the extremal
eigenvalues of $A_{n}$ (see e.g. [9]).

Thus, as a conclusion, we deduce that the preconditioners $C_{S,n}$ (Strang-
Circulant), $\tau_{N,n}$ (Natural-$\tau$), $\tau_{F,n}$ (Frobenius-Optimal
$\tau$), and $\Delta_{n}$ (Discrete FD Laplacian) are all optimal in the sense
that the iteration count is bounded by constants independent of the matrix
size $n$ and the cost per iteration is that of the Fast Fourier Transform,
which amounts to $O(n\log n)$ arithmetic operations (for a formal notion of
optimality see [2]). As the data indicate, the preconditioners $\tau_{N,n}$,
$\tau_{F,n}$ are the best, but also $\Delta_{n}$ is of interest given its
sparsity.

In Table 2-Table 4 the minimal and maximal eigenvalues of the considered
preconditioned matrices are reported to highlight the efficiency of the
proposed preconditioner. For sure, both the data regarding the preconditioners
$\tau_{N,n}$ and $\tau_{F,n}$ deserve further attention. The outliers analysis
reported in Tables 5 and 6, respectively, seems to show a strong cluster at
$1$. For sure, a weak clustering is observed, being the outliers percentage
decreasing as long as $n$ increases: we notice that the weak clustering can be
deduced theoretically using the GLT theory [16], while the strong clustering
is nontrivial given the ill-conditioning of $A_{n}$.

The last choice is a multigrid method designed with the use of the spectral
information available. In fact, the projector and the restriction operators
are those classically used when dealing with the standard discrete Laplacian
$\Delta_{n}=T_{n}(2-2\cos\theta)$. Indeed, even if $A_{n}\over n$ is dense and
seemingly there are no similarities with $\Delta_{n}$, from a spectral point
of view both matrix-sequences have minimal eigenvalue collapsing to zero as
$1\over n^{2}$, up to some moderate constants.

More precisely, the transfer operators are designed, in a very standard way,
as the product of the Toeplitz matrix generated by the symbol $2+2\cos\theta$
and a proper cutting matrix (see e.g. [15, 17]). A pure algebraic multigrid is
considered, so that matrices at coarser levels are obtained by projection via
the transfer operators according to the Galerkin condition. When setting
smoothers, we tested several choices by considering a standard Gauss-Seidel
iteration, or PCG with sine transform preconditioners, according to the
natural approach and the Frobenius optimal choice: $\nu_{\textrm{pre}}$ steps
are applied for the presmoother and $\nu_{\textrm{post}}$ steps for the
postsmoother, respectively.

As it can be seen in Table 7, the numerical results with the multigrid choice
are fully satisfactory, as well: the method is optimal [2] in the sense that
the iteration count is bounded by a constant independent of $n$ and the cost
per iteration is proportional to that of the matrix-vector product.
A special mention deserves a last test in which we set as presmoother the PCG
with Finite Difference discrete Laplacian preconditioner with
$\nu_{\textrm{pre}}=1$, and as postsmoother the PCG with Natural-$\tau$ or
Frobenius-Optimal $\tau$ preconditioner with $\nu_{\textrm{post}}=1$, but only
at the finest level. In all the coarser levels the smoothers are simply chosen
as one iteration of the standard Gauss-Seidel iteration, so further reducing
the computational cost. The number of required iterations is $3$ as for the
case $\gamma$.
Though the efficiency is comparable with those of the preconditioned PCG in
Table 1 in the present unilevel setting, multigrid could become the most
promising choice in the multilevel one, due to the theoretical barriers
studied in [24, 32, 33], regarding the non optimality of the PCG with matrix-
algebra preconditioners in the multilevel context.

Size | $I_{n}$ | $C_{S,n}$ | $C_{F,n}$ | $\tau_{N,n}$ | $\tau_{F,n}$ | $\Delta_{n}$
---|---|---|---|---|---|---
32 | 34 | 7 | 11 | 6 | 5 | 8
64 | 73 | 8 | 14 | 5 | 5 | 9
128 | 154 | 8 | 15 | 5 | 5 | 10
256 | 307 | 8 | 18 | 5 | 5 | 10
512 | 593 | 8 | 22 | 5 | 5 | 11
1024 | 1095 | 8 | 27 | 5 | 5 | 11
2048 | 2112 | 8 | 34 | 5 | 5 | 12
Table 1: Number of PCG’s iterations to reach convergence with respect to scaled residual less than $10^{-7}$ with preconditioner $I_{n}$ (no preconditioning), $C_{S,n}$ (Strang-Circulant), $C_{F,n}$ (Frobenius-Optimal Circulant), $\tau_{N,n}$ (Natural-$\tau$), $\tau_{F,n}$ (Frobenius-Optimal $\tau$), and $\Delta_{n}$ (Discrete FD Laplacian). Size |  $\lambda_{\min}(C_{S,n}^{-1}A_{n})$ |  $\lambda_{\max}(C_{S,n}^{-1}A_{n})$ |  $\lambda_{\min}(C_{F,n}^{-1}A_{n})$ |  $\lambda_{\max}(C_{F,n}^{-1}A_{n})$
---|---|---|---|---
$32$ |  $5.2258\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{-1}$ |  $3.4325\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{1}$ |  $1.7524\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{-1}$ |  $4.1234\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{0}$
$64$ |  $5.1302\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{-1}$ |  $5.4268\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{1}$ |  $1.1034\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{-1}$ |  $5.6639\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{0}$
$128$ |  $5.0743\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{-1}$ |  $9.1152\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{1}$ |  $6.6754\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{-2}$ |  $7.8412\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{0}$
$256$ |  $5.0419\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{-1}$ |  $1.5814\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{2}$ |  $3.9098\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{-2}$ |  $1.0911\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{1}$
$512$ |  $5.0234\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{-1}$ |  $2.8003\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{2}$ |  $2.2338\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{-2}$ |  $1.5232\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{1}$
$1024$ |  $5.0129\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{-1}$ |  $5.0319\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{2}$ |  $1.2526\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{-2}$ |  $2.1315\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{1}$
$2048$ |  $5.0071\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{-1}$ |  $9.1428\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{2}$ |  $6.9259\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{-3}$ |  $2.9876\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{1}$
Table 2: Minimal and maximal eigenvalues of preconditioned matrices with preconditioner $C_{S,n}$ (Strang-Circulant) and $C_{F,n}$ (Frobenius-Optimal Circulant).  Size |  $\lambda_{\min}(\tau_{N,n}^{-1}A_{n})$ |  $\lambda_{\max}(\tau_{N,n}^{-1}A_{n})$ |  $\lambda_{\min}(\tau_{F,n}^{-1}A_{n})$ |  $\lambda_{\max}(\tau_{F,n}^{-1}A_{n})$
---|---|---|---|---
$32$ |  $8.1574\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{-1}$ |  $1.1475\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{0}$ |  $9.3206\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{-1}$ |  $1.1356\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{0}$
$64$ |  $8.0608\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{-1}$ |  $1.1527\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{0}$ |  $9.0938\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{-1}$ |  $1.1463\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{0}$
$128$ |  $7.9895\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{-1}$ |  $1.1584\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{0}$ |  $8.8984\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{-1}$ |  $1.1551\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{0}$
$256$ |  $7.9396\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{-1}$ |  $1.1641\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{0}$ |  $8.7441\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{-1}$ |  $1.1624\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{0}$
$512$ |  $7.9063\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{-1}$ |  $1.1693\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{0}$ |  $8.6271\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{-1}$ |  $1.1685\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{0}$
$1024$ |  $7.8851\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{-1}$ |  $1.1742\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{0}$ |  $8.5401\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{-1}$ |  $1.1737\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{0}$
$2048$ |  $7.8727\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{-1}$ |  $1.1785\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{0}$ |  $8.4763\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{-1}$ |  $1.1783\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{0}$
Table 3: Minimal and maximal eigenvalues of preconditioned matrices with preconditioner $\tau_{N,n}$ (Natural-$\tau$) and $\tau_{F,n}$ (Frobenius-Optimal $\tau$).  Size |  $\lambda_{\min}(\Delta_{n}^{-1}A_{n})$ |  $\lambda_{\max}(\Delta_{n}^{-1}A_{n})$
---|---|---
$32$ |  $3.1941\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{-1}$ |  $5.4802\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{-1}$
$64$ |  $2.6529\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{-1}$ |  $5.3509\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{-1}$
$128$ |  $2.2651\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{-1}$ |  $5.2979\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{-1}$
$256$ |  $1.9750\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{-1}$ |  $5.2727\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{-1}$
$512$ |  $1.7504\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{-1}$ |  $5.2604\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{-1}$
$1024$ |  $1.5713\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{-1}$ |  $5.2543\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{-1}$
$2048$ |  $1.4252\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{-1}$ |  $5.2513\raisebox{0.85358pt}{\mbox{\tiny$\times$}}10^{-1}$
Table 4: Minimal and maximal eigenvalues of preconditioned matrices with preconditioner $\Delta_{n}$ (Discrete FD Laplacian).  Size | $\varepsilon=10^{-1}$ | $\varepsilon=10^{-2}$
---|---|---
$n_{\textrm{out}}^{l}$ | $n_{\textrm{out}}^{r}$ | $\%$ | $n_{\textrm{out}}^{l}$ | $n_{\textrm{out}}^{r}$ | $\%$
32 | 1 | 2 | 9.37% | 3 | 4 | 21.8%
64 | 1 | 2 | 4.69% | 3 | 4 | 10.9%
128 | 1 | 2 | 2.34% | 4 | 4 | 6.25%
256 | 2 | 2 | 1.56% | 5 | 4 | 3.52%
512 | 2 | 2 | 0.78% | 5 | 4 | 1.76%
1024 | 2 | 2 | 0.39% | 5 | 4 | 0.88%
2048 | 2 | 2 | 0.19% | 6 | 4 | 0.49%
Table 5: Number of left ($n_{\textrm{out}}^{l}$) and right ($n_{\textrm{out}}^{r}$) outliers (eigenvalues not belonging to $(1-\varepsilon,1+\varepsilon)$ with $\varepsilon=10^{-1}$ and $10^{-2}$ and their percentage with respect to the dimension in the case of Natural-$\tau$ preconditioner $\tau_{N,n}$. Size | $\varepsilon=10^{-1}$ | $\varepsilon=10^{-2}$
---|---|---
$n_{\textrm{out}}^{l}$ | $n_{\textrm{out}}^{r}$ | $\%$ | $n_{\textrm{out}}^{l}$ | $n_{\textrm{out}}^{r}$ | $\%$
32 | 0 | 2 | 6.25% | 18 | 4 | 68.7%
64 | 0 | 2 | 3.13% | 2 | 6 | 12.5%
128 | 1 | 2 | 2.34% | 2 | 9 | 8.59%
256 | 1 | 2 | 1.17% | 3 | 10 | 5.08%
512 | 2 | 2 | 0.78% | 4 | 10 | 2.73%
1024 | 2 | 2 | 0.39% | 4 | 10 | 1.36%
2048 | 2 | 2 | 0.19% | 4 | 10 | 0.68%
Table 6: Number of left ($n_{\textrm{out}}^{l}$) and right ($n_{\textrm{out}}^{r}$) outliers (eigenvalues not belonging to $(1-\varepsilon,1+\varepsilon)$ with $\varepsilon=10^{-1}$ and $10^{-2}$ and their percentage with respect the dimension in the case of Frobenius-Optimal $\tau$ preconditioner $\tau_{F,n}$. Size |  Case $\alpha$ |  Case $\beta$ |  Case $\gamma$ |  Case $\delta$
---|---|---|---|---
TGM | Vcycle | TGM | Vcycle | TGM | Vcycle | TGM | Vcycle
31 | 9 | 9 | 4 | 4 | 3 | 3 | 2 | 2
63 | 10 | 10 | 4 | 4 | 3 | 3 | 2 | 2
127 | 11 | 11 | 4 | 4 | 3 | 3 | 2 | 2
255 | 11 | 11 | 4 | 4 | 3 | 3 | 2 | 2
511 | 11 | 11 | 3 | 3 | 3 | 3 | 2 | 2
1023 | 11 | 11 | 4 | 4 | 3 | 3 | 2 | 2
2047 | 11 | 11 | 3 | 3 | 3 | 3 | 2 | 2
Table 7: Number of Multigrid iterations to reach convergence with respect to
scaled residual less than $10^{-7}$: Case $\alpha$ = Gauss-Seidel
($\nu_{\textrm{pre}}=1$)/Gauss–Seidel ($\nu_{\textrm{post}}=1$), case $\beta$
= Gauss–Seidel ($\nu_{\textrm{pre}}=1$) /PCG with Natural-$\tau$ or Frobenius-
Optimal $\tau$ preconditioner ($\nu_{\textrm{post}}=1$), case $\gamma$ = PCG
with Discrete FD Laplacian preconditioner ($\nu_{\textrm{pre}}=1$ or $2$) /PCG
with Natural-$\tau$ or Frobenius-Optimal $\tau$ preconditioner
($\nu_{\textrm{post}}=1$), case $\delta$ = PCG with Discrete FD Laplacian
preconditioner ($\nu_{\textrm{pre}}=1$) /PCG with Natural-$\tau$ or Frobenius-
Optimal $\tau$ preconditioner this time ($\nu_{\textrm{post}}=2$).

##  4 Conclusions

In the current note we have considered a type of matrix stemming when
considering the numerical approximations of distributed order FDEs (see [5,
21] for example). The main contribution relies in precise bounds for the
minimal eigenvalue of the involved matrices. In fact the new presented bounds
improve those already present in the literature and give a more accurate
spectral information. The latter knowledge has been used in the design of fast
numerical algorithms for the associated linear systems approximating the given
distributed order FDEs: an interesting challenge is to consider a
$d$-dimensional version of the considered FDE (see [21]), in order to show how
the presented techniques and numerical methods can be adapted and extended in
$d$-level setting with $d\geqslant 2$.

## References

  * [1] Axelsson, O., and Lindskog, G. On the rate of convergence of the preconditioned conjugate gradient method.  Numer. Math. 48, 5 (1986), 499–523.
  * [2] Axelsson, O., and Neytcheva, M. The algebraic multilevel iteration methods-theory and applications.  In Proceedings of the Second International Colloquium on Numerical Analysis (1994), VSP, Utrecht, pp. 13–24.
  * [3] Bogoya, M., Böttcher, A., Grudsky, S.M., and Maximenko, E.A. Eigenvalues of Hermitian Toeplitz matrices with smooth simple-loop symbols.  Oper. Theory: Adv. and Appl. 422 (2015), 1308–1334.
  * [4] Bogoya, M., Böttcher, A., Grudsky, S.M., and Maximenko, E.A. Eigenvectors of Hermitian Toeplitz matrices with smooth simple-loop symbols.  Linear Algebra Appl. 493 (2016), 606–637.
  * [5] Bogoya, M., Grudsky, S.M., Mazza, M., and Serra-Capizzano, S. On the spectrum and asymptotic conditioning of a class of positive definite Toeplitz matrix-sequences, with application to fractional-differential approximations.  arXiv:2112.02685 (2021).
  * [6] Böttcher, A., and Grudsky, S.M. On the condition numbers of large semi-definite Toeplitz matrices.  Linear Algebra Appl. 279, 1/3 (1998), 285–301.
  * [7] Böttcher, A., and Grudsky, S.M. Spectral properties of banded Toeplitz matrices.  Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2005.
  * [8] Chan, R., and Jin, X. An introduction to iterative Toeplitz solvers.  Society for Industrial and Applied Mathematics (SIAM), 2007.
  * [9] Di-Benedetto, F., and Serra-Capizzano, S. A unifying approach to abstract matrix algebra preconditioning.  Numer. Math. 82, 1 (1999), 57–90.
  * [10] Donatelli, M., Garoni, C., Manni, C., Serra-Capizzano, S., and Speleers, H. Robust and optimal multi-iterative techniques for IgA Galerkin linear systems.  Comput. Methods Appl. Mech. Engrg. 284 (2015), 230–264.
  * [11] Donatelli, M., Garoni, C., Manni, C., Serra-Capizzano, S., and Speleers, H. Symbol-based multigrid methods for Galerkin B-spline isogeometric analysis.  SIAM. J. Numer. Anal. 55, 1 (2017), 31–62.
  * [12] Donatelli, M., Mazza, M., and Serra-Capizzano, S. Spectral analysis and structure preserving preconditioners for fractional diffusion equations.  J. Comput. Phys. 307 (2016), 262–279.
  * [13] Donatelli, M., Mazza, M., and Serra-Capizzano, S. Spectral analysis and multigrid methods for finite volume approximations of space-fractional diffusion equations.  SIAM J. Sci. Comput. 40, 6 (2018), A4007–A4039.
  * [14] Ervin, V., and Roop, J. Variational formulation for the stationary fractional advection dispersion equation.  Numer. Methods Partial Differ. Equ. 22, 3 (2006), 558–576.
  * [15] Fiorentino, G., and Serra-Capizzano, S. Multigrid methods for Toeplitz matrices.  Calcolo 28, 3-4 (1992), 283–305.
  * [16] Garoni, C., and Serra-Capizzano, S. Generalized Locally Toeplitz sequences: Theory and applications. Vol. I.  Springer, Cham, 2017.
  * [17] Hackbusch, W. Multigrid methods and applications, vol. 4 of Springer Series in Computational Mathematics.  Springer-Verlag, Berlin, 1985.
  * [18] Hejazi, H., Moroney, T., and Liu, F. Stability and convergence of a finite volume method for the space fractional advection-dispersion equation.  J. Comput. Appl. Math. 255 (2014), 684–697.
  * [19] Lin, Z., and Wang, D. A finite element formulation preserving symmetric and banded diffusion stiffness matrix characteristics for fractional differential equations.  Comput. Mech. 62, 2 (2018), 185–211.
  * [20] Mao, Z., and Karniadakis, G. A spectral method (of exponential convergence) for singular solutions of the diffusion equation with general two-sided fractional derivative.  SIAM J. Numer. Anal. 56, 1 (2018), 24–49.
  * [21] Mazza, M., Serra-Capizzano, S., and Usman, M. Symbol-based preconditioning for Riesz distributed-order space-fractional diffusion equations.  Electr. Trans. Num. Anal. 54 (2021), 499–513.
  * [22] Meerschaert, M., and Tadjeran, C. Finite difference approximations for fractional advection-dispersion flow equations.  J. Comput. Appl. Math. 172, 1 (2004), 65–77.
  * [23] Ng, M. Iterative methods for Toeplitz systems.  Oxford University Press, 2004.
  * [24] Noutsos, D., Serra-Capizzano, S., and Vassalos, P. Matrix algebra preconditioners for multilevel Toeplitz systems do not insure optimal convergence rate.  Theoret. Comput. Sci. 315, 2-3 (2004), 557–579.
  * [25] Pedas, A., and Tamme, E. On the convergence of spline collocation methods for solving fractional differential equations.  J. Comput. Appl. Math. 235, 12 (2011), 3502–3514.
  * [26] Serra-Capizzano, S. New PCG based algorithms for the solution of Hermitian Toeplitz systems.  Calcolo 32 (1995), 53–176.
  * [27] Serra-Capizzano, S. On the extreme eigenvalues of Hermitian (block) Toeplitz matrices.  Linear Algebra Appl. 270 (1998), 109–129.
  * [28] Serra-Capizzano, S. A Korovkin-type theory for finite Toeplitz operators via matrix algebras.  Numer. Math. 82, 1 (1999), 117–142.
  * [29] Serra-Capizzano, S. The rate of convergence of Toeplitz based PCG methods for second order nonlinear boundary value problems.  Numer. Math. 81, 3 (1999), 461–495.
  * [30] Serra-Capizzano, S. Toeplitz preconditioners constructed from linear approximation processes.  SIAM J. Matrix Anal. Appl. 20, 2 (1999), 446–465.
  * [31] Serra-Capizzano, S. Convergence analysis of two-grid methods for elliptic Toeplitz and PDEs matrix-sequences.  Numer. Math. 92, 3 (2002), 433–465.
  * [32] Serra-Capizzano, S. Matrix algebra preconditioners for multilevel Toeplitz matrices are not superlinear. Special issue on structured and infinite systems of linear equations.  Linear Algebra Appl. 343/344 (2002), 303–319.
  * [33] Serra-Capizzano, S., and Tyrtyshnikov, E. Any circulant-like preconditioner for multilevel matrices is not superlinear.  SIAM J. Matrix Anal. Appl. 21, 2 (1999), 431–39.
  * [34] Tian, W., Zhou, H., and Deng, W. A class of second order difference approximations for solving space fractional diffusion equations.  Math. Comput. 84, 294 (2015), 1703–1727.
  * [35] Wang, H., and Du, N. A superfast-preconditioned iterative method for steady-state space-fractional diffusion equations.  J. Comput. Phys. 240 (2013), 49–57.
  * [36] Xu, K., and Darve, E. Isogeometric collocation method for the fractional Laplacian in the 2D bounded domain.  Computer Methods in Applied Mechanics and Engineering 364 (2020), 112936.
  * [37] Zeng, F., Mao, Z., and Karniadakis, G. A generalized spectral collocation method with tunable accuracy for fractional differential equations with end-point singularities.  SIAM J. Sci. Comput. 39, 1 (2017), A360–A383.