# SpeechBlender: Speech Augmentation Framework
for Mispronunciation Data Generation

###### Abstract

One of the biggest challenges in designing mispronunciation detection models
is the unavailability of labeled L2 speech data. To overcome such data
scarcity, we introduce SpeechBlender – a fine-grained data augmentation
pipeline for generating mispronunciation errors. The SpeechBlender utilizes
varieties of masks to target different regions of a phonetic unit, and use the
mixing factors to linearly interpolate raw speech signals while generating
erroneous pronunciation instances. The masks facilitate smooth blending of the
signals, thus generating more effective samples than the ‘Cut/Paste’ method.
We show the effectiveness of our augmentation technique in a phoneme-level
pronunciation quality assessment task, leveraging only a good pronunciation
dataset. With SpeechBlender augmentation, we observed a 3% and 2% increase in
Pearson correlation coefficient (PCC) compared to no-augmentation and goodness
of pronunciation augmentation scenarios respectively for Speechocean762
testset. Moreover, a 2% rise in PCC is observed when comparing our single-task
phoneme-level mispronunciation detection model with a multi-task learning
model using multiple-granularity information.

Index Terms—  Mispronunciation detection, Data augmentation

##  1 Introduction

With the advent of globalization, more people are now inclined towards
multilingualism. Such a shift in interest promotes the demand of the computer-
aided pronunciation training (CAPT) system to language learning application
[1, 2]. The CAPT system potentially exploits the advancement of speech
technology along with curriculum management and learner assessment among
others. Mispronunciation detection (MD) and diagnosis system – a crucial
component of the CAPT – detects the mispronunciation in a L2 learner‘s speech
and highlights the error to give effective feedback.

Over years, researchers have exploited many methods to assess pronunciation
quality. A majority of these methods use pre-trained automatic speech
recognition (ASR) systems to either (i) pinpoint the mismatch between the
aligned ASR output sequence with the reference; or (ii) utilize the log-
posterior probability from the ASR to compute variants of goodness of
pronunciation goodness of pronunciation (GOP) scores [3, 4, 5] to asses
pronunciation. One successful line of research treated the mispronunciation
detection as a regression or a classification task [6]. Using deep learning,
the network is either trained in an end-to-end manner [7, 8, 9, 10] or built
upon the computed GOP features from the pre-trained ASR [11, 12, 13, 14] to
detect errors. These approaches are heavily dependent on the availability of a
large amount of manually annotated datasets for the target demographic L2
speakers. However, such annotated datasets are quite scarce and often contain
uneven positive class labels.111In the case of MD, the positive class
represents – missed, mispronounced, or accented phonemes, and the negative
classes are good pronounced phonemes.

Data augmentation techniques are recently proved to be quite effective to
overcome the data scarcity for MD. Techniques involved (i) modifying the
canonical text with mismatched character pair keeping its original word-level
speech intact [9]; (ii) using Text-To-Speech to synthesize ‘incorrect stress’
samples along with modified lexical stress [15]; (iii) using a mixup
technique, in feature space – based on phone-level GOP pool to design word
level training data [16]; and (iv) using error distance of the clustered self
supervised learning model embedding to replace the phoneme sound with another
close sound unit [17, 18].

The aforementioned techniques simulate mispronunciation errors by finding a
target phoneme unit, in a different – text, feature or audio space, and use it
to replace the current unit. However, the L2 learners’ pronunciation is highly
affected by their native phoneme set [19] and ‘Cut/Paste’ approach reuses the
available information without generating any new feature variation in phoneme
level. Thus making such mixup – Cut/Paste augmentation methods less efficient.

In this study, we introduce SpeechBlender – a simple fine-grained data
augmentation pipeline for generating mispronunciation errors in low-resource
settings. The SpeechBlender framework linearly interpolate raw input signals
utilizing different masks and mixing factors, allowing to modify different
regions of a target phoneme. Moreover, the proposed framework generates new
training samples using only good pronunciation units and thus capable of
generating MD training data from any native/non-native speech datasets.

We show the effectiveness of the proposed augmentation framework over text-
and GOP-based augmentation techniques using a significantly faster non-
autoregressive LSTM [20] model for phoneme-level MD task. We compare the
performance of different masks – such as Smooth Overlay, Smooth Concatenation,
CutMix, and Smooth Gaussian Overlay – in the SpeechBlender and show its
efficacy over the well-known Cut/Paste technique. We present our results using
a well-establish publicly available MD dataset, Speechocean762 [21]. While the
proposed technique is studied for a specific task (MD) in this paper, the
SpeechBlender can be broadly applicable as an augmentation technique for most
sequence labeling such as accented ASR or classification task with
audio/speech as input. Our contributions are:

  1. (a)

Introduce the SpeechBlender, a fine-grained data augmentation pipeline capable
of smoothly modifying regions of phonetic unit to generate erroneous
pronunciation;

  2. (b)

Show the effectiveness of the proposed method over text-based and GOP-based
augmentation;

  3. (c)

Analyse the performance of different masks in the SpeechBlender and compare it
with cut/paste signal approach;

  4. (d)

We make the SpeechBlender augmentation module publicly available.222The code
will be shared upon paper acceptance.

Fig. 1: Overview of the SpeechBlender, a fine-grained data augmentation
pipeline, for zero-shot mispronunciation class generation. $U_{C}$: input
speech utterance for candidate phoneme; $U_{D}$: selected speech utterance for
donor phoneme; $x_{C}(T)$, $x_{D}(T)$ candidate and donor phoneme frames of
length $T$ and $L$ respectively; $\hat{x}$ is the resultant augmented phoneme.

##  2 Methodology

Figure 1 shows an overview of the SpeechBlender – a fine-grained data
augmentation pipeline, for generating mispronunciation classes in a low-
resource settings.333Presence of a few or zero erroneous pronunciation
examples in the data. From the canonical phoneme sequence of an input
utterance $U_{C}$, we select a few good phonetic units as augmentation
candidates. For each candidate unit, we pick a donor phoneme using a close
phoneme-pair dictionary and choose another utterance, $U_{D}$ as a source
signal of the donor phoneme. We obtain the alignments of the utterances
$U_{C}$ and $U_{D}$ using a pre-trained ASR, and get the corresponding phoneme
time boundaries. Using the start and end boundary of the candidate, $x_{C}$,
and donor, $x_{D}$, phoneme, we select the speech frames and apply
augmentation with the SpeechBlender. We then replace the candidate phoneme
$x_{C}$ with the augmented phoneme frames $\hat{x}$ in the candidate utterance
$U_{C}$.

###  2.1 Close Phoneme Pair Dictionary

The ‘close phoneme pair dictionary’ is designed to randomly select a donor
phoneme from a list of confused phonemes associated with the candidate
phoneme. The dictionary is constructed using both the Speechocean762 [21] and
L2-ARCTIC [19] training set. The entries are selected by analysing the
confused phoneme pairs of speaker’s pronunciation using ASR errors like
substitution along with human annotation labels. Some examples of close
phoneme pairs are presented in Table 1.

Close Phoneme Pairs (Confused %)
---
$SH\Leftrightarrow S$ (76%), $V\Leftrightarrow F$ (44%) , $NG\Leftrightarrow
N$ (43%),
$IY\Leftrightarrow IH$ (33%), $Z\Leftrightarrow S$ (77%)

Table 1: Examples of entries in Close Phoneme Pair Dictionary, with their % of
confusion.

###  2.2 Pre-trained ASR and Alignment Module

We align the acoustic frames with the corresponding reference phoneme to
obtain the start and end timestamps. For the force alignment, we employ the
Hidden Markov Model-Time Delay Neural Network (HMM- TDNN)444https://kaldi-
asr.org/models/m13 as the acoustic module. It is trained on Librispeech ASR
task with Kaldi [22].

Fig. 2: Comparison of the different masks used in the SpeechBlender
augmentation framework. $x_{C}$ represent candidate phoneme (with length $T$)
in the good data, augmented with the donor phoneme $x_{D}$ (length =$L$).
$\mu$: the center-point; $R$: the number of regions the mask is divided into;
$W$: list of widths of the regions; $\lambda$: the mixing factor in each
region assigned in different augmentation scenarios given in (a) - (e).

Input: Input phoneme frame sequences $x_{C}$ and $x_{D}$ where
$x_{C}=\\{f^{1}_{C}\cdots f^{T}_{C}\\}$ and $x_{D}=\\{f^{1}_{D}\cdots
f^{L}_{D}\\}$, where ($T$, $L$) are total frames in ($x_{C}$, $x_{D}$)
respectively. $y_{C}$ represent the original output score assigned to the good
phoneme.

Output: Augmented frame output sequences $\hat{x}$ and the corresponding label
$\hat{y}$ of the output.

Function  _Blender(_$x_{C},x_{D},y_{C}$_)_:

$\theta_{\hat{x}}\longleftarrow[]$ $\triangleright$ Store augmented output per
region

$\theta_{\hat{y}}\longleftarrow[]$ $\triangleright$ Store corresponding
partial output score

$mid\longleftarrow rand([1,4])$ $\triangleright$ Mask template id

$\rho\longleftarrow get\\_property(mid,len(x_{C}),len(x_{D}))$
$\triangleright$ Mask property

$W,\Lambda\longleftarrow generate\\_mask(\rho)$

$\triangleright$ $W$: list of regional mask width, $\Lambda$: list of regional
mixing factor values.

foreach  _$w\in W$ and $\lambda\in\Lambda$_ do

$\theta_{\hat{x}}\longleftarrow$ $\lambda\cdot x{{}_{C}^{w}}$ \+
$(1-\lambda)\cdot x{{}_{D}^{w}}$

$\theta_{\hat{y}}\longleftarrow\lambda\cdot y{{}_{C}}$

end foreach

$\hat{x}\longleftarrow Concat(\theta_{\hat{x}})$

$\hat{y}\longleftarrow\left\lfloor\sum\theta_{\hat{y}}/len(\Lambda)\right\rfloor$

return $\hat{x},\hat{y}$

Algorithm 1 SpeechBlender function

###  2.3 SpeechBlender

Given the input frames of good-phoneme candidate $x_{C}$ and the donor-phoneme
$x_{D}$, of length $T$ and $L$ respectively, our approach outputs the
augmented phoneme $\hat{x}$ along with its corresponding output label
$\hat{y}$. We first normalize the energy of $x_{D}$ to match $x_{C}$.
Subsequently, we apply the SpeechBlender function, shown in Algorithm 1, with
a selected mask templates.

#### Mask Templates:

We first randomly select a mask template id, $mid\in[1-4]$. We then generate
the property $\rho$ of the mask based on the $mid$, length of the phonemes $T$
and $L$. The property $\rho\in\\{\mu,R,W,\Lambda\\}$ is parameterized by the
features of the mask including: (a) the center-point $\mu=\frac{min(T,L)}{2}$;
(b) the number of regions $R\in\\{1,3\\}$ the mask is divided into; (c) along
with a list containing all widths of the regions, $W$; and (d) a list
$\Lambda$, with its corresponding mixing factor with each value
$\lambda\in[0-1]$. The width of each region is determined by the minimum frame
length of inputs when both the input pair is considered with $\lambda>0$. A
$\lambda=0$ for a signal represent a mixing weight of zero, meaning the full
signal is blocked out. Details of different mask templates along with its
property, $\rho$, are presented in Figure 2.

#### Blending the Speech signals:

The Given property $\rho$ of the mask, we generated the mask with its $W$ list
and its corresponding $\Lambda$ list containing region-based $\lambda$ values.
Then for each regions in the mask, we multiply the $x_{C}$ and $x_{D}$ with
$\lambda$ and $1-\lambda$ respectively (shown in Algorithm 1) to create the
regional augmented data $\theta_{\hat{x}}(.)$. We create the new augmented
phoneme, $\hat{x}$ by concatenating all the $\theta_{\hat{x}}$ regional data.
The label $\hat{y}\in\\{0,1\\}$ of the augmented phoneme $\hat{x}$ is
generated by averaging the regional label ($\theta_{\hat{y}}=\lambda\cdot y$).
The $y$ is the value assigned to a good pronunciation of a phoneme.555These
values are determined based on the original data annotation scheme. For this
study, $y=2$. The augmented output $\hat{y}=0$ represent the mispronounced or
missing phoneme; whereas $\hat{y}=1$ represent the accented or aberrant
phoneme unit. In our experiment, the regional augmented output is assigned
$\theta_{\hat{y}}=0$ if the mixing factor $\lambda<0.25$, otherwise it is
assigned $\theta_{\hat{y}}=1$.

###  2.4 Baseline Augmentation Methods

We compare the performance of the SpeechBlender with text-level augmentation
and GOP feature-level augmentation.

#### Text-based Augmentation

For text-based phoneme data augmentation, we increase the number of error
phoneme by altering the phoneme embedding from the canonical phoneme sequence.
The candidate phoneme was replaced with a selected donor phoneme, using the
same ‘close phoneme pair dictionary’, when creating the canonical phoneme
embedding space. An output label $\hat{y}=0$ is assigned when the candidate
phoneme is swapped with a distant donor phoneme, for e.g., $SH\Leftrightarrow
V$; whereas $\hat{y}=1$ is assigned for a close phoneme exchange such as
$SH\Leftrightarrow S$.

#### GOP-based Augmentation

For GOP-based data augmentation, we first create a bag of GOP feature
embedding for all good pronounced phonemes and grouped them by phoneme pairs
using the same ‘close phoneme pair dictionary’. We then replaced the candidate
phoneme, in the utterance with a donor phoneme GOP embedding randomly selected
from the bag. We follow the same scoring scheme as text-augmentation for
$\hat{y}$ labels.

##  3 Experiments

In this Section, we describe the dataset and the system architectures used to
show the efficacy of the proposed SpeechBlender augmentation.

###  3.1 Speech Corpus

For the study, we opt for the Speechocean762 [21], a free open-source dataset
designed to assess the pronunciation quality of L2 learners. The dataset
consists of 5000 English utterances from 250 non-native speakers with rich
annotation labels. Along with different utterances and the word level quality
aspects, the dataset include accuracy scores at phoneme-level that are divided
based on three class labels, as shown in Table 2. Each phoneme is annotated by
five experts and the final phoneme accuracy score is calculated as their
average.

Label | Label Description | Train Set (%) | Test Set (%)
---|---|---|---
0 |  | mispronounced or
---
missing phoneme
1991 (4.2%) | 1339 (2.8%)
1 | accented pronunciation | 1967 (4.1%) | 1828 (3.8%)
2 | good pronunciation | 43118 (91.7%) | 44079 (93.4%)

Table 2: Details phoneme-level statistics and score definition used in
Speechocean762 Dataset

###  3.2 Mispronunciation Detection Model

An overview of the mispronunciation detection model architecture is presented
in Figure 3. Given the input audio and the canonical transcription, we first
extract the frame-level phonetic features using pre-trained ASR, described in
Section 2.2. Using the acoustic model, we calculated the 84-dimensional GOP
features for 42 pure phonemes. Each phone GOP representation contain two
different features: (i) log phone posterior (LPP) [6] of a phone $p$ of length
$T$ frames with per frame observation $o_{t}$ see Equation 1; and (ii) the log
posterior ratio (LPR) [6] of phone $p_{j}$ with respect to phone $p_{i}$, as
given in Equation 2

$\small LPP(p)=\frac{1}{T}\sum log\ p(p|o_{t})$ (1) $LPR(p_{j}|p_{i})=log\
p(p_{j}|o_{t})-log\ p(p_{i}|o_{t})$ (2)

These 84-dimensional features are then passed through a feed-forward layer to
project the GOP embeddings in 24 dimensional feature space. Simultaneously, we
inject the textual-content information using the one-hot vector representation
of the canonical phoneme transcription. These vector representations are then
projected to 24-dimensional phoneme embedding and are then added to the GOP
representations.

Fig. 3: Illustration of the proposed LSTM architecture with a sample utterance
“My Name”. The number of encoder layers $N=3$

Consequently, we pass the final representation to the encoder network. As our
encoder, we opt for a long short-term memory (LSTM)666Additional models such
as SVR[21], MLP are tested with SpeechBlender data augmentation. As the
results follow the same trend , for brevity we are reporting the results using
LSTM only. network due to the model’s strength in capturing preceding
contextual information, and its reduced computational cost without losing much
on the model accuracy. For further comparison, we implement the architecture
with similar configuration as [11] with 3 stacked LSTM layers and an output
embedding dimension of 24. We then pass the embedding per phoneme to the
output layer for the regression task.

Model Parameters: We train the mispronunciation detection model using an Adam
optimizer for $20$ epochs with a starting learning rate of $1e-2$ and a batch
size of $32$. The learning rate is then cut in to its half in every 3 epochs,
after 10 initial epochs. As for the loss function, we use the mean squared
error (MSE). We report the results based on the last epoch.

###  3.3 Evaluation Measures

We measure the proposed model performance using standard statistical metric:
MSE and Pearson correlation coefficient (PCC). We base our discussion using
PCC as it is more representative of the performance in such an imbalanced
class scenario.

##  4 Results and Discussion

In this Section, we report the results in Table 3 with different augmentation
settings and compare it with no-augmentation and multi-task learning setup
proposed in [11]. We validate our proposed augmentation pipeline using the
Speechocean762 test data. Furthermore, we show the strength of each mask
proposed with the SpeechBlender and compare it with the Cut/Paste method in
Table 4.

#### Comparison with Baselines:

The reported PCC results in Table 3 indicates the importance of augmentation
in general. Our results indicate both the GOP-augmentation and SpeechBlender
augmentation outperforms the model with only Speechocean762 trainset by $1\%$
and $3\%$ respectively. However, the same is not seen for text-augmentation
baseline. The GOP-augmentation, a feature space replacing technique, performs
comparable with the multi-task learning setup with a PCC of $0.59$. Our
proposed augmentation technique out performs both the GOP-augmentation and the
multi-task setup by 2%. using only phoneme-level information for training.

#### Effectiveness of the Masks:

We assess the effectiveness of each standalone mask incorporated in the
SpeechBlender for creating different mispronunciation instances. From our
results, reported in Table 4, we observe that the smooth blending masks
(Smooth Overlay and Smooth-Gaussian Overlay) create more practical variation
in mispronunciation and accented phoneme than their counter-part. When
compared directly with the ‘Cut/Paste’ augmentation mask, our bending mask
also shows its efficacy.

Training Dataset | PCC | MSE
---|---|---
Baselines
Original | 0.58 | 0.090
Original Multi-task Learning [11] | 0.59 | 0.089
Augmentation Baselines
Original + Text-augmentation | 0.53 | 0.114
Original + GOP-augmentation | 0.59 | 0.090
Our Proposed Augmentation Method
Original + SpeechBlender-Aug | 0.61 | 0.085
Table 3: Reported phoneme-level PCC and MSE on Speechocean762 test data

Training Dataset | PCC | MSE
---|---|---
Mask creating accented phoneme class, $\hat{y}=1$
Original + Smooth Overlay ($\lambda=0.5$) | 0.61 | 0.086
Original + CutMix | 0.59 | 0.088
Original + Smooth Concatenation | 0.59 | 0.089
Original + Smooth-Gausian Overlay | 0.60 | 0.088
Original + Smooth Overlay ($\lambda=0.6$) | 0.60 | 0.087
Mask creating mispronounced phoneme class, $\hat{y}=0$
Original + Cut/Paste | 0.55 | 0.104
Original + Smooth Overlay ($\lambda=0.1$) | 0.56 | 0.096
Original + Smooth Overlay ($\lambda=0.2$) | 0.58 | 0.092

Table 4: Reported phoneme-level PCC and MSE on Speechocean762 test for
different augmentation mask.

##  5 Conclusion

In this paper, we introduce a fine-grained phoneme-level augmentation
technique, SpeechBlender, and show its efficacy in generating new error
classes for a phoneme-level pronunciation quality detection task. The
SpeechBlender linearly interpolate raw input signals using different blending
masks like ‘SmoothOverlay’ and mixing factors $\lambda$. Varieties of masks
included in the proposed method allows smooth modification in different
regions of the phoneme. With SpeechBlender, we observe improved performance
compared to (i) no-augmentation, (ii) multi-task learning scenario, (iii)
text-augmentation, and (iv) GOP-augmentation setting. We noticed the
‘Cut/Paste’ method be it in text-, feature- or raw-audio domain, is not as
efficient as the proposed interpolation technique. Moreover, our results
suggest that with blending the signals, we can generate new positive error
classes instead of reusing samples from the data itself in zero-shot setting.
The different choices of masks allow variations in transition and style of
errors. In future, we will further explore the capability of SpeechBlender for
generating accented/erroneous phoneme candidate for other speech tasks like
ASR, while exploring multiple granularities.

## References

  * [1] Pamela M Rogerson-Revell,  “Computer-assisted pronunciation training (CAPT): Current issues and future directions,”  RELC Journal, vol. 52, no. 1, pp. 189–205, 2021.
  * [2] Ambra Neri, Ornella Mich, Matteo Gerosa, and Diego Giuliani,  “The effectiveness of computer assisted pronunciation training for foreign language learning by children,”  Computer Assisted Language Learning, vol. 21, no. 5, pp. 393–408, 2008.
  * [3] Silke M Witt and Steve J Young,  “Phone-level pronunciation scoring and assessment for interactive language learning,”  Speech communication, vol. 30, no. 2-3, pp. 95–108, 2000.
  * [4] Jiatong Shi, Nan Huo, and Qin Jin,  “Context-aware goodness of pronunciation for computer-assisted pronunciation training,”  arXiv preprint arXiv:2008.08647, 2020.
  * [5] Sweekar Sudhakara, Manoj Kumar Ramanathi, Chiranjeevi Yarra, and Prasanta Kumar Ghosh,  “An improved goodness of pronunciation (GOP) measure for pronunciation evaluation with DNN-HMM system considering HMM transition probabilities.,”  in INTERSPEECH, 2019, pp. 954–958.
  * [6] Binghuai Lin and Liyuan Wang,  “Deep feature transfer learning for automatic pronunciation assessment.,”  in Interspeech, 2021, pp. 4438–4442.
  * [7] Wai-Kim Leung, Xunying Liu, and Helen Meng,  “CNN-RNN-CTC based end-to-end mispronunciation detection and diagnosis,”  in ICASSP 2019-2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2019, pp. 8132–8136.
  * [8] Yiqing Feng, Guanyu Fu, Qingcai Chen, and Kai Chen,  “SED-MDD: Towards sentence dependent end-to-end mispronunciation detection and diagnosis,”  in ICASSP 2020-2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2020, pp. 3492–3496.
  * [9] Kaiqi Fu, Jones Lin, Dengfeng Ke, Yanlu Xie, Jinsong Zhang, and Binghuai Lin,  “A full text-dependent end to end mispronunciation detection and diagnosis with easy data augmentation techniques,”  arXiv preprint arXiv:2104.08428, 2021.
  * [10] Wenping Hu, Yao Qian, and Frank K Soong,  “A new DNN-based high quality pronunciation evaluation for computer-aided language learning (call).,”  in Interspeech, 2013, pp. 1886–1890.
  * [11] Binghuai Lin, Liyuan Wang, Xiaoli Feng, and Jinsong Zhang,  “Automatic scoring at multi-granularity for L2 pronunciation.,”  in Interspeech, 2020, pp. 3022–3026.
  * [12] Yuan Gong, Ziyi Chen, Iek-Heng Chu, Peng Chang, and James Glass,  “Transformer-based multi-aspect multi-granularity non-native English speaker pronunciation assessment,”  in ICASSP 2022-2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2022, pp. 7262–7266.
  * [13] Binghuai Lin, Liyuan Wang, Xiaoli Feng, and Jinsong Zhang,  “Automatic scoring at multi-granularity for L2 pronunciation.,”  in Interspeech, 2020, pp. 3022–3026.
  * [14] Fu-An Chao, Tien-Hong Lo, Tzu-I Wu, Yao-Ting Sung, and Berlin Chen,  “3M: An effective multi-view, multi-granularity, and multi-aspect modeling approach to english pronunciation assessment,”  arXiv preprint arXiv:2208.09110, 2022.
  * [15] Daniel Zhang, Ashwinkumar Ganesan, Sarah Campbell, Daniel Korzekwa, and Amazon Alexa AI,  “L2-GEN: A neural phoneme paraphrasing approach to L2 speech synthesis for mispronunciation diagnosis,”  Proc. Interspeech 2022, pp. 4317–4321, 2022.
  * [16] Kaiqi Fu, Shaojun Gao, Kai Wang, Wei Li, Xiaohai Tian, and Zejun Ma,  “Improving non-native word-level pronunciation scoring with phone-level mixup data augmentation and multi-source information,”  arXiv preprint arXiv:2203.01826, 2022.
  * [17] Zhan Zhang, Yuehai Wang, and Jianyi Yang,  “End-to-end mispronunciation detection with simulated error distance,”  Proc. Interspeech 2022, pp. 4327–4331, 2022.
  * [18] Xiaoshuo Xu, Yueteng Kang, Songjun Cao, Binghuai Lin, and Long Ma,  “Explore wav2vec 2.0 for mispronunciation detection.,”  in Interspeech, 2021, pp. 4428–4432.
  * [19] Guanlong Zhao, Sinem Sonsaat, Alif Silpachai, Ivana Lucic, Evgeny Chukharev-Hudilainen, John Levis, and Ricardo Gutierrez-Osuna,  “L2-ARCTIC: A non-native English speech corpus.,”  in INTERSPEECH, 2018, pp. 2783–2787.
  * [20] Sepp Hochreiter and Jürgen Schmidhuber,  “Long short-term memory,”  Neural computation, vol. 9, no. 8, pp. 1735–1780, 1997.
  * [21] Junbo Zhang, Zhiwen Zhang, Yongqing Wang, Zhiyong Yan, Qiong Song, Yukai Huang, Ke Li, Daniel Povey, and Yujun Wang,  “speechocean762: An open-source non-native English speech corpus for pronunciation assessment,”  arXiv preprint arXiv:2104.01378, 2021.
  * [22] Daniel Povey, Arnab Ghoshal, Gilles Boulianne, Lukas Burget, Ondrej Glembek, Nagendra Goel, Mirko Hannemann, Petr Motlicek, Yanmin Qian, Petr Schwarz, et al.,  “The kaldi speech recognition toolkit,”  in IEEE 2011 workshop on automatic speech recognition and understanding. IEEE Signal Processing Society, 2011, number CONF.


# Domain Selection for Gaussian Process Data: An application to
electrocardiogram signals

Nicolás Hernández
Department of Statistical Science
University College London
London, UK
<EMAIL_ADDRESS>
& Gabriel Martos
Departmento de Matemáticas y Estadística
Universidad Torcuato Di Tella
Buenos Aires, Argentina
<EMAIL_ADDRESS>

###### Abstract

Gaussian Processes and the Kullback-Leibler divergence have been deeply
studied in Statistics and Machine Learning. This paper marries these two
concepts and introduce the local Kullback-Leibler divergence to learn about
intervals where two Gaussian Processes differ the most. We address subtleties
entailed in the estimation of local divergences and the corresponding interval
of local maximum divergence as well. The estimation performance and the
numerical efficiency of the proposed method are showcased via a Monte Carlo
simulation study. In a medical research context, we assess the potential of
the devised tools in the analysis of electrocardiogram signals.

_Keywords_ Domain selection $\cdot$ Gaussian processes $\cdot$ Kullback-
Leibler divergence $\cdot$ intervals of local maximum divergence $\cdot$
electrocardiogram signals

##  1 Introduction

Everyday millions of complex data patterns flow around the world at
unprecedented speed, leading to an explosion on the demand for modelling
random process data, such as time series and functional data.
Electrocardiogram signals (ECG) are an example of such high-dimensional data,
that is usually structured in the form of curves almost continuously recorded
over a grid of discrete time points. From a medical point of view, when a
patient is admitted in the emergency room during a cardiac arrest, one of the
few pieces of information available to make a diagnosis is an ECG. It is
therefore an extremely useful tool for immediate decision making and it can
also be helpful in determining the causes and the gravity of a cardiac
pathology Mullainathan and Obermeyer (2022). Nevertheless, the analysis of ECG
data faces important challenges in practice, in particular due to its high
dimensionality. Therefore, the study of local features of ECG signals play a
key role in Medicine Wang et al. (2013); Rodríguez et al. (2015); Mullainathan
and Obermeyer (2022) for at least two important reasons:

  1. i)

Diagnosis: Identifying time intervals during the cardiac cycle where the ECG
signal present atypical patterns is crucial in order to improve the early
diagnosis of different cardiac diseases, increasing the patient survival
probability during a cardiac episode.

  2. ii)

Causes and effects: Learning time intervals with an anomalous pattern during
the cardiac cycle will help to understand the origins and the consequences of
different heart diseases.

Taking into account (i) and (ii) above, the main goal of this paper is to
learn from ECG data an interval with a certain length where the signals
corresponding to disease and healthy subjects differ the most. Throughout this
paper, we refer to this problem as that of domain selection for ECG signals,
but the method devised in this paper also apply to other random processes data
in a broad sense. We model ECG signals using Gaussian Processes (GPs), a
versatile and flexible tool for modelling complex signal patterns Pérez-Cruz
et al. (2013), and introduce the _local_ Kullback-Leibler (KL) divergence
Kullback and Leibler (1951) as we formally discussed in Section 2, so to learn
intervals of maximum local divergence between GPs. Recent contributions in the
related literature explore domain selection methods–a.k.a. variable selection
in functional data Berrendero et al. (2016); Baíllo et al. (2011); Pini and
Vantini (2017)– to achieve accurate prediction for functional data
classification methods and to assess local differences between functional
means in a two sample problem. Also in the context of statistics in medicine,
the authors in Martos and de Carvalho (2018) propose a Mann-Whitney type of
statistic for functional data to learn about intervals at which two processes
differ the most, based on aspects related with symmetry. Our approach differs
from the ones mentioned above in at least two important ways: (1) Here the
ultimate goal is neither to classify nor to test hypothesis with random
processes data, but rather to learn the interval with a given length where two
random processes, which corresponds to groups of ECG signals, differ the most;
(2) Our approach relies on GP and the KL divergence, whereas the
aforementioned methodologies have mainly been designed in the context of
functional data. As a byproduct, we also contribute on the following points:

  * •

Optimisation: We introduce the interval of local maximum divergence through a
set function optimisation problem. Therefore, the estimation methods proposed
in the paper contributes to the literature on set function optimisation.

  * •

Classification: When the analysis of ECG signals also entails the
discrimination between groups (i.e. healthy vs disease), our method could
benefit other standard functional classifiers if they are applied on a small
interval where the two processes differ the most, rather than treating the
entire signal domain equally.

  * •

Storage efficiency: If only a subset of the entire domain is found to be
relevant in order to assess differences between healthy subjects and disease
patients, then this suggests the potential benefit for collecting and saving
only a smaller subset of ECG signals.

  * •

Miscellany: In the paper we also establish conditions for the existence of an
interval of local maximum divergence, consider subtleties entailed in the
estimation of local KL divergences and the interval of local maximum
divergence; and also discuss variants and extensions around the proposed
methodology.

The remind of the paper is organised as follows: In Section 2 we introduce the
local KL divergence for GP and the interval of local maximum divergence, and
also discuss suitable corresponding estimation methods. In Section 3 we
present Monte Carlo evidence to assess the consistency of our estimator, while
in Section 4 we illustrate the method with an electrocardiogram signals data
application. Finally in Section 5 we discuss the results and conclude our
work.

##  2 Materials and Methods

The goal in this section is to introduce the probabilistic framework to assess
local differences between GPs. To this end, let $X\sim N(\mu_{X},\sigma_{X})$
and $Y\sim N(\mu_{Y},\sigma_{Y})$ be two normally distributed independent
random variables, where $\mu_{\ell}$ and $\sigma_{\ell}$ for
$\ell\in\\{X,Y\\}$ denotes the corresponding mean and variance parameters;
then the KL divergence can be written in a closed from as follows:

$\text{KL}(X||Y)=\frac{1}{2}\left(\frac{\sigma^{2}_{X}}{\sigma^{2}_{Y}}-1+\frac{(\mu_{X}-\mu_{Y})^{2}}{\sigma^{2}_{Y}}+\ln\left(\frac{\sigma^{2}_{Y}}{\sigma^{2}_{X}}\right)\right).$
(1)

The divergence in Eq. (1) is a functional that quantifies the dissimilarity in
the distribution of two Gaussian random variables $X$ and $Y$; and more
importantly Eq. (1) is easy and computationally cheap to evaluate when $X$ and
$Y$ are multivariate normal random vectors. Next we discuss how to extrapolate
Eq. (1) to compute intervals of local maximum divergence for GPs.

###  2.1 Local Kullback-Leibler divergence for Gaussian Processes

Let $X(t)\sim GP(\mu_{X}(t),\sigma_{X}(t,s))$ be a GP with mean function
$\mu_{X}(t)=E\\{X(t)\\}$ and variance function
$\sigma_{X}(t,s)=E\\{(X(t)-\mu_{X}(t))(X(s)-\mu_{X}(s))\\}$, indexed on the
compact set $T\subset\mathbb{R}$; we address the following problem on the
comparison of two GP:

> Learning Problem. Given data drawn from GPs $X(t)\sim
> GP(\mu_{X}(t),\sigma_{X}(t,s))$ and $Y(t)\sim
> GP(\mu_{Y}(t),\sigma_{Y}(t,s))$, indexed on the same compact domain
> $T\subset\mathbb{R}$, learn the interval with a given certain length where
> they statistically differ the most.

Since the mean and variance functions completely determines the law of $X(t)$
and $Y(t)$; hereafter we assume that there exists a compact subset $A\subset
T$ with $\lambda(A)>0$, where $\lambda(A)$ is the Lebesgue measure of the set
$A$, and a positive constant $\nu$, such that one of the following scenarios
of local differences holds:

  1. (A)

For every $t\in A$: $|\mu_{X}(t)-\mu_{Y}(t)|>\nu$, while for all
$t^{\prime}\notin A$ it holds that
$|\mu_{X}(t^{\prime})-\mu_{Y}(t^{\prime})|\leq\nu$.

  2. (B)

For every $(t,s)\in A\times A$: $|\sigma_{X}(t,s)-\sigma_{Y}(t,s)|>\nu$, while
for all $(t^{\prime},s^{\prime})\notin A\times A$ it holds that
$|\sigma_{X}(t^{\prime},s^{\prime})-\sigma_{Y}(t^{\prime},s^{\prime})|\leq\nu$.

  3. (C)

Local mean and variance differences corresponds to scenario (A) and (B)
simultaneously, possibly on different subsets $A_{\mu}$ and $A_{\sigma}$ and
for different thresholds constant $\nu_{\mu}$ and $\nu_{\sigma}$.

In Figure 1–(a) we illustrate these scenario A: The coloured lines represents
the mean functions of two GP; notice that over a relatively small time
interval–around the time point $t=1.5$– the difference
$|\mu_{X}(t)-\mu_{Y}(t)|$ exceeds a certain threshold $\nu_{\mu}$. In Figure
1-(b) we depict scenario B: The coloured surfaces represents the covariance
functions of two GP; notice that over a relatively small time
interval–highlighted with a black square on the top– the difference
$|\sigma_{X}(t,s)-\sigma_{Y}(t,s)|$ is above a certain threshold
$\nu_{\sigma}$. The scenario C corresponds to Figure 1–(a) and (b)
simultaneously.

Figure 1: (a) Mean functions $\mu_{X}(t)$ and $\mu_{Y}(t)$. (b) Variance
functions $\sigma_{X}(t,s)$ and $\sigma_{Y}(t,s)$. The corresponding interval
of local maximum divergence for $c=0.3$ depicted in brackets in (a) and with a
black square in (b).

The goal of this paper is to develop a local KL divergence measure and related
estimation methods to learn about intervals or regions, such as $A^{*}(c)$ in
Figure 1-(a) and (b), where two GP’s differ the most (i.e. where the
functional parameters are the most dissimilar). To this end, it will be
technically convenient and computationally efficient to use a finite-
dimensional (or discrete) representation of GP data. Without loss of
generality, we assume that data is recorded and stored over the same discrete
and equally spaced grid of points $\mathcal{T}=(t_{1},\dots,t_{p})\subset T$,
being $p\gg 0$ (other more general cases are easy to tackle as we discuss in
Section 3). In this setting, GP data corresponds to realisations of
$p$-variate Gaussian random vectors
$\textbf{X}_{\mathcal{T}}\equiv(X(t_{1}),\dots,X(t_{p}))\sim
N_{p}(\boldsymbol{\mu}_{\mathcal{T},X},\boldsymbol{\Sigma}_{\mathcal{T},X})$
and $\textbf{Y}_{\mathcal{T}}\equiv(Y(t_{1}),\dots,Y(t_{p}))\sim
N_{p}(\boldsymbol{\mu}_{\mathcal{T},Y},\boldsymbol{\Sigma}_{\mathcal{T},Y})$
where $\boldsymbol{\mu}_{\mathcal{T},X}=(\mu_{X}(t_{1}),\dots,\mu_{X}(t_{p}))$
and $\boldsymbol{\mu}_{\mathcal{T},Y}=(\mu_{Y}(t_{1}),\dots,\mu_{Y}(t_{p}))$
are the corresponding means, and
$\boldsymbol{\Sigma}_{\mathcal{T},X}\in\mathbb{R}^{p\times p}$ and
$\boldsymbol{\Sigma}_{\mathcal{T},Y}\in\mathbb{R}^{p\times p}$ the
corresponding $p\times p$ variance matrices (i.e.
$[\boldsymbol{\Sigma}_{\mathcal{T},X}]_{ij}=\sigma_{X}(t_{i},t_{j})$ and
$[\boldsymbol{\Sigma}_{\mathcal{T},Y}]_{ij}=\sigma_{Y}(t_{i},t_{j})$ for
$i=1,\dots,p$ and $j=1,\dots,p$).

Under this high dimensional GP’s representation, the KL divergence Kullback
and Leibler (1951) is a natural metric to assess differences in distribution
between two GP’s. The KL divergence between GP $X$ and $Y$ over the grid
$\mathcal{T}$ is computed as follows Pardo (2018):

$\text{KL}_{\mathcal{T}}(X||Y)\equiv{\frac{1}{2}}\left(\operatorname{tr}\left(\boldsymbol{\Sigma}_{\mathcal{T},Y}^{-1}\boldsymbol{\Sigma}_{\mathcal{T},X}-\textbf{I}_{p}\right)+\Delta_{\mathcal{T}}^{\mathsf{T}}\boldsymbol{\Sigma}_{\mathcal{T},Y}^{-1}\Delta_{\mathcal{T}}+\ln\left({\frac{\det\boldsymbol{\Sigma}_{\mathcal{T},Y}}{\det\boldsymbol{\Sigma}_{\mathcal{T},X}}}\right)\right),$
(2)

where
$\Delta_{\mathcal{T}}=\left(\boldsymbol{\mu}_{\mathcal{T},Y}-\boldsymbol{\mu}_{\mathcal{T},X}\right)$,
$\operatorname{tr}(\boldsymbol{\Sigma})$ denote the trace of
$\boldsymbol{\Sigma}$ and $\operatorname{det}(\boldsymbol{\Sigma})$ the
determinant of $\boldsymbol{\Sigma}$. Some comments on the KL divergence are
in order: (i) The expression in Eq. (2) is a generalisation of the
corresponding univariate KL divergence in Eq. (1). (ii) The KL divergence is
not symmetric, nevertheless the symmetrisation is straightforward: consider
for instance
$\text{KL}_{\mathcal{T}}(X||Y)/2+\text{KL}_{\mathcal{T}}(Y||X)/2$. (iii)
Interestingly, the KL divergence considers simultaneously differences in mean
and variance, i.e. both mean and variance appear together in Equation (2).
(iv) In the case of GP with the same covariance function, then
$2\text{KL}_{\mathcal{T}}(X||Y)=\Delta_{\mathcal{T}}^{\mathsf{T}}\boldsymbol{\Sigma}_{\mathcal{T},Y}^{-1}\Delta_{\mathcal{T}}$
corresponds to the squared Mahalanobis distance McLachlan (1999) between the
two GP.

Our goal is to learn subsets of $T$ where the two processes differ the most,
therefore the KL divergence is a suitable “statistical distance" to assess
local differences between GP. For any subset
$\mathcal{A}\subseteq\mathcal{T}$, we define the local-KL divergence as
follows:

$\text{KL}_{\mathcal{A}}(X||Y)\equiv{\frac{1}{2}}\left(\operatorname{tr}\left(\boldsymbol{\Sigma}_{\mathcal{A},Y}^{-1}\boldsymbol{\Sigma}_{\mathcal{A},X}-\textbf{I}_{|\mathcal{A}|}\right)+\Delta_{\mathcal{A}}^{\mathsf{T}}\boldsymbol{\Sigma}_{\mathcal{A},Y}^{-1}\Delta_{\mathcal{A}}+\ln\left({\frac{\det\boldsymbol{\Sigma}_{\mathcal{A},Y}}{\det\boldsymbol{\Sigma}_{\mathcal{A},X}}}\right)\right),$
(3)

where
$\Delta_{\mathcal{A}}=\left(\boldsymbol{\mu}_{\mathcal{A},Y}-\boldsymbol{\mu}_{\mathcal{A},X}\right)$
and $\\{\boldsymbol{\mu}_{\mathcal{A}},\boldsymbol{\Sigma}_{\mathcal{A}}\\}$
denotes suitable partitions of
$\\{\boldsymbol{\mu}_{\mathcal{T}},\boldsymbol{\Sigma}_{\mathcal{T}}\\}$ in
correspondence with the subset $\mathcal{A}\subseteq\mathcal{T}$. Local KL
divergences for GP’s data also have a number of interesting properties which
we summarise next.

######  Proposition 2.1.

The $\text{KL}_{\mathcal{A}}(X||Y)$ is a set function that satisfies the
following properties:

  1. (a)

Non-negative: For fixed GPs $X$ and $Y$, it holds:
$\text{KL}_{\mathcal{A}}(X||Y):\mathcal{P}_{\mathcal{T}}\to\mathbb{R}^{+}_{0}$
where $\mathcal{P_{T}}$ is the power set of $\mathcal{T}$; and
$\text{KL}_{\mathcal{A}}(X||Y)=0$ if and only if $\mu_{X}(t)=\mu_{Y}(t)$ for
all $t\in\mathcal{A}$ and $\sigma_{X}(t,s)=\sigma_{Y}(t,s)$ for all
$(t,s)\in\mathcal{A}\times\mathcal{A}$.

  2. (b)

The local KL divergence is upper bounded (i.e.
$\text{KL}_{\mathcal{T}}(X||Y)<\infty)$ and a monotone set function (i.e. for
$\mathcal{A}^{\prime}\subseteq\mathcal{A}$ it holds:
$\text{KL}_{\mathcal{A}^{\prime}}(X||Y)\leq\text{KL}_{\mathcal{A}}(X||Y)$),
under suitable conditions on the mean and variance functions.

  3. (c)

The local divergence $\text{KL}_{\mathcal{A}}(X||Y)$ is a continuous set
function in $\mathcal{C}_{\mathcal{T}}$, the collection of all contiguous
subsets from the ground set $\mathcal{T}$.

In the Appendix we give formal proofs on previous assertions. Since
$\text{KL}_{\mathcal{A}}(X||Y)$ accounts simultaneously for differences in
mean and variance, we propose this particular metric to define the most
discriminating subset of points in $\mathcal{T}$ as follows:

######  Definition 2.1.

Let $|\cdot|$ be the counting measure on $\mathcal{P}_{\mathcal{T}}$, we
define the subset $\widetilde{\mathcal{A}}^{*}(\tilde{c})\subset\mathcal{T}$,
for any $1\leq\tilde{c}\leq|\mathcal{T}|$, as the variable selection subset
that solves the following set function optimisation problem:

$\max_{\mathcal{A}\subset\mathcal{T}}\text{KL}_{\mathcal{A}}(X||Y),\text{ s.t.
}|\mathcal{A}|\leq\tilde{c}.$ (4)

Interestingly, this definition resembles some developments in the context of
variable selection for functional data classification Berrendero et al.
(2016), where the authors propose a maxima–hunting approach to search for
isolated time points in the domain corresponding to random processes $X(t)$
and $Y(t)$ where the covariance distance Székely et al. (2007) is maximal. Our
approach is somehow similar in the sense that we maximise a metric that
accounts for local differences between the processes, but differs in two
important ways: (i) we work under the GP assumption in order to rely on a fast
and easy metric to compute (the KL divergence); and (ii) our main goal is to
learn about the local maximum divergence interval. Some additional comments
are in order: First, the discrete set $\widetilde{\mathcal{A}}^{*}(\tilde{c})$
does not necessarily corresponds to an interval; and second, Eq. (4) entails a
cumbersome combinatorial set function optimisation problem even for moderate
values of $p$. Since the goal of the paper is to select a subset of the domain
of the GPs instead of isolated points, we put some additional structure on the
shape of our candidate set $\mathcal{A}^{*}$ as in the following definition.

######  Definition 2.2.

Consider $\mathcal{C}_{\mathcal{T}}$ as the collection of all contiguous
subsets from the ground set $\mathcal{T}$, such that for any
$\mathcal{A}\in\mathcal{C}_{\mathcal{T}}$ then
$\mathcal{A}=\\{t_{k},t_{k+1},\dots t_{k+l}\\}$ for some $1\leq k\leq p$ and
$0\leq l\leq p-k$. Then, the interval of local maximum KL divergence of size
$c\in(0,1)$, denoted onward as $\mathcal{A}^{*}(c)$, is defined throughout the
following set function optimisation problem:

$\max_{\mathcal{A}\in\mathcal{C}_{\mathcal{T}}}\text{KL}_{\mathcal{A}}(X||Y),\text{
s.t. }\text{len}(\mathcal{A})\leq c\lambda(T),$ (5)

where
$\text{len}(\mathcal{A})=\max_{t\in\mathcal{A}}(t)-\min_{t\in\mathcal{A}}(t)$
is the length function corresponding to subset $\mathcal{A}$.

The parameter $c\in(0,1)$ in Eq. (5) determines the proportion of the domain
to be selected and avoids the use of a particular length scale as is the case
of $\tilde{c}$ in Eq. (4). In the empirical Section 3 we discuss data driven
strategies to choose the value of $c$. The existence of $\mathcal{A}^{*}(c)$
follows from the result stated in Proposition–2.1 point (b) and the fact that
$\mathcal{P}_{\mathcal{T}}$ is finite (see the Appendix for further details).
However, the interval of local maximum divergence $\mathcal{A}^{*}(c)$ does
not need to be unique, as can be easily seen by considering the limiting case
where GP’s $X(t)$ and $Y(t)$ have the same mean and variance functions. In
such a case every compact subset of $T$ with measure $c$ is a set of maximum
KL local divergence. Moreover, one can learn several and disjoints intervals
of maximum KL local divergence by applying sequential learning. This means
selecting a new interval of maximum KL local divergence, after discarding
previously selected intervals.

###  2.2 Learning intervals with local maximum KL divergence from data

To learn $\mathcal{A}^{*}(c)$ from data we consider samples recorded over the
same discrete grid $\mathcal{T}$:
$\mathcal{D}_{X}=\\{\textbf{x}_{i}\\}_{i=1}^{n}$ and
$\mathcal{D}_{Y}=\\{\textbf{y}_{j}\\}_{j=1}^{m}$ drawn from
$GP(\mu_{X}(t),\sigma_{X}(t,s))$ and $GP(\mu_{Y}(t),\sigma_{Y}(t,s))$
respectively (more general sampling designs are discussed in below). Let
$\boldsymbol{\Psi}_{\mathcal{T},X}=\\{\boldsymbol{\mu}_{\mathcal{T},X},\boldsymbol{\Sigma}_{\mathcal{T},X}\\}$
be the true parameters corresponding to GP $X(t)$, the Maximum Likelihood
estimates are given by:

$\widehat{\boldsymbol{\Psi}}_{\mathcal{T},X}=\left\\{\widehat{\boldsymbol{\mu}}_{\mathcal{T},X}=\frac{1}{n}\sum_{i=1}^{n}\mathbf{x}_{i};\quad\widehat{\boldsymbol{\Sigma}}_{\mathcal{T},X}=\frac{1}{n}\sum_{i=1}^{n}(\mathbf{x}_{i}-\widehat{\boldsymbol{\mu}}_{\mathcal{T},X})(\mathbf{x}_{i}-\widehat{\boldsymbol{\mu}}_{\mathcal{T},X})^{T}\right\\},$

and an analogous expression holds for
$\widehat{\boldsymbol{\Psi}}_{\mathcal{T},Y}$ as well. Plugging suitable
partitions $\widehat{\boldsymbol{\Psi}}_{\mathcal{A},X}$ and
$\widehat{\boldsymbol{\Psi}}_{\mathcal{A},Y}$ corresponding to the subset
$\mathcal{A}$ into Eq. (5) gives us an estimate of
$\widehat{\text{KL}}_{\mathcal{A}}(X||Y)$. Notice that the trace and
determinant are continuous functions in the space of real symmetric matrices,
then $\widehat{\text{KL}}_{\mathcal{A}}(X||Y)$ inherits interesting
statistical properties, in particular for a fixed $p$, consistency, asymptotic
normality and efficiency, providing that $p/n\to 0$, $p/m\to 0$ and
$m/n\to\rho\in(0,1)$ as $n\to\infty$ and $m\to\infty$. Having in hand a
consistent estimator of the local KL divergence, to learn about
$\mathcal{A}^{*}(c)$, for $c\in(0,1)$, we follow the step in provided in
Algorithm 1.

Algorithm 1 Estimating the interval of local maximum KL divergence from GP
data.

Inputs:

$\mathcal{D}_{X}$ and $\mathcal{D}_{Y}$ and a length constrain $0<c\leq 1$.

Step 1:

Parameter estimation.

  1. IF:

$\mathcal{D}_{X}$ and $\mathcal{D}_{Y}$ are recorded over the same grid
$\mathcal{T}$, compute Maximum Likelihood estimates
$\widehat{\Psi}_{\mathcal{T},X}$ and $\widehat{\Psi}_{\mathcal{T},Y}$ from
data.

  2. ELSE:

Use your favourite (semi/non) parametric estimation method to compute
$\widehat{\mu}_{\ell}(t)$ and $\widehat{\sigma}_{\ell}(t,s)$ from
$\mathcal{D}_{\ell}$ for $\ell=\\{X,Y\\}$. Define a suitable common grid
$\mathcal{T}$ and compute Maximum Likelihood estimates
$\widehat{\Psi}_{\mathcal{T},X}$ and $\widehat{\Psi}_{\mathcal{T},Y}$ via the
evaluation of $(\widehat{\mu}_{\ell}(t),\widehat{\sigma}_{\ell}(t,s))$ over
$\mathcal{T}$ for $\ell=\\{X,Y\\}$ respectively.

Step 2:

Exhaustive optimisation: Let $\mathcal{C}_{\mathcal{T}}(c)$ be the set of all
contiguous subsets in $\mathcal{T}$ such as $\text{len}(\mathcal{A})\leq
c\lambda(T)$ if $\mathcal{A}\in\mathcal{C}_{\mathcal{T}}(c)$, then for all
$\mathcal{A}\in\mathcal{C}_{\mathcal{T}}(c)$ compute
$\widehat{\text{KL}}_{\mathcal{A}}(X||Y)$ and return:
$\widehat{\mathcal{A}}^{*}_{c}\equiv\operatorname{arg\,max}_{\mathcal{A}\in\mathcal{C}_{\mathcal{T}}(c)}\widehat{\text{KL}}_{\mathcal{A}}(X||Y).$

Algorithm 1 warrants some comments: Step 1 is carried only once even in the
case of estimating several intervals of local maximum KL divergence for
different values of $c$. In our R implementation, the set function
optimisation problem in Step 2 is solved via the evaluation of
$\widehat{\text{KL}}_{\mathcal{A}}(X||Y)$ over all sets in
$\mathcal{C}_{\mathcal{T}}(c)$; nevertheless other more efficient
derivative–free optimisation approaches Nocedal and Wright (1999) can be
considered as well. In Section 3 we assess the computational efficiency of the
estimation method entailed in Algorithm 1 over different data generating
scenarios.

#### The curse of dimensionality

In the context of GP data, usually $n$ and $m$ are relatively small in
comparison to $p=|\mathcal{T}|$, being often the case where
$p\gg\max\\{n,m\\}$. Since the estimation of $\mathcal{A}^{*}(c)$ depends on
the estimation of two $p\times p$ covariance matrices, some remedy actions are
necessary in order to obtain suitable estimations for
$\boldsymbol{\Sigma}_{\mathcal{T},X}$ and
$\boldsymbol{\Sigma}_{\mathcal{T},Y}$ from data whenever $p$ is relatively
large in comparison to sample sizes $n$ and $m$. A well known strategy is to
impose structure in the GP’s covariance functions, assuming for instance that
both processes are conditionally independent in time (i.e.
$\sigma_{X}(t,s)=\sigma_{Y}(t,s)=0$ for all $t\neq s$), which corresponds to
assume $\boldsymbol{\Sigma}_{\mathcal{T},X}$ and
$\boldsymbol{\Sigma}_{\mathcal{T},Y}$ are diagonal covariance matrices.
Another less contrived approach to circumvent the curse of dimensionality is
to consider a penalised likelihood covariance matrix estimator as follows:

$\widehat{\boldsymbol{\Sigma}}_{\eta}=\eta\,\widehat{\boldsymbol{\Sigma}}+(1-\eta)\,\text{diag}(\widehat{\boldsymbol{\Sigma}}),\text{
for }\eta\in[0,1],$

where $\eta$ is a regularisation parameter that shrinks the Maximum Likelihood
estimator of the variance matrix towards its diagonal Hastie et al. (2009).
The value of $\eta$ is typically determined using cross validation methods.
Nevertheless, other approaches such as Banding, Tapering and alternative
thresholding methods are also available, we refer to Pourahmadi (2013)[Ch. 6]
and references therein for more details.

#### Sampling designs and missalignments

In order to simplify the exposition, in section 2.1, we assume that $X(t)$ and
$Y(t)$ are recorded over the same equally space time point grid $\mathcal{T}$,
but other sampling designs are also frequent in practice. In such cases where
the processes are not recorded over the same grid, to use Algorithm 1, we
first need to define a common auxiliary grid $\mathcal{T}$ and then proceed as
follows:

  1. 1.

Use a suitable smoothing technique such as kernel smoothing, smoothing
Splines, or GP regression among many others; and estimate the functional
parameters $\mu_{\ell}(t)$ and $\sigma_{\ell}(t,s)$ for $\ell=\\{X,Y\\}$ from
data.

  2. 2.

Consider
$\widehat{\boldsymbol{\mu}}_{\ell}=(\widehat{\mu}_{\ell}(t_{1})\dots,\widehat{\mu}_{\ell}(t_{p}))$
and
$[\widehat{\boldsymbol{\Sigma}}]_{i,j}=\widehat{\sigma}_{\ell}(t_{i},t_{j})$,
for $\ell=\\{X,Y\\}$ and $i=\\{1,\dots,p\\}$ and $j=\\{1,\dots,p\\}$ (i.e. the
evaluation of the estimated functional parameters over the grid
$\mathcal{T}$).

After the estimation of mean vectors and covariance matrices over a common
grid $\mathcal{T}$, Step 2 in Algorithm 1 follows straightforwardly. In
principle the auxiliary grid $\mathcal{T}$ contains equally spaced time points
and has length $p$ according to the resolution defined by the user.
Nevertheless, other sampling designs can be considered as well, for instance
sampling $\mathcal{T}$ from a multivariate prior distribution
$\pi_{\mathcal{T}}(t_{1},\dots,t_{p})$. This prior distribution encodes
knowledge about the most local divergent interval in $T$ by putting more
probability mass over a specific interval $A\subset T$. This procedure is a
natural avenue to study a Bayesian extension of the proposed method to select
domain with GP.

GP data are sometimes recorded with different types of random amplitude and
time variations that produce a misalignment in data Ramsay and Li (1998). The
GP data asynchrony may act as a confounding factor when the aim is the
estimation of the interval of local maximum divergence, since the mean and
variance functions estimators defined in § 2.2 are useless. In this case,
before the implementation of Step 1 in Algorithm 1, we suggest to pre-process
GP data using standard alignment or synchronisation tools (a.k.a. as curve
registration in functional data) such as the methods described in Berndt and
Clifford (1994); Kazlauskaite et al. (2019) among others.

#### Smoothing

GP data sometimes present a low signal to noise ratio leading to problems in
the estimation of the interval of local maximum divergence. To alleviate low
signal to noise ratio issues, we recommend to smooth GP data using standard
methods such as Natural Splines, or Kernel Regression, among others Ullah and
Finch (2013). The smoothing process also enables us to use a finer grid of
points (see sampling designs paragraph). Once the smoothed data is obtained,
they can be used as the input in Algorithm 1.

###  2.3 Inference and prediction

#### Inference

The assess the variability of $\widehat{\mathcal{A}}^{*}(c)$ as an estimator
of $\mathcal{A}^{*}(c)$, we resort to non-parametric bootstrap techniques. For
every $c\in(0,1)$, the interval of local maximum divergence and its
corresponding estimator can be parametrised in terms of a ball with centre
$t_{c}$ and radius111Notice that $r_{c}$ is a nuisance parameter that only
depends on $c$. $r_{c}$ as follows:
$B(t_{c}^{*},r_{c})\equiv[\min_{t\in\mathcal{A}^{*}(c)}t,\max_{t\in\mathcal{A}^{*}(c)}t]$
which corresponds to $\mathcal{A}^{*}(c)$; and
$B(\hat{t}_{c}^{*},r_{c})\equiv[\min_{t\in\widehat{\mathcal{A}}^{*}(c)}t,\max_{t\in\widehat{\mathcal{A}}^{*}(c)}t]$
which corresponds to $\widehat{\mathcal{A}}^{*}(c)$. Then a $1-\alpha$
bootstrap confidence interval for $t^{*}_{c}$ is given by
$\text{CI}_{1-\alpha}(t^{*}_{c})\equiv[\hat{t}^{*}_{c,\alpha/2},\hat{t}^{*}_{c,1-\alpha/2}]$,
where $\hat{t}^{*}_{c,\alpha}$ is the $\alpha\in(0,1)$ quantile corresponding
to the empirical distribution $\widehat{F}_{B}(t^{*}_{c})$ obtained using
$B\gg 0$ bootstrap samples from GP data. From the later confidence interval,
we define our $1-\alpha$ confidence set for $\mathcal{A}^{*}(c)$ as
$\text{CS}_{1-\alpha}(\mathcal{A}^{*}(c))\equiv\cup_{t\in\text{CI}_{1-\alpha}(t^{*}_{c})}B(t,r_{c})$.
In the experimental sections we report $\widehat{F}_{B}(t^{*}_{c})$ and
$\text{CI}_{1-\alpha}(t^{*}_{c})$, the computation of
$\text{CS}_{1-\alpha}(\mathcal{A}^{*}(c))$ is straightforward from the later
confidence interval.

Some additional comments on the uncertainty quantification are in order.
Regardless the experimental design, there is a trade-off between the value of
$c$ and the variability of our estimator. For instance, for $p=100$ and
$c=0.1$ there are 91 possible intervals in the domain, while for $c=0.9$ there
are only 11 intervals. Hence, the uncertainty associated to
$\mathcal{\widehat{A}}^{*}(c)$ converges to zero as $c\to 1$–these corresponds
to narrower intervals for $t^{*}_{c}$ as $c\to 1$. Other confidence interval
based on bootstrap procedures can be considered as well, such as the
percentile or the Student-t method. Moreover, the parametric bootstrap is also
another possible approach, taking into account the Gaussian assumptions.
Notice that bootstrap confidence intervals are neither exact nor optimal, but
they are largely used in practice since the method is easy to implement and
its accuracy is near-exact. For a general discussion on the asymptotic
properties regarding the coverage probability of a bootstrap confidence
interval we refer to Efron and Tibshirani (1994).

#### Domain selection and GP classification

Domain selection is an important preliminary step before training a
classification model Berrendero et al. (2016). In the context of GP data, the
Discriminant Analysis (DA) is the Bayes Optimal Classifier Fraley and Raftery
(2002); hence to classify a new unlabelled instance $\mathbf{z}$–a realisation
from $X(t)$ or $Y(t)$ recorded over $\mathcal{T}$–, the DA considers the sign
of the following discriminant function:

$D_{\mathcal{T}}(\mathbf{z})=\frac{1}{2}\left((\mathbf{z}-\boldsymbol{\mu}_{\mathcal{T},Y})^{T}\boldsymbol{\Sigma}_{\mathcal{T},Y}^{-1}(\mathbf{z}-\boldsymbol{\mu}_{\mathcal{T},Y})-(\mathbf{z}-\boldsymbol{\mu}_{\mathcal{T},X})^{T}\boldsymbol{\Sigma}_{\mathcal{T},X}^{-1}(\mathbf{z}-\boldsymbol{\mu}_{\mathcal{T},X})-\ln\bigg{(}\frac{\text{det}\,\boldsymbol{\Sigma}_{\mathcal{T},X}}{\text{det}\,\boldsymbol{\Sigma}_{\mathcal{T},Y}}\Bigg{)}\right)+\ln\bigg{(}\frac{\pi_{X}}{\pi_{Y}}\bigg{)},$
(6)

where $\pi_{X}$ and $\pi_{Y}=1-\pi_{X}$ are the corresponding prior
probabilities. $D_{\mathcal{T}}(\mathbf{z})$ encodes the rule to classify
$\mathbf{z}$ as generated from $X(t)$ or $Y(t)$. In high dimensional contexts
($p\gg 0$), an important drawback of DA is the lack of reliable estimates of
the involved mean vectors and variance matrices. This leads to inaccurate
classification results in case no action is taken on the course of
dimensionality. In this regard, it will be convenient to use only a small
compact subset $A\subset T$ of GP data so as to compute the discrimination
function. This corresponds to replacing $D_{\mathcal{T}}(\mathbf{z})$ by
$D_{\mathcal{A}^{*}(c)}(\mathbf{z})$ for a suitable value $c\in(0,1)$,
reducing the number of parameter entailed in the discrimination analysis.
Choosing a suitable value for $c$ will be crucial to obtain accurate
classification results; in Section 4 we show how to fix $c$ using cross
validation methods to improve the discrimination power between normal and
myocardial infarction heartbeats. It is interesting to mention that a choice
of $c\text{ such that }\tilde{c}=c\lambda(T)=1$, would transform the domain
selection into a variable selection problem, that is, selecting a variable on
the domain of the processes. While this procedure remains valid, it does not
allow for the exploitation of the existing information in the covariance
matrix of the processes by considering only the variability at the selected
variable.

##  3 Simulation Study

In this section we assess the performance of Algorithm 1 in the estimation of
$\mathcal{A}^{*}(c)$ via a Monte Carlo simulation study. To this end, we
consider GP data $\mathcal{D}_{X}=\\{\textbf{x}_{i}\\}_{i=1}^{n}$ and
$\mathcal{D}_{Y}=\\{\textbf{y}_{j}\\}_{j=1}^{m}$ in the interval $T=[0,\pi]$
recorded over a discrete grid $\mathcal{T}=\\{0/p,\dots,(p-1)/p\\}$. Following
the examples in Figure 1, we set 3 different scenarios.

#### Scenario A: Local differences in mean

Gaussian Process data are generated according to the following specification

$X(t)=(\boldsymbol{\beta}_{X}+\boldsymbol{\varepsilon})^{T}\boldsymbol{\Phi}(t),\quad\text{
and }\quad
Y(t)=(\boldsymbol{\beta}_{Y}+\boldsymbol{\varepsilon})^{T}\boldsymbol{\Phi}(t),$

where $\boldsymbol{\beta}_{X}=\\{1,-2,-1,1,2,-1,2,3,-0.5\\}$,
$\boldsymbol{\beta}_{Y}=\\{-1,-2,-1,1,2,-1,2,5,-0.5\\}$,
$\boldsymbol{\Phi}(t)\equiv\\{\phi_{1}(t),\phi_{2}(t),\dots,\phi_{9}(t)\\}$ is
a vector function containing the first 9 Fourier basis functions and
$\boldsymbol{\varepsilon}$ is a normally distributed random vector
$\boldsymbol{\varepsilon}\sim\mathcal{N}_{9}(\textbf{0},0.25\mathbf{I}_{9})$.
For the seek of simplicity, the distribution of $\boldsymbol{\varepsilon}$
remain fixed in all scenarios. In this simulation scenario it holds that
$\sigma_{X}(t,s)=\sigma_{Y}(t,s)$ for all $(t,s)\in T\times T$, nonetheless
$\mu_{X}(t)$ is remarkably different to $\mu_{Y}(t)$ around grid point $50$ as
can be seen in Figure 2–(a).

#### Scenario B: Local differences in variance

Data are generated under the follow specification

$X(t)=\big{(}\boldsymbol{\beta}+\boldsymbol{\varepsilon}+\boldsymbol{\gamma}e^{-(t-3\pi/4)^{2}}\big{)}^{T}\boldsymbol{\Phi}(t),\quad\text{
and }\quad
Y(t)=(\boldsymbol{\beta}+\boldsymbol{\varepsilon})^{T}\boldsymbol{\Phi}(t),$

where $\boldsymbol{\beta}=\\{1,-2,-1,1,2,-1,2,3,-0.5\\}$, and
$\boldsymbol{\gamma}\sim\mathcal{N}(\boldsymbol{0},\tau^{2}\textbf{I}_{\mathbf{9}})$
is a multivariate normal random vector with $\tau^{2}=1$. Some comments on
this scenario are in order: It holds that $\mu_{X}(t)=\mu_{Y}(t)$ for all
$t\in T$, nevertheless $\sigma_{X}(t,s)$ present more differences with
$\sigma_{Y}(t,s)$ around the grid point 75 (which corresponds to time point
$t=3\pi/4$) as can be seen in Figure 2–(b).

#### Scenario C: Local differences in mean and variance

For this scenario we consider

$X(t)=(\boldsymbol{\beta}_{X}+\boldsymbol{\varepsilon}+\boldsymbol{\gamma}e^{-(t-3\pi/4)^{2}}\big{)}^{T}\boldsymbol{\Phi}(t),\quad\text{
and }\quad
Y(t)=(\boldsymbol{\beta}_{Y}+\boldsymbol{\varepsilon})^{T}\boldsymbol{\Phi}(t).$

In this scenario $\boldsymbol{\beta}_{X}$ and $\boldsymbol{\beta}_{Y}$ take
the same values as in scenario A, and $\boldsymbol{\gamma}$ is defined
likewise in Scenario B, therefore it holds that mean and variance functions
differ over different subsets in $T$ as can be seen in Figure 2–(c).

In Figure 2 panels (a) to (c), we depict one shot GP data examples ($n=m=25$
and $p=100$) drawn from scenario A to C respectively; while in panels (d) to
(f) we show the values of the estimated $\text{KL}_{\mathcal{A}(c)}(X||Y)$ as
a function of the central point of each interval corresponding to
$\mathcal{A}(c)$ for $c=\\{0.10,0.20\\}$ (i.e.
$\text{KL}_{\mathcal{A}(c)}(X||Y)=\text{KL}_{t_{c}}(X||Y)$ where $t_{c}$ is
the centre of the ball
$B(t_{c},r_{c})=[\min_{t\in\mathcal{A}(c)}t,\max_{t\in\mathcal{A}(c)}t]$).
Notice that the maximum value of the estimated local divergence in the lower
panels in Figures 2 (d)–(f), corresponds to an estimated interval of local
maximum divergence depicted in the upper panel of Figures 2 (a)–(c). In all
cases of this one shot experiment, the estimated interval of local maximum
divergence is close to the true interval of local maximum divergence. To
validate the accuracy of the estimation method, we run a Monte Carlo
simulation study.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 2: One shot experiment: Panels (a)–(c) illustrate the realisations of
$X(t)$ and $Y(t)$ for scenarios A to C and the corresponding estimated mean
functions (depicted in solid red and blue respectively). Panels (d)–(f) show
the estimated local KL divergence for different interval lengths: $c=0.1$ ( )
and $c=0.2$ ( ); the multiple local maxima in panel (f) corresponds to the
particular casuistry of scenario C (for a sufficiently large $c$, the local
divergence curve will be unimodal as well).

#### Monte-Carlo Results

The Monte Carlo simulation study consider, in each scenario, $M=1000$ data
replicates for sample sizes $m=n\in\\{50,100,250,500,1000\\}$ and grid
resolution levels $p\in\\{50,100,200,500\\}$. To assess the estimation
performance of Algorithm 1, we depict in Figure 3 the empirical distribution
of the Average Integrated Jaccard Distance (AIJD) defined as:

$\text{AIJD}=E\Bigg{\\{}\int_{0}^{1}\bigg{[}1-\frac{|\mathcal{A}^{*}(c)\cap\widehat{\mathcal{A}}^{*}(c)|}{|\mathcal{A}^{*}(c)\cup\widehat{\mathcal{A}}^{*}(c)|}\bigg{]}\mathtt{d}c.\Bigg{\\}}$
(7)

The Jaccard Index Jaccard (1912) is the natural measure to assess similarity
between sets, and we estimate AIJD using the trapezoidal rule over a uniform
grid for $c$. The numerical analysis of our estimation method can be seen in
Figure 3. Some comments about the Monte Carlo results are in order. As sample
size $n$ and $m$ increases, then the AIJD decreases in all scenarios, this
suggest a consistent estimation method. For Scenario A and C, the numerical
experiment show highly accurate results even for low sample sizes (see for
instance $n=m=50$), while in Scenario B, the results present more variability.
Notice also that for fixed $n$ and $m$, the the estimated AIJD increases in
average as $p$ increases (in almost all scenarios). This phenomena is rather
natural, since as data dimension $p$ increases, the troublesome estimation of
$p\times p$ covariance matrices lead to low quality local KL estimations.

(a)

(b)

(c)

Figure 3: Empirical distribution of AIJD for different sample sizes $n$, $m$
and grid resolution levels $p$. Scenarios A–to–C in panels (a)–to–(c)
respectively.

#### Computational complexity

We also study the numerical efficiency of the proposed method using the Monte
Carlo simulations study222The experiment was executed in a computer with 36
cores and 1.5TB per node, paralleling per Monte Carlo iteration, sample size,
grid resolution and Scenario.. In Figure 4–(a) we depict the average time (in
seconds) required to execute Algorithm 1 for different sample sizes and fix
$p=500$ (the largest value in the Monte Carlo simulation) under the 3
scenarios, while and in panel (b) we fix the sample size at its largest value
$n=m=1000$ and consider different grid resolutions levels $p$. As can be seen
in both panels, small values of $c$ (short intervals) involves more
computation work since Algorithm 1 proceed in an exhaustive search for the
maxima. However, even in the case of $c=0.1$, $m=n=1000$ and $p=500$ (the
shortest interval in the experiments and the largest data sets), it takes no
more than 20 seconds (on average) to estimate the interval of local maximum
divergence. In addition, the computational time remains constant on average as
the sample size increases for fixed $p$ and $c$333Notice that the time
required to compute means and variances does not increase in a sensitive way
as $n$ and $m$ increases. Since mean and variances are estimated only once in
Algorithm 1, then increasing $n$ and $m$ has approximately zero impact in the
computational time required to estimate local divergences and the
corresponding interval of local maximum divergence.; meanwhile for fixed $n$,
$m$ and $c$, the computational time increases exponentially in $p$. Redefine
Step 2 in Algorithm 1 in order to search for a local maxima in a more
efficient way will be part of the future research directions in order to
tackle domain selection problems in the context of extremely large values of
$p$.

(a)

(b)

Figure 4: Average computational time required to estimate intervals of local
maximum divergence for different sample sizes in (a) and grid resolution
levels in (b).

##  4 Monitoring electrocardiogram signals

The ECG signal is the visual representation of the heart electrical activity
as a function of time. Learning which part of the signal spectrum is more
relevant to diagnose a cardiac disease, is of fundamental importance in order
to increase the probability of survival during a cardiac episode. In this
section we illustrate the relevance of the proposed method in the analysis of
ECG data.

#### Data and goals

The ECG data set Olszewski (2001) is available in the UEA & UCR repository. It
consists of 200 signals sampled over a grid of 96 equally spaced instances444A
healthy heart runs on a typical rate of 70 to 75 beats per minute and take
about 0.8 second to complete the cardiac cycle.. Each observation represent
the cardiac electrical activity recorded during one heartbeat and there are
two groups of signals: 133 normal heartbeats and 67 myocardial infarction
signals. The inferential and predictive tasks relevant in this section are: i)
Learn about $\mathcal{A}^{*}(c)$ from data and quantify the uncertainty around
the estimation of such interval–we use the bootstrap strategy described in
Subsection 2.3. ii) Although learning about $\mathcal{A}^{*}(c)$ is in
principle unrelated to classification, the ECG data is a popular benchmark for
new classifiers; it may be sensible to ask whether the accuracy of
Discriminant Analysis555DA is the optimal classifier for Gaussian data;
nevertheless other classifiers can be considered as well. The analysis of
other classification methods reach out of the scope of this paper. can be
improved by focusing on $\widehat{\mathcal{A}}^{*}(c)$ rather than treating
the entire time horizon $\mathcal{T}$ equally.

#### Implementation and results

To learn about $\mathcal{A}^{*}(c)$, we consider $c\in\\{0.1,0.2,0.25\\}$ that
corresponds to intervals of 10, 19 and 24 deciseconds (dcs) respectively. In
Figure 5-(a) we display the ECG raw signals along with the corresponding
estimates $\widehat{\mathcal{A}}^{*}(c)$ using brackets. The selected domains
corresponds to grid points between 20 and 55. All in all, for small values of
$c$, the analysis suggests that while normal heartbeats and myocardial
infarction signals have similar ‘peaks’ at the beginning of the sample period
(i.e. they have similar Q waves, in ECG signal analysis terminology),
immediately right after that period (i.e. over their ST segments) they greatly
differ. To quantify the variability of our estimates, we replicate the
estimation of $\widehat{\mathcal{A}}^{*}(c)$ using $B=1000$ bootstrap samples.
Figure 5-(b) display the empirical density of the center of
$\widehat{\mathcal{A}}^{*}(c)$ (we denote this distribution as
$\widehat{F}_{B}(t^{*})$ in § 2.3) obtained via the bootstrap samples and
suggest low variability in our estimator $\widehat{\mathcal{A}}^{*}(c)$ for
$c\in\\{0.1,0.2,0.25\\}$.

To assess the classification performance of Discriminant Analysis over
different slices on the domain, we consider
$c\in\\{0.1,0.2,0.25,0.3,0.4,0.5,0.6,0.7,0.75,0.8,0.9,1.0\\}$ which
corresponds to intervals of $10,19,24,29,38,48,58,67,72,77,86$, and $96$
deciseconds respectively; and randomly split the data into training–testing
samples in a 50%-50% fashion. For each value of $c$, we learn
$\widehat{\mathcal{A}}^{*}(c)$ and the corresponding parameters of the
discriminant function with train data; while test data is used to estimate the
missclassification error rate $\widehat{\text{err}}(c)$. To assess the
estimation variability on $\widehat{\mathcal{A}}^{*}(c)$ and
$\widehat{\text{err}}(c)$, we consider $B=1000$ bootstrap samples with train
and test data respectively. In Figure 6-(a) we show the estimated centres of
$\widehat{\mathcal{A}}^{*}(c)$ (black dots in the vertical box–plots
corresponds to $\hat{t}^{*}_{c}$) and the corresponding 95% bootstrap
confidence interval for such centres as a function of $c$ (in a deciseconds
scale). In Figure 6-(b) we display the estimated $\widehat{\text{err}}(c)$
(black dots in the vertical box–plots) and its corresponding 95% bootstrap
confidence interval as a function of $c$ (in a deciseconds scale). As can be
seen from Figure 5-(b), the discrimination power of DA is significantly larger
if we consider a small interval of local maximum divergence on the ECG
data–say $c\in\\{0.1,0.2,0.25,0.3\\}$– rather than considering the full domain
corresponding to ECG data.

(a)

(b)

Figure 5: (a) ECG signals (lighted coloured curves) alongside the estimated
mean functions –solid curves corresponding to healthy (—) and myocardial
infarction (—) signals. Selected domain for interval lengths $10,19$ and $24$
dcs. displayed with black "[ - ]", gray "[ - ]" and light–gray "[ - ]"
brackets respectively. (b) Bootstrap densities for the interval centre (the
median of each empirical distribution is reported on the horizontal axis)
corresponding to intervals lengths $10,19$ and $24$ dcs.

(a)

(b)

Figure 6: (a) Bootstrap interval centre and 95% bootstrap CI throughout all
the interval lengths considered. (b) Classification error using a Discriminant
model for different selected domains.

##  5 Discussion and closing remarks

From a methodological outlook, the main goal of this article entails the
combination of Kullback–Leibler divergence and Gaussian Processes to develop a
super-fast and easy to implement algorithm for domain selection. We define a
local KL divergence, introduce its fundamental properties, and prove the
existence of an interval of local maximum divergence $\mathcal{A}^{*}(c)$
under suitable conditions on the mean and variance functions. In addition, we
also propose an estimator for $\mathcal{A}^{*}(c)$, devise a nonparametric
approach to assess the estimation uncertainty, and also discuss relevant
variants and extensions.

Through a Monte Carlo simulation study we numerically assess the consistency
of our estimator and demonstrate that even for large $n$ and $p$ (i.e.
$n=m=1000$ and $p=500$, which corresponds to high–dimensional and large data
sets), Algorithm 1 is very efficient. Learning about intervals of local
maximum divergence in the context electrocardiogram data contributes to
improve diagnostic tools, as we demonstrate in the analysis of ECG data. In
addition, we also explore how the discrimination power of DA (healthy vs
disease heartbeats) can be improved by making emphasis on a small interval of
ECG data rather than using the full domain.

Although the proposed method can be used for domain selection with GP, there
are natural opportunities for further analysis: (i) In this paper consistency
is illustrated whit a Monte Carlo simulation. Further theoretical developments
will be conducted in order to establish general conditions in the model in
order to ensure that $\widehat{\mathcal{A}}^{*}(c)$ is a consistent estimator
for $\mathcal{A}^{*}(c)$. (ii) To define local KL divergences for GP we rely
on a discretized version of $X(t)$ and $Y(t)$, therefore a natural point to
address is to extend Eqs. (2) and (3) to the non-discretized case; studying
also the impact on the numerical complexity of adapting Algorithm 1 to such
context. (iii) Another important avenue for future research is the extension
of the proposed domain selection method to address non GP data, perhaps using
alternative metrics to assess differences between non Gaussian processes. (iv)
In the framework of curve classification, is important to balance predictive
classification accuracy against the number of covariates which, in the case of
curves, corresponds to the length of the interval. Therefore a natural follow-
up within the remit of this paper, is the analysis of domain selection as a
shrinkage method for curve classification. (v) We also propose to study
alternative ways to introduce Bayesian tools for domain selection with GP as
we discuss in the paragraph Sampling Designs in § 2.2. (vi) Finally, from a
computational view point, Step 2 in Algorithm 1 needs to be reformulated for
extremely high–dimensional data since its computational complexity grows
exponentially in $p$ as we mention in Section 2 and illustrate numerically in
Section 3.

## Declarations

Conflicts of interest: None of the authors has a conflict of interest.

Ethics approval: Authors have no affiliations with or involvement in any
organisation or entity with any financial interest or non-financial interest
in the subject matter or materials discussed in this manuscript.

Consent for publication: Authors give consent for publication.

Availability of data and code: Data is available on UEA & UCR and source code
to reproduce the results is included as a supplementary file in the
submission.

Authors’ contributions: Authors contributed equally to this work.

###### Proof of Proposition 2.1.

Property (a) follows directly from the properties of KL divergence. To prove
the upper bound in (b), let $\boldsymbol{\Sigma}_{\mathcal{T},X}$ and
$\boldsymbol{\Sigma}_{\mathcal{T},Y}$ be positive definite (PD) matrices, then
$\boldsymbol{\Gamma}_{\mathcal{T}}\equiv\boldsymbol{\Sigma}_{\mathcal{T},Y}^{-1}\boldsymbol{\Sigma}_{\mathcal{T},X}$
is also PD with eigenvalues
$\infty>\gamma_{\mathcal{T},1}\geq\dots\geq\gamma_{\mathcal{T},p}>0$ where
$p=|\mathcal{T}|$; then if $\|\Delta_{\mathcal{T}}\|_{2}<\infty$ it holds from
Eq. (2) that:

$2\text{KL}_{\mathcal{T}}(X||Y)=\sum_{i=1}^{p}\big{\\{}\gamma_{\mathcal{T},i}-\ln(\gamma_{\mathcal{T},i})\big{\\}}-p+\Delta_{\mathcal{T}}^{\mathsf{T}}\boldsymbol{\Sigma}_{\mathcal{T},Y}^{-1}\Delta_{\mathcal{T}}<\infty,$

since $1\leq x-\ln(x)<\infty$, for all $0<x<\infty$. Moreover,
$\text{KL}_{\mathcal{A}}\equiv\text{KL}_{\mathcal{A}}(X||Y)\leq\text{KL}_{\mathcal{A}^{\prime}}$
for all $\mathcal{A}^{\prime}\subseteq\mathcal{A}\subseteq\mathcal{T}$ with
$|\mathcal{A}^{\prime}|\equiv d^{\prime}\leq|\mathcal{A}|\equiv d$ if the
eigenvalues of $\boldsymbol{\Gamma}_{\mathcal{A}^{\prime}}$ are bigger than 1.
To prove the assertion, consider:

$\begin{split}2(\text{KL}_{\mathcal{A}}-\text{KL}_{\mathcal{A}^{\prime}})=\underbrace{\sum_{i=1}^{d^{\prime}}\big{\\{}\gamma_{\mathcal{A},i}-\gamma_{\mathcal{A}^{\prime},i}-\big{(}\ln(\gamma_{\mathcal{A},i})-\ln(\gamma_{\mathcal{A}^{\prime},i})\big{)}\big{\\}}}_{A}+\\\
\underbrace{\sum_{j=d^{\prime}+1}^{d}\big{\\{}\gamma_{\mathcal{A},j}-\ln(\gamma_{\mathcal{A},j})\big{\\}}-(d-d^{\prime})}_{B}+\underbrace{\big{(}\Delta_{\mathcal{A}}^{\mathsf{T}}\boldsymbol{\Sigma}_{\mathcal{A},Y}^{-1}\Delta_{\mathcal{A}}-\Delta_{\mathcal{A}^{\prime}}^{\mathsf{T}}\boldsymbol{\Sigma}_{\mathcal{A}^{\prime},Y}^{-1}\Delta_{\mathcal{A}^{\prime}}\big{)}}_{C}\end{split}$

then $C\geq 0$ since
$\Delta_{\mathcal{A}^{\prime}}^{\mathsf{T}}\boldsymbol{\Sigma}_{\mathcal{A}^{\prime},Y}^{-1}\Delta_{\mathcal{A}^{\prime}}\leq\Delta_{\mathcal{A}}^{\mathsf{T}}\boldsymbol{\Sigma}_{\mathcal{A},Y}^{-1}\Delta_{\mathcal{A}}$
for all $\mathcal{A}^{\prime}\subseteq\mathcal{A}$; $B\geq 0$ since
$x-\ln(x)\geq 1$ if $0<x<\infty$; and $A\geq 0$ by the eigenvalues interlacing
inequality (i.e.
$\gamma_{\mathcal{A},1}\geq\gamma_{\mathcal{A}^{\prime},1}\geq\gamma_{\mathcal{A},2}\geq\gamma_{\mathcal{A}^{\prime},2}\geq\dots\gamma_{\mathcal{A},d^{\prime}}\geq\gamma_{\mathcal{A^{\prime}},d^{\prime}}>0$),
and Napier’s inequality (i.e. $\ln(x)-\ln(y)\leq(x-y)/y$ for $x\geq y>0$):

$A=\sum_{i=1}^{d^{\prime}}\big{\\{}\gamma_{\mathcal{A},i}-\gamma_{\mathcal{A}^{\prime},i}-\big{(}\ln(\gamma_{\mathcal{A},i})-\ln(\gamma_{\mathcal{A}^{\prime},i})\big{)}\big{\\}}\geq\sum_{i=1}^{d^{\prime}}\frac{(\gamma_{\mathcal{A},i}-\gamma_{\mathcal{A}^{\prime},i})(\gamma_{\mathcal{A}^{\prime},i}-1)}{\gamma_{\mathcal{A}^{\prime},i}}\geq
0,$

if
$\infty>\gamma_{\mathcal{A}^{\prime},1}\geq\dots\geq\gamma_{\mathcal{A}^{\prime},d^{\prime}}\geq
1$. Notice that $\mathcal{A}^{*}(c)$, the solution of the set function
optimisation problem stated in Eq. (5), exist under the previous stated
conditions on the mean and variance functions, since
$\mathcal{P}_{\mathcal{T}}$ is a finite collection of sets and
$\text{KL}_{\mathcal{A}}(X||Y)$ is a non-decreasing and bounded set function.
To prove property (c), consider with out loss of generality the following
sequence of non–decreasing sets
$\mathcal{A}^{(k)}\equiv\\{t_{1},t_{2},\dots,t_{k}\\}$ and notice that
$\cup_{k=1}^{\infty}\mathcal{A}^{(k)}=\cup_{k=1}^{p}\mathcal{A}^{(k)}=\mathcal{T}$,
then the left–continuity follows from:

$\lim_{k\to\infty}2\text{KL}_{\mathcal{A}^{(k)}}=\lim_{k\to
p}2\text{KL}_{\mathcal{A}^{(k)}}=\lim_{k\to
p}\sum_{i=1}^{k}\big{\\{}\gamma_{\mathcal{A}^{(k)},i}-\ln(\gamma_{\mathcal{A}^{(k)},i})\big{\\}}-k+\Delta_{\mathcal{A}^{(k)}}^{\mathsf{T}}\boldsymbol{\Sigma}_{\mathcal{A}^{(k)},Y}^{-1}\Delta_{\mathcal{A}^{(k)}}=2\text{KL}_{\mathcal{T}},$

since $\gamma_{\mathcal{A}^{(k)},i}=\gamma_{\mathcal{T},i}$ for all
$i\in\\{1,\dots,p\\}$ and
$\Delta_{\mathcal{A}^{(k)}}^{\mathsf{T}}\boldsymbol{\Sigma}_{\mathcal{A}^{(k)},Y}^{-1}\Delta_{\mathcal{A}^{(k)}}=\Delta_{\mathcal{T}}^{\mathsf{T}}\boldsymbol{\Sigma}_{\mathcal{T},Y}^{-1}\Delta_{\mathcal{T}}$
if $k=p$. To prove the right–continuity consider, without loss of generality,
the sequence of non–increasing sets
$\mathcal{A}_{(k)}\equiv\\{t_{k},t_{k+1},\dots,t_{p}\\}$ and then
$\cap_{k=1}^{\infty}\mathcal{A}_{(k)}=\cap_{k=1}^{p}\mathcal{A}_{(k)}=\\{t_{p}\\}$.
The rest of the proof is tantamount to the previous case. ∎

## References

  * Mullainathan and Obermeyer [2022] Sendhil Mullainathan and Ziad Obermeyer.  Solving medicine’s data bottleneck: Nightingale open science.  _Nature Medicine_ , 28(5):897–899, 2022.
  * Wang et al. [2013] Jin Wang, Mary She, Saeid Nahavandi, and Abbas Kouzani.  Human identification from ecg signals via sparse representation of local segments.  _IEEE Signal Processing Letters_ , 20(10):937–940, 2013.
  * Rodríguez et al. [2015] Ricardo Rodríguez, Adriana Mexicano, Jiri Bila, Salvador Cervantes, and Rafael Ponce.  Feature extraction of electrocardiogram signals by applying adaptive threshold and principal component analysis.  _Journal of applied research and technology_ , 13(2):261–269, 2015.
  * Pérez-Cruz et al. [2013] Fernando Pérez-Cruz, Steven Van Vaerenbergh, Juan José Murillo-Fuentes, Miguel Lázaro-Gredilla, and Ignacio Santamaria.  Gaussian processes for nonlinear signal processing: An overview of recent advances.  _IEEE Signal Processing Magazine_ , 30(4):40–50, 2013.
  * Kullback and Leibler [1951] Solomon Kullback and Richard A Leibler.  On information and sufficiency.  _The annals of mathematical statistics_ , 22(1):79–86, 1951.
  * Berrendero et al. [2016] José R Berrendero, Antonio Cuevas, and José L Torrecilla.  Variable selection in functional data classification: a maxima-hunting proposal.  _Statistica Sinica_ , pages 619–638, 2016.
  * Baíllo et al. [2011] Amparo Baíllo, Antonio Cuevas, and Ricardo Fraiman.  Classification methods for functional data.  _The Oxford handbook of functional data analysis_ , 2011.
  * Pini and Vantini [2017] Alessia Pini and Simone Vantini.  Interval-wise testing for functional data.  _Journal of Nonparametric Statistics_ , 29(2):407–424, 2017.
  * Martos and de Carvalho [2018] G. Martos and M. de Carvalho.  Discrimination surfaces with application to region-specific brain asymmetry analysis.  _Statistics in Medicine_ , 11(37):1859–1873, 2018\.
  * Pardo [2018] Leandro Pardo.  _Statistical inference based on divergence measures_.  CRC press, 2018.
  * McLachlan [1999] Goeffrey J McLachlan.  Mahalanobis distance.  _Resonance_ , 4(6):20–26, 1999.
  * Székely et al. [2007] Gábor J Székely, Maria L Rizzo, and Nail K Bakirov.  Measuring and testing dependence by correlation of distances.  2007\.
  * Nocedal and Wright [1999] Jorge Nocedal and Stephen J Wright.  _Numerical optimization_.  Springer, 1999.
  * Hastie et al. [2009] Trevor Hastie, Robert Tibshirani, Jerome H Friedman, and Jerome H Friedman.  _The elements of statistical learning: data mining, inference, and prediction_ , volume 2.  Springer, 2009.
  * Pourahmadi [2013] Mohsen Pourahmadi.  _High-dimensional covariance estimation: with high-dimensional data_ , volume 882.  John Wiley & Sons, 2013.
  * Ramsay and Li [1998] James O Ramsay and Xiaochun Li.  Curve registration.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 60(2):351–363, 1998.
  * Berndt and Clifford [1994] Donald J Berndt and James Clifford.  Using dynamic time warping to find patterns in time series.  In _KDD workshop_ , volume 10, pages 359–370. Seattle, WA, USA:, 1994\.
  * Kazlauskaite et al. [2019] Ieva Kazlauskaite, Carl Henrik Ek, and Neill Campbell.  Gaussian process latent variable alignment learning.  In _The 22nd International Conference on Artificial Intelligence and Statistics_ , pages 748–757. PMLR, 2019.
  * Ullah and Finch [2013] Shahid Ullah and Caroline F Finch.  Applications of functional data analysis: A systematic review.  _BMC medical research methodology_ , 13:1–12, 2013.
  * Efron and Tibshirani [1994] Bradley Efron and Robert J Tibshirani.  _An introduction to the bootstrap_.  CRC press, 1994.
  * Fraley and Raftery [2002] Chris Fraley and Adrian E Raftery.  Model-based clustering, discriminant analysis, and density estimation.  _Journal of the American statistical Association_ , 97(458):611–631, 2002.
  * Jaccard [1912] Paul Jaccard.  The distribution of the flora in the alpine zone. 1.  _New phytologist_ , 11(2):37–50, 1912.
  * Olszewski [2001] Robert T Olszewski.  Generalized feature extraction for structural pattern recognition in time-series data.  Technical report, CARNEGIE-MELLON UNIV PITTSBURGH PA SCHOOL OF COMPUTER SCIENCE, 2001.