# Stochastic mean-shift clustering

Itshak Lapidot  itshakl@afeka.ac.il$^1,2$ 1Afeka Tel-Aviv Academic College of
Engineering, ACLP, Israel
2Avignon University, LIA, France

###### Abstract

In this paper we presented a stochastic version mean-shift clustering
algorithm. In the stochastic version the data points ”climb” to the modes of
the distribution collectively, while in the deterministic mean-shift, each
datum ”climbs” individually, while all other data points remains in their
original coordinates. Stochastic version of the mean-shift clustering is
comparison with a standard (deterministic) mean-shift clustering on
synthesized 2- and 3-dimensional data distributed between several Gaussian
component. The comparison performed in terms of cluster purity and class data
purity. It was found the the stochastic mean-shift clustering outperformed in
most of the cases the deterministic mean-shift.

###### keywords:

Unsupervised learning , Mean-shift clustering , Stochastic mean-shift , non-
parametric PDF estimation.

##  1 Introduction

Mean-shift clustering algorithm is a non-parametric iterative clustering
algorithm. It estimates the probability density function of a random variable
[13]. The clustering algorithm is applied to a variety of areas, like
segmentation images, [23, 17], particularly medical and satellite images [15,
1, 27, 4], videos [26], and also applied to high dimensional data clustering
[21]. An adapted version of mean-shift clustering was applied to short
segments speaker clustering [19, 18, 20, 10]. This algorithm is deterministic
and in an iterative procedure estimates the multi-modal probability density
function (pdf) via the ”climbing” path of each datum to its mode in a multi-
modal distribution. All the data points that reached the same mode are grouped
to the same cluster. We named it ”deterministic” as the algorithm is
indifferent to the order of the data that are chosen to find their modes, as
at the end of the ”climbing”, each datum is placed back to its original place.
There is a selective version of the mean-shift that chooses only some of the
data points in order to speed up the clustering [22]. It is sub-optimal,
however, although a sub-set of the data points is chosen randomly; after the
end of the ”climbing” each datum is placed at the original place. A version of
collective ”climbing” for Gaussian kernel known as Gaussian blurring mean-
shift, [7]. In the Gaussian blurring mean-shift, the data points do not placed
back to their original places, but collectively ”climb” each time. This
version is also deterministic but each new move depends on the current
position of all the data points and not on the original their places.

In the current work we present a stochastic variant, where all the data points
”climb” together, and after each single shift of one datum, a new datum is
chosen randomly. The data are placed in their new positions and are not placed
back to their original places. We will compare its performance versus the
well-known deterministic mean-shift clustering. This will be done on several
multi-model 2- and 3-dimensional distributions of synthesized data.

In this work is assumed that the problem is well defined - each cluster should
represent a different Gaussian, while at the same time, each Gaussian must be
concentrated in a single cluster. As such, a subjective, external evaluation
is applied, i.e, a criterion that uses the true labels of the data. We applied
two similar, but different criteria. The first is a based on the
quality/purity of a cluster, [16] taking the value of the class with the
maximum intersection with a given cluster, i.e., the winner takes all. We add
in the similar way the quality/purity of the class, i.e. the maximal
intersection between the class with all clusters. The second criterion is more
smooth an represent the contribution of all the classes to each cluster and
vice versa. This criterion is commonly used to validate the speaker clustering
and diarization quality [2, 19, 18, 20, 25, 3, 12] and will be discussed in
subsection 2.3. First the average cluster purity ($ACP$), and the average
speaker purity ($ASP$) are calculated, and then the $K$ value is evaluated as
the geometric mean between the first two measures. In [14] this criterion used
for speaker clustering while compared the adapted versions of deterministic
and stochastic mean-shift algorithms. The stochastic mean-shift always
outperformed the deterministic mean-shift clustering.

The paper is organized as follows. In Section 2, clustering algorithms are
presented, both deterministic and stochastic mean-shifts. Clustering
evaluation criteria are also explained. Synthesized data are described is
described in Section 3. At Section 4, the experiments with synthesized data
are presented. The conclusions and the discussion are in Section 5.

##  2 Mean-shift algorithms

This section describes and evaluates the two mean-shift algorithms applied in
this study. The difference between the deterministic mean-shift and the
stochastic mean-shift is illustrated in Figure 1. All the algorithms are
variants of the mean-shift algorithm that is based on Euclidean distance [11]
(subsection 2.1, Figure 1a). In this case, each time only one datum is shifted
until convergence. Each time the neighbors of the same datum are found, and
the shift of the datum is performed to the mean of the neighbors. This
procedure is repeated until convergence. The final place is recorded and the
datum returned to the original place, and the same procedure repeated with the
next datum. This is done until all the data points are shifted. The stochastic
version is then described in subsection 2.2, Figure 1b. In this case, each
time another datum is randomly chosen. The neighbors are found, and the datum
is shifted. The data points are never placed back in the original place, but
participate in the finding for the neighborhood according to the shifted
places (current places). At each iteration another datum is shifted, until the
collective convergence of all the dataset.

Figure 1:  Comparison between (a) deterministic mean-shift and (b) stochastic
mean-shift.

###  2.1 Mean-shift algorithm

The mean-shift algorithm is a non-parametric iterative algorithm. It estimates
the probability density function (pdf) of a random variable [13]. The
algorithm is inspired by the Parzen window approach to non-parametric density
estimation. The algorithm does not require any prior knowledge regarding the
number of clusters, and assumes that the shape of each cluster is unimodal.
Dense regions in the data space correspond to local maxima or to the pdf
modes. As such, for each data point, in order to reach the local maximum of
the pdf, a gradient ascent on the estimated local density is performed until
convergence is reached. Each stationary point represents a mode of the density
function. Data points that are associated with the same stationary point are
assigned to the same cluster.

The gradient of the density function is required in order to find the above-
mentioned modes. Following the mathematical formulation in [22, 11, 13], the
mean shift vector ${m_{h}}(x)$ expression is derived according to Eq. (1).

${m_{h}}(x)=\frac{{\sum\limits_{{x_{i}}\in{S_{h}}\left(x\right)}{{x_{i}}}}}{\\#S_{h}\left(x\right)}-x$
(1)

where $x$ is the current position of the $d$-dimensional datum; $h$ is the
neighborhood (or the bandwidth) from which the gradient is estimated, $\\#$
donates the cardinality of the set of neighbor points, and
${S_{h}}\left(x\right)$ is the set of data points that are the neighbors of
datum $x$, ${S_{h}}(x)\equiv\left\\{{{x_{i}}:\left\|{x-{x_{i}}}\right\|\leq
h}\right\\}$.

Let ${{\cal X}}=\left\\{{{x_{j}}}\right\\}_{j=1}^{J}$ be $d$-dimensional data
to be clustered, then the mean-shift algorithm is as described in Algorithm 1.
For each $x_{j}$ from the set ${{\cal X}}$, an iterative procedure is applied:
find the neighbors of $x_{j}$; calculate the shift according to Eq. (1); shift
$x_{j}$; repeat the procedure with the shifted datum until convergence.
Convergence in this case is defined as a shift’s norm which is smaller than a
threshold $Th_{1}$. The last step is merging the shifted data points that have
a Euclidean distance smaller than $Th_{2}$ to the same cluster.

Algorithm 1 Mean-shift clustering algorithm

A set of vectors, ${{\cal X}}=\left\\{{{x_{j}}}\right\\}_{j=1}^{J}$
$\triangleright$ $x\in{\mathbb{R}^{d\times 1}}$

Neighborhood size $h$ $\triangleright$ $h\in{\mathbb{R}^{+}}$

A cluster shift threshold $Th_{1}$ $\triangleright$
$Th_{1}\in{\mathbb{R}^{+}}$

A cluster merging threshold $Th_{2}$ $\triangleright$
$Th_{2}\in{\mathbb{R}^{+}}$

for $j:=1$ to $J$ step $1$ do

Set ${\hat{x}_{j}}={x_{j}}$.

repeat

Find the subset $S_{h}\left(\hat{x}_{j}\right)$.

Calculate $m_{h}\left(\hat{x}_{j}\right)$, the shift of the vector
$\hat{x}_{j}$, using Eq. (1).

${\hat{x}_{j}}\leftarrow{\hat{x}_{j}}+{m_{h}}\left({{\hat{x}_{j}}}\right)$

until $\left|m_{h}\left(x\right)\right|>Th_{1}$

Cluster the shifted vectors $\hat{\cal
X}=\left\\{\hat{x}_{j}\right\\}_{j=1}^{J}$ such that the distance between $2$
shifted vectors will be less then $Th_{2}$.

Return: Cluster index of each vector.

###  2.2 Stochastic mean-shift algorithm

In this section, a stochastic version of the mean-shift is presented. Two main
differences exist between the deterministic and the stochastic algorithms:

  1. 1.

At each iteration, a new datum is taken randomly (with uniform probability),
with repetitions, from the pool, while in the deterministic case, the same
datum is iterated until convergence and then a new datum is taken. If the
number of iterations is significantly larger than the size of the dataset,
then with high probability, most (if not all) of the data points will be
chosen several times and as such will probably converge to local optimum. Data
points which will rarely be chosen will probably end as outlier clusters. In
the experimental sections of this paper, it will be shown that the stochastic
mean-shift frequently produces more clusters than the deterministic mean-
shift. This might be the reason for this phenomenon.

  2. 2.

The neighbors that are chosen for the shift are taken from the data points in
their current positions and not the original positions, as in the
deterministic mean-shift. In this way, all the data paints ”climb”
collectively to their modes. With time, when the data points are shifted, the
temporal density of the points should be greater at the modes, and as a
result, should shift each next datum chosen in the more accurate direction.

Algorithm 2 is the stochastic version of the mean-shift of Algorithm 1.

The term ”convergence” is present in both deterministic and stochastic
algorithms. While ”convergence” in the deterministic case is straightforward,
i.e., iterate the same datum until the shift is smaller than a pre-defined
threshold (or pre-defined number of iterations), in the stochastic case it
must be defined more carefully. The use in this paper of the term
”convergence” is: either the shift is smaller than a pre-defined threshold for
a pre-defined number of iterations (as at each iteration a different datum is
shifted); or there is a global pre-defined number of iterations for the entire
procedure. If the number of iterations is sufficiently large and the data
points are randomly taken according to a uniform distribution, then it is
probable that most of the data points are shifted a sufficient number of times
and it reach a place sufficiently close to treated as converged to a local
optimum.

Algorithm 2 Stochastic mean-shift clustering algorithm

A set of vectors, ${{\cal X}}=\left\\{{{x_{j}}}\right\\}_{j=1}^{J}$
$\triangleright$ $x\in{\mathbb{R}^{d\times 1}}$

Neighborhood size $h$ $\triangleright$ $h\in{\mathbb{R}^{+}}$

A cluster shift threshold $Th_{1}$ $\triangleright$
$Th_{1}\in{\mathbb{R}^{+}}$

A cluster merging threshold $Th_{2}$ $\triangleright$
$Th_{2}\in{\mathbb{R}^{+}}$

repeat

Randomly choose ${\hat{x}_{j}}={x_{j}}$.

Find the subset $S_{h}\left(\hat{x}_{j}\right)$.

Calculate $m_{h}\left(\hat{x}_{j}\right)$, the shift of the vector
$\hat{x}_{j}$, using Eq. (1).

${x_{j}}\leftarrow{\hat{x}_{j}}+{m_{h}}\left({{\hat{x}_{j}}}\right)$

until Convergence

Cluster the shifted vectors ${\cal X}=\left\\{x_{j}\right\\}_{j=1}^{J}$ such
that the distance between $2$ shifted vectors will be less then $Th_{2}$.

Return: Cluster index of each vector.

###  2.3 Clustering evaluation criteria

There are two approaches evaluating clustering performance. The first is
objective, internal evaluation, measuring clusters’ compactness versus how far
the clusters are one from the other. Parameters such as single-link, complete-
link and many others are part of this approach, [16, 5, 24, 8]. These criteria
do not assume any true labeling of the data, and as such it is task-
independent. The second approach is task-dependent, external evaluation, and
the clustering is compared with the ground-truth. In this work, external
evaluation is applied.

One common way to evaluate clustering performances versus the ground-truth is
using the purity of the obtained clusters. Evaluation criteria are formulated
according to Eq. (2). The notation is as follows:

  * •

$R$ \- Number of data attributions, different classes,

  * •

$Q$ \- Number of clusters,

  * •

$n_{qr}$ \- Total number of data points in cluster $q$ which are associated
with class $r$,

  * •

$n_{.r}$ \- Total number of data points which are associated with class $r$,

  * •

$n_{q.}$ \- Total number of data points in cluster $q$,

  * •

$N$ \- Total number of data points to be clustered,

  * •

$\mathcal{D}=\left\\{d_{r}\right\\}_{r=1}^{R}$ \- The partitioning of the
ground-truth data attribution,

  * •

$\mathcal{C}=\left\\{c_{q}\right\\}_{q=1}^{Q}$ \- The partitioning of the data
according to the clustering.

$\begin{array}[]{*{20}{c}}{Pur_{\mathcal{C}}=\frac{1}{N}\sum\limits_{q=1}^{Q}\underset{d_{r}\in\mathcal{D}}{\max}\,{\left|c_{q}\cap
d_{r}\right|}}\\\
{Pur_{\mathcal{D}}=\frac{1}{N}\sum\limits_{r=1}^{R}\underset{c_{q}\in\mathcal{C}}{\max}\,{\left|c_{q}\cap
d_{r}\right|}}\\\ {G=\sqrt{Pur_{\mathcal{C}}\cdot
Pur_{\mathcal{D}}}}\end{array}$ (2)

Equation 2 describes, in the upper expression, a well-known way to measure the
quality or purity of clustering, [16]. It measures, in each cluster, which
class is most prevalent and sums it over all the clusters. At the end it
normalizes it by the total number of data points to bound it in the range
$\left[0,1\right]$. However, this can lead to a trivial solution, when the
number of clusters is equal to the size of the data, $N$, and each cluster
contains exactly one datum. For this reason also, we suggest also to measure
the purity of the data classes in the second expression. The logic behind this
is to prefer clustering which clusters all the data points of the same class
to the same cluster. This also can lead to a trivial solution of one cluster
only. The optimal is to combine the two measures and in the current work it is
done as a geometrical mean in the last expression.

Another criterion is the one that is used in speaker diarization and speaker
clustering applications. This criterion was presented in [2], and was used in
many cases [25, 3, 12] for speaker clustering. We applied it for short
segments speaker clustering [19, 18, 20] and for speaker clustering quality
estimation [9, 10]. The idea, similar to the first criterion, is to measure
the impact of the data distribution on each cluster and the impact of
distribution of the class between the clusters. Unlike the previous criterion,
where only the maximal intersection is taken into account, here all the
intersections are taken into account. Evaluation criteria are formulated
according to Eq. (3), First, the average cluster purity ($ACP$) and the
average class purity have to be calculated. In previous works, this was named
as average speaker purity ($ASP$) and we will follow this notation for
consistency. $ACP$ measures how well a cluster is limited, on average, to only
one speaker, while $ASP$ measures how a speaker is limited, on average, to
only one cluster. The ideal case is when both $ACP$ and $ASP$ are equal to
$1.0$. Like in the first case, the geometrical mean of $ACP$ and $ASP$ is
applied as an evaluation criterion, $K$.

$\begin{array}[]{*{20}{c}}{\begin{array}[]{*{20}{c}}{ACP=\frac{1}{Q}\sum\limits_{q=1}^{Q}{\frac{1}{{n_{q.}^{2}}}}\sum\limits_{r=1}^{R}{\left|c_{q}\cap
d_{r}\right|}^{2}=\frac{1}{Q}\sum\limits_{q=1}^{Q}{{p_{q.}}}}&;&{{p_{q.}}=\sum\limits_{r=1}^{R}{\frac{{n_{qr}^{2}}}{{n_{q.}^{2}}}}}\\\
{ASP=\frac{1}{R}\sum\limits_{r=1}^{R}{\frac{1}{{n_{.r}^{2}}}}\sum\limits_{r=1}^{Q}{\left|c_{q}\cap
d_{r}\right|}^{2}=\frac{1}{R}\sum\limits_{r=1}^{R}{{p_{.r}}}}&;&{{p_{.r}}=\sum\limits_{r=1}^{Q}{\frac{{n_{qr}^{2}}}{{n_{.r}^{2}}}}}\end{array}}\\\
{K=\sqrt{ACP\cdot ASP}}\end{array}$ (3)

It is important to notice that $ACP$ is based on cluster purities
$\left\\{{{p_{q.}}}\right\\}_{q=1}^{Q}$ while $ASP$ is based on the speaker
purities $\left\\{{{p_{.r}}}\right\\}_{q=1}^{Q}$. In both cases, the purities
are not probabilities and they do not sum to one. As the criteria $ACP$ and
$ASP$ are averages over purities and the information about the class or
cluster sizes is omitted, it is biased toward the small classes/clusters. We
will maintain this formulation in order to be consistent with previous work.

The values of the cluster purities $\left\\{{{p_{q.}}}\right\\}_{q=1}^{Q}$ and
the speaker purities $\left\\{{{p_{.r}}}\right\\}_{q=1}^{Q}$ are closely
related to the multi-class Gini measure [6], and actually can be seen as
$1-Gini\left(p_{q.}\right)$ and $1-Gini\left(p_{.r}\right)$ respectively.

There are advantages and disadvantages to both criteria. The disadvantage of
one is the advantage of the other and vice versa. The advantage of the
criterion $G$, which is based on $Pur_{\mathcal{C}}$ and $Pur_{\mathcal{D}}$,
is in the fact that it gives a portion of the ”correct” data points from all
the clusters, i.e., at the end it is all divided by the dataset size $N$. For
ACP and ACP, their values depend on averaging the purities of all the
clusters, i.e., all clusters have the same influence, independent of their
sizes. The disadvantage of the $G$ criterion is that it takes into account
only the size of the largest class in the cluster without taking into account
the distribution of other classes in the cluster. It is not the same, if in
addition to the largest class, there is only one additional large class in the
cluster or many small classes. Similarly whether one class splits between two
clusters or many clusters makes a difference. The $K$ criterion takes this
into account, as $ASP$ and $ACP$ are calculated according to the class and
cluster distributions and it is not just the majority which is the only
important quantity in $Pur_{\mathcal{C}}$ and $Pur_{\mathcal{D}}$.

##  3 Synthesized Datasets

In this section, the data designed for comparison between the standard
(deterministic) mean-shift and the stochastic mean-shift algorithms are
presented.

As mean-shift clustering is an estimate of the multi-modal probability density
by iterative ”climbing” of each datum to its mode, several multi-modal
datasets were designed. Synthesized datasets include $2$\- and $3$-dimensional
multi-modal distributions. Datasets are differ by their distributions in
several parameters. Mixture of Gaussian multivariate classes are defer in:

  * •

Data dimension - $d$,

  * •

Number of classes (Gaussian mixtures) - $R$,

  * •

Amount of data per class - $n_{.r}$,

  * •

The mean vector of each class $\eta_{r}$,

  * •

The covariance matrix of each class $\Sigma_{r}$.

Seven datasets were designed:

  1. Set 1:

Three 2 dimensional classes with:

$\begin{array}[]{*{20}{c}}{\eta_{1}=\left[{\begin{array}[]{*{20}{c}}1\\\
1\end{array}}\right]}&{\eta_{2}=\left[{\begin{array}[]{*{20}{c}}-1\\\
-1\end{array}}\right]}&{\eta_{3}=\left[{\begin{array}[]{*{20}{c}}1\\\
-1\end{array}}\right]}&{\Sigma_{r}=\left[{\begin{array}[]{*{20}{c}}0.36&0.00\\\
0.00&0.36\end{array}}\right]}\end{array}$

same isotropic covariance matrix for all classes and $n_{.r}=250$ for all
classes,

  2. Set 2:

same as Set 1, but:

${\Sigma_{r}=\left[{\begin{array}[]{*{20}{c}}0.64&0.00\\\
0.00&0.64\end{array}}\right]}$

Broader Gaussians, while the means are the same, should be harder to cluster.

  3. Set 3:

Same as Set 2, but $n_{.r}=50$, i.e. $5$ times less data,

  4. Set 4:

Same as Set 2, but $n_{.r}=1500$, i.e. $6$ times more data,

  5. Set 5:

Same as Set 2, but
$\left\\{n_{.r}\right\\}=\begin{Bmatrix}100&300&20\end{Bmatrix}$, i.e.
unbalanced data,

  6. Set 6:

The mean vectors are the same, however the covariance matrices are not all
spherical, $n_{.r}=250$ for all classes:

$\begin{array}[]{*{20}{c}}{\Sigma_{1}=\left[{\begin{array}[]{*{20}{c}}0.64&0.00\\\
0.00&0.64\end{array}}\right]}&{\Sigma_{2}=\left[{\begin{array}[]{*{20}{c}}0.73&0.48\\\
0.48&0.73\end{array}}\right]}&{\Sigma_{3}=\left[{\begin{array}[]{*{20}{c}}1.09&-0.60\\\
-0.60&1.09\end{array}}\right]}\end{array}$

  7. Set 7:

Four 3-dimensional classes with
$\left\\{n_{.r}\right\\}=\begin{Bmatrix}250&300&200&200\end{Bmatrix}$:

$\begin{array}[]{*{20}{c}}{\eta_{1}=\left[{\begin{array}[]{*{20}{c}}1\\\ 1\\\
1\end{array}}\right]}&{\eta_{2}=\left[{\begin{array}[]{*{20}{c}}-1\\\ -1\\\
0\end{array}}\right]}&{\eta_{3}=\left[{\begin{array}[]{*{20}{c}}1\\\ -1\\\
0\end{array}}\right]}&{\eta_{4}=\left[{\begin{array}[]{*{20}{c}}-2\\\ 2\\\
2\end{array}}\right]}\end{array}$

$\begin{array}[]{*{20}{c}}{\Sigma_{1}=\left[{\begin{array}[]{*{20}{c}}0.64&0.00&0.00\\\
0.00&0.64&0.00\\\
0.00&0.00&0.64\end{array}}\right]}&{\Sigma_{2}=\left[{\begin{array}[]{*{20}{c}}0.74&0.50&-0.20\\\
0.50&0.77&-0.31\\\ -0.2&-0.31&0.41\end{array}}\right]}\\\
{\Sigma_{3}=\left[{\begin{array}[]{*{20}{c}}1.09&-0.60&-0.24\\\
-0.60&1.73&1.60\\\
-0.24&1.60&1.64\end{array}}\right]}&{\Sigma_{4}=\left[{\begin{array}[]{*{20}{c}}0.77&0.58&0.73\\\
0.58&0.56&0.46\\\ 0.73&0.46&0.78\end{array}}\right]}\end{array}$

##  4 Experiments and results

In this section we present the comparison between deterministic and stochastic
mean-shift algorithms applied to a $2$\- and $3$-dimensional synthesized data
(subsection 4.1).

###  4.1 Multi-modal data clustering

In this section we present the results on the datasets described in section 3.
In Figure 2 the clustering process is presented. Set 3 was taken as a test
case. The true labels data and the unlabeled data are presented at sub-plots
(a) and (b) respectively. Sub-plot (c) shows the paths of the data points
during their shifts when the deterministic algorithm is applied, while sub-
plot (d) shows the paths when the stochastic algorithm is applied. It can be
seen that, in the deterministic case, the directed paths to the modes can be
observed, while in the stochastic case they are much less visible. However,
the high-density areas at the modes are clear. Sub-plots (e) and (f) shows the
found modes of both algorithms. In this example there are far fewer modes in
the stochastic case. This does not mean that deterministic clustering results
in more clusters, as at the last stage of the algorithm, close modes are
merged into the same cluster. Sub-plots (g) and (h) show the final clusters of
the deterministic and the stochastic algorithms respectively. In this example,
deterministic mean-shift clustering ended with $14$ clusters while stochastic
mean-shift clustering provided only $9$ clusters.

Figure 2:  Results of Set 3 (a) data with true labels; (b) unlabeled data; (c)
converging paths of the deterministic mean-shift; (d) converging paths of the
stochastic mean-shift; (e) the found modes of the deterministic mean-shift;
(f) the found modes of the stochastic mean-shift; (g) the clustering of the
deterministic mean-shift; (f) the found clustering of the stochastic mean-
shift.

The results according to the datasets are as follows:

  * •

Set 1: The dataset consists of the Gaussian mixtures with diagonal covariance
matrix. All the variances are $0.36^{2}$ and $250$ data points per class. The
results are summarized in Table1. We can observe that using $K$ as evaluation
criterion, the results of the stochastic algorithm is better than the
deterministic, while according to $G$ as evaluation criterion they are almost
the same. This can be explained as followed: $K$ takes into account the
distribution of all classes in the cluster and not only the one that has the
majority in the cluster to evaluate the cluster purity (ACP versus
$Pur_{\mathcal{C}}$), and the same for class purity (ASP versus
$Pur_{\mathcal{D}}$). In addition, the number of detected clusters is also
provided.

Table 1: The results of Set $1$ (D \- deterministic mean-shift; S \- stochastic mean-shift). | $ACP$ | $ASP$ | $K$ | $Pur_{\mathcal{C}}$ | $Pur_{\mathcal{D}}$ | $G$ | #Clusters
---|---|---|---|---|---|---|---
D | $0.78$ | $0.84$ | $0.81$ | $0.92$ | $0.92$ | $0.92$ | $4$
S | $0.90$ | $0.82$ | $0.86$ | $0.92$ | $0.90$ | $0.91$ | $9$
  * •

Set 2: \- It is similar to Set $1$, but the variances are $0.64^{2}$. As the
means are the same in Set $1$, the clustering problem becomes harder due to
the larger overlap between the clusters. The results are summarized in Table 2
and confirm the increase in clustering dificulty. Comparison of Tables 1 and 2
confirms the degradation, however in terms of $K$, the stochastic algorithm
degrades more, for $G$ a little less.

Table 2: The results of Set $2$ (D \- deterministic mean-shift; S \- stochastic mean-shift). | $ACP$ | $ASP$ | $K$ | $Pur_{\mathcal{C}}$ | $Pur_{\mathcal{D}}$ | $G$ | #Clusters
---|---|---|---|---|---|---|---
D | $0.85$ | $0.70$ | $0.77$ | $0.85$ | $0.82$ | $0.84$ | $8$
S | $0.81$ | $0.72$ | $0.77$ | $0.87$ | $0.84$ | $0.86$ | $10$
  * •

Set 3: In this case we repeat the Set $2$ case but with five times less data
per cluster (Figure 2). The results are summarized in Table 3. In the
deterministic case, the clustering ended with many clusters; each was pure,
however, as the data spread across many clusters, the class purity is very
low. In the stochastic case, there are much fewer clusters, approximately as
for Set $2$ and the results are similar.

Table 3: The results of Set 3 (D \- deterministic mean-shift; S \- stochastic mean-shift). | $ACP$ | $ASP$ | $K$ | $Pur_{\mathcal{C}}$ | $Pur_{\mathcal{D}}$ | $G$ | #Clusters
---|---|---|---|---|---|---|---
D | $0.92$ | $0.36$ | $0.57$ | $0.85$ | $0.53$ | $0.67$ | $14$
S | $0.89$ | $0.65$ | $0.76$ | $0.88$ | $0.79$ | $0.84$ | $9$
  * •

Set 4: Here the data is still balanced, but with much more data points per
class ($1500$). The results are summarized in Table 4. Two things are of great
interest. First, the number of clusters in the deterministic mean-shift is
much higher than in the stochastic mean-shift. We observed this phenomenon
when we repeated this experiment several times. Second, there is a very large
difference between the ACP and the $Pur_{C}$ in the deterministic mean-shift.
This can be explained when looking at the results in Figure 3. The data is
very dense and in the deterministic case, as at the end of the convergence of
one datum it is placed back to the original place before starting with the
next datum, all the modes are very close to each other as can be seen in
subplot 3e. This is not the situation in the stochastic case, as if the data
points shifted and drifted apart, they are not placed back and each datum
continues to move to a different mode (subplot 3f). This causes one large
cluster during the merging and several outliers in the deterministic case
(subplot 3g), versus 3h where we can see three well-separated clusters; one
small cluster at the intersection of the other three, and all the others can
be viewed as outliers. The ACP calculation averages the cluster purity of each
cluster, without taking into account the size of each cluster. This leads to a
relatively high ACP as four clusters of the outliers are pure and one is not.
In contrast, $Pur_{\mathcal{C}}$ counts only the total number of well-
clustered data points. This leads to a low value of $Pur_{\mathcal{C}}$. It is
important to mention that the results of the stochastic mean-shift clustering
almost unchanged as a function of the amount of the data per class, Sets
$2-4$.

Table 4: The results of Set $4$ (D \- deterministic mean-shift; S \- stochastic mean-shift). | $ACP$ | $ASP$ | $K$ | $Pur_{\mathcal{C}}$ | $Pur_{\mathcal{D}}$ | $G$ | #Clusters
---|---|---|---|---|---|---|---
D | $0.87$ | $1.00$ | $0.93$ | $0.33$ | $1.00$ | $0.58$ | $5$
S | $0.88$ | $0.69$ | $0.78$ | $0.85$ | $0.82$ | $0.84$ | $20$

Figure 3:  Results of Set $4$ (a) data with true labels; (b) unlabeled data;
(c) converging paths of the deterministic mean-shift; (d) converging paths of
the stochastic mean-shift; (e) the modes found of the deterministic mean-
shift; (f) the modes found of the stochastic mean-shift; (g) the clustering of
the deterministic mean-shift; (f) the found clustering of the stochastic mean-
shift.

  * •

Set 5: In this case the data was unbalanced, but the number of data points was
not too large. The result are summarized in Table 5. From our observation, it
is not conclusive which algorithm is better, similar to the results we
obtained for Sets $1$ and $2$, Tables 1 and 2 respectively. The results become
constantly better for the stochastic mean-shift when the amount of the data is
increased by about an order: instead of
$\left\\{n_{.r}\right\\}=\begin{Bmatrix}100&300&20\end{Bmatrix}$, cluster
sizes were increased to
$\left\\{n_{.r}\right\\}=\begin{Bmatrix}1250&3750&250\end{Bmatrix}$ ($12.5$
times more). This is due to the same effect as in the previous experiment with
Set $4$.

Table 5: The results of Set $5$ (D \- deterministic mean-shift; S \- stochastic mean-shift). | $ACP$ | $ASP$ | $K$ | $Pur_{\mathcal{C}}$ | $Pur_{\mathcal{D}}$ | $G$ | #Clusters
---|---|---|---|---|---|---|---
D | $0.95$ | $0.59$ | $0.75$ | $0.92$ | $0.82$ | $0.87$ | $9$
S | $0.88$ | $0.69$ | $0.78$ | $0.94$ | $0.92$ | $0.93$ | $9$
  * •

Set 6: In this experiment the data was generated with non-diagonal covariance
matrices with same cluster size, $250$. The results are summarized in Table 6.
The stochastic algorithm ended with $2$ times more clusters and the much
better fit to the same class (both according to ACP and $Pur_{\mathcal{C}}$).
In total, the stochastic mean-shift performed more accurate.

Table 6: The results of Set $6$ (D \- deterministic mean-shift; S \- stochastic mean-shift). | $ACP$ | $ASP$ | $K$ | $Pur_{\mathcal{C}}$ | $Pur_{\mathcal{D}}$ | $G$ | #Clusters
---|---|---|---|---|---|---|---
D | $0.71$ | $0.69$ | $0.70$ | $0.71$ | $0.78$ | $0.74$ | $5$
S | $0.91$ | $0.65$ | $0.77$ | $0.77$ | $0.75$ | $0.76$ | $10$
  * •

Set 7: In the last experiment the data was $3$-dimensional with four clusters
that have similar, but not the same, sizes (all between $200$ and $300$ data
points). The results are summarized in Table 7. The performance of the
stochastic mean-shift are slightly better than of the deterministic mean-
shift. It seems that increasing both, the dimenssionality of the data and the
number of clusters results that the number of ”outlier” clusters is increases.
This can be tuned by increasing the threshold of merging the data point to
clusters.

Table 7: The results of Set $7$ (D \- deterministic mean-shift; S \- stochastic mean-shift). | $ACP$ | $ASP$ | $K$ | $Pur_{\mathcal{C}}$ | $Pur_{\mathcal{D}}$ | $G$ | #Clusters
---|---|---|---|---|---|---|---
D | $0.89$ | $0.47$ | $0.65$ | $0.87$ | $0.62$ | $0.73$ | $32$
S | $0.92$ | $0.52$ | $0.69$ | $0.86$ | $0.68$ | $0.76$ | $24$

The summary of all the experiments is given in Table 8. The best results of
each experiment according to each of the criteria ($K$ or $G$) are appears in
bold. It can be seen that the stochastic mean-shift is mostly performed better
than the deterministic mean-shift. Only for the Set $4$ for $K$ criterion the
deterministic mean-shift gave a significantly higher score. It can be
explained by the fact that the stochastic mean-shift produced $20$ clusters
instead of actual $3$, while the deterministic mean-shift produced only $5$
clusters. In such a case, the ASP for the in the stochastic case was
relatively low. However, at Figure 3 we see the reason of such a difference.
In the deterministic case, almost all the data converged to a single cluster.
It leads to a very high ASP. The ACP is not very low due to the fact that the
outlier clusters are very pure.

Table 8: The results of Set $7$ (D \- deterministic mean-shift; S \- stochastic mean-shift). | Set 1 | Set 2 | Set 3 | Set 4 | Set 5 | Set 6 | Set 7
---|---|---|---|---|---|---|---
| D | S | D | S | D | S | D | S | D | S | D | S | D | S
$K$ | $0.81$ | $0.86$ | $0.77$ | $0.77$ | $0.57$ | $0.76$ | $0.93$ | $0.78$ | $0.75$ | $0.78$ | $0.70$ | $0.77$ | $0.55$ | $0.69$
$G$ | $0.92$ | $0.91$ | $0.84$ | $0.86$ | $0.67$ | $0.84$ | $0.58$ | $0.84$ | $0.87$ | $0.93$ | $0.74$ | $0.76$ | $0.73$ | $0.76$

##  5 Discussion and conclusions

In this work we presented a stochastic version of the mean-shift clustering
algorithm. Originally, the mean-shift algorithm is deterministic. Each vector
in the data ”climbs” separately to one of the pdf modes, while other vectors
remain at their original places. Then, in turn, the same procedure is applied
to another data point, and this procedure is repeated until all the vectors
reach their modes. In the algorithm presented, for each iteration a new datum
is chosen randomly, with replacement for a single shift only. The neighbors of
each datum are found according to data points’ current positions and not
according to the original position (as in the deterministic case). In such a
strategy, the data points converge to their modes jointly rather than one by
one.

The original mean-shift clustering was compared with its stochastic version.
Two- and three-dimensional data sets were synthesized with three and four
Gaussian mixtures. Different scenarios were examined: balanced and unbalanced
data, small ans large amounts of data per class and different covariance
matrices. It can be seen that the stochastic algorithm performed as well, or
better than, the deterministic mean-shift (Tables 1-7). The summary of the
results is presented in Table 8. It seems that the stochastic mean-shift has
an advantage when the amount of data is large and the modes becomes more
blurred (Figure 3 and Table 4, G criterion). Frequently, the stochastic mean-
shift produces more clusters than the deterministic mean-shift. Nevertheless,
the performance is at least as good as in the deterministic case. This is
because the extra clusters are small, and the result of data points that were
rarely chosen (low density) and become outlier clusters.

Another observation relates to the number of estimated clusters. Both
algorithms suffer from under-clustering and produce more clusters than the
true number of classes. However, the stochastic algorithm produces more
clusters than the deterministic mean-shift clustering in most of the
experiments. This can be maintained by tuning the hyper-parameter of $Th_{2}$
which responsible for the merging procedure. Increase of $Th_{2}$ leads to
smaller number of larger clusters., however, it does not lead to a better
clustering.

In summary, in most of the examples, stochastic mean-shift clustering
performed better or at least as good as the deterministic mean-shift
clustering. In the future it will be tested on real data to compare these
clustering algorithms.

## References

  * Ai & Xiong [2014] Ai, X., L.and Gao, & Xiong, J. (2014).  Application of mean-shift clustering to blood oxygen level dependent functional mri activation detection.  BMC Medical Imaging, 14.
  * Ajmera et al. [2002] Ajmera, J., Bourlard, H., Lapidot, I., & McCowan, I. (2002).  Unknown-multiple speaker clustering using HMM.  In In Proceedings of ICSLP-2002 (pp. 573–576).
  * Anguera et al. [2006] Anguera, X., Wooters, C., & Hernando, J. (2006).  Purity algorithms for speaker diarization of meetings data.  In 2006 IEEE International Conference on Acoustics Speech and Signal Processing Proceedings (pp. I–I).  volume 1.
  * Banerjee et al. [2012] Banerjee, B., G., S. V., & Buddhiraju, K. M. (2012).  Satellite image segmentation: A novel adaptive mean-shift clustering based approach.  In 2012 IEEE International Geoscience and Remote Sensing Symposium (pp. 4319–4322).
  * Bolshakova & Azuaje [2003] Bolshakova, N., & Azuaje, F. (2003).  Cluster validation techniques for genome expression data.  Signal Processing, 83, 825 – 833.
  * Breiman et al. [1984] Breiman, L., Friedman, J., Stone, C. J., & Olshen, R. (1984).  Classification and Regression Trees.  Chapman and Hall/CRC.
  * Carreira-Perpiñán [2006] Carreira-Perpiñán, M. Á. (2006).  Fast nonparametric clustering with gaussian blurring mean-shift.  Proceedings of the 23rd international conference on Machine learning, .
  * Chen et al. [2002] Chen, G., Jaradat, S. A., Banerjee, N., Tanaka, T. S., Ko, M. S. H., & Zhang, M. Q. (2002).  Evaluation and comparison of clustering algorithms in analyzing es cell gene expression data.  Statistica Sinica, (pp. 241–262).
  * Cohen & Lapidot [2017] Cohen, Y., & Lapidot, I. (2017).  Estimating speaker clustering quality using logistic regression.  In Proc. Interspeech 2017 (pp. 3577–3581).
  * Cohen & Lapidot [2021] Cohen, Y., & Lapidot, I. (2021).  Speaker clustering quality estimation with logistic regression.  Computer Speech and Language, 65, 101139.
  * Comaniciu & Meer [2002] Comaniciu, D., & Meer, P. (2002).  Mean shift: A robust approach toward feature space analysis.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 24, 603–619.
  * Dabbabi et al. [2017] Dabbabi, K., Hajji, S., & Cherif, A. (2017).  Integration of evolutionary computation algorithms and new AUTO-TLBO technique in the speaker clustering stage for speaker diarization of broadcast news.  EURASIP J. Audio Speech Music. Process., 2017, 21.
  * Fukunaga & Hostetler [1975] Fukunaga, K., & Hostetler, L. (1975).  The estimation of the gradient of a density function, with applications in pattern recognition.  IEEE Transactions on Information Theory, 21, 32–40.
  * Lapidot [2023] Lapidot, I. (2023).  Stochastic mean-shift for speaker clustering.  In 2023 57th Annual Conference on Information Sciences and Systems (CISS) (pp. 1–5).
  * Lu et al. [2011] Lu, J., Chen, J., Zhang, J., & Zou, L. (2011).  Medical image segmentation using mean shift algorithm and general edge detection.  IFAC Proceedings Volumes, 44, 9656–9661.  18th IFAC World Congress.
  * Manning et al. [2008] Manning, C. D., Raghavan, P., & Schütze, H. (2008).  Introduction to Information Retrieval.  Cambridge, UK: Cambridge University Press.
  * Paris & Durand [2007] Paris, S., & Durand, F. (2007).  A topological approach to hierarchical segmentation using mean shift.  In 2007 IEEE Conference on Computer Vision and Pattern Recognition (pp. 1–8).
  * Salmun et al. [2016a] Salmun, I., Opher, I., & Lapidot, I. (2016a).  Improvements to plda i-vector scoring for short segments clustering.  In 2016 IEEE International Conference on the Science of Electrical Engineering (ICSEE) (pp. 1–4).
  * Salmun et al. [2016b] Salmun, I., Opher, I., & Lapidot, I. (2016b).  On the use of plda i-vector scoring for clustering short segments.  In ODYSSEY 2016 -The Speaker and Language Recognition Workshop.
  * Salmun et al. [2017] Salmun, I., Shapiro, I., Opher, I., & Lapidot, I. (2017).  Plda-based mean shift speakers’ short segments clustering.  Computer Speech and Language, 45, 411 – 436.
  * Saptarshi et al. [2021] Saptarshi, C., Debolina, P., & Das, S. (2021).  Automated clustering of high-dimensional data with a feature weighted mean shift algorithm.  In The Thirty-Fifth AAAI Conference on Artificial Intelligence (AAAI-21).
  * Senoussaoui et al. [2014] Senoussaoui, M., Kenny, P., Stafylakis, T., & Dumouchel, P. (2014).  A study of the cosine distance-based mean shift for telephone speech diarization.  IEEE/ACM Transactions on Audio, Speech, and Language Processing, 22, 217–227.
  * Tao et al. [2007] Tao, W., Jin, H., & Zhang, Y. (2007).  Color image segmentation based on mean shift and normalized cuts.  IEEE Transactions on Systems, Man, and Cybernetics, Part B (Cybernetics), 37, 1382–1389.
  * Tibshirani et al. [2001] Tibshirani, R., Walther, G., & Hastie, T. (2001).  Estimating the number of clusters in a data set via the gap statistic.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 63, 411–423.
  * Valente & Wellekens [2004] Valente, F., & Wellekens, C. (2004).  Variational Bayesian speaker clustering.  In Proc. The Speaker and Language Recognition Workshop (Odyssey 2004) (pp. 207–214).
  * Wang et al. [2004] Wang, J., Thiesson, B., Xu, Y., & Cohen, M. (2004).  Image and video segmentation by anisotropic kernel mean shift.  In T. Pajdla, & J. Matas (Eds.), Computer Vision - ECCV 2004 (pp. 238–249).  Berlin, Heidelberg: Springer Berlin Heidelberg.
  * Wu & Luo [2015] Wu, X., G.and Zhao, & Luo, S. (2015).  Histological image segmentation using fast mean shift clustering method.  BioMed Eng OnLine, 14.


# Solving FDE-IVPs by using Fractional HBVMs:
some experiments with the fhbvm code

Luigi Brugnano 111 Dipartimento di Matematica e Informatica “U. Dini”,
Università di Firenze, Italy
<EMAIL_ADDRESS>https://orcid.org/0000-0002-6290-4107 Gianmarco
Gurioli 222 Dipartimento di Matematica e Informatica “U. Dini”, Università di
Firenze, Italy
<EMAIL_ADDRESS>https://orcid.org/0000-0003-0922-8119 Felice
Iavernaro 333 Dipartimento di Matematica, Università di Bari Aldo Moro, Italy
<EMAIL_ADDRESS>https://orcid.org/0000-0002-9716-7370

###### Abstract

In this paper we report a few numerical tests by using a slight extension of
the Matlab${}^{\copyright}\,$code fhbvm in [9], implementing Fractional HBVMs,
a recently introduced class of numerical methods for solving Initial Value
Problems of Fractional Differential Equations (FDE-IVPs). The reported
experiments are aimed to give evidence of its effectiveness.

Keywords: fractional differential equations, fractional integrals, Caputo
derivative, Jacobi polynomials, Fractional HBVMs, FHBVMs.

MSC: 34A08, 65R20, 65-04.

##  1 Introduction

Fractional differential equations have become popular in many applications
(see, e.g., the classical references [15, 19] for an introduction). For this
reason, their numerical solution is an active field of investigation. In
particular, we here report a few numerical tests obtained by using the
Matlab${}^{\copyright}\,$code fhbvm [9], available at the URL [22]. The code
itself is based on a class of methods named Fractional HBVMs (FHBVMs),
recently introduced in [7], as a generalization of Hamiltonian Boundary Value
Methods (HBVMs), a class of low-rank Runge-Kutta methods initally devised for
the efficient numerical solution of Hamiltonian dynamical systems (see, e.g.,
[10, 11]), and, subsequently, generalized along several directions (see, e.g.,
[2, 4, 6, 7, 8, 12]). A main feature of HBVMs is the fact that they can gain
spectral accuracy, when approximating ODE-IVPs [3, 13, 14], and this holds
true also in the FDE case [7]. In this paper we shall use a slight extension
of the code in [9], modified for handling FDE-IVPs in the form:444The code in
[9] was conceived for the case $\ell=1$.

$\displaystyle\alpha$ $\displaystyle\in$ $\displaystyle(\ell-1,\ell),$ (1)
$\displaystyle y^{(\alpha)}(t)$ $\displaystyle=$ $\displaystyle f(y(t)),\qquad
t\in[0,T],\qquad y^{(i)}(0)=y_{0}^{i}\in\mathbb{R}^{m},\quad
i=0,\dots,\ell-1,$

where $y^{(\alpha)}(t)$ is the Caputo fractional derivative and, for the sake
of brevity, we have omitted the argument $t$ for $f$. It is known that under
suitable regularity assumptions (see, e.g. [16, 17]), the solution of (1) is
given by:

$y(t)=\sum_{i=0}^{\ell-1}\frac{t^{i}}{i!}y_{0}^{i}+I^{\alpha}f(y(t))\equiv\sum_{i=0}^{\ell-1}\frac{t^{i}}{i!}y_{0}^{i}+\frac{1}{\Gamma(\alpha)}\int_{0}^{t}(t-x)^{\alpha-1}f(y(x))\mathrm{d}x,\qquad
t\in[0,T],$ (2)

where $I^{\alpha}f(y(t))$ is the Riemann-Liouville integral of $f(y(t))$. With
this premise, the structure of the paper is as follows: in Section 2 we recall
the main facts about the numerical solution of FDE-IVPs proposed in [7]; in
Section 3 we report some numerical experiments with the code, also providing
some comparisons with another existing one, to give evidence of its
effectiveness; at last, a few conclusions are given in Section 4.

##  2 Fractional HBVMs

To begin with, in order to obtain a piecewise approximation to the solution of
the problem, we consider a partition of the integration interval in the form:

$t_{n}=t_{n-1}+h_{n},\qquad h_{n}>0,\qquad n=1,\dots,N,\qquad t_{0}=0,\quad
t_{N}=T.$ (3)

In particular, the code fhbvm implements either a graded mesh,

$h_{n}=rh_{n-1}\equiv r^{n-1}h_{1},\qquad n=1\dots,N,$ (4)

where $r>1$ and $h_{1}>0$ satisfy, by virtue of (3)–(4),
$h_{1}\dfrac{r^{N}-1}{r-1}=T,$ or a uniform mesh,

$h_{n}\equiv h_{1}:=\frac{T}{N},\quad n=1,\dots,N,\quad\Rightarrow\quad
t_{n}=nh_{1},\quad n=0,\dots,N.$ (5)

The former mesh is clearly appropriate when the vector field in (1) is not
very smooth at the origin, whereas the latter mesh is recommended, otherwise.
The code automatically selects the most appropriate mesh and its parameters
[9]. The method is aimed at obtaining a piecewise approximation $\sigma(t)$
defined over the partition (3). In order to derive it, let us set,

$y_{n}(ch_{n}):=y(t_{n-1}+ch_{n}),\quad\sigma_{n}(ch_{n}):=\sigma(t_{n-1}+ch_{n}),\quad
c\in[0,1],\qquad n=1,\dots,N,$ (6)

the restriction of the solution of (1) and of its approximation $\sigma(t)$ on
the interval $[t_{n-1},t_{n}]$, respectively. Consequently, the problem (1)
can be rewritten as:

$\left\\{\begin{array}[]{rcl}y_{n}^{(\alpha)}(ch_{n})&=&f(y_{n}(ch_{n})),\qquad
c\in[0,1],\qquad n=1,\dots,N,\\\\[5.69054pt] y_{1}^{(i)}(0)&=&y_{0}^{i},\qquad
i=0,\dots,\ell-1.\end{array}\right.$ (7)

It is to be noticed that the r.h.s. of each local problem in (7) only depends
on the corresponding restriction of the solution, i.e., $y_{n}$. Further, we
consider the expansion of the local vector fields along the following basis of
orthonormal Jacobi polynomials:555Here, $\alpha>0$ is the same parameter in
(1).

$P_{i}\in\Pi_{i},\qquad\int_{0}^{1}\omega(c)P_{i}(c)P_{j}(c)\mathrm{d}c=\delta_{ij},\qquad
i,j=0,1,\dots,\qquad\omega(c):=\alpha(1-c)^{\alpha-1}.$ (8)

Consequently, one obtains:

$f(y_{n}(ch_{n}))=\sum_{j\geq 0}P_{j}(c)\gamma_{j}(y_{n}),\qquad
c\in[0,1],\qquad\gamma_{j}(y_{n})=\int_{0}^{1}\omega(\tau)P_{j}(\tau)f(y_{n}(\tau
h_{n}))\mathrm{d}\tau,$ (9)

so that (7) can be equivalently rewritten as

$\left\\{\begin{array}[]{rcl}y_{n}^{(\alpha)}(ch_{n})&=&\sum_{j\geq
0}P_{j}(c)\gamma_{j}(y_{n}),\qquad c\in[0,1],\qquad
n=1,\dots,N,\\\\[5.69054pt] y_{1}^{(i)}(0)&=&y_{0}^{i},\qquad
i=0,\dots,\ell-1.\end{array}\right.$ (10)

The approximation $\sigma(t)$ is then defined by truncating the above series
to finite sums with $s\geq 1$ terms:

$\left\\{\begin{array}[]{rcl}\sigma_{n}^{(\alpha)}(ch_{n})&=&\sum_{j=0}^{s-1}P_{j}(c)\gamma_{j}(\sigma_{n}),\qquad
c\in[0,1],\qquad n=1,\dots,N,\\\\[5.69054pt]
\sigma_{1}^{(i)}(0)&=&y_{0}^{i},\qquad i=0,\dots,\ell-1,\end{array}\right.$
(11)

with $\gamma_{j}(\sigma_{n})$ defined according to (9), by formally replacing
$y_{n}$ with $\sigma_{n}$.

For the continuous solution, for $t\,\equiv\,t_{n-1}+ch_{n}$, $c\in[0,1]$, one
then obtains (see (3)):

$\displaystyle y(t)$ $\displaystyle\equiv$ $\displaystyle
y_{n}(ch_{n})~{}=~{}\sum_{i=0}^{\ell-1}\frac{t^{i}}{i!}y_{0}^{i}~{}+~{}I^{\alpha}f(y(t_{n-1}+ch_{n}))$
$\displaystyle=$
$\displaystyle\sum_{i=0}^{\ell-1}\frac{t^{i}}{i!}y_{0}^{i}~{}+~{}\frac{1}{\Gamma(\alpha)}\int_{0}^{t_{n-1}+ch_{n}}(t_{n-1}+ch_{n}-x)^{\alpha-1}f(y(x))\mathrm{d}x$
$\displaystyle=$
$\displaystyle\sum_{i=0}^{\ell-1}\frac{t^{i}}{i!}y_{0}^{i}~{}+~{}\frac{1}{\Gamma(\alpha)}\sum_{\nu=1}^{n-1}\int_{t_{\nu-1}}^{t_{\nu}}(t_{n-1}+ch_{n}-x)^{\alpha-1}f(y(x))\mathrm{d}x$
$\displaystyle\qquad+~{}\frac{1}{\Gamma(\alpha)}\int_{t_{n-1}}^{t_{n-1}+ch_{n}}(t_{n-1}+ch_{n}-x)^{\alpha-1}f(y(x))\mathrm{d}x$
$\displaystyle=$
$\displaystyle\sum_{i=0}^{\ell-1}\frac{t^{i}}{i!}y_{0}^{i}~{}+~{}\frac{1}{\Gamma(\alpha)}\sum_{\nu=1}^{n-1}\int_{0}^{h_{\nu}}(t_{n-1}-t_{\nu-1}+ch_{n}-x)^{\alpha-1}f(y_{\nu}(x))\mathrm{d}x$
$\displaystyle\qquad+~{}\frac{1}{\Gamma(\alpha)}\int_{0}^{ch_{n}}(ch_{n}-x)^{\alpha-1}f(y_{n}(x))\mathrm{d}x$
$\displaystyle=$
$\displaystyle\sum_{i=0}^{\ell-1}\frac{t^{i}}{i!}y_{0}^{i}~{}+~{}\frac{1}{\Gamma(\alpha)}\sum_{\nu=1}^{n-1}h_{\nu}^{\alpha}\int_{0}^{1}\left(\frac{t_{n-1}-t_{\nu-1}}{h_{\nu}}+c\frac{h_{n}}{h_{\nu}}-\tau\right)^{\alpha-1}f(y_{\nu}(\tau
h_{\nu}))\mathrm{d}\tau$
$\displaystyle\qquad+~{}\frac{h_{n}^{\alpha}}{\Gamma(\alpha)}\int_{0}^{c}(c-\tau)^{\alpha-1}f(y_{n}(\tau
h_{n}))\mathrm{d}\tau,\qquad c\in[0,1].$

Consequently, by considering, hereafter, the case of the graded mesh (4), and
taking into account (10), one has:666Similar arguments can be used in the case
of the uniform mesh (5), see [7, 9].

$\displaystyle y_{n}(ch_{n})$ $\displaystyle=$
$\displaystyle\sum_{i=0}^{\ell-1}\frac{t^{i}}{i!}y_{0}^{i}~{}+~{}\frac{1}{\Gamma(\alpha)}\sum_{\nu=1}^{n-1}h_{\nu}^{\alpha}\int_{0}^{1}\left(\frac{r^{n-\nu}-1}{r-1}+cr^{n-\nu}-\tau\right)^{\alpha-1}\sum_{j\geq
0}P_{j}(\tau)\gamma_{j}(y_{\nu})\mathrm{d}\tau$ (12)
$\displaystyle\qquad+~{}\frac{h_{n}^{\alpha}}{\Gamma(\alpha)}\int_{0}^{c}(c-\tau)^{\alpha-1}\sum_{j\geq
0}P_{j}(\tau)\gamma_{j}(y_{n})\mathrm{d}\tau,\qquad c\in[0,1].$

According to (11), the corresponding restriction of the approximation
$\sigma(t)$ on the interval $[t_{n-1},t_{n}]$ is obtained by truncating the
infinite series in (12) to finite sums:

$\displaystyle\sigma_{n}(ch_{n})$ $\displaystyle=$
$\displaystyle\sum_{i=0}^{\ell-1}\frac{t^{i}}{i!}y_{0}^{i}~{}+~{}\frac{1}{\Gamma(\alpha)}\sum_{\nu=1}^{n-1}h_{\nu}^{\alpha}\int_{0}^{1}\left(\frac{r^{n-\nu}-1}{r-1}+cr^{n-\nu}-\tau\right)^{\alpha-1}\sum_{j=0}^{s-1}P_{j}(\tau)\gamma_{j}(\sigma_{\nu})\mathrm{d}\tau$
(13)
$\displaystyle\qquad+~{}\frac{h_{n}^{\alpha}}{\Gamma(\alpha)}\int_{0}^{c}(c-\tau)^{\alpha-1}\sum_{j=0}^{s-1}P_{j}(\tau)\gamma_{j}(\sigma_{n})\mathrm{d}\tau$
$\displaystyle=:$
$\displaystyle\phi_{n-1}^{\alpha,s}(c)~{}+~{}\frac{h_{n}^{\alpha}}{\Gamma(\alpha)}\int_{0}^{c}(c-\tau)^{\alpha-1}\sum_{j=0}^{s-1}P_{j}(\tau)\gamma_{j}(\sigma_{n})\mathrm{d}\tau$
$\displaystyle\equiv$
$\displaystyle\phi_{n-1}^{\alpha,s}(c)~{}+h_{n}^{\alpha}\sum_{j=0}^{s-1}I^{\alpha}P_{j}(c)\gamma_{j}(\sigma_{n}),\qquad
c\in[0,1],$

with $I^{\alpha}P_{j}(c)$ the Riemann-Liouville integral of $P_{j}(c)$ (see
(2)). We observe that in (13) the memory term $\phi_{n-1}^{\alpha,s}(c)$ only
depends on past values of the approximation (in fact, from $\sigma_{\nu}$,
$\nu=1,\dots,n-1$). Consequently, the local problem consists in computing the
Fourier coefficients $\gamma_{j}(\sigma_{n})$, $j=0,\dots,s-1$. These latter
coefficients can be approximated to full machine accuracy by using a Gauss-
Jacobi formula of order $2k$, which is based at the zeros,
$c_{1},\dots,c_{k}$, of $P_{k}(c)$, with corresponding weights (see (8))

$b_{i}=\int_{0}^{1}\omega(c)\ell_{i}(c)\mathrm{d}c,\qquad\ell_{i}(c)=\prod_{j\neq
i}\frac{c-c_{j}}{c_{i}-c_{j}},\qquad i=1,\dots,k,$

by choosing a value of $k$, $k\geq s$, large enough. In other words,

$\gamma_{j}^{n}\,:=\,\sum_{i=1}^{k}b_{i}P_{j}(c_{i})f(\sigma_{n}(c_{i}h_{n}))\,\doteq\,\gamma_{j}(\sigma_{n}),\qquad
j=0,\dots,s-1,$ (14)

where $\,\doteq\,$ means equal within machine precision. Because of this fact,
and for sake of brevity, we shall continue using $\sigma_{n}$ to denote the
fully discrete approximation:

$\sigma_{n}(ch_{n})~{}=~{}\phi_{n-1}^{\alpha,s}(c)~{}+~{}h_{n}^{\alpha}\sum_{j=0}^{s-1}I^{\alpha}P_{j}(c)\gamma_{j}^{n},\qquad
c\in[0,1].$ (15)

We observe that, in order to compute the (discrete) Fourier coefficients, it
is enough evaluating the memory term only at the quadrature abscissae
$c_{1},\dots,c_{k}$. In so doing, from (14)-(15) one obtains a discrete
problem in the form:

$\gamma_{j}^{n}=\sum_{i=1}^{k}b_{i}P_{j}(c_{i})f\left(\phi_{n-1}^{\alpha,s}(c_{i})+h_{n}^{\alpha}\sum_{\ell=0}^{s-1}I^{\alpha}P_{\ell}(c_{i})\gamma_{\ell}^{n}\right),\qquad
j=0,\dots,s-1.$ (16)

Once it has been solved, by taking into account (15) and the fact that
$I^{\alpha}P_{j}(1)=\dfrac{\delta_{j0}}{\Gamma(\alpha+1)}$, the approximation
to $y(t_{n})$ is given by:

$\bar{y}_{n}\,:=\,\sigma_{n}(h_{n})\,=\,\phi_{n-1}^{\alpha,s}(1)\,+\,\frac{h_{n}^{\alpha}}{\Gamma(\alpha+1)}\gamma_{0}^{n}.$
(17)

The equations (16)–(17)) define a fractional HBVM with parameters $(k,s)$. In
short, FHBVM$(k,s)$. We refer to [7] for their analysis, and to [9] for all
implementation details.

##  3 Numerical Tests

In this section we report a few numerical tests using the
Matlab${}^{\copyright}\,$code fhbvm, implementing a FHBVM(22,20) method. The
calling sequence of the code is:

[t,y,stats,err] = fhbvm( fun, y0, T, M )

where:

  * •

fun is the function evaluating the r.h.s. of the equation (in vector mode),
its Jacobian, and the order $\alpha$ of the fractional derivative;

  * •

y0 is a matrix $\ell\times m$ containing the initial conditions in (1);

  * •

T is the final integration time;

  * •

M is an integer parameter such that $h_{\max}=~{}$T/M is approximately the
largest stepsize in the final mesh (it should be as small as possible);

  * •

t,y contain the computed mesh and solution;

  * •

stats (optional) is a vector containing some time statistics;

  * •

err (optional) contains the estimate of the absolute error obtained on a
doubled mesh.

We refer to the inline help (help fhbvm) for further details. The code is
available at the url [22]. All numerical tests have been done on a M2-Silicon
based computer with 16GB of shared memory, using
Matlab${}^{\copyright}\,$R2023b. We shall also make a comparison with the
Matlab${}^{\copyright}\,$code flmm2 [17],777In particular, the BDF2 method is
selected (method=3), with the parameters tol=1e-15 and itmax=1000. in order to
emphasize the potentialities of the new code. The comparisons will be done by
using a so called Work Precision Diagram (WPD), where the execution time (in
sec) is plotted against accuracy. The accuracy, in turn, is measured through
the mixed error significant computed digits (mescd) [20], defined, by using a
Matlab-like notation, as 888This definition corresponds to set atol=rtol in
the definition used in [20].

${\tt
mescd}:=-\log_{10}\max_{i=0,\dots,N}\|(y(t_{i})-\bar{y}_{i})./(1+|y(t_{i})|)\|_{\infty},$

being $t_{i}$, $i=0,\dots,N$, the computational mesh of the considered solver,
and $y(t_{i})$ and $\bar{y}_{i}$ the corresponding values of the solution and
of its approximation. The following parameters are used for all the examples
below, to compute the corresponding WPD:

  * •

flmm2 : $h=\frac{1}{5}2^{-\nu}$, $\nu=1,\dots,20$;

  * •

fhbvm : $M=2,3,4,5$.

### Example 1

The first problem is given by the Example C.1, in [15, page 207], for
$\alpha\in(1,2)$:

$\displaystyle y^{(\alpha)}(t)$ $\displaystyle=$
$\displaystyle-y(t)^{3/2}+\frac{40320}{\Gamma(9-\alpha)}t^{8-\alpha}-3\frac{\Gamma(5+\alpha/2)}{\gamma(5-\alpha/2)}t^{4-\alpha/2}+\left(\frac{3}{2}t^{\alpha/2}-t^{4}\right)^{3}+\frac{9}{4}\Gamma(\alpha+1),$
(18) $\displaystyle t\in[0,1],\qquad y(0)=y^{\prime}(0)=0,$

having solution,

$y(t)=t^{8}-3t^{4+\alpha/2}+\frac{9}{4}t^{\alpha}.$

We consider the value $\alpha=1.3$. Figure 1 contains the obtained WPD, from
which one deduces that flmm2 can reach at most 10 digits of accuracy (further
decreasing of the stepsize only increases the execution time, whereas the
errors grow, due to round-off errors) in about 6 sec. On the other hand, fhbvm
is able to gain full machine accuracy, with execution times of the order of
one tenth of a second.

Figure 1: Work-precision diagram for problem (18), $\alpha=1.3$.

### Example 2

The next problem is taken from [21, Example 2]:

$\displaystyle y^{(3/2)}(t)$ $\displaystyle=$
$\displaystyle\frac{y(t)^{2}-(t^{1.9}-1)^{2}}{2}+\frac{\Gamma(2.9)}{\Gamma(1.4)}t^{0.4},\qquad
t\in[0,1],$ (19) $\displaystyle y(0)=-1,\quad y^{\prime}(0)=0,$

whose solution is $y(t)=t^{1.9}-1$. Figure 2 contains the obtained results,
from which one may see that flmm2 can reach at most 8 significant digits,
whereas fhbvm obtains 11 significant digits with a much smaller time.

Figure 2: Work-precision diagram for problem (19).

### Example 3

The third example is obtained by suitably combining Examples 2 and 3 in [18]:

$\displaystyle y^{(\alpha)}_{1}(t)$ $\displaystyle=$
$\displaystyle\frac{\Gamma(4+\alpha)}{6}t^{3}-t^{8+2\alpha}+y_{2}(t)^{2},$
$\displaystyle y^{(\alpha)}_{2}(t)$ $\displaystyle=$
$\displaystyle\frac{\Gamma(5+\alpha)}{24}t^{4}+t^{3+\alpha}-y_{1}(t),\qquad
t\in[0,1],$ (20) $\displaystyle
y_{1}(0)=y_{1}^{\prime}(0)=y_{2}(0)=y_{2}^{\prime}(0)=0,$

whose solution is given by

$y_{1}(t)=t^{3+\alpha},\qquad y_{2}(t)=t^{4+\alpha}.$

We consider the same value $\alpha=1.25$ considered in [18]. Figure 3
summarizes the obtained results, from which one may see that flmm2 reaches
about 12 significant digits in about 8 sec, while fhbvm obtains about 17
significant digits in less than one tenth of a second.

Figure 3: Work-precision diagram for problem (20), $\alpha=1.25$.

### Example 4

At last, we consider the following stiff problem, which is adapted from
Example 2 in [9]:

$y^{(\alpha)}(t)=\left(\begin{array}[]{rr}-100&0\\\
-99&-1\end{array}\right)y(t),\qquad t\in[0,20],\qquad
y(0)=\left(\begin{array}[]{c}2\\\ 3\end{array}\right).$ (21)

Its solution, for $\alpha\in(0,1)$, is given by:

$y(t)=\left(\begin{array}[]{c}2E_{\alpha}(-100t^{\alpha})\\\
2E_{\alpha}(-100t^{\alpha})+E_{\alpha}(-t^{\alpha})\end{array}\right),$

with $E_{\alpha}$ the one-parameter Mittag-Leffler function. We consider the
value $\alpha=0.25$. Figure 4 summarizes the obtained results, from which one
may see that flmm2 reaches less than 3 significant digits in about 60 sec,
whereas fhbvm reaches over 10 significant digits in less than 1 sec.

Figure 4: Work-precision diagram for problem (21), $\alpha=0.25$.

##  4 Conclusions

In this paper we have recalled the main facts about Fractional HBVMs (FHBVMs),
a class of numerical methods aimed at solving FDE-IVPs of Caputo type. In
particular, the method FHBVM$(22,20)$ has been recently implemented in the
Matlab${}^{\copyright}\,$code fhbvm [9], and is freely available at the URL
[22]. Numerical comparisons with another state-of-art code confirm the
potentialities of the approach and the efficiency of the code fhbvm. In
particular, the present version of the code has been extended to cope with
high-order FED-IVPs, w.r.t. the version presented in [9].

### Acknowledgements.

The authors are members of the Gruppo Nazionale Calcolo Scientifico-Istituto
Nazionale di Alta Matematica (GNCS-INdAM).

## References

  * [1]
  * [2] P. Amodio, L. Brugnano, F. Iavernaro. Spectrally accurate solutions of nonlinear fractional initial value problems. AIP Conf. Proc. 2116 (2019) 140005. https://doi.org/10.1063/1.5114132
  * [3] P. Amodio, L. Brugnano, F. Iavernaro. Analysis of Spectral Hamiltonian Boundary Value Methods (SHBVMs) for the numerical solution of ODE problems. Numer. Algorithms 83 (2020) 1489–1508. https://doi.org/10.1007/s11075-019-00733-7
  * [4] P. Amodio, L. Brugnano, F. Iavernaro. Arbitrarily high-order energy-conserving methods for Poisson problems. Numer. Algoritms 91 (2022) 861–894. https://doi.org/10.1007/s11075-022-01285-z
  * [5] P. Amodio, L. Brugnano, F. Iavernaro. A note on a stable algorithm for computing the fractional integrals of orthogonal polynomials. Applied Mathematics Letters 134 (2022) 108338. https://doi.org/10.1016/j.aml.2022.108338
  * [6] P. Amodio, L. Brugnano, F. Iavernaro. (Spectral) Chebyshev collocation methods for solving differential equations. Numer. Algoritms 93 (2023) 1613–1638. https://doi.org/10.1007/s11075-022-01482-w
  * [7] L. Brugnano, K. Burrage, P. Burrage, F. Iavernaro. A spectrally accurate step-by-step method for the numerical solution of fractional differential equations. J. Sci. Comput. 99 (2024) 48. https://doi.org/10.1007/s10915-024-02517-1
  * [8] L. Brugnano, G. Frasca-Caccia, F. Iavernaro, V. Vespri. A new framework for polynomial approximation to differential equations. Adv. Comput. Math. 48 (2022) 76 https://doi.org/10.1007/s10444-022-09992-w
  * [9] L. Brugnano, G. Gurioli, F. Iavernaro. Numerical solution of FDE-IVPs by using Fractional HBVMs: the fhbvm code. Numer. Algorithms (2024). https://doi.org/10.1007/s11075-024-01884-y
  * [10] L. Brugnano, F. Iavernaro. Line Integral Methods for Conservative Problems. Chapman et Hall/CRC, Boca Raton, FL, USA, 2016.
  * [11] L. Brugnano, F. Iavernaro. Line Integral Solution of Differential Problems. Axioms 7(2) (2018) 36. https://doi.org/10.3390/axioms7020036
  * [12] L. Brugnano, F. Iavernaro. A general framework for solving differential equations. Ann. Univ. Ferrara Sez. VII Sci. Mat. 68 (2022) 243–258. https://doi.org/10.1007/s11565-022-00409-6
  * [13] L. Brugnano, J.I. Montijano, F. Iavernaro, L. Randéz. Spectrally accurate space-time solution of Hamiltonian PDEs. Numer. Algorithms 81 (2019) 1183–1202. https://doi.org/10.1007/s11075-018-0586-z
  * [14] L. Brugnano, J.I. Montijano, L. Randéz. On the effectiveness of spectral methods for the numerical solution of multi-frequency highly-oscillatory Hamiltonian problems. Numer. Algorithms 81 (2019) 345–376. https://doi.org/10.1007/s11075-018-0552-9
  * [15] K. Diethelm. The analysis of fractional differential equations. An application-oriented exposition using differential operators of Caputo type. Lecture Notes in Math., 2004. Springer-Verlag, Berlin, 2010.
  * [16] K. Diethelm, R. Garrappa, M. Stynes. Good (and Not So Good) Practices in Computational Methods for Fractional Calculus. Mathematics 8(3) (2020) 324. https://doi.org/10.3390/math8030324
  * [17] R. Garrappa. Numerical solution of fractional differential equations: a survey and a software tutorial. Mathematics 6(2) (2018) 16. http://doi.org/10.3390/math6020016
  * [18] T.B. Nguyen, B. Jang. A high-order predictor-corrector method for solving nonlinear differential equations of fractional order. Fractional Calculus and Applied Analysis 20(2) (2017) 447–476. https://doi.org/10.1515/fca-2017-0023
  * [19] I. Podlubny. Fractional differential equations. An introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications. Academic Press, Inc., San Diego, CA, 1999.
  * [20] F. Mazzia, C. Magherini. Test Set for Initial Value Problem Solvers. Release 2.4, February 2008, Department of Mathematics, University of Bari and INdAM, Research Unit of Bari, Italy, available at: https://archimede.uniba.it/~testset/testsetivpsolvers/
  * [21] Z. Satmari. Iterative Bernstein spline technique applied to fractional order differential equations. Mathematical Foundations of Computing 6 (2023) 41–53. https://doi.org/10.3934/mfc.2021039
  * [22] https://people.dimai.unifi.it/brugnano/fhbvm/