# A Stochastic Method for Solving Time-Fractional Differential Equations

Nicolas L. Guidotti<EMAIL_ADDRESS>Juan A. Acebrón
<EMAIL_ADDRESS>José Monteiro<EMAIL_ADDRESS>INESC-ID,
Instituto Superior Técnico, Universidade de Lisboa, Portugal  Department of
Mathematics, Carlos III University of Madrid, Spain

###### Abstract

We present a stochastic method for efficiently computing the solution of time-
fractional partial differential equations (fPDEs) that model anomalous
diffusion problems of the subdiffusive type. After discretizing the fPDE in
space, the ensuing system of fractional linear equations is solved resorting
to a Monte Carlo evaluation of the corresponding Mittag-Leffler matrix
function. This is accomplished through the approximation of the expected value
of a suitable multiplicative functional of a stochastic process, which
consists of a Markov chain whose sojourn times in every state are Mittag-
Leffler distributed. The resulting algorithm is able to calculate the solution
at conveniently chosen points in the domain with high efficiency. In addition,
we present how to generalize this algorithm in order to compute the complete
solution. For several large-scale numerical problems, our method showed
remarkable performance in both shared-memory and distributed-memory systems,
achieving nearly perfect scalability up to $16,384$ CPU cores.

###### keywords:

Monte Carlo method , time-fractional differential equations , anomalous
diffusion , Mittag-Leffler function , matrix functions , parallel algorithms ,
high performance computing

††journal: Journal of Computational Physics

##  1 Introduction

Among the many non-conventional statistical phenomena observed experimentally
in the last few years and characterized typically by unusual scaling
exponents, the anomalous diffusion stands out for its major impact in a
variety of scientific disciplines [1]. To cite just a few relevant examples,
it has been already observed in turbulent atmosphere, transport in amorphous
solids, diffusion in cell nucleus, magnetic resonance imaging, search
behavior, and in finance for an accurate modeling of the stock price
fluctuations.

In the anomalous diffusion, the mean squared displacement (MSD) vary
nonlinearly in time, i.e., $\langle x^{2}(t)\rangle\sim t^{\alpha}$ with
$\alpha\neq 1$, during the entire process. Here $x(t)$ is the relative
position at time $t$ of a particle with respect to a given reference point.
This contrasts with the classical diffusion problems, modeled typically by a
Brownian motion, where it is well known that the MSD grows linearly in time.
In practice, there are two different types of regimes according to the value
of the scaling exponent $\alpha$. The phenomenon is termed as subdiffusion
when $\alpha<1$ and, superdiffusion when $\alpha>1$. It has been proposed in
the literature several pathways leading to anomalous diffusion. The most
commonly found are the presence of long-range correlations, non-identical
displacements, or non-finite mean or variance of the probability density
function for the trajectory displacements.

Within this context, fractional calculus has been proven to be extremely
useful in order to provide realistic models for many real-life processes and
phenomena [2, 3]. The main reason for this is due to fractional derivative
operators being in practice non-local operators, and therefore they are
specially suited for describing the long-time memory and spatial heterogeneity
effects typically found in any anomalous diffusion problem.

The formal solution to the simplest dynamical problem defined by a system of
linear fractional rate equations for an $n$-dimensional process
$\mathbf{y}(t)$ subject to a given initial condition $\mathbf{y}(0)={\bf
y_{0}}$ is

$\mathbf{y}(t)=E_{\alpha}(\mathbf{A}\,t^{\alpha})\,{\bf y_{0}},$ (1)

where $\mathbf{A}$ is an $n\times n$ matrix of constant coefficients. When
$\alpha=1$, the Mittag-Leffler (ML) function
$E_{\alpha}(\mathbf{A}\,t^{\alpha})$ reduces to the matrix exponential
$e^{\mathbf{A}\,t}$. The ML function also appears in other scientific domains,
such as biology [4, 5], physics [6, 7, 8], control systems [7, 9] and complex
network analysis [10, 11]. It is important to remark here that this system of
linear rate equations serve as building-blocks for more complex systems, and
therefore becomes of paramount importance investigating algorithms capable of
computing it efficiently.

However, an efficient algorithm for solving large-scale problems is still
missing. Therefore, the main contribution of this paper is to propose a
probabilistic method for the efficient computation of the ML matrix function.
This method is based on random walks that change states according to a
suitable continuous-time Markov chain, whose sojourn times follows a Mittag-
Leffler distribution. Such method allows us to compute the ML matrix function
for large-scale matrices, focusing primarily here on the action of the
function over a vector. Our second contribution is an extensive analysis of
the performance and convergence of our method using a few relevant numerical
examples, including a formal description of the variance and a performance
comparison against classical algorithms. Our third and last contribution was
to parallelize and analyze the scalability of the stochastic algorithm for a
large number of processors (up to $16,384$ CPU cores) in the Karolina
Supercomputer located in IT4Innovations National Supercomputing Centre.

The paper is organized as follows. In the next section we give an overview of
the state of the art in this area. Section 3 describes the probabilistic
method and its practical implementation. This section also shows how a program
can generate random numbers from the Mittag-Leffler probability distribution.
In Section 4, several numerical experiments are presented in order to
characterize the performance and numerical errors of the stochastic method
compared to classical approaches. Finally, in Section 5, we highlight the main
results and conclude the paper.

##  2 Background

One of the major drawbacks of describing a problem in terms of fractional
operators is the lack of uniqueness in the mathematical formulation. This
contrasts sharply with the classical calculus. In fact, in fractional calculus
there exists a variety of different definitions of fractional operators
(differentials and integrals), and they are chosen according to a given set of
assumptions imposed to satisfy different needs and constraints of the
particular problem. There are already in the literature a number of excellent
texts that review the mathematics of fractional calculus, e.g. [8]. Among the
different definitions of fractional operators, the Caputo fractional
derivative appears as the most commonly used for modeling temporal evolution
phenomena, which accounts for past interactions and also non-local properties.
Moreover, mathematically the initial conditions imposed for fractional
differential equations with Caputo derivatives coincides with the integer-
order differential equations. Such derivative is defined as follows:

${}^{C}D_{t}^{\alpha}f(t)=\int_{0}^{t}\frac{f^{(n)}(\tau)}{(t-\tau)^{\alpha+1-n}}d\tau$
(2)

with $n$ an integer such that $n-1<\alpha<n$. Note that if $0<\alpha<1$, $n$
turns out to be 1.

The simplest dynamical process described in terms of fractional operators
consists in an initial value problem for a linear fractional rate equation,
which for a given variable $y(t)$ is given by

${}^{C}D_{t}^{\alpha}y(t)=-\lambda y(t),\quad y(0)=y_{0}.$ (3)

The solution of this equation can not be expressed in a closed-form, being
rather an infinite series. It was first obtained by the mathematician Mittag-
Leffler, and is formally written as

$y(t)=y_{0}\,E_{\alpha}(-\lambda t^{\alpha}),$ (4)

where $E_{\alpha}(t)$ is the function defined as

$E_{\alpha}(-\lambda\,t^{\alpha})=\sum_{k=0}^{\infty}\frac{(-\lambda\,t^{\alpha})^{k}}{\Gamma(k\alpha+1)},$
(5)

and carries his name. Here $\Gamma(x)$ denotes the Euclidean Gamma function.
Note that when $\alpha=1$, it reduces to the exponential function.
Essentially, the Mittag-Leffler function generalizes the exponential function,
and has been considered for many as the Queen function of the fractional
calculus [12, 13, 14]. For complex parameters $\alpha$ and $\beta$ with
$\mathbb{R}\\{\alpha\\}>0$, the Mittag-Leffler function (ML function) can be
further generalized, and was defined as

$E_{\alpha,\beta}(z)=\sum_{k=0}^{\infty}{\frac{z^{k}}{\Gamma(k\alpha+\beta)}},z\in\mathbb{C}.$
(6)

Here $E_{\alpha,\beta}(z)$ is an entire function of order
$\rho=1/\mathbb{R}\\{\alpha\\}$ and type 1. The exponential function can be
recovered by setting $\alpha=\beta=1$ since $\Gamma(k~{}+~{}1)~{}=~{}k!$ for
$k\in\mathbb{N}$. When $\beta=1$, the notation of the ML function can be
simplified to $E_{\alpha}(z)=E_{\alpha,1}(z)$.

In this paper we consider the $n$-dimensional case of a system described by a
set time-fractional linear differential equations, whose solution is given by
(1). We therefore need to evaluate the ML function over a potentially very
large system’s matrix.

A naive approach to evaluate the ML function comes directly from its
definition in (6). However, this is ill-advised for most applications since
the convergence of the series is very slow if either the modulus $|z|$ is
large or the value of $\alpha$ is small. Moreover, the terms of the series may
grow very large before decreasing, which can cause overflows and catastrophic
numerical cancellations when using standard double-precision arithmetic (IEEE
754). For this reason, the scalar ML function is usually evaluated through
alternative methods [15, 16, 17, 18], such as the Taylor series, inversion of
the Laplace transformation, integral representation and other techniques.

To evaluate the ML matrix function, Garrappa [19] recently proposed a Schur-
Parlett algorithm [20], which in practice requires computing its derivatives.
This method, however, is only suited for small matrices since the Schur-
Parlett algorithm scales with $O(n^{3})$, but it can go up to $O(n^{4})$
depending on the distribution of the eigenvalues. Moreover, if the eigenvalues
are highly clustered, the algorithm may require high-order derivatives which
are more expensive to compute, while also being less accurate. There is also a
Krylov-based method for computing the action of a Mittag-Leffler function over
a vector [21], however, the code is not publicly available.

As an alternative to deterministic methods, Monte Carlo algorithms for
approximating matrix functions have already been proposed in the past [22, 23,
24, 25, 26], primarily for solving linear systems. The main ideia of these
algorithms consists in generating random paths corresponding to the
realizations of a discrete Markov chain. Those random paths evolves through
the different indices of the matrix $\mathbf{A}$, which is used as the
generator. Essentially, the method corresponds to a Monte Carlo sampling of
the Neumann series [27] of the matrix. The convergence of this method was
rigorously established in [26, 25, 28]. It is broadly accepted that these
stochastic methods offer interesting features from a computational point of
view, such as being easily parallelizable, fault-tolerant, and more suited to
heterogeneous architectures (which are extremely important attributes in view
of the current high performance computers). However, the truth is that they
also exhibit some significant weaknesses, namely a very slow convergence rate
to the solution. This can cause the underlying algorithms to be highly
demanding computationally, especially when high accuracy solutions are
required.

The extension of the Monte Carlo method for other matrix functions, such as
the matrix exponential, was only accomplished recently in [29]. The method
consists in generating continuous-time random walks over the matrix
$\mathbf{A}$, whose holding time follows an exponential distribution. Later,
[30] applied a multilevel technique to improve the performance of this method.

For the specific case of the one-dimensional fractional diffusion equation,
[31] introduced a Monte Carlo method for simulating discrete time random walks
which are used for obtaining numerical solutions of the fractional equation,
but not for computing the ML matrix function.

##  3 Description of the Probabilistic Method

In this section, we describe the probabilistic representation of the solution
vector as well as a practical algorithm for approximating the action of the
Mittag-Leffler function over a vector using this representation. With slight
modifications, the algorithm can either calculate a single entry of the
solution vector or the full solution. We also present a way to sample random
numbers from the Mittag-Leffler distribution.

###  3.1 Mathematical description

Consider: $\mathbf{A}=(a_{ij})$, a general n-by-n matrix with $a_{ii}<0$
$\forall i$; $\mathbf{u}$, a given n-dimensional vector; and $\mathbf{y}$, an
n-dimensional vector. This last vector corresponds to the solution vector
after computing the action of the Mittag-Leffler function over the vector
$\mathbf{u}$, that is
$\mathbf{y}=E_{\alpha,\beta}(\mathbf{A}\,t^{\alpha})\,{\bf u}$ with
$0<\alpha\leq 1$, $\beta>0$ and $t\geq 0$. Here $t$ denotes the value of time
when the solution is computed.

Let us define the following matrices: a diagonal matrix $\mathbf{D}$,
represented hereafter as a vector ${\bf d}$, with entries $d_{ij}=0$ for
$\forall i\neq j$ and $d_{ii}=d_{i}=a_{ii}$ for $i=1,\ldots,n$; $\mathbf{M}$,
a matrix obtained as $\mathbf{M}=\mathbf{A}-\mathbf{D}$, and hence with zero
diagonal entries; $\mathbf{Q}$, the matrix with entries $q_{ij}$ given by

$q_{ij}=\begin{cases}0,&\text{if}\ i=j\\\
\frac{|m_{ij}|}{\sum_{j=1}^{n}|m_{ij}|},&\text{otherwise;}\end{cases}$ (7)

$\mathbf{G}$, a matrix with entries $g_{ij}$ taking values $1$ when
$m_{ij}\geq 0$, and $-1$ otherwise; and finally $\mathbf{R}$, a matrix with
entries

$r_{ij}=\begin{cases}r_{i}=\frac{\sum_{j=1}^{n}|m_{ij}|}{d_{i}},&\text{if}\
i=j\\\ 0,\quad\text{otherwise;}\end{cases}$ (8)

Note that $\sum_{j=1}^{n}q_{ij}=1\;\;\forall i$, and therefore matrix
$\mathbf{Q}$ can be used as the transition rate matrix of a given Markov
chain. Furthermore, it holds that
$\mathbf{M}=\mathbf{D}\,\mathbf{R}\,(\mathbf{G}*\mathbf{Q})$, where $*$
denotes element-wise matrix multiplication.

######  Theorem 1.

Let $\\{X_{t}:t\geq 0\\}$ be a stochastic process with finite state space
$\Omega~{}=~{}\\{1,2,\dots,n\\}$ given by

$X_{t}=\sum_{m=1}^{\infty}Z_{m-1}\mathds{1}_{[T_{m-1}\leq t\leq T_{m}]}.$ (9)

Such process changes states according to a Markov chain
$Z=(Z_{m})_{m\in\mathbb{N}}$, which takes values in $\Omega$ and $\mathbf{Q}$
is the corresponding transition matrix. Here, $T_{k}$ is the time of the
$k$-th event, and $\mathds{1}_{E}$ denotes the indicator function, being $1$
or $0$ depending on whether the event $E$ occurs. The sojourn times in the
$i$-th state follow the Mittag-Leffler distribution $ML_{\alpha}(d_{i})$.
Then, we have that any entry $i$ of the solution vector $\mathbf{y}$ can be
represented probabilistically as

$y_{i}=\mathds{E}\left[u_{X_{0}}\mathds{1}_{[T_{0}>t]}+\left(\prod_{j=1}^{\eta}r_{X_{T_{j-1}}}\,g_{X_{T_{j-1}},X_{T_{j}}}\right)u_{X_{T_{\eta}}}\mathds{1}_{[T_{0}\leq
t]}\right],$ (10)

where $X_{0}=i$. Here $\mathds{E}$ is the expected value with respect to the
joint distribution function of the random variables $T$ and $\eta$, where
$\eta$ is the number of events occurring between $0$ and $t$.

###### Proof.

Our goal is to prove that the expected value of the functional of the
stochastic process $X_{t}$ in (10) coincides with the vector solution
$\mathbf{y}$. As a first step in the proof, it is required to compute the
joint distribution function of the random variables $T$ and $\eta$. This
function depends on the probability of jumping from state $i$ to state $j$ at
time $t$, $\mathds{P}(X_{t}=j|X_{0}=i)$. Let us define $\mathbf{P}$ as the
probability matrix such that
$\mathbf{P}(t)=(p_{ij})=\mathds{P}(X_{t}=j|X_{0}=i)$, then it holds that

$p_{ij}(t)=\delta_{ij}\,\mathds{P}(t_{i}>t)+\int_{0}^{t}ds\,f_{i}(s)\sum_{l\neq
i}\mathds{P}(X_{\tau}=l,X_{t}=j|\tau=s),$ (11)

where $t_{i}$ corresponds to the instant of time of a first event when the
state is $i$. Here the function
$f_{i}(s)=-d_{i}\,s^{\alpha-1}\,E_{\alpha,\alpha}(d_{i}\,s^{\alpha})$ is a
density function. In fact, it satisfies that

$\int_{0}^{\infty}ds\,f_{i}(s)=1,$ (12)

and is called the Mittag-Leffler density function [13]. Then the probability
of no events taking place in the interval of time $[0,t]$,
$\mathds{P}(t_{i}>t)$, is given by

$\mathds{P}(t_{i}>t)=\int_{t}^{\infty}ds\,f_{i}(s)=E_{\alpha}(d_{i}\,t^{\alpha}).$
(13)

Therefore, it follows that

$p_{ij}(t)=\delta_{ij}\,E_{\alpha}(d_{i}\,t^{\alpha})+\int_{0}^{t}ds\,f_{i}(s)\sum_{l\neq
i}q_{il}\,p_{lj}(t-s).$ (14)

This is an integral equation for $p_{ij}(t)$ that can be solved recursively,
and the solution can be rewritten formally as

$p_{ij}(t)=\sum_{\eta=0}^{\infty}p_{ij}^{(\eta)}(t),$ (15)

where $p_{ij}^{(\eta)}(t)$ denotes the probability of having a transition from
$i$ to $j$ when $\eta$ events occur during the time interval $[0,t]$. Let
$\nu_{i}^{(\eta)}$ be the corresponding contribution to the solution $y_{i}$
for $\eta$ events. If no events occur, the contribution to the solution
$\nu_{i}^{(0)}$ can be readily obtained as
$E_{\alpha}(d_{i}\,t^{\alpha})\,u_{i}$. The remaining contributions can be
obtained as follows. For one event, from (14) the probability is given by

$p_{ij}^{(1)}(t)=\int_{0}^{t}ds\,f_{i}(s)\sum_{l\neq
i}q_{il}\,\delta_{lj}\,E_{\alpha}(d_{l}\,(t-s)^{\alpha}),$ (16)

and thus the corresponding contribution to the solution is

$\nu_{i}^{(1)}=\sum_{j=1}^{n}\int_{0}^{t}ds\,f_{i}(s)\sum_{l\neq
i}r_{i}\,g_{il}\,q_{il}\,\delta_{lj}E_{\alpha}(d_{l}\,(t-s)^{\alpha})\,u_{j}.$
(17)

Note that the equation above can be rewritten equivalently in matrix form as

$\boldsymbol{\nu}^{(1)}=\int_{0}^{t}ds\,{\bf
f}(s)\mathbf{R}\,(\mathbf{G}*\mathbf{Q})\,E_{\alpha}(\mathbf{D}\,(t-s)^{\alpha})\,{\bf
u},$ (18)

where $\boldsymbol{\nu}^{(1)}$ and ${\bf f}(s)$ are vectors with components
$\nu_{i}^{(1)}$ and $f_{i}(s)$, respectively. Since
$\mathbf{R}\,(\mathbf{G}*\mathbf{Q})=\mathbf{D}^{-1}\,\mathbf{M}$, then

$\boldsymbol{\nu}^{(1)}=-\int_{0}^{t}ds\,s^{\alpha-1}\,E_{\alpha,\alpha}(\mathbf{D}\,s^{\alpha})\mathbf{M}\,E_{\alpha}(\mathbf{D}\,(t-s)^{\alpha})\,{\bf
u}.$ (19)

The contribution corresponding to two events is

$\boldsymbol{\nu}^{(2)}=\int_{0}^{t}ds_{1}{\bf
f}(s_{1})\,\mathbf{D}^{-1}\mathbf{M}\,\int_{0}^{t-s_{1}}ds_{2}{\bf
f}(s_{2})\,\mathbf{D}^{-1}\mathbf{M}\,E_{\alpha}(\mathbf{D}\,(t-s_{1}-s_{2})^{\alpha})\,{\bf
u}.$ (20)

For arbitrary number of events $\eta$, the contribution can be readily
computed, and written compactly as follows

$\boldsymbol{\nu}^{(\eta)}=(-1)^{\eta}\,\mathbf{H}^{\eta}E_{\alpha}(\mathbf{D}\,t)\,{\bf
u},\quad\forall k\geq 0.$ (21)

Here $\mathbf{H}$ is an integral operator $\mathbf{H}$ defined as

$\mathbf{H}\psi(t):=\int_{0}^{t}ds\,{\bf
f}(s)\mathbf{D}^{-1}\mathbf{M}\psi(t-s).$ (22)

The solution $\mathbf{y}$ is finally obtained summing all those partial
contributions, and yields

$\mathbf{y}=\sum_{\eta=0}^{\infty}\boldsymbol{\nu}^{(\eta)}=(\mathbf{I}+\mathbf{H})^{-1}\,E_{\alpha}(\mathbf{D}\,t)\,{\bf
u},$ (23)

which can be expressed equivalently as an integral equation in the form

$\mathbf{y}=E_{\alpha}(\mathbf{D}\,t)\,{\bf u}-\int_{0}^{t}ds\,{\bf
f}(s)\mathbf{D}^{-1}\mathbf{M}\,{\mathbf{y}}(t-s).$ (24)

Taking the Laplace transform, and using that [13]

$\mathscr{L}\\{t^{\beta-1}E_{\alpha,\beta}(c\,t^{\alpha})\\}=s^{-\beta}(1-c\,s^{-\alpha})^{-1},$
(25)

we obtain

${\hat{\mathbf{y}}}(s)=\frac{s^{-1}}{\mathbf{I}-\mathbf{D}\,s^{-\alpha}}\,{\bf
u}+\frac{\mathbf{D}\,s^{-\alpha}}{\mathbf{I}-\mathbf{D}\,s^{-\alpha}}\,\mathbf{D}^{-1}\mathbf{M}\,{\hat{\mathbf{y}}},$
(26)

where ${\hat{\mathbf{y}}}(s)=\mathscr{L}\\{y(t)\\}$. The solution is

${\hat{\mathbf{y}}}=\frac{s^{-1}}{\mathbf{I}-\mathbf{D}\,s^{-\alpha}-\mathbf{D}\,s^{-\alpha}\,\mathbf{D}^{-1}\mathbf{M}}\,{\bf
u}.$ (27)

After applying the inverse Laplace transform and using again (25), we obtain

${\mathbf{y}}=E_{\alpha}((\mathbf{D}+\mathbf{M})t^{\alpha})\,{\bf
u}=E_{\alpha}(\mathbf{A}\,t^{\alpha})\,{\bf u},$ (28)

which concludes the proof. ∎

###  3.2 Implementation

Algorithm 1 Calculates the $i$-th entry of
$\mathbf{y}=E_{\alpha}(\mathbf{A}\,t^{\alpha})\,\mathbf{u}$ based on Theorem
1. $N_{p}$ indicates the number of random paths.

1:function MCML_Single($\mathbf{A}$, $\mathbf{u}$, $\alpha$, $i$, $N_{p}$,
$t$)

2: $y_{i}=0$

3: for each random path do

4: $X_{0}=i;\omega=1$; $k=0$

5: Generate $Z_{\alpha}\sim ML_{\alpha}(d_{X_{0}})$

6: $\tau=Z_{\alpha}$

7: while $\tau<t$ do

8: $X_{k+1}=$ SelectNextState($X_{k}$, $\mathbf{Q}$)

9: $\omega=\omega\times r_{X_{k}}g_{X_{k},X_{k+1}}$

10: Generate $Z_{\alpha}\sim ML_{\alpha}(d_{X_{k}})$

11: $\tau=\tau+Z_{\alpha}$

12: $k=k+1$

13: end while

14: $y_{i}=y_{i}+\omega\,u_{X_{k}}$

15: end for

16: $y_{i}=y_{i}/N_{p}$

17: return $y_{i}$

18:end function

According to Theorem 1, the value of the $i$-th entry of the solution vector
$\mathbf{y}$ can be estimated through the simulation of the stochastic process
$X_{t}$, which consists in generating random paths from the Markov chain $Z$
and then computing the realization of a random variable $\omega$ over these
paths. Each random path starts at state $X_{0}$ and time $\tau=0$ and then
jumps from one state to the next until it reaches the time $\tau=t$. The next
state is always selected at random based on the transition matrix $\mathbf{Q}$
and the sojourn time in each state follows an ML distribution. For a sequence
of states $X_{0},X_{1},\ldots,X_{\eta}$, the value of $\omega$ associated with
this random path can be calculated as

$\omega=\left(\prod_{j=1}^{\eta}r_{X_{j-1}}\,g_{X_{j-1},X_{j}}\right)u_{X_{\eta}}$

The procedure for simulating the process $X_{t}$ is described in Algorithm 1.
In order for the computation to be practical, Algorithm 1 can only generate a
finite number $N_{p}$ of random paths and approximate the value of $y_{i}$ as

$y_{i}\approx\frac{1}{N_{p}}\sum_{j=1}^{N_{p}}{\omega_{j}}$ (29)

where $\omega_{j}$ is the realization of $\omega$ over the random path $j$.
Naturally, replacing the expected value in (10) by the arithmetic mean over a
finite sample size $N_{p}$ in (29) leads to a numerical error $\epsilon$ of
order of $O(N_{p}^{-1/2})$. In fact, it is well-known that the arithmetic mean
provides the best unbiased estimator for the expected value, and, for large
number $N_{p}$ of random paths, the error $\epsilon$ is approximately a random
variable distributed according to a normal distribution with standard
deviation $\epsilon\approx\sigma N_{p}^{-1/2}$ [32]. Here $\sigma$ denotes the
standard deviation of the random variable $\omega$. Note that there is no
other source of error in Algorithm 1.

The main limitation of Algorithm 1 is that it can only estimate one entry of
$\mathbf{y}$ at a time, requiring a separated set of random paths to calculate
the solution in each point in the domain. Therefore, it is better suited for
estimating the solution locally (i.e., for a specific set of number of points)
than the complete solution. A more efficient method for computing the solution
for the entire domain is presented next.

###  3.3 Computing the full solution vector

Algorithm 2 Calculates
$\mathbf{y}=E_{\alpha}(\mathbf{A}\,t^{\alpha})\,\mathbf{u}$ based on its full
probabilistic representation. $N_{p}$ indicate the number of random paths.

1:function MCML_Full($\mathbf{A}$, $\mathbf{u}$, $\alpha$, $N_{p}$, $t$)

2: $\mathbf{y}=0$

3:
$\mathbf{p}^{(0)}=\left\\{p_{i}^{(0)}=\cfrac{|u_{i}|}{\|\mathbf{u}\|_{1}}\right\\}$

4: for each random path do

5: $k=0$

6: $X_{0}=$ SelectInitialState($\mathbf{p}^{(0)}$)

7: $\omega=\cfrac{u_{X_{0}}}{|u_{X_{0}}|}\|\mathbf{u}\|_{1}$

8: Generate $Z_{\alpha}\sim ML_{\alpha}(d_{X_{0}})$

9: $\tau=Z_{\alpha}$

10: while $\tau<t$ do

11: $X_{k+1}=$ SelectNextState($X_{k}$, $\mathbf{Q}^{\intercal}$)

12: $\omega=\omega\times r_{X_{k}}g_{X_{k},X_{k+1}}$

13: Generate $Z_{\alpha}\sim ML_{\alpha}(d_{X_{k}})$

14: $\tau=\tau+Z_{\alpha}$

15: $k=k+1$

16: end while

17: $y_{X_{0}}=y_{X_{0}}+\omega$

18: end for

19: $\mathbf{y}=\mathbf{y}/N_{p}$

20: return $\mathbf{y}$

21:end function

Theorem 1 can be conveniently modified to represent the complete solution
vector $\mathbf{y}$ rather than a single entry. The procedure is identical to
the case of the matrix exponential (namely $\alpha=\beta=1$), which is
described fully in [30]. For completeness, we present here the main details.
The representation requires generating random paths, which now start at a
randomly chosen state $j$ and time $t=0$ according to a suitable distribution
probability $p_{j}^{(0)}=\mathds{P}(X_{0}=j)$, ending at state $i$ at time
$t$. Those paths evolve in time changing states governed now by
$\mathbf{Q}^{\intercal}$ as the corresponding transition matrix of the Markov
chain. This can be explained resorting to the Bayes’ theorem. In fact, it
holds that

$\displaystyle p_{ji}$ $\displaystyle=\mathds{P}(X_{t}=i|X_{0}=j)$
$\displaystyle=\mathds{P}(X_{t}=j|X_{0}=i)\frac{\mathds{P}(X_{t}=i)}{\mathds{P}(X_{0}=j)}$
$\displaystyle=p_{ij}\frac{p_{i}}{p_{j}^{(0)}}$
$\displaystyle=(\mathbf{P}^{\intercal})_{ji}\frac{p_{i}}{p_{j}^{(0)}}$

where $p_{i}=\mathds{P}(X_{t}=i)$. Note that the matrix probability
$(\mathbf{P}^{\intercal})$ corresponds to a Markov chain governed now by the
transition matrix $\mathbf{Q}^{\intercal}$. As it was explained in [30], the
distribution function $p_{j}^{(0)}$ can be arbitrarily chosen. However, the
choice may have a direct impact on the variance, and, in turn, on the
performance of the algorithm. For simplicity, we use here the more reasonable
choice corresponding to $p_{j}^{(0)}$ proportional to $|u_{j}|$, since this
resembles the well known importance sampling method for variance reduction,
where the sampling is done according to the importance of the data.

Algorithm 2 describes the stochastic method based on this alternative
representation. The main advantage of this method is the ability to estimate
the entire solution vector with a single set of random paths, naturally
distributing the contribution of each random path over all points in the
domain. Due to this property, Algorithm 2 is often more efficient than
Algorithm 1 for estimating the entire solution.

###  3.4 Generating random numbers from a Mittag-Leffler distribution

For every state in the random path, Algorithms 1 and 2 generate a random
number $Z_{\alpha}$ from the ML probability distribution to determine the time
spent in this state. The value of $Z_{\alpha}$ is then later added to the
total time $\tau$ of the random path. The random path ends when $\tau\geq t$.

Considering that an ML random variable follows the Geometric Stable Laws [33,
34], it can be represented as a combination of an exponential and a stable
random variable [35, 36]. Using this representation, Kozubowski [36] deducted
the following formula for generating $Z_{\alpha}\sim ML_{\alpha}(\gamma)$ for
$0~{}<~{}\alpha~{}\leq~{}1$:

$Z_{\alpha}=-|\gamma|^{-\frac{1}{\alpha}}\,\ln{(U)}\,[\sin(\alpha\pi)\cot(\alpha\pi
V)-\cos(\alpha\pi)]^{\frac{1}{\alpha}}$ (30)

where $U,V$ are random numbers from the uniform distribution over the $[0,1]$
interval and $\gamma$ is the rate parameter. For $\alpha=1$, (30) reduces to
the inversion formula for the exponential distribution:
$Z_{1}=-|\gamma|^{-1}\ln{U}$. Note that for $\alpha>1$, the ML function is not
monotonic and, thus, cannot be considered a probability distribution [37]. It
is worth mentioning that the rejection sampling technique [38, 39, 40] can
also be used for generating the ML random numbers, however it requires the
construction of a pointwise representation directly from (6), which can be
quite expensive due to the slow convergence of the series [41].

##  4 Numerical Results

In this section, we evaluate the performance of our Monte Carlo algorithm by
solving time-fractional partial differential equations (fPDEs) through the
method of lines. The method of lines [42, 43] consists in discretizing the
spatial variables of the initial problem to obtain a system of coupled
fractional ordinary differential equations with time as the independent
variable. The system can then be solved by computing the action of the Mittag-
Leffler function over the discretized initial values. Using this method, we
solve the Dirichlet boundary-value problem for both a 2D time-fractional
diffusion equation and a 3D time-fractional convection-diffusion equation.

###  4.1 Setup

All simulations regarding the shared-memory architecture were executed on a
commodity workstation with an AMD Ryzen 5800X 8C @4.7GHz and 32GB of RAM,
running Arch Linux. The MCML_Full (Algorithm 2) was implemented in C++ with
OpenMP and uses PCG64DXSM [44] as its random number generator. PCG64DXSM is an
improved version of the default generator of the popular NumPy [45] module for
Python. The code was compiled with the Clang/LLVM v14.0.0 with the -O3 and
-march=znver3 flags. We will refer to the Monte Carlo method that uses dense
and sparse matrices as mc_dense and mc_sparse, respectively.

To the best of our knowledge, the only freely available code capable of
computing the Mittag-Leffler function for matrices was proposed by Garrappa
and Popolizio [19] and it is written in MATLAB 2021a. This code is based on
the Schur-Parlett algorithm [20, 27] that was implemented in MATLAB as the
funm command. The Schur-Parlett algorithm consists in decomposing the matrix
$\mathbf{A}$ as $\mathbf{V}\mathbf{T}\mathbf{V}^{*}$ (Schur decomposition
[46]) and then computing the Mittag-Leffler function on the upper triangular
$\mathbf{T}$ matrix using the Taylor series and the Parlett recurrence. The
derivatives for the Taylor series are calculated through a combination of
series expansion, numerical inversion of the Laplace transform and summation
formula. We will refer to Garrappa and Popolizio’s code as matlab.

Note that the MATLAB backend is written in C/C++ and calls the Intel Math
Kernel Library [47] for many matrix and vector operations, which exploits very
efficiently the hardware resources of modern CPUs, including multi-threading
and SIMD units. It is worth remarking that we only compare our code against a
MATLAB implementation due to the lack of any parallel code freely available in
C/C++.

Nevertheless, there are some noticeable differences between the algorithms.
First of all, matlab only supports dense matrices and has a fixed precision
(IEEE 754 double-precision standard), while ours supports both dense and
sparse matrices as well as an user-defined precision. Moreover, a Monte Carlo
method is not only fully parallelizable but also allows the computation of
single entries of the vector solution, which is not possible with the
classical deterministic algorithm. Due to the random nature of the Monte Carlo
algorithms, all results reported in this section correspond to the mean values
obtained after several runs. With the exception of the variance analysis, the
program calculates the solution for all points in domain in all simulations.

###  4.2 Time-fractional diffusion equation

The first example we consider consists in solving the 2D time-fractional
diffusion equation,

$D^{\alpha}_{t}u(\mathbf{x},t)=\nabla^{2}u(\mathbf{x},t),$ (31)

in a domain $\Omega=[-\mu,\mu]^{2}$, with time $t>0$, a space variable
$\mathbf{x}~{}=~{}(x,y)~{}\in~{}\mathbb{R}^{2}$, an initial condition
$u(\mathbf{x},0)=u_{0}(\mathbf{x})$ and a Dirichlet boundary condition
$u(\mathbf{x},t)\lvert_{\partial\Omega}\ =0$. Here $D^{\alpha}_{t}$ is the
Caputo’s fractional derivative of order $\alpha$ with $0<\alpha\leq 1$.
Considering a discrete mesh with $m$ cells in each dimension, such that
$\Delta x=\Delta y=2\mu/m$, and the standard 5-point stencil finite difference
approximation, the approximated solution $\hat{u}(\mathbf{x},t)$ for (31) can
be written as

$\hat{u}(\mathbf{x},t)=E_{\alpha}\left(\frac{m^{2}}{4\mu^{2}}\,\hat{\mathbf{L}}\,t^{\alpha}\right)\,\hat{u}_{0}(\mathbf{x})=E_{\alpha}(\mathbf{A}\,t^{\alpha})\,\hat{u}_{0}(\mathbf{x})$
(32)

where $\hat{\mathbf{L}}$ denotes the discrete Laplacian operator. If
$\mathbf{I}$ is the $m\times m$ identity matrix and $\mathbf{C}$ is an
$m\times m$ matrix defined as follows, then matrix $\hat{\mathbf{L}}$ is given
as

$\hat{\mathbf{L}}=\begin{bmatrix}~{}\mathbf{C}&~{}\mathbf{I}&~{}0&~{}0&\cdots&~{}0\\\
~{}\mathbf{I}&~{}\mathbf{C}&~{}\mathbf{I}&~{}0&\cdots&~{}0\\\
~{}0&~{}\mathbf{I}&~{}\mathbf{C}&~{}\mathbf{I}&\cdots&~{}0\\\
\vdots&\ddots&\ddots&\ddots&\ddots&\vdots\\\
~{}0&\cdots&~{}0&~{}\mathbf{I}&~{}\mathbf{C}&~{}\mathbf{I}\\\
~{}0&\cdots&\cdots&~{}0&~{}\mathbf{I}&~{}\mathbf{C}\end{bmatrix}\qquad\mathbf{C}=\begin{bmatrix}-4&~{}1&~{}0&~{}0&\cdots&~{}0\\\
~{}1&-4&~{}1&~{}0&\cdots&~{}0\\\ ~{}0&~{}1&-4&~{}1&\cdots&~{}0\\\
\vdots&\ddots&\ddots&\ddots&\ddots&\vdots\\\ ~{}0&\cdots&~{}0&~{}1&-4&~{}1\\\
~{}0&\cdots&\cdots&~{}0&~{}1&-4\end{bmatrix}$

This example is particularly suited for analyzing the numerical errors since
the eigenvalues of the matrix $\mathbf{A}$ can be calculated analytically [48]
and, consequently, the corresponding solution of the equation. In fact, for a
point $\mathbf{x}=(x,y)$ in the discrete domain, the corresponding eigenvalue
$\lambda(\mathbf{x})$ and eigenvector $\mathbf{v}(\mathbf{x})$ are known to be
equal to

$\displaystyle\lambda(\mathbf{x})$
$\displaystyle=\frac{m^{2}}{4\mu^{2}}\left[2\cos\left(\frac{x\pi}{m+1}\right)+2\cos\left(\frac{y\pi}{m+1}\right)-4\right]$
$\displaystyle\mathbf{v}(\mathbf{x})$
$\displaystyle=\left\\{v_{(i,j)}=\sin\left(\frac{ix\pi}{m+1}\right)\sin\left(\frac{jy\pi}{m+1}\right)\text{
for }i,j=1,2,\ldots,m\right\\}$

After organizing the eigenvalues and eigenvectors of $\mathbf{A}$ into the
matrices $\mathbf{D}$ and $\mathbf{V}$, respectively, the solution for the 2D
diffusion problem can be calculated as

$\hat{u}(\mathbf{x},t)=\mathbf{V}\,E_{\alpha}(\mathbf{D}\,t^{\alpha})\,\mathbf{V}^{\intercal}\,\hat{u}_{0}(\mathbf{x})$
(33)

Moreover, we can determine the stiffness ratio $r=\lambda_{max}/\lambda_{min}$
[49] for this problem, which is given asymptotically by $r\sim 4m^{2}/\pi^{2}$
for large values of $m$. In particular, the diffusion problem can be
considered stiff for the range of values of $m$ used in this section and,
thus, is a suitable example for testing the performance of our method in
solving this class of problems.

For all the results shown in this section, we chose the domain to be
$\Omega=[-1,1]^{2}$ and the initial condition
$\hat{u}_{0}(\mathbf{x})=c\,m^{2}\,\delta(\mathbf{x}-\mathbf{x}_{c})$, which
consists on a discrete impulse located at the center of the computational
grid, $\mathbf{x}_{c}~{}=~{}(m/2,m/2)$, and strength $c=1/4096$. Since the
matlab code only works with dense matrices, it requires around $25$GB of
memory when calculating the solution for a discrete mesh with $m=160$. In
comparison, mc_dense uses around $10$GB of memory for the same mesh, while
mc_sparse, only a couple of megabytes. For this reason, the maximum number of
cells $m$ in any experiment will be equal to $160$. Regarding the time
variable $t$, we choose $t=0.1$ where it is possible to clearly distinguish
the diffusion process for different values of $\alpha$.

Figure 1: Maximum absolute error as a function of the number of random paths
when solving the 2D diffusion equation for $t=0.1$ and $m=80$.

Figure 2: A histogram of the maximum absolute error for $8000$ runs, each one
starting with a different random seed. In each run, we solved the 2D diffusion
equation with $t=0.1$, $m=80$ and $N_{p}=2.5\times 10^{5}$.

Fig. 2 shows the maximum absolute error as function of the number of random
paths $N_{p}$. The error has been computed using the analytical solution in
(33). Recall from Section 3.2 that the numerical error of the Monte Carlo
method approximates a Gaussian random variable with standard deviation
determined by $\sigma\,N_{p}^{-1/2}$. This relation is confirmed numerically
by the trend line in the graph, which has a slope of approximately $-1/2$ in
the logarithmic scale. Repeating the Monte Carlo simulation several times with
different initial seeds, we observe in Fig. 2 that the numerical errors are
indeed distributed according to a normal distribution.

Figure 3: Variance $\sigma^{2}$ of the Monte Carlo algorithm as a function of
the matrix size $N=m^{2}$ and $\alpha$ when solving the 2D diffusion equation
at a single point located at the centre of the mesh for $t=0.1$. The number of
random paths was kept fixed to $10^{10}$.

An estimation of the variance can be obtained for this simple initial
condition, and in the limiting case when $m\gg 1$. Note that for this case
there are only two possible outcomes for the random variable $\omega$ in (29),
$0$ or $c\,m^{2}$, with probabilities $1-\xi$ and $\xi$, respectively. The
outcome $c\,m^{2}$ is obtained when the ending state of the random path
coincides with the node of the computational mesh corresponding to the point
$\mathbf{x}_{c}$, being $0$ otherwise. For simplicity, let us assume that we
choose the point where the solution is computed to be the same
$\mathbf{x}_{c}$. Then, the outcome will be different from $0$ only when,
after a random number of jumps $n$ in a random path, the final state is the
same as the initial one. Estimating the probability of $n$ is straightforward,
since this problem is equivalent to the problem of computing the probability
of returning to the initial point after $n$ steps for a simple symmetric 2D
random walk [50], and is given by

$h_{n}=\left[\binom{n}{n/2}\frac{1}{2n}\right]^{2},\quad n=0,2,4,\ldots.$ (34)

Asymptotically when $n\gg 1$, $h_{n}\sim 2/\pi(n+1)$. Concerning the number of
transitions, this is a random variable governed by a probability distribution
$P_{n}$, and depends on the number of events occurring during the time
interval $[0,t]$. Moreover, since the sojourn time in every state is Mittag-
Leffler distributed, then the probability can be estimated, and it turns out
to be the so-called fractional Poisson distribution [51], where the mean
number of events $\bar{n}$ is

$\bar{n}=\frac{m^{2}t^{\alpha}}{\Gamma{(\alpha+1)}},$ (35)

and the probability of no events occurring during this time interval is
$P_{0}=E_{\alpha}(m^{2}\,t^{\alpha})$. When $m\to\infty$, a simple piecewise
constant function can be used as an approximation, and yields

$P_{n}\approx\begin{cases}P_{0},&n\leq\frac{1}{P_{0}}\\\
0,&n>\frac{1}{P_{0}}.\end{cases}$ (36)

Then, it follows that the probability $\xi$ is given by

$\xi=\sum_{n=0}^{\infty}h_{2n}P_{2n}\approx
P_{0}+P_{0}\sum_{n=1}^{1/P_{0}}h_{2n}.$ (37)

Using the known difference equation
$\sum_{k=0}^{n-1}(k+x)^{-1}=\psi(x+n)-\psi(x)$, where $\psi(x)$ is the digamma
function, we obtain

$\xi\approx
P_{0}+\frac{P_{0}}{\pi}\left[\psi\left(\frac{3}{2}+\frac{1}{P_{0}}\right)-\psi\left(\frac{5}{2}\right)\right].$
(38)

When $m\to\infty$, $P_{0}\sim[m^{2}\,t^{\alpha}\,\Gamma(1-\alpha)]^{-1}$, and
then using the asymptotic expansion properties of the digamma function it
holds that $\xi\sim[m^{2}\,t^{\alpha}\,\Gamma(1~{}-~{}\alpha)]^{-1}$. Thus,
asymptotically an estimation of the variance $\sigma^{2}$ is given by

$\sigma^{2}\sim(c\,m^{2})^{2}\xi=c^{2}\,m^{2}\frac{t^{-\alpha}}{\Gamma(1-\alpha)}=c^{2}\,N\frac{t^{-\alpha}}{\Gamma(1-\alpha)}.$
(39)

This relation can be seen in Fig. 3. For a fixed $\alpha$ and $t$, the
variance $\sigma^{2}$ scales linearly with the matrix size $N=m^{2}$.
Likewise, for a fixed $N$, $t$, the variance $\sigma^{2}$ is proportional to
$[t^{\alpha}\,\Gamma(1~{}-~{}\alpha)]^{-1}$. Note that the points in the graph
do not perfectly match the curve from (39) due to the low number of cells $m$
used in the simulation. Still, it provides the best approximation for the
variance $\sigma^{2}$ under these conditions.

Figure 4: Elapsed time for solving the 2D diffusion equation as a function of
the time $t$ and the parameter $\alpha$ for $m=80$. The number of random paths
was kept fixed at $10^{5}$.

Figure 5: Elapsed time for solving the 2D diffusion equation as a function of
the matrix size $N=m^{2}$ and the parameter $\alpha$ for $t~{}=~{}0.1$. The
number of random paths was kept fixed at $10^{6}$.

Concerning the computational cost of the Monte Carlo algorithm, this depends
on the number of random paths and the time for generating each path. Even
though this time is random, we can readily assume to be proportional to the
mean time of events in (35). Therefore, the computational cost of the
algorithm is of the order of $O(N_{p}\,N\,t^{\alpha})$, which is in good
agreement with the results shown in Fig. 5 and 5. For a fixed number of random
paths $N_{p}$ and matrix size $N$, the execution time is proportional to
$t^{\alpha}$ as shown in Fig. 5. Note that the scale is $log-log$ and the
slope of the curves coincides with $\alpha$. Similarly, when we kept fixed
both the time $t$ and number of random paths $N_{p}$, the execution time grows
linearly with the matrix size $N$ as shown in Fig. 5. It is worth mentioning
that exponential random numbers are much cheaper to produce than ML random
numbers, resulting in significantly lower execution times when $\alpha=1$.

To compare the performance and scaling of our methods against a classical
approach, we measure the serial and parallel execution times for different
matrix sizes $N=m^{2}$ as shown in Fig. LABEL:fig:poisson_scale. In terms of
execution time, the Monte Carlo algorithm is several times faster than the
matlab code, especially when fully exploiting the sparse nature of the
problem, while maintaining a reasonable precision (i.e., the maximum error in
any simulation is below $3\times 10^{-4}$). When using 8 threads and
considering a discrete mesh with $m=128$, both matlab and mc_dense shows a
speedup of $6.5$ (or an efficiency of $81.25\%$), while mc_sparse presents a
speedup of $7.5$ ($95\%$ efficiency).

Since matrix $\mathbf{A}$ is real and symmetric, its Schur decomposition can
be written as $\mathbf{A}=\mathbf{V}\mathbf{D}\mathbf{V}^{\intercal}$, where
$\mathbf{D}$ is a diagonal matrix containing the eigenvalues of $\mathbf{A}$,
and the columns of the matrix $\mathbf{V}$ the corresponding eigenvectors. In
this case, the matrix function $f(\mathbf{A})$ can be simply evaluated as
$f(\mathbf{A})=\mathbf{V}f(\mathbf{D})\mathbf{V}^{\intercal}$ [20]. MATLAB
uses a QR algorithm to compute the Schur decomposition, directly calling the
corresponding routine from the Intel Math Kernel Library (e.g., dsteqr for
real symmetric matrices). The QR algorithm has a well-known computational
complexity of order $O(N^{3})$.

In comparison, the Monte Carlo algorithm requires an order of $O(N^{2})$
operations to solve the diffusion equation, as shown in Fig.
LABEL:fig:poisson_scale. Considering that the variance $\sigma^{2}$ and,
consequently, the statistical error grow with the matrix size $N$, the number
of random paths must increase proportionally to the matrix size in order to
keep the same level of accuracy for the numerical solution.

###  4.3 Time-fractional convection-diffusion equation

The second example we discuss in this paper consists in solving a time-
fractional convection-diffusion equation:

$D^{\alpha}_{t}u(\mathbf{x},t)=c\nabla^{2}u(\mathbf{x},t)+\nu\nabla
u(\mathbf{x},t)$ (40)

with boundary and initial conditions, respectively

$u(\mathbf{x},t)\lvert_{\partial\Omega}\ =g(\mathbf{x},t)\qquad
u(\mathbf{x},0)=f(\mathbf{x})$ (41)

where $c$ is the diffusion coefficient and $\nu$ the velocity field. After
applying the standard Galerkin finite element method [52] to the discretized
nodes $x_{i}$ for $i=0,1,\dots,n$, we obtain the following system of
equations:

$\mathbf{B}D^{\alpha}_{t}\mathbf{u}=\mathbf{K}\mathbf{u}+\mathbf{h},\qquad\mathbf{u}(0)=\mathbf{u}_{0}$
(42)

where $\mathbf{u}=\\{u(x_{1},t),u(x_{2},t),\dots,u(x_{n},t)\\}$, $\mathbf{B}$
is the assembled mass matrix, $\mathbf{K}$ is the corresponding assembled
stiffness matrix and $\mathbf{h}$ is the load vector. The diffusion
coefficient $c$, the velocity field $\nu$, and the boundary conditions are
already included in these matrices and vector. To simplify the computation,
the mass matrix was lumped [52], transforming the mass matrix into a diagonal
matrix.

The solution for the inhomogeneous systems of fPDEs in (42) can be written in
terms of the Mittag-Leffler function as follows

$\mathbf{u}(\mathbf{x},t)=E_{\alpha}(\mathbf{A}\,t^{\alpha})\
\mathbf{u}_{0}(\mathbf{x})+\int_{0}^{t}ds(t-s)^{\alpha-1}E_{\alpha}(\mathbf{A}\,(t-s)^{\alpha})\mathbf{v}$
(43)

with $\mathbf{A}=\mathbf{B}^{-1}\mathbf{K}$ and
$\mathbf{v}=\mathbf{B}^{-1}\mathbf{h}$. Note that the second term in the right
side depends on the boundary data, and in particular for numerical purpose it
will require to approximate a definite integral. This will entail another
source of error along with the aforementioned statistical error. Since the
goal of this manuscript is to focus solely on computing numerically the ML
function, and analyzing the associated statistical error, here we set a null
value for the Dirichlet boundary conditions $f(\mathbf{x})=0$, and hence
$\mathbf{h}=0$. Clearly the general case deserves further investigation, a
detailed analysis is left for a future manuscript.

Figure 6: Computational mesh of the domain where the fractional convection-
diffusion equation is solved.

Figure 7: Histogram of the eigenvalues of
$\mathbf{A}=\mathbf{B}^{-1}\mathbf{K}$ from (43). Here the scale parameter was
set to $s=30$, which results in a matrix of size $N~{}=~{}453$.

For the domain, we consider here a block of aluminium of $2s\times 2s\times
20s$ ($W\times H\times L$, in cm) with a scale parameter $s$, an isotropic
diffusion coefficient $c~{}=~{}9.4\times 10^{-5}m^{2}/s$ [53], a field
velocity $\nu=0$ and a Dirichlet boundary condition $u=0$ on the surface of
the block. The finite-element mesh in Fig. 7 was generated using the
scientific application tool COMSOL [54] considering a maximum and minimum
element size equal to $20$ cm and $0.1$ cm, respectively. Based on this mesh,
we generated the corresponding finite-element method (FEM) matrices and
vectors. For the initial conditions we use again a discrete impulse located at
the position $(s,s,10s)$. It is worth mentioning that the resulting matrix is
asymmetric, having entries entirely arbitrary in both sign and magnitude, and
in this sense this can be considered as a much more complex example than the
case analysed in the previous section.

Table 1: The execution time of the matlab code for computing the solution of the 3D convection-diffusion equation for different values of $t$ and $\alpha$ considering a scale parameter $s=30$ (the total number of nodes in the computational mesh was $453$). The parameters of funm were kept at their default values. The number in parenthesis indicates the order of the derivative of the largest block. | $\alpha=1$ | $\alpha=0.9$ | $\alpha=0.7$ | $\alpha=0.5$
---|---|---|---|---
$t=20$ | $6.43\ (19)$ | FAIL | FAIL | $919\ (21)$
$t=40$ | $5.81\ (22)$ | FAIL | FAIL | $965\ (24)$
$t=60$ | $5.55\ (24)$ | FAIL | FAIL | $967\ (24)$
$t=80$ | $6.15\ (27)$ | FAIL | FAIL | $971\ (25)$
$t=100$ | $6.78\ (30)$ | FAIL | FAIL | $982\ (27)$

Similarly to the previous example, we examine the eigenvalues of matrix
$\mathbf{A}$ to gain some insight about its numerical properties. In this
case, the eigenvalues were calculated through the command eigs of MATLAB,
which corresponds in essence to the Arnoldi method. Fig. 7 depicts the
histogram for the eigenvalues of the matrix $\mathbf{A}$. The eigenvalues are
not only small but are also highly clustered between $-0.02$ and $-0.05$,
which is a problem for the matlab code of [19]. On one hand, this leads to
very large blocks during the partitioning of the Schur form of $\mathbf{A}$ in
the Schur-Parlett algorithm [20, 19]. The evaluation of the Taylor series is
particularly expensive for large blocks since it requires order $O(m^{4})$
operations, where $m$ is the size of the block. On the other hand, the
numerical evaluation of high-order derivatives of the Mittag-Leffler function
often entails very large numbers even for small input values due to the
rapidly increasing factorial $(x)_{k}=x(x-1)\dots(x-k+1)$:

$\frac{d^{k}}{dz^{k}}E_{\alpha}(z)=\sum_{j=k}^{\infty}{\frac{(j)_{k}}{\Gamma(\alpha
j+1)}z^{j-k}},\quad k\in\mathbb{N}$ (44)

Indeed, MATLAB reports an infinite derivative when trying to compute the ML
function for $\alpha=0.7$ and the algorithm fails to converge to a solution
for $\alpha=0.9$ even after computing $250$ terms of the Taylor series. This
is shown in Table 1. The matlab algorithm was able to reach a solution for
$\alpha=0.5$ and $\alpha=1$.

In the Schur-Parlett algorithm, the Schur form of $\mathbf{A}$ are divided
into blocks according to its eigenvalues. In particular, if the absolute
difference between two eigenvalues is less than a tolerance $\delta$, they
will be assigned to the same block. With the default tolerance $\delta=0.1$,
almost all eigenvalues of $\mathbf{A}$ will be grouped into a single, large
block. Even after decreasing the tolerance to $\delta=0.01$, most eigenvalues
are still concentrated into a single block. Note that the value of $\delta$
cannot be too small, otherwise the two distinct blocks may contain close
eigenvalues, causing the Parlett’s recurrence to break down due to numerical
cancellations [20, 55]. It is worth mentioning that the derivatives of the ML
function are calculated serially and, thus, the matlab code shows very little
speedup when using multiple threads.

Concerning the computational cost of the Monte Carlo method for solving this
problem, recall that this depends on the number of random paths and the mean
number of events occurred during the prescribed time interval. The former is
chosen according to the desired level of accuracy for the numerical solution.
As it was already explained for the previous example, the numerical error
depends not only on the number of random paths but also on the variance, which
depends on the specific entries of the matrix $\mathbf{A}$ and the input
vector $\mathbf{u}$. The more heterogeneous they are, the larger the variance
will be, and consequently more random paths will be required for attaining the
prescribed accuracy. Regarding the mean number of events, this depends on the
specific values of the diagonal entries of $\mathbf{A}$, being larger for
larger entries.

In contrast to the matlab code which depends strongly on the distribution of
eigenvalues, the Monte Carlo method do not show explicit dependency on them.
For this reason, the stochastic method was able to compute the solution
regardless of the values of $t$ or $\alpha$. The execution time of mc_dense
and mc_sparse for different values of $t$ and $\alpha$ are shown in Fig.
LABEL:fig:fem_time_scale. Regarding the parallel performance, mc_dense and
mc_sparse usually achieve a speedup between $5$ and $6$ when using $8$
threads.

Figure 8: Parallel execution time of the Monte Carlo method as a function of
the number of nodes $N$ in the computational mesh when computing the solution
of the 3D convection-diffusion equation for $t=40$. The number of random paths
was kept fixed at $10^{9}$, resulting in an error $\epsilon\approx 2\times
10^{-5}$.

Figures LABEL:fig:fem_err and LABEL:fig:fem_scale_err show the maximum
absolute error as a function of the number of random paths and the matrix size
$N$, respectively. Note that the numerical error of the method $\epsilon$
decreases with the square root of the number of random paths as theoretically
expected.

In the previous example, the domain was fixed to $\Omega~{}=~{}[-1,1]^{2}$,
and the matrix size was enlarged increasing the number of grid points inside
the domain. As a result the magnitude of the diagonal entries of the matrix
increases accordingly, and so do the mean number of events and the variance.
We now consider a case where the distance between two nodes in the finite-
element mesh was kept constant when the domain is conveniently scaled in order
to increase the matrix size. This leads to a matrix with diagonal entries
independent of the matrix size, and therefore, the execution time and error of
the Monte Carlo algorithm only depend on the number of random paths. This is
shown in Fig. LABEL:fig:fem_scale_err and 8. Note however that the algorithm
needs to perform a binary search for selecting the next state of the random
path, which causes the execution time to slightly grow as the matrix size
increases. This has a negligible effect on the sparse implementation due to a
low number of nonzero entries per matrix row.

###  4.4 Distributed-Memory Performance

Figure 9: Weak scaling of mc_sparse in the Karolina supercomputer. This test
solves a 2D diffusion equation for $t=0.1$, $\mu=1$ and $m=1024$
($N=m^{2}=1,048,576$). The number of random paths was chosen based on the
number of cores, beginning with $2\times 10^{6}$.

Figure 10: Weak scaling of mc_sparse in the Karolina supercomputer. This test
solves a 3D convection-diffusion equation for $t=40$ considering a scalar
parameter $s=300$ (the total number of nodes in the computational mesh was
$753,675$). The number of random paths was chosen based on the number of
cores, beginning with $10^{11}$.

The distributed-memory tests were carried out at the Karolina Supercomputer
located in IT4Innovations National Supercomputing Centre. Each computational
node in the cluster contains two AMD 7H12 64C @2.6GHz CPUs and 256GB of RAM,
running CentOS 7. For all tests, two MPI processes were launched on each node
(one per processor) with as many threads as the number of physical cores (64).
Since there is no freely available code capable of computing the action of a
Mittag-Leffler function over a vector that is suitable for distributed-memory
systems, this section only contains the results for the Monte Carlo method
(mc_sparse). The C++ code was compiled with GCC 12.1.0 and OpenMPI 4.1.4 with
the -O3 and -march=native flags.

Figures 10 and 10 show the weak scaling of mc_sparse between 1 and 128 nodes
($16,384$ cores), considering a single computational node as the baseline.
While the scaling of the domain size is typically used to increase the
workload, for the specific FEM problem in (42), as it was shown in the
previous section, this does not affect the computational time, and therefore
does not provide any useful insight about the scalability of the algorithm.
Instead, we increased the workload by scaling the number of random paths
according with the number of nodes used. In doing so the statistical error
reduces, improving the accuracy of the computed solution.

Since all random paths in the Monte Carlo method are independent, the task of
generating them can be done in parallel and the parallel code only needs to
communicate at the end of the simulation to combine the results. This is the
reason that the mc_sparse shows perfect scalability in both examples.

##  5 Conclusion

In this paper we propose a novel stochastic method for solving time fractional
partial differential equations. These equations are already being used for
modeling natural phenomenon subject to memory effects, and microscopically are
typically described by non-Markovian processes. Fractional equations are
capable of capturing such effects due to the inherent non-locality of the
operator. As a consequence, the classical numerical schemes often based on
time-stepping suffer from heavy memory storage requirements since the
numerical solution of the fPDEs at current time depends on all preceding time
instances. This can be even worse when dealing with high-dimensional problems,
degrading significantly the performance of the corresponding numerical
algorithms.

The main advantage of the proposed method rests on the fact that it allows for
the computation of the solution at single point of a given domain, or
equivalently a single entry of the corresponding vector solution. Moreover,
the numerical algorithm is not based on any time-stepping scheme, and
therefore the solution is obtained without the need of storing previous
results. Rather, the solution is computed through an expected value of a
functional of random processes which resembles the non-Markovian process found
in the microscopic description of the phenomenon, and hence exploits somehow
naturally the non-locality of the fractional operators.

Furthermore, since the method is based on Monte Carlo it inherits all the
known advantages from a computational point of view, such as the comparative
ease of implementation in parallel, fault-tolerance, and in general of being
well suited for heterogeneous architectures. In fact, our parallel
implementation was able to solve large-scale problems efficiently in both
shared-memory and distributed-memory systems, demonstrating the versatility
and scalability of the probabilistic method compared to the classical
numerical schemes.

## Acknowledgments

This work was supported by Fundação para a Ciência e a Tecnologia under grant
2022.11506.BD and projects UIDB/50021/2020, 2022.08838.PTDC and PTDC/CCI-
INF/6762/2020. JA was funded by Ministerio de Universidades and specifically
the requalification program of the Spanish University System 2021-2023 at the
Carlos III University.

## References

  * [1] R. Metzler, J. Klafter, The random walk’s guide to anomalous diffusion: A fractional dynamics approach, Physics Reports 339 (1) (2000) 1–77.  doi:10.1016/S0370-1573(00)00070-3.
  * [2] B. J. West, Fractional Calculus View of Complexity: Tomorrow’s Science, CRC Press, Boca Raton, FL, 2016.
  * [3] A. Lischke, G. Pang, M. Gulian, F. Song, C. Glusa, X. Zheng, Z. Mao, W. Cai, M. M. Meerschaert, M. Ainsworth, G. E. Karniadakis, What is the fractional Laplacian? A comparative review with new results, Journal of Computational Physics 404 (2020) 109009.  doi:10.1016/j.jcp.2019.109009.
  * [4] E. Estrada, Fractional diffusion on the human proteome as an alternative to the multi-organ damage of SARS-CoV-2, Chaos: An Interdisciplinary Journal of Nonlinear Science 30 (8) (2020) 081104.  doi:10.1063/5.0015626.
  * [5] R. Magin, Fractional Calculus in Bioengineering, Part 1, Critical Reviews™ in Biomedical Engineering 32 (1) (2004).  doi:10.1615/CritRevBiomedEng.v32.i1.10.
  * [6] R. Hilfer, Applications Of Fractional Calculus In Physics, World Scientific, Singapore, 2000.
  * [7] J. Sabatier, O. P. Agrawal, J. A. Tenreiro Machado (Eds.), Advances in Fractional Calculus: Theoretical Developments and Applications in Physics and Engineering, Springer, Dordrecht, 2007.
  * [8] F. Mainardi, Fractional Calculus And Waves In Linear Viscoelasticity: An Introduction To Mathematical Models, World Scientific, London, 2010.
  * [9] Y. Chen, I. Petras, D. Xue, Fractional order control - A tutorial, in: 2009 American Control Conference, 2009, pp. 1397–1411.  doi:10.1109/ACC.2009.5160719.
  * [10] F. Arrigo, F. Durastante, Mittag–Leffler Functions and their Applications in Network Science, SIAM Journal on Matrix Analysis and Applications 42 (4) (2021) 1581–1601.  doi:10.1137/21M1407276.
  * [11] J. H. Martínez, S. Romero, J. J. Ramasco, E. Estrada, The world-wide waste web, Nature Communications 13 (1) (2022) 1615.  doi:10.1038/s41467-022-28810-x.
  * [12] F. Mainardi, Why the Mittag-Leffler Function Can Be Considered the Queen Function of the Fractional Calculus?, Entropy 22 (12) (2020) 1359.  doi:10.3390/e22121359.
  * [13] R. Gorenflo, A. A. Kilbas, F. Mainardi, S. Rogosin, Mittag-Leffler Functions, Related Topics and Applications, Springer Monographs in Mathematics, Springer Berlin Heidelberg, Berlin, Heidelberg, 2020.  doi:10.1007/978-3-662-61550-8.
  * [14] R. Gorenflo, F. Mainardi, Fractional Calculus, in: A. Carpinteri, F. Mainardi (Eds.), Fractals and Fractional Calculus in Continuum Mechanics, International Centre for Mechanical Sciences, Springer, Vienna, 1997, pp. 223–276.
  * [15] R. Hilfer, H. J. Seybold, Computation of the generalized Mittag-Leffler function and its inverse in the complex plane, Integral Transforms and Special Functions 17 (9) (2006) 637–652.  doi:10.1080/10652460600725341.
  * [16] R. Garrappa, Numerical Evaluation of Two and Three Parameter Mittag-Leffler Functions, SIAM Journal on Numerical Analysis 53 (3) (2015) 1350–1369.  doi:10.1137/140971191.
  * [17] H. Seybold, R. Hilfer, Numerical Algorithm for Calculating the Generalized Mittag-Leffler Function, SIAM Journal on Numerical Analysis 47 (1) (2009) 69–88.  doi:10.1137/070700280.
  * [18] R. Gorenflo, J. Loutchko, Y. Luchko, Computation of the Mittag-Leffler function and its derivatives, Fractional Calculus & Applied Analysis (FCAA) 5 (2002) 491–518.
  * [19] R. Garrappa, M. Popolizio, Computing the matrix Mittag-Leffler function with applications to fractional calculus, Journal of Scientific Computing 77 (1) (2018) 129–153.  arXiv:1804.04883, doi:10.1007/s10915-018-0699-5.
  * [20] P. I. Davies, N. J. Higham, A Schur-Parlett Algorithm for Computing Matrix Functions, SIAM Journal On Matrix Analysis and Applications 25 (2) (2003) 464–485.  doi:10.1137/S0895479802410815.
  * [21] I. Moret, P. Novati, On the Convergence of Krylov Subspace Methods for Matrix Mittag–Leffler Functions, SIAM Journal on Numerical Analysis 49 (5) (2011) 2144–2164.  doi:10.1137/080738374.
  * [22] G. E. Forsythe, R. A. Leibler, Matrix Inversion by a Monte Carlo Method, Mathematical Tables and Other Aids to Computation 4 (31) (1950) 127–129.  arXiv:2002508, doi:10.2307/2002508.
  * [23] I. Dimov, Monte Carlo Methods for Applied Scientists, World Scientific, Singapore, 2008.
  * [24] I. Dimov, V. Alexandrov, A. Karaivanova, Parallel resolvent Monte Carlo algorithms for linear algebra problems, Mathematics and Computers in Simulation 55 (1-3) (2001) 25–35.  doi:10.1016/S0378-4754(00)00243-3.
  * [25] H. Ji, M. Mascagni, Y. Li, Convergence Analysis of Markov Chain Monte Carlo Linear Solvers Using Ulam-Von Neumann Algorithm, SIAM Journal on Numerical Analysis 51 (4) (2013) 2107–2122.  arXiv:42004066.
  * [26] M. Benzi, T. M. Evans, S. P. Hamilton, M. Lupo Pasini, S. R. Slattery, Analysis of Monte Carlo accelerated iterative methods for sparse linear systems, Numerical Linear Algebra with Applications 24 (3) (2017).  doi:10.1002/nla.2088.
  * [27] N. J. Higham, Functions of Matrices, Other Titles in Applied Mathematics, Society for Industrial and Applied Mathematics, Philadelphia, PA, 2008.  doi:10.1137/1.9780898717778.
  * [28] I. Dimov, S. Maire, J. M. Sellier, A new Walk on Equations Monte Carlo method for solving systems of linear algebraic equations, Applied Mathematical Modelling 39 (15) (2015) 4494–4510.  doi:10.1016/j.apm.2014.12.018.
  * [29] J. Acebrón, A Monte Carlo method for computing the action of a matrix exponential on a vector, Applied Mathematics and Computation 362 (2019) 124545\.  doi:10.1016/j.amc.2019.06.059.
  * [30] J. A. Acebrón, J. R. Herrero, J. Monteiro, A highly parallel algorithm for computing the action of a matrix exponential on a vector based on a multilevel Monte Carlo method, Computers & Mathematics with Applications 79 (12) (2020) 3495–3515.  doi:10.1016/j.camwa.2020.02.013.
  * [31] J. A. Nichols, B. I. Henry, C. N. Angstmann, Subdiffusive discrete time random walks via Monte Carlo and subordination, Journal of Computational Physics 372 (2018) 373–384.  doi:10.1016/j.jcp.2018.06.044.
  * [32] A. C. Berry, The Accuracy of the Gaussian Approximation to the Sum of Independent Variates, Transactions of the American Mathematical Society 49 (1) (1941) 122–136.  arXiv:1990053, doi:10.2307/1990053.
  * [33] T. J. Kozubowski, S. T. Rachev, Univariate Geometric Stable Laws, Journal of Computational Analysis and Applications 1 (2) (1999) 177–217.  doi:10.1023/A:1022629726024.
  * [34] T. J. Kozubowski, Fractional moment estimation of linnik and mittag-leffler parameters, Mathematical and Computer Modelling 34 (9) (2001) 1023–1035.  doi:10.1016/S0895-7177(01)00115-7.
  * [35] T. J. Kozubowski, Exponential Mixture Representation of Geometric Stable Distributions, Annals of the Institute of Statistical Mathematics 52 (2) (2000) 231–238.  doi:10.1023/A:1004157620644.
  * [36] T. J. Kozubowski, Computer simulation of geometric stable distributions, Journal of Computational and Applied Mathematics 116 (2) (2000) 221–229.  doi:10.1016/S0377-0427(99)00318-0.
  * [37] R. N. Pillai, On Mittag-Leffler functions and related distributions, Annals of the Institute of Statistical Mathematics 42 (1) (1990) 157–161.  doi:10.1007/BF00050786.
  * [38] L. Devroye, Non-Uniform Random Variate Generation, 1st Edition, Springer, New York, 1986.
  * [39] G. Marsaglia, W. W. Tsang, The Ziggurat Method for Generating Random Variables, Journal of Statistical Software 5 (8) (2000).  doi:10.18637/jss.v005.i08.
  * [40] H. Rubin, B. C. Johnson, Efficient generation of exponential and normal deviates, Journal of Statistical Computation and Simulation 76 (6) (2006) 509–518.  doi:10.1080/10629360500108004.
  * [41] D. Fulger, E. Scalas, G. Germano, Random numbers from the tails of probability distributions using the transformation method, Fractional Calculus and Applied Analysis 16 (2) (2013) 332–353.  doi:10.2478/s13540-013-0021-z.
  * [42] W. E. Schiesser, The Numerical Method of Lines: Integration of Partial Differential Equations, 1st Edition, Elsevier, San Diego, CA, 1991\.
  * [43] B. Guo, X. Pu, F. Huang, Fractional Partial Differential Equations And Their Numerical Solutions, World Scientific, Singapore, 2015.
  * [44] M. E. O’Neill, PCG: A Family of Simple Fast Space-Efficient Statistically Good Algorithms for Random Number Generation, Tech. Rep. HMC-CS-2014-0905, Harvey Mudd College, Claremont, CA (Sep. 2014).
  * [45] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, T. E. Oliphant, Array programming with NumPy, Nature 585 (7825) (2020) 357–362.  doi:10.1038/s41586-020-2649-2.
  * [46] G. H. Golub, C. F. Van Loan, Matrix Computations, 4th Edition, Johns Hopkins Studies in the Mathematical Sciences, The Johns Hopkins University Press, Baltimore, MD, 2013.
  * [47] E. Wang, Q. Zhang, B. Shen, G. Zhang, X. Lu, Q. Wu, Y. Wang, Intel Math Kernel Library, in: E. Wang, Q. Zhang, B. Shen, G. Zhang, X. Lu, Q. Wu, Y. Wang (Eds.), High-Performance Computing on the Intel® Xeon Phi™: How to Fully Exploit MIC Architectures, Springer International Publishing, Cham, 2014, pp. 167–188.
  * [48] G. Strang, Computational Science and Engineering, Wellesley-Cambridge Press, Wellesley, MA, 2007.
  * [49] J. D. Lambert, Computational Methods in Ordinary Differential Equations, Vol. Vol. 5, John Wiley & Sons Incorporated, New York, 1973\.
  * [50] O. C. Ibe, Elements of Random Walk and Diffusion Processes, John Wiley & Sons, 2013.
  * [51] N. Laskin, Fractional Poisson process, Communications in Nonlinear Science and Numerical Simulation 8 (3-4) (2003) 201–213.  doi:10.1016/S1007-5704(03)00037-6.
  * [52] O. C. Zienkiewicz, R. L. Taylor, J. Z. Zhu, The Finite Element Method: Its Basis and Fundamentals, seventh Edition, Elsevier, Kidlington, Oxfor, 2013.
  * [53] W. J. Parker, R. J. Jenkins, C. P. Butler, G. L. Abbott, Flash Method of Determining Thermal Diffusivity, Heat Capacity, and Thermal Conductivity, Journal of Applied Physics 32 (9) (1961) 1679–1684.  doi:10.1063/1.1728417.
  * [54] COMSOL: Multiphysics Software for Optimizing Designs, https://www.comsol.com/.
  * [55] N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd Edition, Other Titles in Applied Mathematics, Society for Industrial and Applied Mathematics, Philadelphia, PA, 2002.