Muon puzzle in inclined muon bundles
detected by NEVOD-DECOR

A.G. Bogdanov1$\star$, N.S. Barbashina1, S.S. Khokhlov1, V.V. Kindin1, R.P.
Kokoulin1,

K.G. Kompaniets1, A.Yu. Konovalova1, G. Mannocchi2, A.A. Petrukhin1, V.V.
Shutenko1,

G. Trinchero2, V.S. Vorobev1, I.I. Yashin1, E.A. Yurina1, and E.A. Zadeba1

1 National Research Nuclear University MEPhI (Moscow Engineering Physics
Institute), Russia

2 Osservatorio Astrofisico di Torino – INAF, Italy

*<EMAIL_ADDRESS>
21st International Symposium on Very High Energy Cosmic Ray Interactions
(ISVHECRI 2022)
Online, 23-27 May 2022 10.21468/SciPostPhysProc.?

## Abstract

The data of cosmic ray NEVOD-DECOR experiment on the investigation of inclined
muon bundles for a long time period (May 2012 – March 2021) are presented. The
analysis showed that the observed intensity of muon bundles at primary cosmic
ray energies of about 1 EeV and higher can be compatible with the expectation
in frame of widely used hadron interaction models only under the assumption of
an extremely heavy mass composition. This conclusion is consistent with data
of several experiments on investigations of muon content in air showers, but
contradicts the available fluorescence data on $X_{\text{max}}$ which favor a
light mass composition at these energies. In order to clarify the nature of
the “muon puzzle”, investigations of the muon bundle energy deposit in the
detector material were carried out. For the first time, experimental estimates
of the average energy of muons in the bundles of inclined air showers
initiated by primary particles with energies from 10 to 1000 PeV have been
obtained.

###### Contents

  1. 1 Introduction
  2. 2 NEVOD-DECOR setup and experimental data
  3. 3 Local muon density spectra
  4. 4 Energy characteristics of muons in the bundles
  5. 5 Conclusions

##  1 Introduction

The muon component of extensive air showers (EAS) is formed mainly as a result
of the decays of pions and kaons generated in hadronic interactions in the
atmosphere, and is often used to estimate the mass composition of primary
cosmic rays (PCR) and to verify models of hadron interactions at high
energies. Of particular interest are events with muon bundles, which represent
a simultaneous (within tens of nanoseconds) passage through the detector of
several penetrating particles with almost parallel tracks.

To date, NEVOD-DECOR is the only experiment in which long-term systematic
studies of muon bundles in a wide range of zenith angles are conducted. This
article provides a brief overview of the results of these studies obtained so
far.

##  2 NEVOD-DECOR setup and experimental data

The NEVOD-DECOR setup is located at MEPhI and includes a Cherenkov water
calorimeter NEVOD [1] with a volume of about 2000 m3 and a coordinate-tracking
detector DECOR [2] with an area of about 70 m2, united by a trigger system.

The Cherenkov water detector NEVOD consists of 91 quasi-spherical measuring
modules (QSMs) located at the nodes of the volume-centered spatial lattice. In
fact, the lattice is formed by 25 vertical strings of 3 or 4 QSMs. Each QSM
consists of 6 FEU-200 photomultipliers with a flat photocathode of 15 cm in
diameter, oriented along the axes of the orthogonal coordinate system. Such
design provides almost the same detection efficiency of Cherenkov radiation
coming from any direction. A wide dynamic range of the measurements of each
PMT (1 – 105 photoelectrons) is due to two-dynode signal readout. This allows
calorimetric studies, in particular, measurements of the energy deposit of
muon bundles.

The coordinate-tracking detector DECOR consists of 8 vertical supermodules
(SMs) located in the galleries of the experimental building on three sides of
the water volume of the NEVOD detector. The effective area of one SM is
3.1$\times$2.7 m2. Each SM represents 8 parallel planes consisting of plastic
chambers of gas-discharge streamer tubes. The planes are equipped with a
system of external strips for signal readout in two coordinates (along and
across the chambers). Accuracy of localization of tracks of charged particles
in one SM is better than 1 cm, and the angular accuracy of the reconstruction
of tracks crossing the SM is better than 1∘.

In the present work, we used experimental data on muon bundles accumulated
over a long period from May 2012 to March 2021 (“live” observation time is
equal to 58.3 thousand hours). About 99.6 thousand events with muon
multiplicity $m\geq 5$ and zenith angles $\theta\geq 55^{\circ}$ were
selected; in addition 30.4 thousand events with smaller zenith angles from 40
to 55∘ (“live” time is 6.3 thousand hours) were sampled. In order to improve
the identification of muon tracks, the events were selected in two 60∘-wide
azimuth angle sectors where 6 of the 8 DECOR SMs are shielded by the NEVOD
detector water tank. At that, the average threshold energy of muons is about 2
GeV.

The method of identifying muon bundles of atmospheric origin in a coordinate
detector is based on the parallelism of particle tracks (within 5∘ cone)
recorded by the setup. The event selection procedure consists of several
stages: hardware level (3-fold coincidence of the signals from different SMs
within the time gate of 250 ns); program reconstruction and selection; final
classification of events and counting of tracks by several operators.

An example of the spatial reconstruction of the muon bundle ($m$ = 8 and
$\theta$ = 57∘) registered in the NEVOD-DECOR setup, as well as the main parts
of the experimental complex are shown in Figure 1. Thin lines show the
reconstruction of muon tracks according to the DECOR data, the small circles
show the hit PMTs in the NEVOD calorimeter (colors reflect signal amplitudes),
the large rectangles show the SMs of the detector DECOR. Multimuon events have
a bright signature in the coordinate detector, and their interpretation is
almost unambiguous. Obviously, for each such event, it is possible to find the
number of tracks and determine the direction.

Figure 1: An example of geometrical reconstruction of muon bundle event
detected by the NEVOD-DECOR setup

During the whole time of the experiment, an express analysis of the NEVOD and
DECOR data is carried out in order to quickly identify and eliminate equipment
faults and to prepare experimental data for further physical and methodical
processing.

##  3 Local muon density spectra

A new approach to the study of the EAS, the method of local muon density
spectra (LMDS), was developed for the physical analysis of the data of the
first experiment on the study of muon bundles conducted at detector DECOR in
2002-2007 [3]. The typical size of the EAS muon component ($\sim$ km)
significantly exceeds the size of the NEVOD-DECOR setup (tens of meters), so
the detector can be considered as a point-like one. In an individual event
with a muon bundle the local muon density $D$ at the observation point is
estimated, in the first approximation: $D\sim m/S_{\text{det}}$, where $m$ is
the multiplicity of muons in the bundle (the number of muons that hit the
detector), and $S_{\text{det}}$ is the area of the DECOR detector for a given
direction. The distribution of events in the local muon density estimate $D$
forms the LMDS. The spectrum of events by local muon density has a nearly
power-law form, with a slope $\beta$ slightly steeper than the slope of the
PCR spectrum.

It is important that, at the same muon density, different zenith angles
correspond to substantially (by orders of magnitude) different characteristic
energies of primary particles contributing to the selected events. The fact is
that with an increase in the zenith angle, the spread of muons in bundles
increases, mainly due to transverse momenta during the generation and decay of
mother hadrons, deviation in the Earth’s magnetic field and multiple
scattering, since muons travel a longer path in the atmosphere. Thus,
measurements of the LMDS at different zenith angles make it possible to study
a wide energy range of PCRs by means of a relatively small setup. In this
case, the event collection area is determined not by the size of the detector,
but by the cross section of the air shower, which in the muon component near
the horizon reaches several km2, that is sufficient to reach energies of 1018
eV and even higher.

The local muon density spectra are used as a frontier between the experiment
and calculations. The main stages of data analysis on the local muon density
are described in detail in [3, 4, 5].

The experimental spectra are obtained by means of reconstructing
$dF(D,\theta)/dD$ (in a detector-independent form) from the experimental
distributions of event characteristics $N(m,\theta,\phi)$. The reconstruction
procedure takes into account geometric factors, Poisson fluctuations in the
number of muons that hit the detector, the detection efficiency of streamer
tube chambers, trigger and event selection conditions, etc.

The expected spectra $dF(D,\theta)/dD$ are obtained as a result of convolution
of the lateral distribution functions (LDFs) of muons with certain models of
the PCR spectrum. Two-dimensional LDFs are calculated on the basis of EAS
simulations using the CORSIKA program [6] for a set of fixed zenith angles and
PCR energies, also taking into account such factors as the threshold energy of
secondary particles, the altitude of the observation level, the atmosphere,
and the Earth’s magnetic field. Naturally, the shape of the LDF depends on the
chosen model of hadron interactions and assumptions about the PCR mass
composition.

Comparison of experimental and calculated LMDS for 9 intervals of zenith
angles is shown in Figure 2. The NEVOD-DECOR data are marked by symbols. The
solid, dashed, and dash-dotted curves represent the results of calculations
based on the CORSIKA (v. 7.69) program for the hadronic interaction models
tuned according to the LHC data: QGSJET-II-04, SIBYLL-2.3c, and EPOS-LHC,
respectively. Two limiting cases of the PCR mass composition were used in the
simulations: only protons (p, three lower curves for each zenith angle) and
only iron nuclei (Fe, upper curves). The absolute intensity increases for
showers initiated by heavier nuclei, because at a fixed nucleus energy
$E_{0}$, the muon density rises with increasing atomic mass. As a model of the
PCR energy spectrum, the approximation of the data of various experiments from
the PDG review [7] proposed by us in [4] was used. The arrows in the figure
indicate the typical (mean logarithmic) PCR energies, which give the main
contribution to the events sampled by the muon density. In this case, the PCR
energy distributions are relatively narrow ($\sigma_{\text{lgE}}$ $\approx$
0.4 for $D$ = 0.2, which corresponds to 5 – 7 muons that passed through the
detector DECOR). This is due to a rapid decrease in the cosmic ray flux with
increasing energy, although the common integral is given by air showers from
PCRs with different energies, registered at different (random) distances from
the EAS axis.

Figure 2: Differential local muon density spectra for different zenith angles

It can be seen from the figure that at moderate zenith angles (42 – 52∘),
which correspond to PCR energies $\sim$ 1016 eV, the experimental data are
close to the results of calculations for a light mass composition. Further, at
higher energies, there is a relative increase of the intensity of muon
bundles, which, in principle, can be interpreted as a trend toward heavier PCR
mass composition. At large zenith angles (77 – 82∘), corresponding to PCR
energies $\sim$ 1018 eV, the NEVOD-DECOR data are compatible with calculations
only under the assumption of an extremely heavy mass composition of PCRs (iron
nuclei), if we remain within the framework of existing hadron interaction
models. However, this contradicts the data on the measurements of the EAS
development maximum $X_{\text{max}}$ by the fluorescence method (mainly the
electron-photon component) [8, 9, 10], which indicate a light mass composition
of PCRs at such energies.

For the comparison of LMDS at different zenith angles, we used the
$z$-parameter, which was proposed by the Working group on Hadronic
Interactions and Shower Physics (WHISP) [11] to compare data of various
experiments on the study of EAS muon component, including NEVOD-DECOR:

$z=\left(\text{ln}\,N_{\mu}^{\text{obs}}-\text{ln}\,N_{\mu}^{\text{p
sim}}\right)\Big{/}\left(\text{ln}\,N_{\mu}^{\text{Fe
sim}}-\text{ln}\,N_{\mu}^{\text{p sim}}\right),$ (1)

where $N_{\mu}^{\text{obs}}$ is the observed value (muon density, number of
muons, intensity of muon bundles, etc.), and $N_{\mu}^{\text{p sim}}$ and
$N_{\mu}^{\text{Fe sim}}$ are the calculated estimates of this value for EASs
formed by primary protons and iron nuclei, then $z$ = 0 means that PCRs
consist only of protons, and $z$ = 1 corresponds to iron nuclei.

The dependences of the $z$-parameter on the PCR energy for 9 intervals of
zenith angles and three modern hadron interaction models: QGSJET-II-04,
SIBYLL-2.3c, and EPOS-LHC are shown in Figure 3,a-c. In this case, the same
model of the PCR energy spectrum was used to construct the LMDS [4]. It
follows from the figures that the values of the $z$-parameter obtained in
different ranges of zenith angles from 40 to 90∘ overlap and are in good
agreement with each other within the errors. At the same time, steady trend of
growth of $z$-parameter at the energies above 1017 eV is observed. Comparison
of the figures 3,a-c shows the influence of the choice of models of hadronic
interactions on this dependence. Thus, for the EPOS-LHC model in comparison
with the QGSJET-II-04 and SIBYLL-2.3c models, the $z$-parameter is shifted
towards a heavier mass composition over the entire PCR energy range.

Figure 3: Comparison of LMDS at different zenith angles in terms of $z$-scale
(a,b,c – effect of the hadron interaction model; c,d – influence of the
assumed model of primary cosmic ray energy spectrum)

As mentioned above, in order to obtain the calculated LMDS, it is necessary to
make an assumption about the shape of the PCR energy spectrum. In this regard,
it is interesting to see how different versions of the spectrum affect the
dependence of the $z$-parameter on the PCR energy. For comparison with our
approximation from [4], we chose the GSF (Global Spline Fit) spectrum model
proposed in [12], which is widely used by the WHISP group. It is a
parametrization of data from various experiments, including direct
observations, on the PCR flux and mass composition.

A comparison of the dependences of the $z$-parameter on the PCR energy for
these two variants of the spectrum and 9 intervals of zenith angles is shown
in Figure 3c,d, in both cases the same model of hadronic interactions EPOS-LHC
was used. Obviously, the choice of the PCR energy spectrum model affects both
the absolute values of the $z$-parameter and its energy dependence (the slope
changes). Both PCR spectra are shown in Figure 4. In general, their difference
is not too large, but still quite noticeable (about 20%) in the energy range
from 10 to 100 PeV.

Figure 4: PCR energy spectrum models used in [4] (NEVOD-DECOR approximation
based on PDG data) and [9] (GSF)

Generally, the choice of the hadron interaction model and the model of the PCR
energy spectrum will affect the estimate of the mass composition of PCRs when
interpreting the experimental data. However, the main conclusion that the
NEVOD-DECOR data on the intensity of muon bundles can be consistent with
calculations only assuming a heavy PCR mass composition in the energy region
$\sim$ 1018 eV remains.

##  4 Energy characteristics of muons in the bundles

A drastic contradiction between the experimental data on measuring the
density/number of muons in air showers and the $X_{\text{max}}$ by the
fluorescence method at ultra-high energies gives reason to believe that a new
physical features may emerge, which is not taken into account in modern hadron
interaction models. For a deeper understanding of the problem, it is necessary
to study the energy characteristics of the EAS muon component, in particular,
the muon bundles. Thus, the appearance of an excess of ultra-high energy muons
should unambiguously indicate the inclusion of new physical processes (the
formation of a new state of matter) [13].

One of the possible approaches to studying the energy characteristics of the
EAS muon component is to measure the energy deposit of muon bundles in the
detector material (in our case, in the Cherenkov water calorimeter, see e.g.
[14, 15, 16]), since muon energy loss in the matter almost linearly depends on
their energy: $dE/dX\sim a+bE$.

As a measure of the energy deposit of muon bundle in the detector NEVOD, the
sum of the all PMTs signals $\Sigma$ (in units of photoelectrons) is used. The
assumption that the total yield of Cherenkov light is proportional to the
energy losses of muons, including secondary particles and cascades from them,
is confirmed by the results of simulating the detector response.

Since, in the first approximation, the total energy deposit $\Sigma$ is
proportional to the muon density $D$ in the event (according to the DECOR
data), it is more reasonable to consider the specific energy deposit
$\Sigma/D$, and we have used this value some time [14, 15].

However, the specific energy deposit depends on the characteristics of the
detector. In order to move to traditionally used physical quantities, for
example, muon energy losses, it is necessary to find their relationship with
the yield of Cherenkov light (signals of the PMTs). For this purpose, a
mathematical model of the NEVOD-DECOR setup was developed based on the Geant4
simulation toolkit (v. 10.7) [17]. The calibration of the model was carried
out by comparing the results of calculating the response of the NEVOD detector
to the passage of near-horizontal muons with experimental data [16]. Near-
horizontal muons were sampled using SMs of the detector DECOR located in
opposite short galleries of the experimental building. The average energy of
such muons is estimated to be about 100 GeV. So, the difference between the
average values of the total number of photoelectrons of all PMTs $\Sigma$ in
the experiment and simulation is less than 1.5%, at that the distribution
shapes are in good agreement with each other, too. An analysis of the average
response of PMTs and QSMs depending on the distance to the muon track,
distributions in the number of hit (fired) PMTs and QSMs and other
phenomenological characteristics confirmed that the model well (within a few
percent) describes the experimental data.

Next, the response of the NEVOD-DECOR setup to artificial muon bundles with a
fixed energy of 100 GeV was simulated using a Monte Carlo technique [16]. Both
the physical features of the detectors and the conditions for selecting events
with muon bundles used in the experiment were taken into account. The
distributions of simulated events in muon multiplicity in the bundles (the
number of muons that hit the detector), zenith and azimuth angles are almost
identical to the experimental ones obtained from the DECOR data.

Special attention was paid to identifying possible systematic distortions that
could affect the measurement results. The influence of such factors as
residual contribution of the electron-photon and hadron components of the EAS
to the response of the detector NEVOD (about 5% at zenith angles of 55 – 60∘
and less than 1% at $\theta=$ 60 – 65∘), under-estimation of the response due
to the digitization threshold, some changes of the registration conditions in
different measurement series, efficiency of the DECOR response, masking of
tracks, etc. were taken into account.

The desired experimental dependences of the average muon energy in the bundles
on the zenith angle and on the local muon density were obtained (see [16] in
detail) from the ratio of the measured and calculated specific energy deposits
in the NEVOD-DECOR setup. At that, the dependence of the average losses in
water on the muon energy normalized to losses at 100 GeV [18] which is
practically linear in the region of hundreds of GeV was used.

The expected dependences of the average muon energy in the bundles on the
zenith angle and on the local muon density were obtained based on the results
of simulations of the EAS muon component (using the CORSIKA program) in a wide
range of zenith angles and PCR energies for two extreme (protons and iron
nuclei) assumptions about the mass composition and three models of high-energy
hadron interactions: QGSJET-II-04, SIBYLL-2.3c and EPOS-LHC. The features of
the LMDS method [4] were taken into account: the contribution of the area
element $dS$ of the air shower transverse cross section (in the vicinity of
the point r) to the flux of detected events is proportional to
$\rho(\textbf{r})^{\beta}dS$, where $\rho(\textbf{r})$ is the muon density,
$\beta$ is the LMDS slope index. In this case, $\beta=2.2$ was assumed to be
constant; in future, it is supposed to take into account its weak change with
increasing PCR energy.

Figure 5 shows for the first time the dependence of the average muon energy in
the bundles on the zenith angle. The experimental data of NEVOD-DECOR setup
are denoted by symbols; the expected dependences are shown by solid, dashed,
and dash-dotted curves for three models of high-energy hadron interactions:
QGSJET-II-04, SIBYLL-2.3c, and EPOS-LHC, respectively. The lower curves
correspond to calculations for the pure proton PCR mass composition, and the
upper ones for iron nuclei. The arrows indicate the typical PCR energies. As
can be seen from the figure, there is a growth of the average muon energy in
the bundles with an increase in the zenith angle, as well as good agreement
between the experiment and theoretical dependences.

Figure 6 shows the dependence of the average muon energy in the bundles on the
local muon density (in fact, on the PCR energy) obtained for the zenith angle
interval 65 – 75∘. The symbols show the experimental data; the notations of
curves and arrows are the same as in Figure 5. It follows from the figure that
the average muon energies in the bundles for the QGSJET-II-04 and SIBYLL-2.3c
models only slightly differ from each other, but for the EPOS-LHC model they
are approximately 10% lower. Experimental data indicate an increase in the
average muon energy in the bundles at high muon densities corresponding to PCR
energies above 1017 eV, while the expected muon energies show a tendency to a
gradual decrease. The excess of the experiment over the calculations is
3.8$\sigma$ for primary protons and 2.6$\sigma$ for iron nuclei (estimates are
given for the QGSJET-II-04 model and the last four experimental points). The
deviation of the experimental points from the calculated curves may indicate
the inclusion of a new mechanism for the generation of high-energy muons (in
the final state) at ultra high energies of PCR. Accumulation of the
experimental data and analysis of possible systematic errors are being
continued.

Figure 5: Dependences of the average energy of muons in the bundles on the
zenith angle

Figure 6: Dependences of the average energy of muons in the bundles on the
local muon density

##  5 Conclusions

The present data on the intensity of the muon bundles observed in the NEVOD-
DECOR experiment are only compatible with calculations if we assume extremely
heavy (iron group nuclei) mass composition at the primary energies around
$10^{18}$ eV. However, such an assumption contradicts the available
fluorescence data on $X_{\text{max}}$ which favor a light (proton, He) mass
composition at these energies. It seems that the solution of the “muon puzzle”
requires the introduction of serious changes in existing hadron interaction
models, which could change the balance between muon and electron-photon EAS
components.

For the first time, the dependences of the average energy of muons in the
bundles on the zenith angle and local density in the region corresponding to
the energies of primary particles 10 – 1000 PeV have been obtained. An
increase in the average energy of muons in the bundles is observed in
comparison with the expected one for primary energies above $10^{17}$ eV. This
may indicate the inclusion of a new mechanism for the generation of high-
energy muons.

## Acknowledgements

The work was performed at the Unique Scientific Facility “Experimental complex
NEVOD” with the financial support provided by the Russian Ministry of Science
and Higher Education, project “Fundamental problems of cosmic rays and dark
matter”, No. 0723-2020-0040. Simulations were carried out using the resources
of the MEPhI high-performance computing center.

## References

  * [1] A.A. Petrukhin, Cherenkov water detector NEVOD, Phys. Usp. 58 486 (2015), 10.3367/UFNe.0185.201505i.0521.
  * [2] N.S. Barbashina et al., A Coordinate detector for studying horizontal fluxes of cosmic rays, Instrum. Experim. Techniques 43 743 (2000), 10.1023/A:1026603430476.
  * [3] A.G. Bogdanov et al., Investigation of the properties of the flux and interaction of ultrahigh-energy cosmic rays by the method of local-muon-density spectra, Phys. Atom. Nucl. 73 1852 (2010), 10.1134/S1063778810110074.
  * [4] A.G. Bogdanov et al., Investigation of very high energy cosmic rays by means of inclined muon bundles, Astropart. Phys. 98 13 (2018), 10.1016/j.astropartphys.2018.01.003.
  * [5] R.P. Kokoulin et al., Muon excess in ultra-high energy inclined EAS according to the NEVOD-DECOR data, PoS ICRC2021 381 (2021), 10.22323/1.395.0381.
  * [6] D. Heck et al., CORSIKA: A Monte Carlo Code to Simulate Extensive Air Showers, Forschungszentrum Karlsruhe Report FZKA 6019 (1998).
  * [7] C. Patrignani et al. (Particle Data Group), Review of Particle Physics, Rev. Particle Phys. Chin. Phys. C 40 100001 (2016), 10.1088/1674-1137/40/10/100001.
  * [8] R.U. Abbasi et al., Indications of Proton-Dominated Cosmic-Ray Composition above 1.6 EeV, Phys. Rev. Lett. 104 161101 (2010), 10.1103/PhysRevLett.104.161101.
  * [9] A. Aab et al., Depth of maximum of air-shower profiles at the Pierre Auger Observatory. I. Measurements at energies above 1017.8 eV, Phys. Rev. D 90 122005 (2014), 10.1103/PhysRevD.90.122005.
  * [10] R.U. Abbasi et al., Study of Ultra-High Energy Cosmic Ray composition using Telescope Array’s Middle Drum detector and surface array in hybrid mode, Astropart. Phys. 64 49 (2015), 10.1016/j.astropartphys.2014.11.004.
  * [11] H.P. Dembinski et al., Report on Tests and Measurements of Hadronic Interaction Properties with Air Showers, EPJ Web of Conf. 210 02004 (2019), 10.1051/epjconf/201921002004 .
  * [12] H.P. Dembinski et al., Data-driven model of the cosmic-ray flux and mass composition from 10 GeV to 1011 GeV, PoS ICRC2017 533 (2017), 10.22323/1.301.0533.
  * [13] A.A. Petrukhin, Muon puzzle in cosmic ray experiments and its possible solution, Nucl. Instrum. and Meth. in Phys. Res. A 742 228 (2014), 10.1016/j.nima.2013.12.011.
  * [14] R.P. Kokoulin et al., Investigating the Energy Characteristics of Muon Bundles in Inclined EAS, Bull. Russ. Acad. Sci.: Phys. 79 365 (2015), 10.3103/S1062873815030314.
  * [15] R.P. Kokoulin et al., Measuring the Energy Deposited by Muon Bundles of Inclined EAS in the NEVOD-DECOR Experiment, Phys. Part. Nucl. 49 101 (2018), 10.1134/S1063779618010239.
  * [16] E.A. Yurina et al., Measurements of the average muon energy in inclined muon bundles in the NEVOD-DECOR experiment, PoS ICRC2021 383 (2021), 10.22323/1.395.0383.
  * [17] S. Agostinelli et al., Geant4 - a simulation toolkit, Nucl. Instrum. and Meth. in Phys. Res. A 506 250 (2003), 10.1016/S0168-9002(03)01368-8.
  * [18] D.E. Groom, N.V. Mokhov, S.I. Striganov, Muon stopping power and range tables 10-MeV to 100-TeV, Atom. Data and Nucl. Data Tables 78 183 (2001), 10.1006/adnd.2001.0861.


MIDLMedical Imaging with Deep Learning 2021 Short Paper – MIDL 2021 submission
– Under Review Hoang C. Nguyen1<EMAIL_ADDRESS>
Tung T. Le1<EMAIL_ADDRESS>
Hieu H. Pham1,2<EMAIL_ADDRESS>
Ha Q. Nguyen1,2<EMAIL_ADDRESS>
1 Medical Imaging Department, Vingroup Big Data Institute, Hanoi, Vietnam
2 College of Engineering and Computer Science, VinUniversity, Hanoi, Vietnam

# VinDr-RibCXR: A Benchmark Dataset for Automatic Segmentation and Labeling of
Individual Ribs on Chest X-rays

###### Abstract

We introduce a new benchmark dataset, namely VinDr-RibCXR, for automatic
segmentation and labeling of individual ribs from chest X-ray (CXR) scans. The
VinDr-RibCXR contains 245 CXRs with corresponding ground truth annotations
provided by human experts. A set of state-of-the-art segmentation models are
trained on 196 images from the VinDr-RibCXR to segment and label 20 individual
ribs. Our best performing model obtains a Dice score of 0.834 (95% CI,
0.810–0.853) on an independent test set of 49 images. Our study, therefore,
serves as a proof of concept and baseline performance for future research.

###### keywords:

Rib segmentation, CXR, Benchmark dataset, Deep learning.

††editors: Under Review for MIDL 2021

##  1 Introduction

A wide range of diagnostic tasks can benefit from an automatic system that is
able to segment and label individual ribs on chest X-ray (CXR) images. To this
end, traditional approaches [Candemir et al.(2016)Candemir, Jaeger, Antani,
Bagci, Folio, Xu, and Thoma] exploited hand-crafted features to identify the
ribs, but failed with anterior ribs. Recently, deep learning (DL) has shown
superior performance to other methods in the segmentation and labeling of
individual ribs [Wessel et al.(2019)Wessel, Heinrich, von Berg, Franz, and
Saalbach]. However, developing DL algorithms for this task requires annotated
images for each rib structure at pixel-level. To the best of our knowledge,
there exists no such benchmark datasets and protocols. Hence, we present
VinDr-RibCXR – a benchmark dataset for the automatic segmentation and labeling
of individual ribs on CXRs. This work also reports performance of several
state-of-the-art DL-based segmentation models on the VinDr-RibCXR dataset. The
dataset111https://vindr.ai/datasets/ribcxr and
codes222https://github.com/vinbigdata-medical/MIDL2021-VinDr-RibCXR are made
publicly available to encourage new advances.

##  2 Material and Method

###  2.1 Dataset

We built a dataset, namely VinDr-RibCXR, of 245 AP/PA CXR images for
segmentation and labeling of individual ribs. The raw images in DICOM format
were sourced from VinDr-CXR dataset [Nguyen et al.(2020)Nguyen, Lam, Le, Pham,
Tran, Nguyen, Le, Pham, Tong, Dinh, et al.], for which all scans have been de-
identified to protect patient privacy. We then designed an in-house web-based
labeling tool called VinDr Lab (https://vindr.ai/vindr-lab) that allows to
segment individual ribs at pixel-level on DICOM scans. Each image was assigned
to an expert, who manually segmented and annotated each of 20 ribs, denoted as
L1$\to$L10 (left ribs) and R1$\to$R10 (right ribs). The masks of ribs (see
Figure LABEL:fig:demo) were then stored in a JSON file that can later be used
for training instance segmentation models. To the best of our knowledge, the
VinDr-RibCXR is the first publicly released dataset that includes segmentation
annotations of the individual ribs, and for both anterior and posterior ribs.
To develop and evaluate segmentation algorithms, we divided the whole dataset
into a training set of 196 images and a validation set of 49 images.

###  2.2 Model development

We deployed various state-of-the-art DL-based segmentation models for the task
of individual rib segmentation and labeling, including U-Net [Ronneberger et
al.(2015)Ronneberger, Fischer, and Brox], U-Net with EfficientNet-B0 [Tan and
Le(2019)], Feature Pyramid Network (FPN) [Lin et al.(2017)Lin, Dollár,
Girshick, He, Hariharan, and Belongie] with EfficientNet-B0, and U-Net++ with
EfficientNet-B0. These models are well-known to be effective for many medical
image segmentation tasks. To train the networks, we resized all training
images to 512$\times$512 pixels. Dice loss and Adam optimizer were used to
optimize the network’s parameters. All models were trained for 200 epochs with
a learning rate of $1e-3$. Some augmentation techniques (i.e. horizontal flip,
shift, scale, rotate, and random brightness contrast) were applied to avoid
overfitting. The models were implemented using Pytorch (v1.7.1) and trained on
a NVIDIA RTX 2080 Ti GPU.

##  3 Evaluation and Result

The effectiveness of the segmentation models was evaluated by Dice score, 95%
Hausdorff distance (95% HD), sensitivity, and specificity. Experimental
results reported by the proposed segmentation models are shown in Table 1. Our
best performing model, which was Nested U-Net with EfficientNet-B0 encoder,
achieved a Dice score of 0.834 (95% CI, 0.810–0.853), a 95% HD of 15.453 (95%
CI, 13.340–17.450), a sensitivity of 0.841 (95% CI, 0.812–0.858), and a
specificity 0.998 (95% CI, 0.997–0.998). Figure LABEL:fig:demo shows
segmentation and labeling results of the Nested U-Net on some representative
cases from the validation set of VinDr-RibCXR and external datasets including
JSRT and Shenzen. We observe that the algorithm correctly segmented the
individual ribs on both normal and abnormal CXR scans, for both the VinDr-
RibCXR and two external datasets. Although the abnormal CXR contains ribs that
are obscured by lung abnormalities, the proposed algorithm still produced a
considerably accurate segmentation mask for all ribs.

Table 1: Segmentation performance on the validation set of VinDr-RibCXR. Model | Dice | 95% HD | Sensitivity | Specificity
---|---|---|---|---
U-Net | .765 (.737–.788) | 28.038 (22.449–34.604) | .773 (.738–.791) | .996 (.996–.997)
U-Net w. EfficientNet-B0 | .829 (.808–.847) | 16.807 (14.372–19.539) | .844 (.818–.858) | .998 (.997–.998)
FPN w. EfficientNet-B0 | .807 (.773–.824) | 15.049 (13.190–16.953) | .808 (.773–.828) | .997 (.997–.998)
U-Net++ w.  EfficientNet-B0 | .834 (.810–.853) | 15.453 (13.340–17.450) | .841 (.812–.858) | .998 (.997–.998)

##  4 Conclusion

In this work, we introduced a new benchmark dataset and baseline predictive
models for automatically segmenting and labeling individual ribs from chest
radiographs. Our experimental results have shown the effectiveness of the DL-
based segmentation models for this task. This study serves as a proof of
concept and baseline performance for future research. Future works include
investigating more powerful DL models and extending the dataset to enhance the
segmentation performance.

fig:demo

Figure 1: Segmentation results on our validation set and two external
datasets.

##  5 Acknowledgements

This research was supported by the Vingroup Big Data Institute. We are
especially thankful to Tien D. Phan, Dung T. Le, and Chau T.B. Pham for their
help during the data annotation process.

## References

  * [Candemir et al.(2016)Candemir, Jaeger, Antani, Bagci, Folio, Xu, and Thoma] Sema Candemir, Stefan Jaeger, Sameer Antani, Ulas Bagci, Les R Folio, Ziyue Xu, and George Thoma.  Atlas-based rib-bone detection in chest X-rays.  _Computerized Medical Imaging and Graphics_ , 51:32–39, 2016.
  * [Lin et al.(2017)Lin, Dollár, Girshick, He, Hariharan, and Belongie] Tsung-Yi Lin, Piotr Dollár, Ross Girshick, Kaiming He, Bharath Hariharan, and Serge Belongie.  Feature pyramid networks for object detection.  In _IEEE Conference on Computer Vision and Pattern Recognition_ , pages 2117–2125, 2017.
  * [Nguyen et al.(2020)Nguyen, Lam, Le, Pham, Tran, Nguyen, Le, Pham, Tong, Dinh, et al.] Ha Q Nguyen, Khanh Lam, Linh T Le, Hieu H Pham, Dat Q Tran, Dung B Nguyen, Dung D Le, Chi M Pham, Hang TT Tong, Diep H Dinh, et al.  VinDr-CXR: An open dataset of chest X-rays with radiologist’s annotations.  _arXiv preprint arXiv:2012.15029_ , 2020.
  * [Ronneberger et al.(2015)Ronneberger, Fischer, and Brox] Olaf Ronneberger, Philipp Fischer, and Thomas Brox.  U-net: Convolutional networks for biomedical image segmentation.  In _International Conference on Medical Image Computing and Computer-assisted Intervention_ , pages 234–241. Springer, 2015.
  * [Tan and Le(2019)] Mingxing Tan and Quoc Le.  Efficientnet: Rethinking model scaling for convolutional neural networks.  In _International Conference on Machine Learning_ , pages 6105–6114, 2019.
  * [Wessel et al.(2019)Wessel, Heinrich, von Berg, Franz, and Saalbach] Jöran Wessel, Mattias P Heinrich, Jens von Berg, Astrid Franz, and Axel Saalbach.  Sequential rib labeling and segmentation in chest X-ray using Mask R-CNN.  _arXiv preprint arXiv:1908.08329_ , 2019.


# Unsupervised Change Point Detection for heterogeneous sensor signals

Mario Krause
Friedrich-Alexander Universität Erlangen-Nürnberg

###### Abstract

Change point detection is a crucial aspect of analyzing time series data, as
the presence of a change point indicates an abrupt and significant change in
the process generating the data. While many algorithms for the problem of
change point detection have been developed over time, it can be challenging to
select the appropriate algorithm for a specific problem. The choice of the
algorithm heavily depends on the nature of the problem and the underlying data
source. In this paper, we will exclusively examine unsupervised techniques due
to their flexibility in the application to various data sources without the
requirement for abundant annotated training data and the re-calibration of the
model. The examined methods will be introduced and evaluated based on several
criteria to compare the algorithms.

###### Index Terms:

Change point detection, Segmentation, Time series data

##  I Introduction

In the modern world a large amount of data is generated by all the
interconnected devices and networks and collected via sensors [21]. These
sequences of measurements from the sensors are indexed with a time stamp and
called time series data. This massive amount of data can not be monitored
without automation, so time series analysis has become increasingly important.
In the context of change point detection, time series analysis is used to
detect changes in the behaviour of the system, which can be due to internal or
external properties. The task of identifying a specific point in time where
the statistical properties of the underlying model of a signal or time series
change is called change point detection. [2]

Change point detection is an essential topic in time series analysis with a
wide range of applications that require the detection of abrupt changes in the
data, for example to indicate a transition of the system from one state to
another or to indirectly detect the state of a system by estimating that the
state changed [21]. Here are some examples of applications that use change
point detection:

Climate change detection
Climate analysis and climate monitoring have become very important over the
last few decades due to the possible occurrence of climate change caused by
the increase in greenhouse gases in the atmosphere. In this context change
point detection is used to discover climatic discontinuities and changes in
the temperature. [6, 22, 10]

Medical condition monitoring
Change point detection methods can be used to monitor patients health data and
quickly detect abrupt changes of the medical conditions. These algorithms can
be used to analyze changes in physical activity during the strength training
program [14] or to monitor the heart rate of patients undergoing anesthesia
and surgery to indicate disturbances to the cardiovascular system [26].

Stock Market Analysis
In the stock market the fluctuation of any stock price is normal according to
economic theory, but there are some shifts that are abnormal and worth the
investors special attention [4]. In momentum strategies it is necessary to
identify momentum turning points, when a trend reverses from an uptrend to a
downtrend such as in the 2020 market crash due to the covid outbreak [25].

This article presents an overview and comparison of algorithms commonly used
for detecting change points in time series data. The focus is on unsupervised
change point detection, which involves segmenting the data without relying on
large amounts of annotated training data or the need to re-calibrate the model
for each data source. The goal of this article is to help choosing the right
detection method for a particular application, with an emphasis on practical
aspects like the implementation and the calibration of the parameters. Our
selection of methods aims for a good general performance for different data
sources without fine tuning the algorithm. In practice fine tuning a method
for each sensor will most probably yield better performance, but adds
significant overhead if the system has multiple heterogeneous signals with the
number of signals possibly growing in the future. We are especially focusing
on methods that can be applied to different heterogeneous sensor signals, like
temperature or pressure measurements, without the requirement to tune the
parameters separately for each data source. The methods should be able to
detect change points in different 1-dimensional sensor signals from a complex
heterogeneous system to discover dependencies between the signals.

##  II Background

A change point occurs when the statistical properties of the time series data
change abruptly. This change of measurements can be caused, for example, by a
transition in the state of the underlying data generating system. These
changes may be indicated by a simple shift in mean or variance, or more
complex shifts such as in the frequency domain. An abrupt change refers to any
alteration in the system’s parameters that happens either instantly or very
quickly relative to the measurement sampling period. It’s important to note
that abrupt changes don’t necessarily indicate large shifts. In fact,
detecting small changes is often the main challenge in most applications. [3]

In time series analysis mainly two concepts are used to detect deviations from
the normal behavior in a data set: outlier detection and change point
detection. The outlier detection is a method used to identify data points that
are significantly different from the rest of the data, but are mainly very
short shifts in the system or measurement errors. In contrast to that, change
point detection aims to discover points in time where the state of the data
generator changes for a subsequent and longer period. Therefore they are
related, but not entirely similar. Figure 1 illustrates the difference between
these concepts by showing an example time series that contains several change
points, where the mean and variance of the time series switches abruptly, and
outliers, where an individual observation is far from the majority of the
data. The main underlying assumption when investigating change points is that
the properties or parameters describing the data are either constant or slowly
changing over time in each state of the system [3].

Figure 1: Time series with change points and outliers

###  II-A Offline and online algorithms

Change point detection algorithms are traditionally divided into two main
branches: online and offline methods. Although both classes of algorithms aim
to detect a change point in time series data, their methods of operation
differ substantially. The online change point detection algorithms operate in
real-time and run concurrently with the process they are monitoring,
processing each data point as it becomes available. The constraint for these
algorithms is that the processing should be completed before the arrival of
the next data point [5]. These methods are typically used in embedded
applications such as monitoring industrial processes. Online change point
detection methods must be able to process data quickly and make decisions in
near real time. However, it is worth noting that online methods can also be
applied in an offline setting, by analyzing stored data in a sequence.

In contrast to that are the offline algorithms, which retrospectively detect
changes and considers the entire data set. These methods can also be used in
an online setting by collecting a small batch of data and detecting the change
points in this batch. In general the concept of online and offline methods is
a spectrum describing the number of data points the algorithm have to look at
ahead of a potential change point, with complete offline methods considering
the entire data set and complete online methods operating in perfect real-
time. The offline methods allow the detection of multiple change points at
once and focus more on the accuracy of the detection, while the online methods
try to detect the most recent change point as fast as possible. Offline change
point algorithms are typically more computationally intensive than online
algorithms, as they have access to the entire data set and can use more
complex methods to detect change points [2]. The choice between online and
offline methods should be based on the specific requirements of the
application, with the most important consideration being the real-time
constraint. Since the motivation for this paper is an application in an
offline setting with the entire data set available, we will cover online and
offline methods.

###  II-B Supervised and unsupervised algorithms

A variety of supervised and unsupervised methods have been developed and
improved over time to tackle the challenge of change point detection. The
supervised methods utilize machine learning algorithms that are trained on
labeled data, where the change point locations are already known. The
objective of these supervised algorithms is to build a model that can learn
from the limited labeled data and accurately predict the change points for the
entire data set. A variety of supervised methods can be used for this learning
problem, such as decision tree, naive Bayes, support vector machines, logistic
regression and nearest neighbor methods [2]. One of the main drawbacks of
supervised change point detection methods is the need for labeled data, as
well as the need to retrain the model for each individual data source. This
can make the process of applying these supervised methods to different data
sources very costly and time consuming.

Unsupervised learning algorithms are better suited for this as they are used
to identify patterns in unlabeled data. In the context of change point
detection, these algorithms can divide time series data into segments, with
one segment before the change point and one after. These algorithms utilise
statistical aspects of the data to locate change points in the time series.
Unsupervised change point detection is attractive because it can handle a
variety of different data sets without requiring prior training for each data
set [2]. Both types of methods have been shown to be effective, but for the
purposes of this paper, we will focus solely on unsupervised methods. This
decision is based on the underlying application of this paper that requires
methods that are easily applicable to different data sources without
retraining the algorithm. For an overview on supervised methods, we refer
readers to [21, 2].

###  II-C Stability

An important criteria for our application is the stability of the algorithms
in terms of their ability to detect change points in different data sources
without requiring prior fine tuning of the parameters for each data source. We
are comparing the number of parameters and evaluate the effort of calibrating
the algorithm. The parameters should be robust to a wide variation of values
to minimize the effort of the parameter tuning. An important aspect for the
stability of the algorithms is also the assumptions the algorithm makes about
the underlying data source.

In order to increase the stability of the algorithms, we are especially
looking for algorithms that make little to no assumptions, therefor we are
comparing the change point detection algorithms by determining if they are
parametric or non parametric. It is important to make this distinction,
because parametric models make assumptions about the underlying distribution
of the data, such as for example the in the real world often appearing and
therefor commonly used normal distribution. The distribution is described by a
set of parameters, which have to be estimated from the data. The main issue of
parametric models is that the assumptions about the underlying distribution
may not hold for all data sources, which can lead to inaccurate estimates of
the change points. The parametric nature of the algorithm can also be a
benefit and lead to good results if the underlying distribution of the data is
known. In contrast, the non-parametric models do not make assumptions about
the underlying distribution of the data, which makes them more robust to
different data sources. They skip the step of parameter estimation by using a
non parametric approach by for example directly estimating the ratio between
to densities without calculating the densities itself. [2]

###  II-D Algorithm constraints

Change point detection techniques can be differentiated based on the
requirements imposed on the input data and the algorithm itself. These
limitations play a crucial role in determining the most suitable method for
detecting change points in a specific time series. Especially in our
application it is important to have as little constraints as possible to
ensure the algorithm can be applied to multiple data sources. Some algorithms
for example are constraint to stationary or independent and identically
distributed (i.i.d.) data sets. Another important constraint is the number of
change points an algorithm can detect and if the number of change points has
to be known prior to the run of the algorithm. Based on that we are
classifying the algorithms into single and multiple change point detection
algorithms. Algorithms that require the number of change points to be
specified before the run will not be considered, since this is a major
drawback when analysing different time series with an unknown amount of change
points. The algorithms also differ in their outputs, with some algorithms
providing the probability of a change point at a specific time point and
others only returning the change points. The different outputs determine how
easy the results can be interpreted and can for example give insights in the
confidence of the detected change point.

###  II-E Scalability

The amount of time series data generated in the world is rapidly increasing,
both in terms of the number of data points and dimensions. To effectively
handle this massive amount of data, change point detection methods must be
designed to be computationally efficient. Therefore, it is crucial to compare
the computational cost of different change point detection algorithms and
determine which methods can reach an optimal or approximate solution as
quickly as possible. The problem of detecting change points in high
dimensional data sets is not discussed in this paper, since our application is
dealing with univariant time series. The computational costs of the algorithms
are compared based on the information provided by the authors or estimated
based on the algorithmic description.

##  III Review

This section provides an overview of commonly used change point detection
algorithms. We will outline the basic principles of the algorithm and evaluate
them based on the in section II defined criterion. These algorithms include
both online and offline methods and where chosen based on their stability,
algorithm constraints and scalability. All of the algorithms presented are
unsupervised learning methods.

###  III-A Likelihood ratio methods

A typical statistical formulation of change point detection is to compare the
probability distributions $p_{\theta_{0}}$ and $p_{\theta_{1}}$ of the data
before and after a potential change point $x_{t}$. If these distributions are
significantly different, the time point is considered a change point. One
common method involves monitoring the logarithm of the likelihood ratio
between consecutive intervals in the time series data. [12]

$\displaystyle s(x_{i})\text{ =
}\ln\frac{p_{\theta_{1}}(x_{i})}{p_{\theta_{0}}(x_{i})}$

We will discuss two different approaches that utilize the likelihood ratio in
change point detection. The first approach is called the cumulative sum
algorithm and relies on predesigned parametric models, where the probability
density of the two consecutive intervals is calculated separately from the
underlying data and the density ratio is calculated. The second approach is a
direct density ratio estimation method, which is a non parametric algorithm
where the ratio of the probability densities is directly estimated, without
the requirement to perform density estimation for the individual segments.

####  III-A1 CUSUM

On of the most familiar change point detection algorithms is the cumulative
sum algorithm (CUSUM). Page [20] was the first to suggest the use of a
cumulative sum to find changes in a parameter of interest. The CUSUM algorithm
monitors the cumulative sum of the differences between successive measurements
and a reference value or target. If this cumulative sum exceeds a
predetermined threshold, it indicates that there has been a change point in
the process. The log likelihood ratio is a measure of the difference between
two distributions, and can be used to determine whether a change point has
occurred.

The typical behavior of the log likelihood ratio corresponds to a negative
drift before change, and a positive drift after change. Therefore, the
relevant information for change point detection lies in the difference between
the value of the log-likelihood ratio and its current minimum value [3].
Assume we have a sequence ${x_{1},x_{2},...}$ of time series variables with a
unknown change point $x_{t}t$ and the distribution $p_{\theta_{0}}(x_{t})$
before and $p_{\theta_{1}}(x_{t})$ after the change point $x_{t}$. Then the
corresponding decision rule for the CUSUM algorithm is:

$S_{n}\text{ =
}\sum_{i=1}^{n}\ln\frac{p_{\theta_{1}}(x_{i})}{p_{\theta_{0}}(x_{i})}-\min_{k\leq
n}\sum_{i=1}^{k}\ln\frac{p_{\theta_{1}}(x_{i})}{p_{\theta_{0}}(x_{i})}>L\\\ $
(1)

The CUSUM algorithms detects a change point when the CUSUM statistic $S_{n}$
is larger then the predefined threshold $L$. Another common approach, that is
computationally more efficient, is to calculate the CUSUM statistic $S_{n}$
recursively in the following way:

$S_{n}\text{ =
}\max\left(0,S_{n-1}+\ln\frac{p_{\theta_{1}}(x_{i})}{p_{\theta_{0}}(x_{i})}\right)$
(2)

In practice, the CUSUM algorithm is implemented by first setting an initial
value for the CUSUM statistic $S_{0}$, usually set to zero. Then the algorithm
iterates through the observations, updating the CUSUM statistic at each time
step. The algorithm continues until a change point is detected or the end of
the data is reached.

Stability: The CUSUM algorithm is parametric and makes assumptions about the
underlying distribution of the data. In each iteration of the algorithm the
set of parameters describing the distribution before a potential change point
$p_{\theta_{0}}$ and after the change point $p_{\theta_{1}}$ have to be
estimated from the data. A common application is when the process is normally
distributed with the pre-change and post-change distributions having the same
known variance $\sigma$. Then the interest centers on detecting shifts away
from the pre-shift mean $\mu_{0}$ to $\mu_{1}$. The algorithm can handle
different types of distributions, but if the assumption of the underlying
distribution is not correct the result of the algorithm may vary. [11]

The algorithm has only two parameters, with the initial value of the CUSUM
statistic usually set to zero and the threshold value $L$ which is more
difficult to calibrate. The threshold value depends on the specific data set
and the goals of the analysis, for example it can be set such that the false
positive rate or the false negative rate is minimized or determined using
cross-validation techniques. Different methods to calculate values of $L$ have
been developed with a constant or adaptive threshold. A comparison of these
methods is given in [23].

Algorithm constraints: The standard CUSUM algorithm is an online algorithm for
single change point detection, but many variations of this concept were
developed over time for online and offline settings, as well as adaptions to
multiple change point detection by running the algorithm with a sliding
window. An overview of four of these adaptions is provided in [3]. The
algorithm has no constraints or limitations to the data it can be applied to
wide variety of problems. The algorithm does not provide any information about
the confidence of the output, it stops when a single change point is detected,
indicated by the CUSUM statistic exceeding the threshold, and directly outputs
the change point.

Scalability: The main advantage of CUSUM is the simplicity of the underlying
concept, which makes it easy to implement. The computational complexity of the
CUSUM algorithm is generally considered to be low, as it only involves simple
mathematical operations such as addition. The complexity estimated based on
the algorithm is $O(n)$ for a sequence of n points.

####  III-A2 KLIEP

The Kullback–Leibler importance estimation procedure (KLIEP) is the non
parametric counter part to the CUSUM algorithm and also uses the likelihood
ratio to detect change points. KLIEP avoids the problems that come with the
parametric model of the CUSUM algorithm by using a more flexible non
parametric model. This grants the advantage that it does not rely on strong
model assumptions. A naive method would be to use a non parametric approach to
estimate the densities separately and use them to calculate the ratio, but
this approach is ineffective due to the challenges of non parametric density
estimation. Instead of individually estimating each density, we can directly
estimate the density ratio. Different online methods have been developed that
use the idea of direct density ratio estimation [12, 16, 8]. In these methods
the density ratio between two consequent intervals $X$ and $X^{\prime}$ is
modeled by a non parametric Gaussian kernel:

$\hat{w}(X)\text{ = }\frac{p(X)}{p(X^{\prime})}\text{ =
}\sum_{l=1}^{n}\alpha_{l}K_{\sigma}(X,X_{l})$ (3)

where $p(X)$ is the probability distribution of the interval $X$ and the
parameters $\alpha$ are learned from data samples. $K_{\sigma}(Y,Y^{\prime})$
is the Gaussian kernel function with mean $Y^{\prime}$ and standard deviation
$\sigma$. The parameters $\alpha$ are determined in a training phase by
minimizing a dissimilarity measure. From the density ratio estimator
$\hat{w}(X)$ an approximation of the dissimilarity measure between two samples
is calculated and the higher the dissimilarity measure is, the more likely the
point is a change point [8]. One of these methods is the Kullback-Leibler
importance estimation procedure (KLIEP) that estimates the density ratio by
using Kullback-Leibler (KL) divergence:

$\displaystyle KL[p(x)||p^{\prime}(x)]\text{ = }$ $\displaystyle-\int
p^{\prime}(x)\log\frac{p(x)}{p^{\prime}(x)}dx$

The parameters $\alpha$ in the formula 3 are determined so that the empirical
Kullback-Leibler divergence is minimized. The solution of this problem can be
obtained by solving a convex optimization problem. With the estimated
parameters, the logarithm of the likelihood ratio is evaluated as the change
detection score and a change point is detected if the score is beyond a given
threshold $L$

$S\text{ = }\sum_{i=1}^{n}\ln\hat{w}(X_{i})>L$ (4)

Stability: The KLIEP algorithm has several tuning parameters, such as the
kernel width $\sigma$, the lengths of the two intervals $X$ and $X^{\prime}$
and the threshold for the change point score $L$. The main advantage of the
method described in [12] is that it is equipped with a natural cross
validation procedure for tuning the kernel width parameter $\sigma$. The
interval length has an impact on the accuracy of the estimation of the density
ratio, with larger intervals leading to better estimates but also to possible
issues in data sources where change points occur very frequently. The
detection of the change points is sensitive to the threshold parameter $L$,
with possible methods to calibrate $L$ already discussed in chapter 1) CUSUM.

The main advantage in terms of the stability of KLIEP compared to the CUSUM
algorithm is the robustness to different types of data distributions. The
KLIEP algorithm is non parametric and doesn’t make assumptions about the
underlying distribution of the data. This is especially useful if the
underlying distribution of the data is unknown or the algorithm is applied to
different data sources that do not necessarily have the same distribution.

Algorithm constraints: The KLIEP and CUSUM algorithm share a lot of
properties, with their sequential behaviour both are online methods use for
single change point detection and can be adapted for multiple change point
problems with small adjustments. The algorithm has no constraints or
limitations to the data and directly returns the change point.

Scalability: Compared with other approaches like the CUSUM algorithm, methods
based on density ratio estimation tend to be computationally more expensive
because of the cross-validation procedure for model selection [16]. After the
model selection the calculation of the change point score for a given density
ratio is computationally very efficient and involves only simple mathematical
operations that are performed in linear time.

###  III-B Bayesian online change point detection

A common bayesian method for detecting changes, described by Adam and MacKay
[1], involves estimating the posterior distribution of the run length $r_{t}$,
which describes the time that passed since the last change point. The
algorithm is an online method, so in each iteration a new data point is
considered and the run length can either increase by 1 or drop down to 0. The
set of observation $x_{t}^{(r)}$ is associated with the run length $r_{t}$ and
includes all the observations of the current run, with new observations added
to the set until a change point is found, which sets $r_{t}$ to 0 and
$x_{t}^{(r)}$ empty. The algorithm determines the probability
$P(r_{t}=0|x_{1:t})$ of a change point occurring at time $t$ by calculating
the probability of the run length of 0 at time $t$ given the set of
observations $x_{1:t}={x_{1},...,x_{t}}$. The functionality of the algorithm
is visualized in figure 2 with the run length and the change point probability
represented as a grey scale, with darker pixels indicating a higher
probability. The plot shows that the run length drops to zero when the
probability of a change point gets close to one.

Figure 2: Example of bayesian online change point detection from [1]

This probability of a change point can be calculated based on Bayes theorem
as:

$\displaystyle
P(r_{t}|x_{1:t})=\frac{\sum_{r_{t-1}}P(r_{t}|r_{t-1})P(x_{t}|r_{t-1},x_{t}^{r})P(r_{t-1},x_{1:t-1})}{P(x_{1:t})}$

Since $P(r_{t-1}|x_{t-1})$ is a recursive component, it is known from the
previous step and we only need to calculate the conditional run length
probability $P(r_{t}|r_{t-1})$ and evaluate the predictive distribution
$P(x_{t}|r_{t-1},x_{t}^{r})$. The conditional prior on the change point
$P(r_{t}|r_{t-1})$ is only nonzero at the two outcomes $r_{t}$=0 or
$r_{t}$=$r_{t-1}$, which gives the algorithms its computational efficiency.
[1]

$P(r_{t}|r_{r-1})=\left\\{\begin{array}[]{ll}1-H(r_{t-1}+1),&\hbox{if
}r_{t}=r_{t-1}+1\\\ H(r_{t-1}+1),&\hbox{if }r_{t}=0\\\ 0,&\hbox{otherwise}\\\
\end{array}\right.$

The function H(r) is called the hazard function and describes how likely a
change point is occurring at a run length r. A common approach is to make this
process memoryless by setting H(r) = $\frac{1}{\lambda}$ with a timescale
parameter $\lambda$. Another option would be to make H(r) increasing over the
run to penalize longer run lengths. The predictive distribution
$P(x_{t}|r_{t-1},x_{t}^{r})$ represents the probability that the most recent
observation belongs to current run and is the most challenging to calculate.
In [1] it was proposed to use a conjugate exponential model with parameters
$\nu$ and $\chi$ to make the algorithm computationally efficient.

Stability The bayesian online change point detection algorithm is a parametric
method. The method allows us to encode knowledge about the world into the
algorithm by setting the two conjugate priors $\chi$ and $\nu$. The conjugate
priors $\chi$ and $\nu$ are the only two parameters of the algorithm and are
set based on prior knowledge. This prior knowledge could be the mean and
variance of the data, which is estimated roughly based on the first few
observations. A possible option for further tuning the algorithm it by setting
the hazard function H(r), but the common memoryless approach is sufficient for
most applications.

Algorithm constraints: The algorithm is a online multiple change point
detection method that assumes the sequence of observations can be segmented
into non-overlapping partitions and the data within each partition p is i.i.d.
from some probability distribution $P(x_{t}|\mu_{p})$ [1]. The restriction to
i.i.d. time series is a major drawback for this algorithm when applying it to
multiple different data sources. The algorithm directly returns the
probability for a change point occurring at the specific run lengths, which
makes the results easier to interpret compared to methods like CUSUM or KLIEP
that only return the change point itself.

Scalability: The implementation of the algorithm proposed by Adam and MacKay
[1] is quadratic in space and time complexity in the number of data points $n$
so far observed. Clearly this is problematic when analysing long time series.
To overcome this issue a simple approximation was introduced in [17] that
reduced the complexity to $O(n)$.

###  III-C Singular spectrum transformation

The singular spectrum transformation (SST) algorithm first proposed by
Moskvina and Zhigljavsky [19] uses singular spectrum analysis and adapts it
for change point detection. The singular spectrum transform is a non
parametric approach and can be used to analyze time series with complex
structure [19]. This technique involves transforming the original time series
into a new series of change point scores. The resulting time series can be
interpreted as the probability distribution that some change point occurs at
time $t$. In figure 3 such a change point score with the corresponding time
series is shown with a clear peak in the change point score when the frequency
of the time series changed.

Figure 3: Example of a time series and the corresponding change point score
from [9]

The underlying idea is to compute for each time point $x_{t}$ the difference
between a representative pattern of a few time points before and a few points
after $x_{t}$. The dynamics of the time series are represented using a Hankel
matrix. We call the Hankel matrix, representing the change patterns within the
past $w$ points, trajectory matrix $H(t)$ and the Hankel matrix for
representing the future change patterns the test matrix $G(t)$. The
representative pattern of $H(t)$ and $G(t)$ can be extracted by performing a
singular value decomposition on both matrices. The $l<w$ left singular vectors
$U_{i_{1}}(t),...,U_{i_{l}}(t)$ of $H(t)$ with the largest singular values
represent the past change pattern as a hyperplane and build the matrix
$U(t)\text{ = }U_{i_{1}}(t),...,U_{i_{l}}(t)$. The matrix $U(t)$ encodes the
major direction of change in the past signal and the parameter $l$ is the
number of representative patterns that are considered. The importance of each
representative pattern $U_{i}(t)$ is given by the corresponding singular
value, with the most dominant pattern corresponding to the largest singular
value. The direction of maximum change in the future of the signal is given by
the left singular vector $\beta(t)$ of $G(t)$ with the largest singular value.
[9]

To estimate the difference between the past and the future patterns,
$\beta(t)$ will be projected onto $U(t)$ and normalized to calculate a change
point score. If there is no change in the dynamics of the signal, it is
expected that $\beta(t)$ will lie in or very near to the hyperplane
represented by $U_{l}$. The change point score ranges from zero to one, with a
high likelihood of a change point occurring when the score is close to one.
The score can be calculated at any time point $t$ by finding representative
patterns in both the trajectory and test matrices. This can be seen as a
transformation from the original time series $T$ to a new time-series $T_{c}$.
This demonstrates that SST can make variables of different types comparable by
converting a heterogeneous system into a homogeneous one [9], i.e.

$T\rightarrow T_{c}(w,l,g,m,n)$ (5)

Stability: The main problem of the SST algorithm is the need to specify five
parameters: the length of the column vectors of the Hankel matrix $w$, the
number of columns of H(t) $n$, the number of singular vectors $l$, the shift
of the starting point for the future signal $g$ and the number of columns of
G(t) $m$. In [9, 18] it was shown that SST is usually robust to a wide range
in different values of $w$ and $n$ and domain knowledge or visualization can
help finding an appropriate value for both. The problem lies in choosing the
other three parameters, where domain knowledge is not very useful. Despite the
problem of specifying five different parameters it was shown in [9] that for a
wide range of $w$ (6 < $w$ < 40), the results of the algorithms with
$w=-g=m=n$ and $l=3$ are quite robust and the essential features remain
unchanged.

Algorithm constraints: The SST is based on the singular value decomposition of
two Hankel matrix, so it is a non parametric method, as it does not rely on
any specific assumptions about the distribution of the data, making it
suitable for a wide range of time series data [19]. The advantage of the SST
algorithms lies in its comparability of different time series. The algorithms
returns a change point score for each time point, which allows to detect
multiple change points in one run. The change point score can be interpreted
as the probability that a change point occurred at a specific time point,
which makes it easy to interpret the output of the algorithm and discover
dependencies within multiple different heterogeneous time series. [9]

Scalability: The complexity of the singular spectrum transformation is linear
in the length of the time series [18]. In each of the $n$ time steps a
singular value decomposition of the fixed size Hankel matrices $H(t)$ and
$G(t)$ have to be computed, which makes SST less computational efficient as
other change point detection algorithms with linear complexity.

###  III-D Binary Segmentation

Binary segmentation is a common technique for offline change point detection
in time series data due to its simplicity both in the underlying concept and
the implementation of the algorithm. A implementation of the algorithm
available in the ruptures python library for offline change point detection
and is called BinSeg [24]. The binary segmentation procedure is one of the
standard methods for detecting multiple change points by using a test for
single change point detection [15].

The algorithm is a recursive technique for multiple change point detection in
which initially the entire data set is searched for one change point. Once a
change point is detected, the data set is split into two subsegments, defined
by the detected change point. The algorithm is then performed on each of the
two subsegments resulting in further splits. This procedure continuous in each
new split until a stopping criteria is satisfied or no change point is
detected. The indices of the segment boundaries are the change points. This
process of recursively splitting the segments into two smaller subsegments is
shown in the schematic view of the algorithm in figure 4.

Figure 4: Schematic example of binary segmentation from [24]

This technique can be used with a variety of time series and is able to find
changes in the mean, variance, or distribution of the data. Binary
segmentation is considered a greedy algorithm, because it is performed
sequentially, with each stage only visited once and depending on the previous
ones. [7]

The algorithm starts by initializing the starting and ending indices $s$ and
$e$ of the time series sequence, for the first run start $s=1$ and end $e=n$,
where $n$ is the number of data points. Then a cost function $c(x_{a...b})$ is
defined to measure the cost of a segment $x_{a...b}$, with
$x_{a...b}={x_{i}|a<=i<=b}$. The cost function should be chosen based on the
nature of the data and the goal of the analysis. An overview and comparison of
common cost functions is given in [24]. The algorithm calculates the costs for
each possible split by summing the cost function for the two resulting segment
and returns the change point $\hat{k}$ that minimizes the costs:

$\hat{k}=\arg\min_{s<k<e}c(x_{s...k})+c(x_{k...e})$ (6)

The algorithm then recursively calls the binary segmentation algorithm with
the new start and end indices: BinSeg($s$, $k$) and BinSeg($k$, $e$). This
process is repeated until a stopping criterion is met, such as a minimum
number of data points in a segment, the maximum number of change points or a
maximum number of iterations. Once the stopping criterion is met, the change
points are assigned as the indices of the segment boundaries.

Stability: The algorithm is very simple and has only the cost function and a
stopping criteria as inputs. It’s important to note that the choice of the
cost function $c()$ and the stopping criterion are crucial to the performance
of the algorithm. The stopping criterion can be chosen according to the number
of change points expected in the data, the granularity of the segments, or a
threshold of the cost function for a segment. In general, the choice of cost
function and stopping criterion depends on the nature of the data and the goal
of the analysis.

The algorithm can be used in a parametric and non parametric setting by
selecting the respective cost function, depending on the prior knowledge
before performing the change point detection. The choice of the cost function
is determining the assumption the algorithm makes about the data. In [24] an
overview of possible parametric and non parametric cost functions is given
with their underlying assumptions about the data.

Algorithm constraints: Binary segmentation is a simple and versatile offline
algorithm for change point detection, that can be applied to any time series
without any limitations. The algorithm can be used in situations where the
number of change points is known prior to the run but also for an unknown
number of change points. However, the algorithm has a limitation in accurately
detecting change points that are too close together and result in small
segments. This can lead to inaccuracies in the results for time series, that
frequently change between different states resulting in close change points.
The binary segmentation algorithm directly returns the detected change points,
with the most significant change point appearing first. This feature allows
the user to stop the algorithm at any point in time and still obtain a
meaningful result, making it an efficient tool for change point detection.

Scalability: The benefits from the simplicity of binary segmentation includes
the low complexity of the algorithm of $O(Cnlogn)$, where $n$ is the number of
samples and $C$ the complexity of calling the considered cost function on one
subsegment [24].

###  III-E Bottom up Segmentation

Bottom-up segmentation is the natural counterpart of binary segmentation. In
contrast to binary segmentation, bottom up segmentation starts by dividing the
time series into individual observations, each treated as a separate segment,
followed by sequentially merging adjacent segments based on a discrepancy
criterion. A schematic overview of the bottom up segmentation is provided in
figure 5, showing the different steps of the algorithm starting with a
segmented time series in a grid and merging these segments together until a
stopping criterion is met and the approximation of the change points is
returned.

Figure 5: Schematic example of bottom up segmentation from [24]

Bottom up segmentation is like its counterpart binary segmentation available
in the ruptures python library [24]. It is also an offline method with a very
simple underlying concept and low computational complexity.

The algorithm for bottom up segmentation starts by creating the finest
possible representation of the time series of size $n$ by dividing it into
$\frac{n}{2}$ segments. The indices of the segment boundaries are the
potential change points. For each potential change point $t$ a discrepancy
measure $d(x_{a...t},x_{t...b})$ is calculated between the segments
$x_{a...t}$ and $x_{t...b}$ separated by $t$, with these segments defined as
$x_{a...b}=\\{x_{i}|a\leq i\leq b\\}$. For a given cost function $c()$, the
discrepancy measure between two adjacent segments is given by:

$d(x_{a...t},x_{t...b})=c(x_{a...b})-c(x_{a...t})-c(x_{t...b})$ (7)

All potential change points are ranked by the discrepancy measure with the
lowest discrepancy change point getting deleted, meaning that the
corresponding segments are merged. When two adjacent segments $t$ and $t+1$
are merged, the algorithm performs some bookkeeping tasks and calculates the
cost of merging the newly formed segment with its right neighbor segment.
Additionally, the costs for merging the left neighbor with the newly formed
larger segment have to be recalculated. These bookkeeping steps are necessary
to ensure the accuracy of the bottom up segmentation process and its ability
to detect change points in the time series. [13]

The process of merging and evaluating segments is repeated until no further
improvements can be made. This results in a set of segments that represent the
best approximation of the underlying process changes in the time series. The
final step is to identify the change points, which correspond to the
boundaries between the segments and represent the points where the underlying
process changes.

Stability: Due to the simple underlying concept of the bottom up algorithm the
only two inputs are the cost function, used in the discrepancy measure, and a
stopping criteria. The choice of the two inputs have a significant impact on
the performance of the algorithm. The stopping criterion can be chosen
according to the number of change points expected in the data, the size of the
segments, or a threshold of the discrepancy measure for two segments. In
general, the cost function and stopping criterion should be chosen based on
the nature of the time series and the goal of the analysis.

Similar to binary segmentation this algorithm can be parametric and non
parametric, depending on the selected cost function. The cost function should
be chosen based on the prior knowledge of the data source before performing
the change point detection. The choice of the cost function is determining the
assumption the algorithm makes about the data. In [24] an overview of possible
parametric and non parametric cost functions is given with their underlying
assumptions.

Algorithm constraints: The algorithm has no limitations in the time series it
can be applied to and can be used for a known and unknown number of change
points. Due to the very fine grid in the initialization the first iterations
of merging procedures can be unstable because of the small segments they are
performed on, for which the statistical significance is smaller. This problem
could be overcome by starting with larger segments. However if a true change
point does not belong to the original set of boarders of these segments, it
would never be considers for a change point, resulting in an inaccurate
detection of the change points. [24]

Scalability: The algorithm for bottom up segmentation has a similar complexity
as their counterpart binary segmentation with $O(C$ $n$ $log(n))$, where $n$
is the number of samples in the time series and $C$ the complexity of calling
the discrepancy measure for two adjacent subsegments [13].

##  IV Discussion and comparison

The previous sections present an overview of change point detection algorithms
that are commonly used in the literature and have a good general performance
for different data sources without fine tuning each algorithm for the
particular data source. The task of selecting the algorithm best suited for a
particular application can be challenging and depends on which criterion is
most important for the application. We compare change point detection methods
based on their stability, algorithm constraints and their scalability to help
with this task.

###  IV-A Stability

When evaluating the stability of algorithms, one important factor to consider
is the number of parameters and the robustness of these parameters to various
changes. While all algorithms except singular spectrum transformation have a
small number of parameters, the choice of these parameters can have a
significant impact on the performance of the algorithms. For instance, in the
case of the CUSUM or KLIEP algorithm, only one threshold parameter $L$ needs
to be assigned. However, the threshold value in CUSUM or KLIEP and the other
stopping criterion in binary segmentation or bottom up segmentation can be
very sensitive, and choosing an incorrect threshold value can result in either
an early or a late stopping of the algorithm, which leads to inaccurate
results. In contrast, singular spectrum transformation requires five different
parameters to be assigned, but it has been shown to be effective in detecting
change points in a variety of data sources without adjusting the parameters
[9]. Similarly, bayesian online change point detection algorithm only requires
the specification of the conjugate priors as input, which is essentially a
problem of estimating the distribution of the data. The absence of a stopping
criteria in both singular spectrum transformation and bayesian online change
point detection makes these algorithms relatively more stable compared to
other algorithms.

To enhance the stability of a change point detection method, it is important
to consider the assumptions the algorithm makes about the data. In general,
non parametric change point detection methods tend to be more robust than
parametric methods, as the latter heavily depend on the choice of parameters
used to model the distribution of the data [2]. A naive but common approach is
to use a parametric model and always assume that the different data sources
have a normal distribution. Since this assumption is true for a lot of real
world applications this approach may produce good results for data sets that
actually follow a normal distribution, but it can lead to poor performance
when the data deviates from a normal distribution. In contrast, non-parametric
methods, such as KLIEP, binary segmentation, or bottom-up segmentation, do not
make any assumptions about the underlying distribution of the data and are
therefore more robust. A comprehensive overview of the algorithms and their
parameters, along with a distinction between parametric and non-parametric
methods, can be found in Table I.

Table I: Comparison of change point detection algorithms stability Method |  Parametric/
non parametric |  Parameters
---|---|---
Cusum |  Parametric |  threshold value $L$
KLIEP |  Non parametric |  kernel width $\sigma$
lengths of the two intervals
threshold value $L$
BOCPD |  Parametric |  conjugate priors $\chi$ and $\nu$
SST |  Nonparametric |  length of the representative pattern w
number of columns of H(t)
number of singular vectors l
shift of the starting point g
number of columns of G(t)
BinSeg |  Parametric/
non parametric |  cost function
stopping criteria
Bottom Up |  Parametric/
non parametric |  cost function
stopping criteria

###  IV-B Algorithm Constraints

When choosing the appropriate algorithm for a particular application, it is
important to take into account the constraints of the algorithm. This review
focuses on three main constraints: the restrictions on the time series, the
number of change points that can be detected, and the output of the algorithm.
The Likelihood ratio methods have no limitations on the time series, but they
are only capable of detecting a single change point, so multiple runs of the
algorithm may be required to detect multiple change points. This problem is
overcome by directly using methods such as binary or bottom-up segmentation,
which are designed for multiple change point detection and also do not have
any limitations on the time series. Bayesian online change point detection has
the only restriction on the time series in that it requires the time series to
be independently and identically distributed. Singular Spectrum
Transformation, on the other hand, requires the time series to be stationary.
It is important to consider these constraints when selecting the right
algorithm for a given application to ensure that the algorithm is capable of
meeting the specific needs and requirements of the data.

The various algorithms also have different outputs, with some only providing
the change points and others providing a change point probability or score.
The bayesian online change point detection and singular spectrum
transformation algorithms provide a change point probability or score, but
this information requires a separate post processing step to assign the actual
change points. However, this additional information provides insight into the
confidence of the algorithm in the identified change points. On the other
hand, the likelihood ratio methods and binary and bottom up segmentation
algorithms directly return the change points, avoiding the need for a separate
post processing step. A comprehensive overview of the algorithms and their
constraints can be found in table II.

Table II: Comparison of change point detection algorithms constraints Method |  Time series
limitations |  No. change points detected |  Output
---|---|---|---
Cusum |  No limitation |  Single |  Change point
KLIEP |  No limitation |  Single |  Change point
BOCPD |  i.i.d. time series |  Multiple |  Change point probability
SST |  Time series should be stationary |  Multiple |  Change point score
BinSeg |  No limitation |  Multiple |  Change point
Bottom Up |  No limitation |  Multiple |  Change point

###  IV-C Scalability

Another important criteria to consider are the computational costs of change
point detection algorithms. Table III presents a comparison of the
computational cost of the algorithms surveyed. For some algorithms the
computational costs where provided by the authors. In cases where the authors
have not provided this information, the comparison was conducted based on the
descriptions of the algorithms.

The computational cost of CUSUM is relatively low, as it requires only a few
calculations for each time step. The time complexity of CUSUM is linear in the
number of data points $n$, making it a suitable option for large data sets.
The KLIEP algorithm, on the other hand is a slightly more complex algorithm
with a higher time complexity as CUSUM due to the model selection required.

The computational cost of bayesian online change point detection depends on
the choice of priors and the complexity of the model used. In general, it is
computationally more expensive and has a quadratic time complexity $O(n^{2})$,
but this can be reduced to a linear complexity with the help of a simple
approximation.

The singular spectrum transform has also a linear time complexity but is less
efficient as CUSUM or KLIEP as it requires the calculation of a singular value
decomposition in each time step, which has a high time complexity.

The computational complexity of binary segmentation and bottom up segmentation
depends on the cost function used in the algorithm. Both have a complexity of
$O(C$ $n$ $log(n))$, where $n$ is the number of samples in the time series and
$C$ the complexity of calling the cost function.

This review covered both online and offline change point detection methods. In
general the online methods are computationally less expensive, because they
focus on detecting the most recent change point as fast as possible. This was
also confirmed by the algorithms in this review with two offline methods
binary segmentation and bottom up segmentation having the highest
computational complexity.

Table III: Comparison of the algorithms computational complexity

Method |  Online / Offline method |  Computational complexity
---|---|---
Cusum |  online |  O(n)*
KLIEP |  online |  Cusum < KLIEP*
BOCPD |  online |  O(n)
SST |  online |  KLIEP < SST*
BinSeg |  offline |  $O(n\log n)$
Bottom Up |  offline |  $O(n\log n)$

*Estimated based on the algorithm

##  V Conclusion

In this review, we have analyzed various techniques for change point detection
and organized them under a unified framework. Our emphasis was on unsupervised
methods that have a good general performance for diverse data sources. To
assess the methods, we compared them based on three criteria: stability,
algorithmic constraints, and scalability. The aim is to provide a framework
that can be utilized to assess future algorithms in a systematic manner. The
focus of this review was chosen based on a particular application that was the
motivation for this paper. It is worth noting that in practice fine tuning a
specific algorithm for a single time series will most probably yield better
performance, but significantly increases the time and effort needed to detect
change points in a large system of heterogeneous signals. The comparison of
the methods was based solely on the description of the algorithm and leaves
room for further research by for example comparing the performance of the
methods on different data sources and evaluating the results.

## References

  * [1] Ryan Prescott Adams and David J. C. MacKay.  Bayesian online changepoint detection.  2007\.
  * [2] Samaneh Aminikhanghahi and Diane J Cook.  A survey of methods for time series change point detection.  Knowledge and information systems, 51(2):339–367, 2017.
  * [3] Michèle Basseville and Igor V. Nikiforov.  Detection of Abrupt Changes - Theory and Application.  Prentice Hall, Inc. - http://people.irisa.fr/Michele.Basseville/kniga/, 1993.
  * [4] Jie Chen and Arjun K Gupta.  Parametric statistical change point analysis: with applications to genetics, medicine, and finance.  2012\.
  * [5] Allen B. Downey.  A novel changepoint detection algorithm.  2008\.
  * [6] Jean-François Ducré-Robitaille, Lucie A. Vincent, and Gilles Boulet.  Comparison of techniques for detection of discontinuities in temperature series.  International Journal of Climatology, 23(9):1087–1101, 2003.
  * [7] Piotr Fryzlewicz.  Wild binary segmentation for multiple change-point detection.  The Annals of Statistics, 42(6):2243 – 2281, 2014.
  * [8] Mikhail Hushchyn and Andrey Ustyuzhanin.  Generalization of change-point detection in time series data based on direct density ratio estimation.  J. Comput. Sci., 53:101385, 2020.
  * [9] Tsuyoshi Ide and Keisuke Inoue.  Knowledge discovery from heterogeneous dynamic systems using change-point correlations.  Proceedings of the 2005 SIAM International Conference on Data Mining, SDM 2005, 04 2005.
  * [10] Naoki Itoh and Jürgen Kurths.  Change-point detection of climate time series by nonparametric method.  1:445–448, 2010.
  * [11] Daniel R. Jeske, Veronica Montes De Oca, Wolfgang Bischoff, and Mazda Marvasti.  Cusum techniques for timeslot sequences with applications to network surveillance.  Computational Statistics and Data Analysis, 53(12):4332–4344, 2009\.
  * [12] Yoshinobu Kawahara and Masashi Sugiyama.  Sequential change-point detection based on direct density-ratio estimation.  Statistical Analysis and Data Mining: The ASA Data Science Journal, 5(2):114–127, 2012.
  * [13] Eamonn Keogh, Selina Chu, David Hart, and Michael Pazzani.  An online algorithm for segmenting time series.  In Proceedings 2001 IEEE international conference on data mining, pages 289–296. IEEE, 2001.
  * [14] Patrick Lapointe, Kévin Chapron, Isabelle Lessard, Kevin Bouchard, Mélissa Lavoie, Cynthia Gagnon, Elise Duchesne, and Sébastien Gaboury.  Monitoring changes in physical activity data during strength training of people with myotonic dystrophy type 1.  Procedia Computer Science, 198:51–58, 2022.  12th International Conference on Emerging Ubiquitous Systems and Pervasive Networks / 11th International Conference on Current and Future Trends of Information and Communication Technologies in Healthcare.
  * [15] Marc Lavielle and Gilles Teyssiere.  Adaptive detection of multiple change-points in asset price volatility.  Long memory in economics, pages 129–156, 2007.
  * [16] Song Liu, Makoto Yamada, Nigel Collier, and Masashi Sugiyama.  Change-point detection in time-series data by relative density-ratio estimation.  Neural Networks, 43:72–83, 2013.
  * [17] Rakesh Malladi, Giridhar P Kalamangalam, and Behnaam Aazhang.  Online bayesian change point detection algorithms for segmentation of epileptic activity.  pages 1833–1837, 2013.
  * [18] Yasser Mohammad and Toyoaki Nishida.  Robust singular spectrum transform.  pages 123–132, 2009.
  * [19] Valentina Moskvina and Anatoly Zhigljavsky.  An algorithm based on singular spectrum analysis for change-point detection.  Communications in Statistics-simulation and Computation - COMMUN STATIST-SIMULAT COMPUT, 32:319–352, 01 2003.
  * [20] E. S. Page.  Continuous inspection schemes.  Biometrika, 41(1/2):100–115, 1954.
  * [21] Aditya Pushkar, Muktesh Gupta, Rajesh Wadhvani, and Manasi Gyanchandani.  A comparative study on change-point detection methods in time series data.  pages 1–7, 2022.
  * [22] Jaxk Reeves, Jien Chen, Xiaolan L. Wang, Robert Lund, and Qi Qi Lu.  A review and comparison of changepoint detection techniques for climate data.  Journal of Applied Meteorology and Climatology, 46(6):900 – 915, 2007.
  * [23] Nassim Sahki, Anne Gégout-Petit, and Sophie Wantz-Mézières.  Performance study of change-point detection thresholds for cumulative sum statistic in a sequential context.  Quality and Reliability Engineering International, 36(8):2699–2719, 2020.
  * [24] Charles Truong, Laurent Oudre, and Nicolas Vayatis.  Selective review of offline change point detection methods.  Signal Processing, 167:107299, 2020.
  * [25] Kieran Wood, Stephen J. Roberts, and Stefan Zohren.  Slow momentum with fast reversion: A trading strategy using deep learning and changepoint detection.  2021\.
  * [26] P. Yang, G. Dumont, and J.M. Ansermino.  Adaptive change detection in heart rate trend monitoring in anesthetized children.  IEEE Transactions on Biomedical Engineering, 53(11):2211–2219, 2006\.