11institutetext: Research Center for Information Technology Innovation,
Academia Sinica, Taipei, Taiwan 11email:
<EMAIL_ADDRESS>Department of Statistics, National
Cheng-Kung University, Tainan, Taiwan 33institutetext: Department of Applied
Mathematics National Sun Yat-sen University,
Kaohsiung, Taiwan 44institutetext: Department of Nuclear Medicine, Kaohsiung
Chang Gung Memorial Hospital,
Chang Gung University College of Medicine, Taiwan

# Dopamine Transporter SPECT Image Classification for Neurodegenerative
Parkinsonism via Diffusion Maps and Machine Learning Classifiers

Jun-En Ding(🖂) 11 Chi-Hsiang Chu 22 Mong-Na Lo Huang 33
Chien-Ching Hsu 44

###### Abstract

Neurodegenerative parkinsonism can be assessed by dopamine transporter single
photon emission computed tomography (DaT-SPECT). Although generating images is
time-consuming, these images can show interobserver variability and they have
been visually interprete by nuclear medicine physicians to date. Accordingly,
this study aims to provide an automatic and robust method based on Diffusion
Maps and machine learning classifiers to classify the SPECT images into two
types, namely Normal and Abnormal DaT-SPECT image groups. In comparison with
deep learning methods, our contribution is to propose an explainable diagnosis
process with high prediction accuracy. In the proposed method, the 3D images
of $N$ patients are mapped to an $N\times N$ pairwise distance matrix and are
visualized in Diffusion Maps coordinates. The images of the training set are
embedded into a low-dimensional space by using diffusion maps. Moreover, we
use Nystr$\ddot{o}$m’s out-of-sample extension, which embeds new sample points
as the testing set in the reduced space. Testing samples in the embedded space
are then classified into two types through the ensemble classifier with Linear
Discriminant Analysis (LDA) and voting procedure through twenty-five-fold
cross-validation results. The feasibility of the method is demonstrated via
Parkinsonism Progression Markers Initiative (PPMI) dataset of 1097 subjects
and a clinical cohort from Kaohsiung Chang Gung Memorial Hospital (KCGMH-TW)
of 630 patients. We compare performances using Diffusion Maps with those of
three alternative manifold methods for dimension reduction, namely Locally
Linear Embedding (LLE), Isomorphic Mapping Algorithm (Isomap), and Kernel
Principal Component Analysis (Kernel PCA). We also compare results using
through 2D and 3D CNN methods. The diffusion maps method has an average
accuracy of 98% from the PPMI and 90% from the KCGMH-TW dataset with twenty-
five fold cross-validation results. It outperforms the other three methods
concerning the overall accuracy and the robustness in the training and testing
samples.

###### Keywords:

Diffusion distance Diffusion maps Linear discriminant analysis Manifold
learning Nonlinear dimensionality reduction Parkinson’s disease

##  1 Introduction

Parkinson’s disease (PD) is a neurodegenerative disorder, and its pathological
feature is the loss of dopaminergic neurons in the substantia nigra [1].
Currently, the diagnosis of PD is mainly based on the clinical symptoms
(tremor at rest, rigidity, bradykinesia, gait disturbance) and the response to
medication (levodopa or dopamine agonist). Parkinsonian syndromes refer to
some diseases with clinical symptoms which is similar to PD. The etiologies of
Parkinsonian syndromes may also be related to the degeneration of dopamine
neurons (such as multiple system atrophy, progressive supranuclear palsy,
dementia with Lewy bodies, and corticobasal degeneration) or non-
neurodegenerative (such as essential tremor, secondary parkinsonism related to
hydrocephalus, stroke, drugs, toxins, trauma, brain tumor, or infection) [2,
3].

Dopamine transporters (DaT) are located on the presynaptic dopaminergic nerve
terminal and play an important role in regulating extracellular dopamine level
via reuptake dopamine into the nerve terminal. The loss of dopaminergic
neurons in neurodegenerative parkinsonism leads to the reduction of DaT. Many
radiotracers for single photon emission computed tomography (SPECT) such as
[99mTc]TRODAT-1 [4], 123I-beta-CIT, 123I FP-CIT, which can bind specifically
with DaT, were developed to diagnose PD and neurogenerative parkinsonism.

There are studies for diagnosis of Parkinson’s disease based on features
extracted from SPECT images, such as shape and volume etc., combining with
statistical test or machine learning classifiers [5, 10]. Although, features
as the shape of the abnormal striatum may be quite different and in certain
cases the features cannot be captured correctly, especially when most images
require high-dimensional analysis, the proposed analysis can be time consuming
and may have a large estimation error due to noisy image.

In recent years, due to rapid advancements in data storage and hardware
technology, neural network methods have become as a poweful technique for
prediction and image feature space extraction. One of the most well known
neural networks is the convolutional neural network (CNN). CNNs have been
widely applied for medical image analysis, including MRI and fMRI, well-known
CNN models using transfer learning architecture such as AlexNet[20],
VGG-16[21], VGG-19 and Deep Convolution Network (DCNN). The series of VGGNet
was first proposed VGG-16 in 2014-ILSVRC competition followed by VGG-19 as two
successful architectures on ImageNet. VGG-16 and VGG-19 use different
frameworks, and their models make an improvement on AlexNet by replacing large
kernel-sized filters with multiple small kernel-sized filters resulting in 13
and 16 convolution layers for VGG-16 and VGG-19 respectively. In addition,
more and more deep learning classifications of 123I SPECT scans use transfer
learning from deep neural networks pretrained on nonmedical images [23].

However, their results are often intuitive and require a large training set to
improve their classification accuracy. In general, biomedical images have a
high dimnesion feature space, which further complicates the classification
process. Consequently, it is desirable to perform dimension reduction to
improve the efficiency of the training and testing process. Many well-known
non-linear dimension reduction methods such that Kernel PCA, multidimensional
scaling (MDS), Isometric feature mapping, and locally linear embedding were
used to image classification and improve speed [6].

In comparison with previous studies, we use a manifold learning methodology,
namely the Diffusion Maps (DM) [7] to perform dimension reduction of the SPECT
image for classification purposes. Diffusion Maps do not require complicated
parameter adjustment and therefore have a shorter training time and require
fewer samples to capture more important information in low-dimensional space.
In addition to classification or early diagnosis and visualization, a general
standard specification for non-linear dimensionality reduction methods is
spectral decomposition. Through mapping the images into a low-dimensional
space with new coordinates, corresponding to the eigenspace associated with
the largest few eigenvalues, it is hoped that the points that have a
relationship between samples can be as close as possible after dimensionality
reduction, and keep the initial data structure. The method which refers to
Diffuison Maps hopes to find the geometric description of the corresponding
low-dimensional data through the diffusion process (Diffusion process). In
summary, we make three major contributions in this paper:

  * •

We propose a robust a method for images embedding associated with their main
features and together with the ensemble classifier.

  * •

We provide a diagnosis procedure for every patient twenty-five voting to
predict and construct an interpretable two-model confusion matrix.

  * •

We conduct extensive experiments on two real-world datasets. One is obtained
from the benchmark of Parkinson’s Progression Markers Initiative (PPMI)
database [9, 11], and the other is Taiwan clinical database. In our method, we
can use low-cost computing time than deep learning method and get better and
more robust results.

##  2 Datasets of PPMI and Clinical Cohort

###  2.1 PPMI Dataset

Data for this study were obtained from the PPMI database, a longitudinal,
multicentre study to assess the progression of clinical features, imaging, and
biologic markers on PD patients and healthy controls (HC). All the PD subjects
were at an early stage of Hoehn and Yahr stage I or II at baseline. The
diagnosis of PD was confirmed from the PPMI imaging core that the screening
DaT-SPECT (123I FP-CIT) is consistent with DaT deficit. Further details of the
PPMI database can be found at (http:// www.ppmi-info.org). The PPMI dataset in
our study contained 876 PD and 414 HC subjects.

We divide our KCGMH-TW samples into two classes (normal and abnormal). We
perform analyses with data based on a single ($1\times 128\times 128\times 3$)
(i.e., center image) or three ($3\times 128\times 128\times 3$) (i.e., center
three images) KCGMH-TW in our investigation. The PPMI ($1\times 109\times
91\times 3$) dataset has been divided into two classes (HC and PD).

###  2.2 Clinical Dataset of KCGMH-TW

Clinical dataset of KCGMH-TW in clinical dataset, we enrolled 730 patients who
underwent [99mTc] TRODAT-1 brain SPECT between January 2017 and June 2019 in
the Kaohsiung Chang Gung Memorial Hospital, Taiwan (KCGMH-TW). The Chang Gung
Medical Foundation Institutional Review Board approved this retrospective
study and waived the requirement for obtaining informed consent from the
patients. Each patient was intravenously injected with a 925-MBq dose of
[99mTc] TRODAT-1 (Institute of Nuclear Energy Research, Taiwan). [99mTc]
TRODAT SPECT images were acquired using a hybrid SPECT/CT system (Symbia T;
Siemens Medical Solution). SPECT images were obtained with 30s per step
acquiring 120 projections over a circular 360 degree rotation using low-
energy, high-resolution parallel-hole collimators. A $128\times 128$ matrix
and a $\times 1.45$ zoom were used. The CT images were acquired without
contrast medium; they used the following parameters: 130 kV; 45mAs (Image
Quality Reference mAs, CARE Dose 4D; Siemens Medical Solutions); rotation
time, 1.5s; collimation, $2\times 2.5$ mm. CT images were reconstructed to a
$512\times 512$ image matrix with a very smooth kernel, H08s (Siemens Medical
Solutions) for SPECT attenuation correction. Raw SPECT data were reconstructed
into transaxial slices using flash 3D (OSEM reconstruction method with 3D
collimator beam modeling) with 8 subsets and 8 iterations and corrected with
the H08s CT attenuation map. Images were smoothed using a 3D spatial Gaussian
filter (fullwidth at half maximum, 6mm). The reconstructed transaxial slice
thickness was 3.3 mm.

We select nine consecutive SPECT transaxial images showing the whole striatal
radioactivity. Each patient has a separate 9 slices SPECT images, where it is
more visible as shown in Fig. 1. The pixel size of each cell is $128\times
128\times 3$ (RGB), so we take the middle best three slices of striatal images
by nuclear medicine physicians [10] and combining into a single image from the
nine-grid SPECT as our research target.

Figure 1: Nine consecutive transaxial images of dopamine transporter single
photon emission computed tomography (DaT-SPECT) showing the whole striatal
radioactivity were displayed in a $3\times 3$ grid for visually interpreted.
(Left) The normal images of DaT-SPECT show symmetrical, well-delineated comma-
shaped radioactivity in the bilateral striata. (Right) The abnormal images of
DaT-SPECT show reduced radioactivity of the bilateral striata (nearly equal to
the background radioactivity of the brain). The putamen is usually more
severely affected than caudate nucleus resulting in a circular or oval shape.

###  2.3 Labeling Criterion

All DaT-SPECT images were visually interpreted by three experienced board-
certified nuclear medicine physicians according to Society of Nuclear Medicine
practice guideline [12]. The labeling criteria were established manually, (1)
Normal DaT-SPECT images: the normal striata on transaxial images should look
crescent- or comma-shaped and should have symmetric well-delineated borders.
(2) Abnormal DaT-SPECT images: the abnormal striata have reduced intensity on
one or both sides, often shrinking to a circular or oval shape. The putamen is
usually more severely affected than the caudate nucleus.

Blinded to patients’ clinical information except age and gender, three nuclear
medicine physicians visually interpreted the DaT-SPECT images independently.
The final consensus result of normal or abnormal image was assigned at least 2
physicians achieved an agreement. The PPMI dataset images ($1\times 109\times
91\times 3$) were divided into two classes (HC and PD). The clinical dataset
images of KCGMH-TW were divided into two classes (normal and abnormal). The
train/validation/testing dataset were summarized in Table. 1. We analyzed
KCGMH-TW dataset images based on single image ($1\times 128\times 128\times
3$) (i.e. central image in Fig. 1 and three images ($3\times 128\times
128\times 3$) (i.e. three images of middle row in Fig. 1).

Table 1: Dataset statistics | Dataset
---
(train/ validation /testing)
Two classes | Mean age
| PPMI
---
(872/225/193)
| HC (n=414)
---
PD (n=876)
-
| KCGMH-TW
---
(504/126/100)
| Normal (n=353)
---
Abnormal (n=377)
| 68.4
---
68.2

##  3 Methodology

In this study, we firstly propose the use of diffusion maps for dimension
reduction of the data, and then find the corresponding classification methods
for disease diagonsis. Accuracy of the proposed methodology is evaluated
through a cross-validation procedure for the training samples and later for
testing samples. For making a cross-validation study, we split a total of 630
original data samples to normal and abnormal group. Then we randomly divide
each group into five folds and cross combine them into 25 folds of paired
data. Later each fold of partition to five folds according to individual
normal and abnormal group. The divided data will be training cross combination
and validate until each combination sample has been trained. Our proposed
method applies each train set and conducts Nystrom’s out-of-sample expansion
[14] new sample to project training space Fig. 3, which reduced space will be
classified by optimal classifier in validation processing. Finally, we will
add a new sample and fixed test set (n=100) for prediction, and each patient
in this test set will have a total of twenty-five votes.

###  3.1 Training Sample Reduction via Diffusion Maps

Our framework uses a graph model treating data points (samples) as nodes
connected by edges with distances, defined by a weighting scheme with $w_{ij}$
denoting distance of node $i$ to node $j$, for all $i,j=1,\cdots,n$. Given $n$
data point set $\left\\{X_{i}\right\\}_{i=1}^{n}$, where
$X_{i}\in\mathbb{R}^{N\times N}$, $i=1,2,...,n$. The data point is embedded to
a manifold surfaces $\mathcal{M}$ in high-dimensional space and later is
mapped into a lower dimensional space through a diffusion maps. For
illustration let an undirected graph $G(X,W)$, in which $V$ and $W$ are the
set of nodes and weighted edges, respectively, be defined that each node is
connected by the adjacency matrix $W=(w_{ij})$, where $w_{ij}$ is the weight
of the edge connecting, i.e. $i$ and $j$.

Figure 2: Compute pairwise of SPECT features distance and then mapping to low-
dimensional space for image embedding.

In diffuison maps, the weights between node $i$ and $j$ is usually defined
through a kernel function $\mathcal{K}(X_{i},X_{j})$ such as the Gaussian [15]
kernel with $w_{ij}=\mathcal{K}(X_{i},X_{j})$ ,where $\mathcal{K}$ is a
positive and symmetric kernel matrix using Euclidean distance
$\left\|X_{i}-X_{j}\right\|^{2}$ for the $X_{i}$ and $X_{j}$ metric

$\mathcal{K}_{ij}=\mathcal{K}(X_{i},X_{j})=\exp(-\frac{\|X_{i}-X_{j}\|^{2}}{\alpha}).$
(1)

The $\alpha$ is a scale parameter which measures the $X_{i}$ and $X_{j}$ in
manifold neighborhood distance, we can construct a matrix with normalized rows
with unit length as

$P=D^{-1}\mathcal{K}.$ (2)

where $i$th element of degree matrix $D$ is computed by
$d_{ii}=\sum_{j=1}^{n}\mathcal{K}_{ij}$ and denoted as $D=$
$\textbf{diag}{(d_{11},d_{22},...,d_{nn})}$. The resulting matrix $P$ is
actually a Markov transition matrix with every entry to be nonnegative and has
all row sums to be equal to 1. The position of the data point $X_{i}$ and
$X_{j}$ is connected in manifold $\mathcal{M}$ by weights through where
$p(X_{i},X_{j})$ is normalized by $\mathcal{K}_{ij}/d(x_{i})$, i.e.

$P=(p(X_{i},X_{j})).$ (3)

If points $X_{i}$ and $X_{j}$ are similar, the distance between the two points
is closer, and there are higher probability jumping from node $i$ to nearby
node $j$. The similarity of the two data points can be evaluated by the
following distance measure

$\displaystyle D^{2}(X_{i},X_{j})$ $\displaystyle=\sum_{u\in
X}\|p(X_{i},u)-p(X_{j},u)\|^{2}$ (4)
$\displaystyle=\sum_{k}\|p(X_{i},u)-p(X_{j},u)\|^{2}$
$\displaystyle=\sum_{k}\|P_{ik}-P_{kj}\|^{2}.$

And let $Y_{i}$ be the $n\times 1$ distance vector with elements composed of
distances from node $i$ to all nodes $j=1,\dots,n$ (including itself), mapped
from the original data point $X_{i},i=1,\dots,n$, namely

$Y_{i}:=\begin{bmatrix}p(X_{i},X_{1})\\\ p(X_{i},X_{2})\\\ \vdots\\\
p(X_{i},X_{n})\end{bmatrix}=P_{i\cdot}^{T\quad}$ (5)

Note that the Euclidean distance between $Y_{i}$ and $Y_{j}$ can be expressed
as

$\displaystyle\|Y_{i}-Y_{j}\|_{E}$
$\displaystyle=\sum_{k}\|p(X_{i},u)-p(X_{j},u)\|^{2}$
$\displaystyle=\sum_{k}\|P_{ik}-P_{kj}\|^{2}$
$\displaystyle=D(X_{i},X_{j})^{2}.$

###  3.2 Nystr$\ddot{o}$m’s Out-Of-Sample Extension

Manifold learning methods usually require recalculate the kernel matrix with
new added samples, and the above steps are repeated for the entire dataset. It
becomes more difficult when new sample points are added sequentially.
Nystr$\ddot{o}$m’s out-of-sample extension [16] allows the original sample to
be extended by adding new samples and embedding them into the existing low-
dimensional space to form the manifold geometry of the new sample points. Fig.
3 demonstrates the procedure of using out-of-sample expansion, especially when
it comes to do cross-validation analysis. The benefits of Nystrom’s out-of-
sample procedure extending the new samples are illustrated as follows: (1) The
compliance of machine learning is divided into training set and test set
standard verification steps; (2) When new samples are added, it is not
necessary to repeatedly calculate diffusion distance of whole Markov Matrix
also solves the problem of computation time; (3) Maintain the training set
geometric structure from the original samples.

Figure 3: Out-of-sample in high dimensional data projection.

We can extend new sample point $X_{new}\in\mathbb{R}^{N\times N}$ from the
validation set, and then recalculate the Euclidean distance and Kernel matrix,
the new element of Kernel matrix has an augmented vector
$\mathcal{K}_{new}^{T}=(\mathcal{K}_{new,1},\dots,\mathcal{K}_{new,n})$ with

$\mathcal{K}_{new,j}=\mathcal{K}(X_{new},X_{j})=\exp(-\frac{\|X_{new}-X_{j}\|^{2}}{\alpha}).$
(6)

The augmented vector ${P}_{new}^{T}=D_{new}^{-1}\mathcal{K}_{new}^{T}$ is
normalized vector of $\mathcal{K}_{new}^{T}$ by dividing its row sum
$D_{new}=\sum_{j}^{n}\mathcal{K}_{new,j}$. So the augmented Markov transition
matrix $P_{n+1}$ can be rewritten as

$P_{n+1}=\begin{pmatrix}P_{n}\\\
P_{new}^{T}\end{pmatrix},\mathcal{K}_{n+1}=\begin{pmatrix}\mathcal{K}_{n}\\\
\mathcal{K}_{new}^{T}\end{pmatrix}$ (7)

Finally, project a new sample to the diffusion space of chosen $k$ eigenvector
above by

$\psi_{l}(X_{new})=\frac{1}{\lambda_{l}}\sum_{j=1}^{n}p(X_{new},X_{j})\psi_{l}(X_{j}),l=1,...,k.$
(8)

Figure 4: The process of the high-dimensional SPECT data reduction uses out-
of-sample porjects for classifications.

Each point for each patient is projected onto the original trained low-
dimensional diffusion space, and it goes on to use the selected classifier to
carry on the classification as shown in Fig. 4. The algorithm for constructing
the corresponding diffusion maps and out-of-sample extension is presented in
Algorithm. 1 below.

Algorithm 1 Diffusion Maps and out-of-sample extension algorithm.

0: Data of images $\\{X_{i}\\}_{i=1}^{n}\in\mathbb{R}^{N\times N}$;

0: Projected new coordinate vector $Y^{\prime}_{new}\in\mathbb{R}^{N\times
k}$;

1: Normalized data: $\tilde{X}_{i},i=1,2,..,n.$
$\tilde{X}_{i}=(X_{i}-\bar{X})/S(X_{i})$, which
$\bar{X}_{i}=\frac{1}{N}\sum_{j,k}X_{i}(j,k)$$S({X}_{i})=\frac{1}{N^{2}-1}\sum_{j,k}(X_{i}(j,k)-\bar{X}_{i}(j,k))^{2}$;

2: Construct the Kernel
matrix:$K=(K_{i,j}),K_{i,j}=K(\tilde{X}_{i},\tilde{X}_{j})$;

3: Build the Markov transitoin
matrix:$P=D^{-1}K$,$D=\textbf{diag}{(d_{11},d_{22},...,d_{nn})}$, where
$d_{ii}=\sum_{j=1}^{n}K_{ij}$;

4: Spectral decomposition $P$ matrix with corresponding eigenpairs
$\\{\lambda_{j},\psi_{j}\\}_{j=1}^{n}$ as coordinate:$Y_{i}^{\prime}$ T
$:=(\lambda_{1}\psi_{1}(i),\lambda_{2}\psi_{2}(i)...,\lambda_{k}\psi_{k}(i))$

5: Out-of-sample extension

6: Compute new extension vectors:$P_{new}^{T}=d_{new}^{-1}K_{new}^{T}$, where
$K_{new}^{T}$= $(K_{new,1},...,K_{new,n})$;

7: Repeat Step 1 to Step 4 and compute extented Markov
matrix:$P_{n+1}^{T}=(P_{n},P_{new}^{T})$, where
$K_{n+1}^{T}=(K_{n},K_{new}^{T})$;

8: Eigenvecotrs of the coordinates projected from original
space:$\psi_{l}(X_{new})=\frac{1}{\lambda_{l}}\sum_{j=1}^{k}p(X_{new},X_{j})\psi_{l}(X_{j}),l=1,...,k$;

9: return $\psi_{l}(X_{new}),l=1,...,k$;

Figure 5: The workflow of the two steps model diagnosis architecture with
twenty-five folds cross-validation. (A) Procedure for nonlinear dimension
reduction through Manifold method and diagnosis with given classifier. (B)
Procedure for new test samples embedding and diagnoses with DM and LLE as well
as the corresponding two model confusion matrix.

##  4 Experiments

In this work, we implement the same data preprocessing procedure in our
clinical data KCGMH-TW and PPMI dataset. Our approach mainly consists of two
steps as demonstrated in Fig. 5:

  * •

First step: several manifold learning methods are used for dimension reduction
on the training sets for low-dimensional embedding, and then out-of-sample
extension is applied to the new samples of the validation sets[19]. Next,
classify the embedded samples and aggregate the ensemble classifiers with the
twenty-five voting results to making prediction on the corresponding test set.

  * •

Second step: according to the results of the first step, we use the DM and the
classifier with the best performance to predict new test samples. Each patient
is diagonalized through the twenty-five classification voting procedure and
the prediction results with different manifold learning methods can be
compared through the two model confusion matrices.

In the first step, three classifiers are compared, namely the: LDA, SVM,
Logistic regression and the classifier with the best performance is chosen to
enter into the second step process. In this work, the DM parameter is set to
be $\alpha=8$ and time step $t=1$ in Markov transition matrix $P^{(t)}$. In
out-of-sample step, there is no need to recompute the entire diffusion matrix,
simply project the new samples to the original chosen DM space which will be
used as the feature space for diagonsis analysis.

###  4.1 Two Steps Model Ensemble and Classifer Selection

In order to compare diagnosis results using DM with those by manifold learning
methods such as LLE, ISOMP, KPCA, and classifiers such as LDA, SVM, Logistic
Regression in another scenario. We compare the performances of DM, LLE,
ISOMAP, and KPCA diagonalized results based on the twenty-five folds
classifications under four low-dimensionality (30, 100, 200 and 300) with
different classifiers. Overall accuracy averages based on the twenty-five
folds predictions after voting with the three classifiers are presented in
Table. 2 . It can be seen that (DM, LDA) pair with 200 dimensionality seems to
be the best. An interesting observation on the accuracy of the four methods,
is that after dimension reduction, the linear discriminant classifier (LDA)
works better than the other two non-linear classifiers. The performances of
the DM and LLE are superior to the other two methods, therefore in the
following we use DM and LLE with LDA for comparisons of the performances for
testing the new 100 samples.

Table 2: Performance comparisons of the four manifold learning methodology and three ensemble classifiers on KCGMH-TW  Method | LDA | SVM | Logistic
---|---|---|---
Dimension | Acc. | Sens. | Spec. | Acc. | Sens. | Spec. | Acc. | Sens. | Spec.
DM | 30 | 0.88 | 0.93 | 0.82 | 0.86 | 0.93 | 0.79 | 0.85 | 0.90 | 0.80
100 | 0.89 | 0.94 | 0.84 | 0.86 | 0.89 | 0.83 | 0.88 | 0.90 | 0.86
200 | 0.90 | 0.91 | 0.91 | 0.85 | 0.86 | 0.84 | 0.87 | 0.88 | 0.86
300 | 0.91 | 0.91 | 0.91 | 0.87 | 0.92 | 0.82 | 0.88 | 0.92 | 0.84
LLE | 30 | 0.84 | 0.87 | 0.81 | 0.83 | 0.88 | 0.77 | 0.81 | 0.84 | 0.79
100 | 0.89 | 0.96 | 0.81 | 0.83 | 0.88 | 0.77 | 0.81 | 0.84 | 0.79
200 | 0.89 | 0.93 | 0.83 | 0.87 | 0.93 | 0.81 | 0.85 | 0.87 | 0.82
300 | 0.89 | 0.93 | 0.83 | 0.87 | 0.81 | 0.93 | 0.85 | 0.87 | 0.82
ISOMAP | 30 | 0.80 | 0.81 | 0.79 | 0.81 | 0.83 | 0.79 | 0.79 | 0.77 | 0.81
100 | 0.86 | 0.89 | 0.84 | 0.85 | 0.86 | 0.84 | 0.85 | 0.84 | 0.85
200 | 0.87 | 0.89 | 0.85 | 0.86 | 0.86 | 0.85 | 0.87 | 0.87 | 0.87
300 | 0.87 | 0.89 | 0.85 | 0.86 | 0.86 | 0.85 | 0.87 | 0.87 | 0.87
KPCA | 30 | 0.79 | 0.80 | 0.79 | 0.76 | 0.80 | 0.71 | 0.75 | 0.85 | 0.70
100 | 0.83 | 0.83 | 0.82 | 0.77 | 0.82 | 0.71 | 0.77 | 0.83 | 0.70
200 | 0.86 | 0.90 | 0.81 | 0.79 | 0.85 | 0.72 | 0.78 | 0.84 | 0.71
300 | 0.87 | 0.93 | 0.81 | 0.80 | 0.89 | 0.73 | 0.78 | 0.84 | 0.71

For each test patient, based on the two-step procedure above, there are
twenty-five ensemble votes for prediction, then we can calculate the
proportion of machine predicted abnormal probability:
$p_{k}=\frac{V_{k}}{\sum_{k=1}V_{k}}$, with $V_{k}$ denoting the $k$th voting
results, where $V_{k}=0$ or $1$ representing normal or abnormal respectively,
$k=1,\ldots,25$. So we build a threshold to (0,1)
_e.g.,_($\mathbbm{1}_{\\{p_{k}>0.5\\}}$) for testing the voting prediction.

###  4.2 Classification

We compare the performances of the two-step manifold approach with 2D-CNN,
3D-CNN on KCGMH-TW and PPMI datasets in twenty-five folds predict. And then we
choose the best performing 200 dimensions as a baseline and compare different
well-known CNN models such as AlexNet[20], VGG-16[21], VGG-19 and Deep
Convolution Network (DCNN). We implement the KCGMH-TWdataset with three SPECT
images from $9\times 9$ grids as a high-order brain image tensor as input to
the 3D-Convolutional Neural Network model[22]. The volumes of SPECT have a
size of $3\times 128\times 128\times 3$ voxels. For the CNN approach, we set
193 images as the test set and 100 images in the manifold method test in
Table. LABEL:Group_table.

Evaluate out-of-sample extension to new samples, we use Frobenius norm
distance to measure testing sample point embedding on low-dimensional space.
In KCGMH-TW (three images) case, three images show that between the training
set over to embed testing set significant effect classification accuracy like
a Table. LABEL:Group_table.

We consider the early abnormal in our KCGMH-TW for training and testing
(testing = 100) and implement our procedure in the PPMI dataset (testing =
197). The classification results reported in Table. LABEL:Group_table. DM+LDA
method generally achieves the best performance on two datasets with single and
three types. Compare to our approach methods, 2D-CNN often have excellent
prediction than 3D-DCNN such that VGG-16 in PPMI 93% and DCNN in KCGMH-TW
single 87% cases but also has higher variability of 0.13 and 0.22. The
comparison of our method DM+LDA in PPMI single image and KCGMH-TW three images
classification has average accuracy 98% and 90% performance.

Table 3: Binary classification in different datasets Datasets | Single image | Model | Acc. | Prec. | Senc. | Spec. |
---|---|---|---|---|---|---|---
PPMI | HC vs PD | DCNN |  0.89 ($\pm 0.54$) | 0.92 | 0.85 | 0.90 |
|  | AlexNet |  0.88 ($\pm 0.27$) | 0.89 | 0.78 | 0.92 |
|  | VGG-16 |  0.93 ($\pm 0.13$) | 0.94 | 0.95 | 0.92 |
|  | VGG-19 |  0.89 ($\pm 0.21$) | 0.91 | 0.93 | 0.86 |
|  | Diffusion Maps + LDA |  0.98 ($\pm 0.02$) | 0.98 | 0.96 | 0.97 |
KCGMH-TW | Normal vs Abnormal | DCNN |  0.87 ($\pm 0.22$) | 0.88 | 0.89 | 0.84 |
|  | AlexNet |  0.83 ($\pm 0.35$) | 0.86 | 0.84 | 0.80 |
|  | VGG 16 |  0.85 ($\pm 0.44$) | 0.84 | 0.89 | 0.76 |
|  | VGG 19 |  0.87 ($\pm 0.11$) | 0.88 | 0.92 | 0.81 |
|  | Diffusion Maps + LDA |  0.86 ($\pm 0.04$) | 0.85 | 0.93 | 0.77 |
Datasets | Three images | Model | Acc. | Prec. | Senc. | Spec. |
KCGMH-TW | Normal vs Abnormal | 3D-DCNN |  0.82 ($\pm 0.10$) | 0.83 | 0.84 | 0.79 |
|  | Diffusion Maps + LDA |  0.90 ($\pm 0.05$) | 0.88 | 0.95 | 0.84 |

From the overall results in Fig. 6, we draw DM+LDA and CNN methods in twenty-
five folds situation. In Fig. 6(a)(b), display that we only use a single PPMI
image to predict 193 patients, DM+LDA has robust performance than VGG-16,
VGG-19, AlexNet, and DCNN. The PPMI dataset has high image quality compared to
our clinical KCGMH-TW dataset. To reduce effects of the noise in our KCGMH-TW,
we also use three images for DM+LDA for comparisons. We find that DM+LDA
procedure has more robust results with three images, see Fig. 6(c) and Fig.
6(d). It may be due to that three images have included the most obvious
symmetrical strata in high order tensor space.

(a) PPMI (single image) (b) KCGMH-TW (single image)

Figure 6: The performace of DM and CNN comparison on PPMI and KCGMH-TW
dataset. Figure 7: Performance of DM in single image, three images and 3D-DCNN
methods

##  5 Diagnosis and Discussion

In this section, we focus on machine prediction on the 100 test patients and
diagnosis of misclassification samples. Our diagnostic framework as exhibited
in Fig. 5 (B), provides a general interpretable misclassification table for
diagnosis in the final step with the confusion table. It is helpful to see
whether an incorrect diagnosis is dut to the patient’s age or not.

###  5.1 Two Model Confusion Matrices

From Section 4.1 , it is obvious that DM and LDA classifier have high
performance in classification after the twenty-five folds voting. As LLE also
has quite good performance with the twenty-five folds procedure, we examine
the diagnosis results based on the two models and form a two-model confusion
matrix. First, we set the doctor label as true label (i.e., two classes) and
the ensemble DM and LLE diagnoses in two classes as in Fig. 8(a) and then
check how many misclassifications under the true label normal or abnrmal.

The diagonal elements of the two model confusion matrix table are the numbers
of correct predictions (positive is with label 1, negative is with 0) for the
DM and LLE ensembles respectively. As for the numbers of misclassified samples
are shown in the non-diagonal sub-table. For example in the left table in Fig.
8(b)  left column is ture label of normal, DM and LLE predict ID:678 patient
is early abnormal. In the right column Fig. 8(b) , ID:677 is predicted by LLE
and abnormal is predicted incorrectly as normal by LLE and correctly as
abnormal by DM, while ID:678 is predicted incorrectly as normal by both
models. In the left column of the table Fig. 8(c)  as abnormal, three of the
patients with true label as normal are diagnosed correctly as normal (ID: 657,
674, 714 ) by DM and incorrectly as abnormal.

Figure 8: Two-model confusion matrices

Diagnosis Table: Sometimes we care more about false negative situation as the
seven cases as indicated in Fig. 8(d). For these misclassifications, we can
provide more detailed information for the doctor’s diagnosis table and the
probability of LDA matching as demonstrated in Table. 9, in order to identify
early abnormal in the true designation of three classes. We can also check if
our model predicts an error due to high age or shrinking striatum abnormal.
For example, in the case of the misclassification ID:709, whose age is 74.4
higher than the average age of our dataset of 68.3.

Figure 9: The diagnosis table for the seven false negative samples on their
age and proportion of the twenty-fold voting results.

###  5.2 Visualization

Traditional analysis of the SPECT image is for the physician to visualize the
symmetry of the left and right striatum. However, the early abnormal symptoms
is more controversial. To find mild changes of Parkinson’s disease, we use DM
for dimension reduction on 630 original SPECT images in low-dimensional space,
and visualize the embedding test samples in the corresponding three-
dimensional manifold with the first three eigenvector
$\psi_{1},\psi_{2},\psi_{3}$ with the largest three eigenvalues. The are shown
to be like an U curve with abnormal samples appeared on the right side of the
curve as shown in Fig. 10. The color from dark purple to dark red indicates
the distance between sample points of each patient. Moreover, we label the
true diagnosis of each patient (two classes) to examine the spatial status.
The data point (ID:342) with the largest distance from the origin on the right
extreme is the most serious abnormal case. In contrast, the point (ID:250) on
the left extreme is normal, and the early abnormal appears in the middle
position, as ID: 18.

Figure 10: Visualize of SPECT in different position and annotate the true
label on each data points, the color present the distance of DM sorted by
eigenvectors mapping.

##  6 Conclusion

In addition to maintaining excellent geometric structure in low dimensions, DM
is computationally less expensive than deep learning methods and does not
require too many parameter adjustments. The average accuracy of DM +LDA
classification with significant effect is up to 90% for normal vs. abnormal
(three images) and up to 98% for PPMI dataset (single). The overall
performance on KCGMH-TW is 90% and the variation is lower for our methods.
Using the confusion matrix of the two models, we can detect the early abnormal
samples with high accuracy and therefore have provided a complete workflow for
diagnosis.

## References

  * [1] Kalia, Lorraine V. and Lang, Anthony E.,:Parkinson’s disease, The Lancet , 9996, 896–912, Elsevier Ltd, 386 (2015). https://doi.org/10.1016/S0140-6736(14)61393-3
  * [2] Keener AM, Bordelon YM.: Parkinsonism. Semin Neurol. 36(4), 330-4 (2016). https://doi.org/10.1055/s-0036-1585097
  * [3] Hayes MT.: Parkinson’s Disease and Parkinsonism. Am J Med. 132(7), 802-807 (2019). https://doi.org/10.1016/j.amjmed.2019.03.001
  * [4] Kung, Hank F. and Kung, Mei Ping and Wey, Shiaw Pyng and Lin, Kun Ju and Yen, Tzu Chen.: Clinical acceptance of a molecular imaging agent: a long march with [99mTc]TRODAT, Nuclear Medicine and Biology. 34(7), 787-9 (2007) 132(7), 802-807 (2019). https://doi.org/10.1016/j.nucmedbio.2007.03.010
  * [5] Prashanth, R and Dutta Roy, Sumantra and Mandal, Pravat K and Ghosh, Shantanu.: High-Accuracy Classification of Parkinson’s Disease Through Shape Analysis and Surface Fitting in 123I-Ioflupane SPECT Imaging. IEEE Journal of Biomedical and Health Informatics. 21, 794–802 (2016) 132(7), 802-807 (2019). https://doi.org/10.1109/jbhi.2016.2547901
  * [6] Faaeq, Ainuddin and Guruler, Huseyin and Peker, Musa.: Image classification using manifold learning based non-linear dimensionality reduction. 26th IEEE Signal Processing and Communications Applications Conference. 1-4 (2018) 132(7), 802-807 (2019). https://doi.org/10.1109/SIU.2018.8404441
  * [7] Coifman, R R and Lafon, S and Lee, A B and Maggioni, M and Warner, F and Zucker, S.: Geometric diffusions as a tool for harmonic analysis and structure definition of data: Diffusion maps, Proceedings of the National Academy of Sciences. 102(21), 7426-7431, (2005) https://doi.org/10.1073/pnas.0500334102
  * [8] Belkin, Mikhail and Niyogi, Partha.: Laplacian eigenmaps for dimensionality reduction and data representation. Neural Computation. 15, 1373 (2003) https://doi.org/10.1162/089976603321780317
  * [9] Cummings, Jeffrey L. and HencNormalliffe, Claire and Schaier, Sharon and Simuni, Tanya and Waxman, Alan and Kemp, Paul.: The role of dopaminergic imaging in patients with symptoms of dopaminergic system neurodegeneration. Brain 134(11), 3146-66 (2011) https://doi.org/10.1093/brain/awr177
  * [10] A. Faro and D. Giordano and C. Spampinato and S. Ullo and A. Di Stefano.: Basal Ganglia Activity Measurement by Automatic 3-D Striatum Segmentation in SPECT Images. IEEE Transactions on Instrumentation and Measurement. 60(10), 3269-3280 (2011) https://doi.org/10.1109/TIM.2011.2159315
  * [11] Quan, Justin and Xu, Lin and Xu, Rene and Tong, Tyrael and Su, Jean.: DaTscan SPECT image classification for Parkinson’s disease. arXiv. 1-9 (2019) https://doi.org/arXiv:1909.04142v1
  * [12] Djang, David S.W. and Janssen, Marcel J.R. and Bohnen, Nicolaas and Booij, Jan and Henderson, Theodore A. and Herholz, Karl and Minoshima, Satoshi and Rowe, Christopher C. and Sabri, Osama and Seibyl, John and Van Berckel.: SNM practice guideline for dopamine transporter imaging with 123I-ioflupane SPECT 1.0. Journal of Nuclear Medicine. 53, 154-163 (2012) https://doi.org/10.2967/jnumed.111.100784
  * [13] Shimizu, Soichiro and Namioka, Nayuta and Hirose, Daisuke and Kanetaka, Hidekazu and Hirao, Kentaro and Hatanaka, Hirokuni and Takenoshita, Naoto and Kaneko, Yoshitsugu and Ogawa, Yusuke and Tsugawa, Akito and Umahara, Takahiko and Sakurai, Hirofumi and Hanyu, Haruo.: Comparison of diagnostic utility of semi-quantitative analysis for DAT-SPECT for distinguishing DLB from AD, Journal of the Neurological Sciences. 377, 50-54 (2017) https://doi.org/10.1016/j.jns.2017.03.040
  * [14] ScNormallar, Alon and Rokach, Lior and Amit, Amir.: Diffusion ensemble classifiers, IJCCI 2012 - Proceedings of the 4th International Joint Conference on Computational Intelligence. 443-450 (2012) https://doi.org/10.5220/0004102804430450
  * [15] Porte, J De La and Herbst, B M and Hereman, W and Walt, S J Van Der.: An introduction to diffusion maps, In The 19th Symposium of the Pattern Recognition Association of South Africa. (2008)
  * [16] Bengio, Yoshua and Paiement, Jean-François and Vincent, Pascal and Delalleau, Olivier and Roux, Nicolas Le and Ouimet, Marie.: Out-of-Sample Extensions for LLE, Isomap, MDS, Eigenmaps, and Spectral Clustering. In Advances in Neural Information Processing Systems.
  * [17] Fernández, Ángela and González, Ana M. and Díaz, Julia and Dorronsoro, José R.: Diffusion Maps for dimensionality reduction and visualization of meteorological data. Neurocomputing. 163, 25-37 (2015) https://doi.org/10.1016/j.neucom.2014.08.090
  * [18] Prashanth, R. and Dutta Roy, Sumantra and Mandal, Pravat K. and Ghosh, Shantanu.: Automatic classification and prediction models for early Parkinson’s disease diagnosis from SPECT imaging. Expert Systems with Applications. 41, 3333-3342 (2014) https://doi.org/10.1016/j.eswa.2013.11.031
  * [19] Raeper, Rory and Lisowska, Anna and Rekik, Islem.: Joint correlational and discriminative ensemble classifier learning for dementia stratification using shallow brain multiplexes. MICCAI, 11070 LNCS, 599-607, (2018)
  * [20] Krizhevsky.: Imagenet classification with deep convolutional neural networks, Alex and Sutskever, Ilya and Hinton, Geoffrey E, Neural Information Processing Systems. 25 2012 https://doi.org/10.1145/3065386
  * [21] Simonyan, Karen and Zisserman, Andrew.: Very deep convolutional networks for large-scale image recognition. ICLR. 1-14, 2015 https://arxiv.org/abs/1409.1556
  * [22] Esmaeilzadeh, S and Yang, Y and Adeli, E.: End-to-End Parkinson Disease Diagnosis using Brain MR-Images by 3D-CNN, Arxiv, 2018 https://arxiv.org/abs/1806.05233
  * [23] Kim DH, Wit H, Thurston M.: Artificial intelligence in the diagnosis of Parkinson’s disease from ioflupane-123 single-photon emission computed tomography dopamine transporter scans using transfer learning. Nucl Med Commun. 39(10), 887-893 (2018) https://doi.org/10.1097/MNM.0000000000000890


# Evaluating the reliability of acoustic speech embeddings

###### Abstract

Speech embeddings are fixed-size acoustic representations of variable-length
speech sequences. They are increasingly used for a variety of tasks ranging
from information retrieval to unsupervised term discovery and speech
segmentation. However, there is currently no clear methodology to compare or
optimise the quality of these embeddings in a task-neutral way. Here, we
systematically compare two popular metrics, ABX discrimination and Mean
Average Precision (MAP), on 5 languages across 17 embedding methods, ranging
from supervised to fully unsupervised, and using different loss functions
(autoencoders, correspondence autoencoders, siamese). Then we use the ABX and
MAP to predict performances on a new downstream task: the unsupervised
estimation of the frequencies of speech segments in a given corpus. We find
that overall, ABX and MAP correlate with one another and with frequency
estimation. However, substantial discrepancies appear in the fine-grained
distinctions across languages and/or embedding methods. This makes it
unrealistic at present to propose a task-independent silver bullet method for
computing the intrinsic quality of speech embeddings. There is a need for more
detailed analysis of the metrics currently used to evaluate such embeddings.

Index Terms: unsupervised speech processing, speech embeddings, frequency
estimation, evaluation metrics, representation learning, $k$-nearest
neighbours

##  1 Introduction

Unsupervised representation learning is the area of research that aims to
extract units from unlabelled speech that are consistent with the phonemic
transcription [1, 2, 3]. As opposed to text, speech is subject to large
variability. Two speech sequences with the same transcription can have
significantly different raw speech signals. In order to work on speech
sequences in an unsupervised way, there is a need for robust acoustic
representations. To address that challenge, recent methods use speech
embeddings, i.e. fixed-size representations of variable-length speech
sequences [4, 5, 6, 7, 8, 9, 10, 11].111A speech sequence is a non-silent part
of the speech signal (not necessarily a word). It can be transcribed into a
phoneme $n$-gram.

Speech embeddings can be used in many applications, such as key-word
spotting[12, 13, 14], spoken term discovery[15, 16, 17], and segmentation of
speech into words [18, 19, 20]. It is convenient to evaluate the reliability
of speech embeddings without being tied to a particular downstream task. One
way to do that is to compute the intrinsic quality of speech embeddings. The
basic idea is that a reliable speech embedding should maximise the information
relevant to its type and minimise irrelevant token-specific information. Two
popular metrics have been used: the mean average precision (MAP) [21] and the
ABX discrimination score [22].

ABX and MAP are mathematically distinct yet they are expected to correlate
well with each other as they both evaluate the discriminability of speech
embeddings in terms of their transcription. However, [5] revealed a surprising
result: the best model according to the ABX, is also the worst one according
to the MAP. Following [5]’s results, we observed that this kind of
discrepancies is much more common than we had expected. If a model performs
well according to the MAP and bad according to the ABX, which metric should be
trusted? For research in this field to go forward, there is a need to quantify
the correlation of these two metrics.

In this paper, we wanted to go further and check that MAP and ABX can also
predict performances on a downstream task. Such tasks are numerous, but one of
them has not yet received enough interest: the unsupervised frequency
estimation. We define the frequency of a speech sequence as the number of
times the phonetic transcription of this sequence appears in the corpus. When
dealing with text corpora, frequencies can be computed exactly with a lookup
table and are used in many NLP applications. In the absence of labels,
deriving the frequency of a speech sequence becomes a problem of density
estimation. Estimated frequencies can be useful in representation learning by
enabling efficient sampling of tokens in a speech database [7]. Also,
frequencies could be used for the unsupervised word segmentation using
algorithms similar to those used in text [18].

In Section 2, we present the range of embedding models that can be grouped in
five categories of increasing expected reliability: hand-crafted,
unsupervised, self-supervised, supervised plus a top-line embedding. In
Section 3, we present the MAP and ABX metrics and introduce our frequency
estimation task. In Section 4, we present results on the five speech datasets
from the ZeroSpeech Challenge [2, 3]. From these results, we draw guidelines
for future improvements in the field of acoustic speech embeddings.

##  2 Embedding Methods

###  2.1 Acoustic features

Neural networks learn representations on top of input features. Therefore we
used two types of acoustic features known as the log-MEL filterbanks
(Mel-F)[23] and the Perceptual Linear Prediction (PLP)[24]. These two features
can be considered as two levels of phonetic abstraction: a high-level one
(PLP) and a low-level one (Mel-F). Formally, let us define a speech sequence
$s_{t}$ by $x_{1}$,$x_{2}$,…,$x_{T}$, where $x_{i}\in\mathbb{R}^{n}$ is called
a frame of the acoustic features. $T$ is the number frames in the sequence
$s_{t}$. In our setting, these frames are spaced out every 10 ms each
representing a 25 ms span of the raw signal.

###  2.2 Hand-crafted model: Gaussian downsampling

Holzenberger and al. [5] described a method to create fixed-size embedding
vectors that requires no training of neural networks: the Gaussian down-
sampling (GD). Given a sequence $s_{t}$, $l$ equidistant frames are sampled
and a Gaussian average is computed around each sample. It returns an embedding
vector $e_{t}$ of size $l\times n$ for any size $T$ of input sequences.
Therefore, given our two acoustic features, two baselines model are derived:
the Gaussian-down-sampling-PLP (GD-PLP) and the Gaussian-down-sampling-Mel-F
(GD-Mel-F).

Similarly, we derived a simple top-line model. Instead of using hand-crafted
features, we can use the transcription of a given random segment. Each frame
$x_{i}$ in a sequence $s_{t}$ will be assigned a 1-hot vector referring
directly to the phoneme being said. This model goes through the same Gaussian
averaging process to form the Gaussian-down-sampling-1hot (GD-1hot) model.
This model is almost the true labels notwithstanding the information loss due
to compression.

###  2.3 Unsupervised model: RNNAE

A more elaborate way to create speech embeddings is to learn them on top of
acoustic features using neural networks. Specifically, recurrent neural
networks (RNN) can be trained with back-propagation in an auto-encoding (AE)
objective: the RNNAE [5, 9]. Formally, the model is composed of an encoder
network, a decoder network and a speaker encoder network. The encoder maps
$s_{t}$ to $e_{t}$, a fixed-size vector. The speaker encoder maps the speaker
identity to a fixed size vector $spk_{t}$. Then, the decoder concatenate
$e_{t}$ and $spk_{t}$ and maps them to $\hat{s}_{t}$, a reconstruction of
$s_{t}$. The three networks are trained jointly to minimise the Mean Square
Error between $\hat{s}_{t}$ and $s_{t}$.

###  2.4 Self-supervised and supervised models: CAE, Siamese and CAE-Siamese

####  2.4.1 Advanced training objectives

We consider two popular embedding models. They are also encoder-decoders but
they use additional side information. One is trained according to the Siamese
objective [7, 25, 6] the other is a correspondence auto-encoder (CAE)
objective [11]. Both models assume a set of pairs of sequences from the
training corpus. Positive pair are assumed to have the same transcription,
negative pairs, different transcriptions. Let $p_{t}=(s_{t},s_{t^{\prime}},y)$
where $(s_{t},s_{t^{\prime}})$ is a pair of sequences of lengths $T$ and
$T^{\prime}$. A binary value $y$ indicates the positive or negative nature of
the pair. We will see how to find such pairs in the next sub-section.

The CAE objective uses only positive pairs. The auto-encoder is asked to
encode $s_{t}$ into $e_{t}$ and decode it into $\hat{s}_{t}$. The speaker
encoder network is used similarly as for the RNNAE. To satisfy the CAE
objective, $\hat{s}_{t}$ has to minimise the Mean Square Error between
$\hat{s}_{t}$ and $s_{t}^{\prime}$. It forces the auto-encoder to learn a
common representation for $s_{t}$ and $s_{t^{\prime}}$ .

The Siamese objective does not need the decoder network. It encodes both
$s_{t}$ and $s_{t^{\prime}}$ and forces the encoder to learn a similar or
different representation depending on whether the pair is positive or
negative.

$L_{s}(e_{t},e_{t^{\prime}},y)=y\cos(e_{t},e_{t^{\prime}})\,-\,(1-y)\max(0,\cos(e_{t},e_{t^{\prime}})-\gamma)$

where $\cos$ is the cosine similarity and $\gamma$ is a margin. This latter
accounts for negative pairs whose transcriptions have phonemes in common.
These pairs should not have embeddings ’too’ far away from each other. The CAE
and Siamese objective can also be combined into a CAE-Siamese loss by a
weighted sum of their respective loss function [26].

####  2.4.2 Finding and choosing pairs of speech embeddings

Finding positive pairs of speech sequences is an area of research called
unsupervised term discovery (UTD) [15, 16, 17, 27]. Such UTD systems can be
DTW alignment based [15] or involve a k-Nearest-Neighbours search [27]. We
opted for the latter, as it is both scalable and among the state-of-the-art
methods. It encodes exhaustively all possible speech sequences with an
embedding model, and used optimised $k$-NN search [28] to retrieve
acoustically similar pairs of speech sequences (see the details in [27]). In
our experiments, we used the pairs retrieved by $k$-NN on GD-PLP encoded
sequences to train our self-supervised models (CAE,Siamese, CAE-Siamese).

As a supervised alternative, it is possible to sample ‘gold’ pairs, i.e. pairs
of elements that have the exact same transcription. These ‘gold’ pairs are
given to the CAE, Siamese and CAE-Siamese to train supervised models. These
supervised models indicate how good these self-supervised models could be if
we enhanced the UTD system.

##  3 Evaluation metrics and frequency estimation

###  3.1 Intrinsic quality metrics: ABX and MAP

The intrinsic quality of an acoustic speech embedding can be measured using
two types of discrimination tasks: the MAP (also called same-different)[21]
and ABX tasks [22]. Let us consider a set of $n$ acoustic speech embeddings:
$((e_{1},t_{1}),(e_{2},t_{2}),...(e_{n},t_{n}))$ where $e_{i}$ are the
embeddings and $t_{i}$ the transcriptions. The ABX task creates all possible
triplets ($e_{a}$,$e_{b}$,$e_{x}$) such that: $t_{a}=t_{x}$ and $t_{b}\neq
t_{x}$. The model is asked to predict 1 or 0 to indicate if $e_{x}$ is of type
$t_{a}$ or $t_{b}$. Such triplets are instances of the phonetic contrast
between $t_{a}$ and $t_{b}$. Formally for a given a triplet, the task is to
predict:

$y(e_{x},e_{a},e_{b})=\mathds{1}_{d(e_{a},e_{x})\leq d(e_{b},e_{x})}$

The error rate on this classification task is the ABX score. It is first
averaged by type of contrast (all triplets having the same $t_{a}$ and
$t_{b}$) then average over all contrasts.
The MAP task forms a list of all possible pairs of embeddings
($e_{a}$,$e_{x}$). The model is asked to predict 1 or 0 to indicate if $e_{x}$
and $e_{a}$ have the same type, i.e the same transcription, or not. Formally
for a given pair, the model predicts:

$y(e_{a},e_{x},\theta)=\mathds{1}_{d(e_{a},e_{x})\leq\theta}$

The precision and recall on this classification task are computed for various
values of $\theta$. The final score or the MAP task is obtained by integrating
over the precision-recall curve.

###  3.2 Downstream task: unsupervised frequency estimation

####  3.2.1 The R2 metric

Here, we introduce the novel task of frequency estimation as the assignment,
for each speech sequence, of a positive real value that correlate with how
frequent the transcription of this sequence is in a given reference
corpus222This estimation could be up to a scaling coefficient; the task of
finding exact count estimates is a harder task, not tackled in this paper.. To
evaluate the quality of frequency estimates, we use the correlation
determinant $R^{2}$ between estimation and true frequencies. We compute this
number in log space, to take into account the power-law distribution of
frequencies in natural languages [29]. This coefficient is between 0 and 1 and
tells what percentage of variance in the true frequencies can be explained by
the estimated frequencies.

####  3.2.2 $k$-NN and density estimation

We propose to estimate frequencies using density estimation, also called the
Parzen-Rosenblatt window method [30]. Let $N$ be the number of speech sequence
embeddings. First, these $N$ embeddings are indexed into a $k$-NN graph, noted
$G$, where all distances between embeddings are computed. Then, for each
embedding, we search for the $k$ closest embeddings in $G$. Formally, given an
embedding $e_{t}$ from the $k$-NN graph $G$, we compute its $k$ distances to
its $k$ closest neighbours ($d_{n_{1}}$,…$d_{n_{k}}$). The frequency
estimation is a density estimation function $\kappa$ of the $k$-NN graph $G$
that has three parameter: a Gaussian kernel $\beta$, the number of neighbours
$k$ and the embedding $e_{t}$.

$\kappa_{G}(e_{t},\beta,k)=\sum_{i=1}^{k}\exp^{-\beta d_{n_{i}}^{2}}$

This density estimation yields a real number in $[1,k]$, which we take as our
frequency estimation. We set $k$ to $2000$, the maximal frequency that should
be predicted using the transcription of our training corpus (the Buckeye, see
section 4.1). Then, we must tune $\beta$, the dilation of the space of a given
embedding model. For each model, we choose $\beta$ such that it maximises the
variance of the estimated log frequencies, thereby covering the whole spectrum
of possible log frequencies, in our case $[0,\log(k)]$, which is beneficial
for power-law types of distribution. Note that the $\beta$ kernel cannot be
too large (resp. small), as it would predict only high (resp. low) values.

####  3.2.3 Density estimation versus clustering

Models/methods | K-means | HC-K-means |  $k$-NN
---|---|---|---
GD-1hot | 0.67 | 0.73 | 0.74
RNNAE Mel-F | 0.30 | 0.35 | 0.41
CAE Siamese Mel-F | 0.26 | 0.37 | 0.43

Table 1: Frequency estimations using K-means, HC-K-means and $k$-NN density
estimation on a subset of the Buckeye corpus

We compared density estimation with an alternative method: the clustering of
speech embeddings. Jansen et al.[31] did a thorough benchmark of clustering
methods on the task of clustering speech embeddings. Across all their metrics,
the model that performs best is Hierarchical-K-means (HC-K-means), an improved
version of K-means for a higher computational cost. In particular HC-K-means
performs better than GMM HC-K-means is not scalable to our data sets, so we
extracted 1% of the Buckeye corpus in order to compare it with our method. A
similar size of corpus is used by Jansen et al.[31].

We applied $k$-NN, K-means and HC-K-means from the python library scikit-learn
[32] on three of our models on this subset. For K-means and HC-K-means, we
used the hyper-parameters that gave the best scores in [31], namely k-means++
initialisation and average linkage function for HC-K-means. On our subset, the
ground truth number of clusters is $K=33000$. Yet, we did a grid-search on the
value of $k$ that maximises the $R^{2}$ score for frequency estimation. We
found that K-means and HC-K-means perform better for $K=20000$. It shows these
algorithms are not tuned to handle data distributed according to the Zipf’s
law. Indeed K-means is subject the so-called ‘uniform effect’ and tends to
find clusters of uniform sizes [33]. Table 1 shows that even by optimising the
number of clusters $K$, the $k$-NN method outperforms K-means and HC-K-means.

##  4 Experiments

###  4.1 Data sets

Five data sets at our disposal from the ZeroSpeech challenge: Conversational
English (a sample of the Buckeye [34] corpus), English, French, Xitsonga and
Mandarin [2, 3]. These are multi-speaker non-overlapping (i.e one speaker per
file) recordings of speech. All silences were removed using voice activity
detection and corrected manually.

Each corpus was split into all possible segmentations to produce random speech
sequences as described in [27]. Random speech sequences span from $70ms$ to
$1s$. Shorter than $70ms$ sequences may contain less than one phoneme or be
ill-pronounced phonemes. Therefore we removed very short sequences to avoid
issues that are out of scope of this study.

The Buckeye sample corpus contains 12 speakers and 5 hours of speech. The
French and English corpora being much larger, we reduced their number of
speech sequences and speakers to the size of the Buckeye. Mandarin and
Xitsonga are smaller data sets and were left untouched.

Figure 1: Value of metrics and the downstream task across models, corpora. The
average column is the average score over all corpora

###  4.2 Training and hyperparameters

Our encoder-decoder network is a specific use of a three-layers bi-directional
LSTM as described by Holzenberger et al. [5] with hyper-parameters selected to
minimise the ABX error on the Buckeye corpus. The speaker embedding network is
a single fully connected layer with fifteen neurons. Our UTD system [27] uses
the embeddings of the GD-PLP model. A set of speech pairs is returned, sorted
by cosine similarity. We selected the pairs that have a cosine similarity
above $0.85$ as it seemed to be optimal on the Buckeye corpus according to the
ABX metric. In comparison, we trained our supervised models with ‘gold’ pairs,
i.e pairs with the exact same transcription.

Each corpus was randomly split into train (90%), dev(5%) and test (5%). Neural
networks were trained on the train set, early stopping was done using the
development set and metrics computed on the test set. Specifically, we trained
each model on the five training sets using the Buckeye’s hyper-parameters.

MAP was ABX were computed on the test sets. Frequency estimation was computed
by indexing the five training sets and building $k$-NN graphs. For each
element of a given test set, we searched neighbours and estimated frequencies
using the $k$-NN graphs. We used the FAISS [28] library that provides an
optimised $k$-NN implementation.

###  4.3 Results

####  4.3.1 Across models

The results of the two metrics and downstream task are shown in Figure 1 and
the following broad trends can be observed.

  * •

Supervised models yield substantially lower performance than the ground truth
1-hot encodings, on all metric and all languages. These supervised models have
a margin for improvement as they do not learn optimal embeddings despite
having access to ground truth labels.

  * •

Supervised models outperform their corresponding self-supervised model, in
almost all metrics and for all languages. It means that self-supervision has
also a margin for improvement given better pairs from the UTD systems.

  * •

Among self-supervised and supervised models, the CAE-Siamese Mel-F takes the
pole position. This model seems to be able to combine the advantages of both
training objectives. A result already claimed by [26].

  * •

(self) supervised neural neural networks trained on low-level acoustic
features (Mel-F) performs better or equally well as high-level acoustic
features (PLP). This shows that neural networks can learn their own high-level
acoustic features from low-level information.

  * •

Self-supervised models are expected to outperform unsupervised models because
they use side information. Yet many configurations do not show this
consistently. Only the Buckeye data set seems consistent, but this dataset is
the one on which pairs were selected through a grid-search to minimise the ABX
error. This may be due to the variable quality of the pairs found by UTD;
better UTD is therefore needed to help self-supervised models.

  * •

Unsupervised models are supposed to be better than hand-crafted models because
they can adjust by learning from the dataset. Yet, this is not consistently
found. Hand crafted models are worse than unsupervised models for ABX and
frequency estimation but not for MAP.

  * •

In detail, which model is best in a particular language depends on the metric.

####  4.3.2 Across metrics and frequency estimation

In Table 2, we quantified the possibility to observe the discrepancies that we
have just discussed . We computed the correlation $R^{2}$ across the three
‘average’ columns. Cross correlation scores range from $R^{2}=0.33$ to $0.53$;
the top-line model is not included when computing these scores.

R2 | Frequency est. | MAP | ABX
---|---|---|---
Frequency est. | 1.0 | 0.34 | 0.53
MAP | 0.34 | 1.0 | 0.45
ABX | 0.53 | 0.45 | 1.0

Table 2: Correlation $R^{2}$ across the ’average’ column of MAP, ABX and
frequency estimation

These correlations are low enough to permit sizeable discrepancies across
metrics and the downstream task. One of our model, the RNNAE Mel-F, epitomises
the problem. This model is comparatively bad according to the MAP but good
according to ABX and the frequency estimation. It means that MAP and ABX
reveal different aspect of the reliability of embedding models. Therefore,
only large progress according to one metric assures a progress according to an
other metric. It shows the limit of intrinsic evaluation of speech embeddings.
Moderate variations on a intrinsic metric cannot guarantee a progress on a
given downstream task.

ABX and MAP scores are averages over multiple phonetic contrasts. These
contrasts could be clustered based on their phonetic frequencies, average
lengths or number of phonemes in common. Such fined-grained analyses can
sometimes give understanding divergences across metrics. However, we have been
unable to find a categorisation of results that make sense of Figure 1 as a
whole. There are currently no fully reliable metrics to assess the intrinsic
quality of speech embeddings.

##  5 Conclusion

We quantified the correlation across two intrinsic metrics (MAP and ABX) and a
novel downstream task: frequency estimation. Although MAP and ABX agree on
general categories (like supervised versus unsupervised embeddings), we also
found large discrepancies when it comes to select a particular model
highlighting the limits of these intrinsic quality metrics. However convenient
intrinsic metrics may be, they only show partial views of the overall
reliability of a model. We showed using frequency estimation that variations
on intrinsic quality metrics should not be accounted for certain progress on
downstream tasks. More attention should be brought on downstream tasks that
have the credit to answer practical problems.

##  6 Acknowledgements

We thank Matthijs Douze for useful comments on density estimation. We also
thank IDRIS from CNRS for offering GPU resources on the supercomputer Jean
Zay. This work was funded in part by the Agence Nationale pour la Recherche
(ANR-17-EURE-0017 Frontcog, ANR-10-IDEX-0001-02 PSL*, ANR-19-P3IA-0001 PRAIRIE
3IA Institute), CIFAR, and a research gift by Facebook.

## References

  * [1] A. Jansen, E. Dupoux, S. Goldwater, M. Johnson, S. Khudanpur, K. Church, N. Feldman, H. Hermansky, F. Metze, R. Rose, M. Seltzer, P. Clark, I. McGraw, B. Varadarajan, E. Bennett, B. Borschinger, J. Chiu, E. Dunbar, A. Fourtassi, D. Harwath, C. Lee, K. Levin, A. Norouzian, V. Peddinti, R. Richardson, T. Schatz, and S. Thomas, “A summary of the 2012 jhu clsp workshop on zero resource speech technologies and models of early language acquisition,” in _2013 IEEE International Conference on Acoustics, Speech and Signal Processing_ , 2013, pp. 8111–8115.
  * [2] M. Versteegh, R. Thiollière, T. Schatz, X.-N. Cao Kam, X. Anguera, A. Jansen, and E. Dupoux, “The zero resource speech challenge 2015,” 09 2015.
  * [3] E. Dunbar, X. Cao, J. Benjumea, J. Karadayi, M. Bernard, L. Besacier, X. Anguera, and E. Dupoux, “The zero resource speech challenge 2017,” _CoRR_ , vol. abs/1712.04313, 2017. [Online]. Available: http://arxiv.org/abs/1712.04313
  * [4] H. Kamper, “Truly unsupervised acoustic word embeddings using weak top-down constraints in encoder-decoder models,” _CoRR_ , vol. abs/1811.00403, 2018\. [Online]. Available: http://arxiv.org/abs/1811.00403
  * [5] N. Holzenberger, M. Du, J. Karadayi, R. Riad, and E. Dupoux, “Learning word embeddings: Unsupervised methods for fixed-size representations of variable-length speech segments,” 09 2018, pp. 2683–2687.
  * [6] S. Settle and K. Livescu, “Discriminative acoustic word embeddings: Recurrent neural network-based approaches,” _CoRR_ , vol. abs/1611.02550, 2016. [Online]. Available: http://arxiv.org/abs/1611.02550
  * [7] R. Riad, C. Dancette, J. Karadayi, N. Zeghidour, T. Schatz, and E. Dupoux, “Sampling strategies in siamese networks for unsupervised speech representation learning,” _CoRR_ , vol. abs/1804.11297, 2018. [Online]. Available: http://arxiv.org/abs/1804.11297
  * [8] K. Levin, K. Henry, A. Jansen, and K. Livescu, “Fixed-dimensional acoustic embeddings of variable-length segments in low-resource settings,” 12 2013, pp. 410–415.
  * [9] Y.-A. Chung, C.-C. Wu, C.-H. Shen, H.-y. Lee, and L.-S. Lee, “Audio word2vec: Unsupervised learning of audio segment representations using sequence-to-sequence autoencoder,” 03 2016.
  * [10] Y. Chung, W. Weng, S. Tong, and J. R. Glass, “Unsupervised cross-modal alignment of speech and text embedding spaces,” _CoRR_ , vol. abs/1805.07467, 2018. [Online]. Available: http://arxiv.org/abs/1805.07467
  * [11] D. Renshaw, H. Kamper, A. Jansen, and S. Goldwater, “A comparison of neural network methods for unsupervised representation learning on the zero resource speech challenge,” _Annual Conference of the International Speech Communication Association_ , 2015.
  * [12] T. J. Hazen, W. Shen, and C. White, “Query-by-example spoken term detection using phonetic posteriorgram templates,” in _2009 IEEE Workshop on Automatic Speech Recognition Understanding_ , 2009, pp. 421–426.
  * [13] K. Levin, A. Jansen, and B. Van Durme, “Segmental acoustic indexing for zero resource keyword search,” in _2015 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 2015, pp. 5828–5832.
  * [14] Y. Wang, H. Lee, and L. Lee, “Segmental audio word2vec: Representing utterances as sequences of vectors with applications in spoken term detection,” _CoRR_ , vol. abs/1808.02228, 2018. [Online]. Available: http://arxiv.org/abs/1808.02228
  * [15] A. S. Park and J. R. Glass, “Unsupervised pattern discovery in speech,” _IEEE Transactions on Audio, Speech, and Language Processing_ , vol. 16, no. 1, pp. 186–197, 2008.
  * [16] C.-y. Lee, T. J. O’Donnell, and J. Glass, “Unsupervised lexicon discovery from acoustic input,” _Transactions of the Association for Computational Linguistics_ , vol. 3, pp. 389–403, 2015. [Online]. Available: https://www.aclweb.org/anthology/Q15-1028
  * [17] O. Räsänen, G. Doyle, and M. Frank, “Unsupervised word discovery from speech using automatic segmentation into syllable-like units,” 09 2015.
  * [18] S. Goldwater, T. Griffiths, and M. Johnson, “A bayesian framework for word segmentation: Exploring the effects of context,” _Cognition_ , vol. 112, pp. 21–54, 04 2009.
  * [19] H. Kamper, A. Jansen, and S. Goldwater, “A segmental framework for fully-unsupervised large-vocabulary speech recognition,” _Computer Speech and Language_ , 2017.
  * [20] K. Kawakami, C. Dyer, and P. Blunsom, “Unsupervised word discovery with segmental neural language models,” _CoRR_ , vol. abs/1811.09353, 2018. [Online]. Available: http://arxiv.org/abs/1811.09353
  * [21] M. A. Carlin, S. Thomas, A. Jansen, and H. Hermansky, “Rapid evaluation of speech representations for spoken term discovery,” in _INTERSPEECH_ , 2011\.
  * [22] T. Schatz, V. Peddinti, F. Bach, A. Jansen, H. Hermansky, and E. Dupoux, “Evaluating speech features with the minimal-pair abx task: Analysis of the classical mfc/plp pipeline,” _INTERSPEECH 2013: 14th Annual Conference of the International Speech Communication Association_ , pp. 1–5, 01 2013.
  * [23] S. B. Davis and P. Mermelstein, “Comparison of parametric representations for monosyllabic word recognition in continuously spoken sentences,” _ACOUSTICS, SPEECH AND SIGNAL PROCESSING, IEEE TRANSACTIONS ON_ , pp. 357–366, 1980.
  * [24] H. Hermansky, “Perceptual linear predictive (plp) analysis of speech.” _The Journal of the Acoustical Society of America_ , vol. 87 4, pp. 1738–52, 1990.
  * [25] R. Thiollière, E. Dunbar, G. Synnaeve, M. Versteegh, and E. Dupoux, “A hybrid dynamic time warping-deep neural network architecture for unsupervised acoustic modeling,” in _INTERSPEECH_ , 2015.
  * [26] P. Last, H. A. Engelbrecht, and H. Kamper, “Unsupervised feature learning for speech using correspondence and siamese networks,” _IEEE Signal Processing Letters_ , vol. 27, pp. 421–425, 2020.
  * [27] A. Thual, C. Dancette, J. Karadayi, J. Benjumea, and E. Dupoux, “A k-nearest neighbours approach to unsupervised spoken term discovery,” in _2018 IEEE Spoken Language Technology Workshop (SLT)_ , 2018, pp. 491–497.
  * [28] J. Johnson, M. Douze, and H. Jégou, “Billion-scale similarity search with gpus,” _arXiv preprint arXiv:1702.08734_ , 2017.
  * [29] G. K. Zipf, _The psycho-biology of language_ , 1935.
  * [30] E. Parzen, “On estimation of a probability density function and mode,” _Ann. Math. Statist._ , vol. 33, no. 3, pp. 1065–1076, 09 1962. [Online]. Available: https://doi.org/10.1214/aoms/1177704472
  * [31] H. Kamper, A. Jansen, S. King, and S. Goldwater, “Unsupervised lexical clustering of speech segments using fixed-dimensional acoustic embeddings,” in _2014 IEEE Spoken Language Technology Workshop (SLT)_ , 2014, pp. 100–105.
  * [32] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay, “Scikit-learn: Machine learning in Python,” _Journal of Machine Learning Research_ , vol. 12, pp. 2825–2830, 2011.
  * [33] J. Wu, _The Uniform Effect of K-means Clustering_ , 07 2012, pp. 17–35.
  * [34] M. A. Pitt, K. Johnson, E. Hume, S. Kiesling, and W. Raymond, “The buckeye corpus of conversational speech: Labeling conventions and a test of transcriber reliability,” _Speech Communication_ , vol. 45, no. 1, pp. 89–95, 2005.


* [39] Mohammad Atif, Praveen Kumar Kolluru, Chakradhar Thantanapally, and Santosh Ansumali.  Essentially entropic lattice Boltzmann model.  Physical Review Lett., 119:240602, Dec 2017.
  * [40] Jonas Latt, Christophe Coreixas, Joël Beny, and Andrea Parmigiani.  Efficient supersonic flow simulations using lattice Boltzmann methods based on numerical equilibria.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 378(2175):20190559, jul 2020.
  * [41] Dominique D’Humières.  Generalized Lattice-Boltzmann Equations.  Rarefied Gas Dynamics: Theory and Simulations, 159:450–458, jan 1994.
  * [42] Pierre Lallemand and Li-Shi Luo.  Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, Galilean invariance, and stability.  Physical Review E, 61(6):6546–6562, 2000.
  * [43] Dominique D’Humières, Irina Ginzburg, Manfred Krafczyk, Pierre Lallemand, and Li-Shi Luo.  Multiple-relaxation-time lattice Boltzmann models in three dimensions.  Phil. Trans. R. Soc. A, 360(1792):437–451, 2002.
  * [44] R. Courant, K Friedrichs, and H. Lewy.  On the Partial Difference Equations of Mathematical Physics.  IBM Journal of Research and Development, 11(2):215–234, mar 1967\.
  * [45] Amiram Harten.  On the symmetric form of systems of conservation laws with entropy.  Journal of Computational Physics, 49(1):151–164, 1983.
  * [46] T.J.R Hughes, L.P. Franca, and M. Mallet.  A new finite element formulation for computational fluid dynamics: I symmetric forms of the compressible Euler and Navier-Stokes equations and the second law of thermodynamics.  Computer Methods in Applied Mechanics and Engeneering, 54:223–234, 1986.
  * [47] D. Aregba-Driollet and R. Natalini.  Discrete kinetic schemes for multi-dimensional systems of conservation laws.  SIAM J. Numer. Anal., 37(6):1971–2004, 2000.
  * [48] R. Abgrall and D. Torlo.  Asymptotic preserving deferred correction residual distribution schemes.  https://arxiv.org/abs/1811.09284, November 2018.
  * [49] A. V. Bobylev.  Instabilities in the Chapman-Enskog expansion and hyperbolic burnett equations.  Journal of Statistical Physics, 124(2-4):371–399, 2006.
  * [50] R. Natalini.  A discrete kinetic approximation of entropy solution to multi-dimensional scalar conservation laws.  Journal of differential equations, 148:292–317, 1998.
  * [51] A.Iserles.  Order stars and saturation theorem for first-order hyperbolics.  IMA J. Numer. Anal., 2:49–61, 1982.
  * [52] E. Hairer and G. Wanner.  Solving ordinary differential equations. II, volume 14 of Springer Series in Computational Mathematics.  Springer-Verlag, Berlin, 2010.  Stiff and differential-algebraic problems, Second revised edition, paperback.
  * [53] Maria Han Veiga, Lorenzo Micalizzi, and Davide Torlo.  On improving the efficiency of ader methods, 2023.
  * [54] Remi Abgrall and Davide Torlo.  Some preliminary results on a high order asymptotic preserving computationally explicit kinetic scheme, 2021.
  * [55] Lorenzo Micalizzi and Davide Torlo.  A new efficient explicit deferred correction framework: analysis and applications to hyperbolic pdes and adaptivity, 2023.
  * [56] Edward R. Benton and George W. Platzman.  A table of solutions of the one-dimensional Burgers equation.  Quarterly of Applied Mathematics, 30(2):195–212, 1972.
  * [57] Ya. B. Zeldovich.  Physics of Shock Waves and High-Temperature Hydrodynamic Phenomena.  Elsevier, 1967.