# Generalized multiple change-point detection in the structure of
multivariate, possibly high-dimensional, data sequences

Andreas Anastasiou111Department of Mathematics and Statistics, University of
Cyprus, P.O. Box: 20537, 1678, Nicosia, Cyprus<EMAIL_ADDRESS>and Angelos Papanastasiou222Department of Mathematics and Statistics,
University of Cyprus, P.O. Box: 20537, 1678, Nicosia, Cyprus,
<EMAIL_ADDRESS>
###### Abstract

The extensive emergence of big data techniques has led to an increasing
interest in the development of change-point detection algorithms that can
perform well in a multivariate, possibly high-dimensional setting. In the
current paper, we propose a new method for the consistent estimation of the
number and location of multiple generalized change-points in multivariate,
possibly high-dimensional, noisy data sequences. The number of change-points
is allowed to increase with the sample size and the dimensionality of the
given data sequence. Having a number of univariate signals, which constitute
the unknown multivariate signal, our algorithm can deal with general
structural changes; we focus on changes in the mean vector of a multivariate
piecewise-constant signal, as well as changes in the linear trend of any of
the univariate component signals. Our proposed algorithm, labeled Multivariate
Isolate-Detect (MID), allows for consistent change-point detection in the
presence of frequent changes of possibly small magnitudes in a computationally
fast way.

Keywords: Change-point detection; high-dimensional setting; piecewise-constant
structure; piecewise-linear structure.

##  1 Introduction

Change-point detection algorithms have been actively developed and
investigated from the scientific community. Their ability to segment data into
smaller, homogeneous parts have helped researchers and practitioners to
develop flexible statistical models that can adapt in non-stationary
environments. Due to the natural data heterogeneity in many real problems,
such algorithms have been applied in a wide range of application areas, such
as bioinformatics (Picard et al. (2011); Hocking et al. (2013)), cyber
security (Siris and Papagalou (2004)), or finance (Lavielle and Teyssiere
(2007); Schröder and Fryzlewicz (2013)). The advantages of detecting changes
in the behaviour of the data fall into two main categories; interpretation and
forecasting. Interpretation comes naturally since the detected changes are
usually connected with life events that took place near the estimation time.
Associating the changes with such real-life phenomena can lead to a better
understanding and quantification of the effect that these events had on the
behaviour of the stochastic process. With respect to forecasting, the role of
the final segment is important because it allows for a more accurate
prediction of the future values of the data sequence at hand.

Based on whether we have full knowledge of the data to be analysed, change-
point detection is split into two main categories; offline detection, where
the data are already obtained, and online detection, in which the observations
arrive sequentially at present. With respect to the dimensionality of the
data, change-point detection can be further separated into algorithms that act
only on univariate data and to those that are suitable for change-point
detection in multivariate data sequences. In this paper, we focus on
multivariate, possibly high-dimensional, offline settings; our aim is to
estimate the number and locations of certain structural changes in the
behaviour of given multivariate data. The model is

$\bm{X_{t}}=\bm{f_{t}}+\Sigma\bm{\epsilon_{t}},\hskip 28.45274ptt=1,\ldots,T,$
(1.1)

where $\bm{X_{t}}\in\mathbb{R}^{d\times 1}$ are the observed data and
$\bm{f_{t}}\in\mathbb{R}^{d\times 1}$ is the $d$-dimensional deterministic
signal with structural changes at certain points. The signals that we treat in
the current manuscript are those that changes appear in the mean structure or
in the vector of the first order derivatives. The diagonal matrix
$\Sigma\in\mathbb{R}^{d\times d}$ has diagonal elements denoted by
$\sigma_{1},\ldots,\sigma_{d}$, while for any
$t\in\left\\{1,\ldots,T\right\\}$, the noise terms
$\bm{\epsilon_{t}}\in\mathbb{R}^{d\times 1}$ are random vectors with mean the
zero vector and covariance the identity matrix. The elements
$\sigma_{j},j=1,\ldots,d$ of the diagonal matrix $\Sigma$ might be unknown. In
such cases, $\sigma_{j}$ are estimated using the Median Absolute Deviation
method explained in Hampel (1974). The true number, $N$, of the change-points,
as well as their locations $r_{1},\ldots,r_{N}$, are unknown and our aim is to
estimate them; $N$ is free to grow with the sample size $T$ and the
dimensionality $d$.

The initial purpose of change-point detection algorithms has been to detect a
single change in the mean structure of a univariate signal under the setting
of Gaussian noise, but much progress has since been made. Researchers have
heavily focused on the detection of multiple change-points in the mean
structure of a univariate data sequence. Towards this purpose, optimization-
based methods have been developed, in which the estimated signal is chosen
based on its fit to the data, penalized by a complexity rule. To solve the
implied penalization problem, dynamic programming approaches, such as the
Segment Neighborhood (SN) and Optimal Partitioning (OP) methods of Auger and
Lawrence (1989) and Jackson et al. (2005), have been developed. Killick et al.
(2012) and Rigaill (2015) introduce improvements over the classical OP and SN
algorithms, respectively. In the context of regression problems, Frick et al.
(2014) introduced the simultaneous multiscale change-point estimator (SMUCE)
for change-point detection in exponential family regression. Apart from
optimization based algorithms, a popular method in the literature is binary
segmentation where changes are detected one at a time through an iterative
binary splitting of the data. Recent variants of binary segmentation with
improved performance are the Wild Binary Segmentation (WBS) of Fryzlewicz
(2014) and its recently developed second version of Fryzlewicz (2020), the
Narrowest-Over-Threshold (NOT) method of Baranowski et al. (2019), and the
Seeded Binary Segmentation of Kovács et al. (2020). The Isolate-Detect (ID)
algorithm has been developed in Anastasiou and Fryzlewicz (2022) to detect,
one by one, structural changes in a data sequence. It is based on an isolation
technique, which leads to very good accuracy on the estimated number and
locations of the change-points, particularly in scenarios with many frequent
change-points. Our proposed method is partly based on ID; therefore we
elaborate on its important parts later. For a more thorough review of the
literature on the detection of multiple change-points in the mean of
univariate data sequences, see Cho and Kirch (2020) and Yu (2020). Apart from
changes in the mean of a univariate data sequence, research has also been done
for the detection of change-points under more complex scenarios, such as
detection of changes in the slope for piecewise-linear models (Anastasiou and
Fryzlewicz (2022); Baranowski et al. (2019); Fearnhead et al. (2019); Maeng
and Fryzlewicz (2019); Tibshirani (2014)) changes in the variance (Inclán and
Tiao (1994)), as well as for distributional changes under a non-parametric
setting (Matteson and James (2014); Zou et al. (2014)).

Even though there is an extensive literature on change-point detection for
univariate data sequences, the multivariate, possibly high-dimensional
setting, which is the focus of this paper, has not been investigated in such
degree. Working under the model in (1.1), Vert and Bleakley (2010) proposes a
method for approximating the signal $\bm{f_{t}}$ as the solution of a convex
optimization problem. In order to achieve this, the problem is first
reformulated to a group LASSO one and then the group least-angle regression
(LARS) procedure explained in Yuan and Lin (2006) is employed. Another
interesting approach with very good behaviour has been introduced in Wang and
Samworth (2018). The algorithm, called inspect, estimates the number and
locations of the change-points in the mean structure of $\bm{f_{t}}$ as in
(1.1). Firstly, inspect applies a cumulative sum (CUSUM) transformation to the
original data matrix. Secondly, a projection direction of the transformed
matrix is computed as its leading left singular vector, and, finally, a
univariate change-point detection algorithm is applied to the projected
series. It is among the many methods that employ CUSUM-type statistics for
change-point detection in the multivariate setting. In general, methods that
belong to this category, mainly either use CUSUM aggregations of the $d$
component data sequences in order to test the obtained values against a
threshold or to construct alternative test-statistics. For instance, Groen et
al. (2013) examines the asymptotic behaviour of the maximum absolute and
average CUSUM and gives finite-sample performance results. Focusing on testing
for the existence of a change-point in the mean structure of the multivariate
signal, Enikeeva and Harchaoui (2019) and Horváth and Hušková (2012) propose
tests based on the $\ell_{2}$ aggregation of the CUSUM statistics for each
univariate component, while Jirak (2015) employs the $\ell_{\infty}$
aggregation of the aforementioned values. In Cho (2016), the Double CUSUM (DC)
operator is introduced, which takes as input the ordered absolute CUSUM values
of each individual component and performs a weighted $\ell_{1}$ aggregation to
construct the DC statistic which is then compared against a test criterion.
Departing from the detection of changes in the mean structure of a
multivariate signal, Cho and Fryzlewicz (2015) propose the Sparsified Binary
Segmentation algorithm (SBS) for the detection of multiple change-points in
the second-order structure of a multivariate data sequence. SBS is based on a
first, “sparsifying” step which is used to exclude individual component data
sequences from an $\ell_{1}$ aggregation; a pre-defined threshold is used for
the exclusion. The recent work of Anastasiou et al. (2022) introduces Cross-
covariance isolate detect (CCID), which, motivated from the necessity of
estimating changes in time-varying functional connectivity networks, detects
multiple change-points in the second-order (cross-covariance or network)
structure of multivariate, possibly high-dimensional time series. Ombao et al.
(2005) investigates the application of smooth localized complex exponentials
(SLEX) waveforms, to the detection of changes in spectral characteristics of
EEG data. Lavielle and Teyssiere (2006) detects changes in the covariance
structure of i.i.d. multivariate time series based on the minimization of a
penalized Gaussian log likelihood, while Bücher et al. (2014) uses a test
statistic based on sequential empirical copula processes to detect changes in
the cross-covariance structure. For a survey of various offline change-point
detection algorithms on multivariate time series see Truong et al. (2020).

In this paper, we propose a method called Multivariate Isolate-Detect (MID)
for the consistent estimation of multiple change-points under the
multivariate, possibly high-dimensional, structure of the model in (1.1). Our
method builds on the foundations of the ID algorithm developed in Anastasiou
and Fryzlewicz (2022); that is, we first isolate each true change-point within
subintervals of the domain $[1,\ldots,T]$ and then we proceed to detect them.
Isolation enhances detection power, especially in frameworks with frequently
occurring change-points. MID is explained in detail in Section 2; here, we
only give a brief description of its important steps. The main idea is that
for the observed data sequences $x_{t,j}\hskip 5.69046ptt=1,\ldots,T,\quad
j=1,\ldots,d,$ and with $\lambda_{T}$ a positive constant, playing the role of
an expansion step as in Anastasiou and Fryzlewicz (2022), our method first
creates two ordered sets of $K=\lceil T/\lambda_{T}\rceil$ right- and left-
expanding intervals. For $i=1,\ldots,K$, the $i^{th}$ right expanding interval
is $R_{i}=[1,\min\left\\{i\lambda_{T},T\right\\}]$ while the $i^{th}$ left-
expanding interval is $L_{i}=[\max\left\\{1,T-i\lambda_{T}+1\right\\},T].$ We
collect these intervals in the ordered set
$S_{RL}=\\{R_{1},L_{1},R_{2},L_{2},\ldots,R_{K},L_{K}\\}$. The algorithm first
acts on the interval $R_{1}=[1,\lambda_{T}]$ by calculating, for every
univariate component data sequence, the contrast function value for the $Q$
possible candidates in this interval (details are given in Section 2). This
process will return $Q$ vectors $\bm{y_{j}},j=1,\ldots,Q$ of length $d$ each;
for example, the elements of $\bm{y_{1}}\in\mathbb{R}^{d}$ will be the
contrast function values related to the first change-point candidate in
$R_{1}$, for each of the $d$ component data sequences, the elements of
$\bm{y_{2}}\in\mathbb{R}^{d}$ will be the relevant values for the second
candidate in $R_{1}$, and so on. The next step is to apply to each
$\bm{y_{j}}$ a mean-dominant norm $L:\mathbb{R}^{d}\rightarrow\mathbb{R}$. The
definition of such norms can be found in Section 2 of Carlstein (1988) and
examples include

$\displaystyle
L_{2}:=L_{2}(\bm{y_{j}})=\frac{1}{\sqrt{d}}\sqrt{\sum_{i=1}^{d}y_{j,i}^{2}}$
$\displaystyle
L_{\infty}:=L_{\infty}(\bm{y_{j}})=\sup_{i=1,\ldots,d}\left\\{y_{j,i}\right\\},$
(1.2)

where $y_{j,i}\geq 0,\;\forall i,j$. Applying $L(\cdot)$ to each $\bm{y_{j}}$,
will return a vector $\bm{v}$ of length $Q$. We identify $\tilde{b}_{R_{1}}$
:= ${\rm argmax}_{j}\left\\{v_{j}\right\\}$. If $v_{\tilde{b}_{R_{1}}}$
exceeds a certain threshold, denoted by $\zeta_{T,d}$, then
$\tilde{b}_{R_{1}}$ is taken as a change-point. If not, then the process tests
the next interval in $S_{RL}$. Upon detection, the algorithm makes a new start
from the end-point (respectively, start-point) of the right- (respectively,
left-) expanding interval where the detection occurred. Upon correct choice of
$\zeta_{T,d}$, MID ensures that we work on intervals with at most one change-
point.

The rest of the paper is structured as follows. Section 2 describes our
proposed methodology and its associated theory. Useful variants of our
algorithm and the choice of important parameter values are explained in
Section 3. In Section 4, we perform a thorough simulation study to compare our
algorithm with state-of-the-art methods, while Section 5 illustrates the
behaviour of MID on two examples of real data; the monthly percentage changes
in the UK house price index over a period of twenty two years in twenty London
Boroughs, and the daily number of new COVID-19 cases in the four constituent
countries of the United Kingdom; England, Northern Ireland, Scotland, and
Wales. We conclude the paper in Section 6 with some general remarks and
reflections on our proposed algorithm. The proofs of Theorem 2.1 and Theorem
2.2 are given in the Appendix.

##  2 Methodology and Theory

###  2.1 Methodology

We work under the model given in (1.1). The objective is to estimate both the
number, $N$, and the locations $r_{1},\ldots,r_{N}$ where the multivariate
deterministic signal $\bm{f_{t}}$ exhibits structural changes. We note that
$N$ can possibly grow with the sample size $T$ and with the dimensionality
$d$. In addition, a change-point does not necessarily appear in all univariate
component signals. Before providing a full, step-by-step explanation of our
algorithm, two simple examples are given to assist in ease of understanding.
In Figure 1, we have a three dimensional data sequence of length $T=200$ with
three change-points in the mean vector at locations $r_{1}=27,r_{2}=73$ and
$r_{3}=165$. To be more precise,

$\displaystyle f_{t,1}=\begin{cases}0,&t=1,\ldots,27\\\ 6,&t=28,\ldots,165\\\
0,&t=166,\ldots,200\end{cases},$ $\displaystyle
f_{t,2}=\begin{cases}0,&t=1,\ldots,73\\\ -6,&t=74,\ldots,165\\\
0,&t=166,\ldots,200\end{cases},$ (2.3) $\displaystyle
f_{t,3}=0,t=1,\ldots,200.$

We take $\sigma_{1}=3,\sigma_{2}=1,$ and $\sigma_{3}=2$, while the random
variables $\epsilon_{t,i}$ follow the standard Gaussian distribution for
$i=1,2,3$. The component data sequences $X_{t,1}$ and $X_{t,2}$ share a common
change-point at $t=165$, while they also have their own change-points at
$t=27$ and $t=73$, respectively. There are no change-points in $X_{t,3}$.

Figure 1: An example of a three dimensional data sequence that undergoes three
changes in its mean structure at locations $r_{1}=27,r_{2}=73$ and
$r_{3}=165$.

Let us denote by $r_{0}=0$, $r_{N+1}=T$. In this toy example, we take the
expansion parameter $\lambda_{T}=10$, while $\zeta_{T,d}$ is a well-chosen
predefined threshold. The choice of the aforementioned parameters is discussed
in detail in Section 3.2; special attention will be given to the
dimensionality of the data sequence in order to make a robust threshold
choice. Figure 2 shows the steps of MID until all change-points are detected.
We will be referring to Phases 1, 2 and 3 involving five, six, and five
intervals, respectively; these are clearly indicated in Figure 2. These phases
are only related to this specific example of three change-points; in cases
with a different number of change-points we would have a different number of
such phases. At the beginning of the detection process, we have that $s=1$ and
$e=T=200$. As already mentioned in Section 1, the proposed algorithm acts both
sequentially and interchangeably on subintervals of the full domain; this
being $[1,\ldots,200]$ for this example. For a well-chosen threshold
$\zeta_{T,d}$, then $r_{1}=27$ is the first change-point to be detected; this
occurs in the interval $[1,30]$ as shown in Phase 1 of Figure 2. We now
briefly explain how the detection occurred for $r_{1}$. Let $A$ be the
$30\times 3$ matrix, with each column being the first 30 observations (since
we are working in the interval $[1,30]$) of each of the three univariate data
sequences. The next step is to compute the contrast function (in this specific
case of piecewise-constant signals, the function is the absolute value of the
widely used CUSUM statistic as given in (2.4)) values for each candidate point
and for all three component data sequences. We end up with a matrix
$B\in\mathbb{R}^{29\times 3}$, with $B_{i,j}$ being the value of the contrast
function for the $i^{th}$ data point of the $j^{th}$ data sequence when we
work in the interval $[1,30]$; the last point of the interval is not among the
change-point candidates. Applying a mean-dominant norm to each row of $B$
gives us a vector of length 29. Figure 3 (Detection 1), illustrates these
values, when we employed the $L_{2}$ and the $L_{\infty}$ mean-dominant norms
as defined in (1). In Figure 3, we see that for both employed norms, $t=27$
has the highest value, which exceeds the predefined threshold value obtained
as in Section 3.2. Therefore, $\hat{r}_{1}=27$ is assigned as the estimated
location for $r_{1}$. After the detection, $s$ is updated as the end-point of
the (right-expanding) interval where the detection occurred; therefore $s=30$,
and MID is, in Phase 2, applied in the interval $[30,200]$. Then, $r_{3}=165$
gets detected at the sixth step of Phase 2 in the interval $[161,200]$. After
this second detection, MID proceeds to Phase 3, where it is applied in the
interval $[30,161]$ and $r_{2}$ gets isolated (for the first time) and
detected in the interval $[30,80]$ as shown in Figures 2 and 3. In the end,
MID is applied in the interval $[80,161]$, where there will be no expanding
interval that contains a point with an aggregated CUSUM value that surpasses
the threshold $\zeta_{T,d}$; therefore, the process will terminate after
scanning all the data.

In Figure 4, we graphically provide an example of a three dimensional data
sequence of length $T=200$ with three change-points in the slope at locations
$r_{1}=53,r_{2}=100$ and $r_{3}=124$. To be more precise,

$\displaystyle f_{t,1}=\begin{cases}-t+1,&t=1,\ldots,53\\\
2t-158,&t=54,\ldots,124\\\ -t+214,&t=125,\ldots,200\end{cases},$
$\displaystyle f_{t,2}=\begin{cases}-t+1,&t=1,\ldots,100\\\
2t-299,&t=101,\ldots,124\\\ -t+73,&t=125,\ldots,200\end{cases},$
$\displaystyle f_{t,3}=t,t=1,\ldots,200.$

We take $\sigma_{1}=\sigma_{2}=\sigma_{3}=7$, while $\epsilon_{t,i}$ follow
the standard Gaussian distribution for $i=1,2,3$. The first two component data
sequences have two change-points each; for $X_{t,1}$ at locations $t=53$ and
$t=124$, while for $X_{t,2}$ at locations $t=100$ and $t=124$. There are no
change-points in $X_{t,3}$.

Figure 2: An example with three change-points in the mean structure;
$r_{1}=27,r_{2}=73$ and $r_{3}=165$. The dashed line is the interval in which
the detection took place in each phase.

Figure 3: CUSUM values in the relevant detection intervals explained in Figure
2. The vertical dashed line indicates the time point with the highest value in
the corresponding interval, while the horizontal dashed line is the optimal
threshold value. On the left column you can see the results when the $L_{2}$
aggregation method was used, while for the right column we employed the
$L_{\infty}$-based aggregation approach.

Figure 4: An example of a three dimensional data sequence with piecewise-
linear structure, that undergoes three changes in its first derivative at
locations $r_{1}=53,r_{2}=100$ and $r_{3}=124$.

After giving two examples of structures that MID can be employed to, our
method can now be described in a general framework. Our proposed algorithm is
based on the same isolation technique as that of the univariate change-point
detection method ID and therefore, extensive details of how this isolation is
achieved are avoided and can be found in Section 3.1 of Anastasiou and
Fryzlewicz (2022). As shown in Figure 2 for the specific example in Figure 1,
our method is looking for change-points interchangeably in right- and left-
expanding intervals. In each of them, MID acts in the same way: Let $I$ be one
of these intervals, with cardinality $\mid I\mid$, and
$A^{I}\in\mathbb{R}^{\mid I\mid\times d}$ is such that $A^{I}_{i,j}=X_{i,j}$
for $i\in I$ and $j\in\left\\{1,\ldots,d\right\\}$. The contrast function
value is then separately applied to the elements of each of the columns of
$A^{I}$; this will result in the matrix $B^{I}\in\mathbb{R}^{J\times d}$,
whose element at position $(i,j)$ corresponds to the contrast function value
of the $i^{th}$ data point in $I$ for the $j^{th}$ data sequence. The number
$J$ of rows of the matrix $B^{I}$ (with the contrast function values) depends
on the type of changes that our algorithm aims to find; for example, it has
already been discussed in the literature that in the case of changes in the
mean structure $J=\mid I\mid-1$, while if we are looking for changes in the
slope of a piecewise-linear signal, then $J=\mid I\mid-2$.

The next step is to apply to each row of $B^{I}$, a mean-dominant norm as
those in (1) which will lead to a vector $v\in\mathbb{R}^{J}$. We identify
$\hat{b}_{I}={\rm argmax}_{t}\\{v_{t}\\}$ and $v_{\hat{b}_{I}}$ is tested
against the predefined threshold $\zeta_{T,d}$; heuristically,
$v_{\hat{b}_{I}}$ will be small if $b$ is not a change-point and large in the
opposite scenario. If $v_{\hat{b}_{I}}>\zeta_{T,d}$, then $\hat{b}_{I}$ is
taken as a change-point and MID restarts and looks for change-points in non-
previously tested intervals. The algorithm stops when there are no more
intervals to be tested.

###  2.2 Theory

We work under the setting in (1.1). Denoting by $r_{0}=0$ and $r_{N+1}=T$, the
illustration scenarios are:
(S1) Changes in the mean structure: For $k=1,\ldots,N+1$,
$\bm{f_{t}}=\bm{\mu_{k}}\in\mathbb{R}^{d}$ for $t=r_{k-1}+1,\ldots,r_{k}$. In
this case, the univariate component signals $f_{t,j}$, for $j=1,\ldots,d$ are
piecewise-constant.
(S2) Changes in the first derivative: For $k=1,\ldots,N+1$,
$\bm{f_{t}}=\bm{\mu_{1,k}}+\bm{\mu_{2,k}}t$ for $t=r_{k-1}+1,\ldots,r_{k}$,
where $\bm{\mu_{1,k}}$ and $\bm{\mu_{2,k}}$ are vectors in $\mathbb{R}^{d}$.
In addition, we require that for $j=1,\ldots,N$, the equality
$\bm{\mu_{1,j}}+\bm{\mu_{2,j}}r_{j}=\bm{\mu_{1,j+1}}+\bm{\mu_{2,j+1}}r_{j}$ is
satisfied. Under this framework, the change-points, $r_{k}$, satisfy that
$\bm{f_{r_{k}-1}}+\bm{f_{r_{k}+1}}\neq 2\bm{f_{r_{k}}}$. Therefore, the
univariate component signals $f_{t,j}$, for $j=1,\ldots,d$ are continuous and
piecewise-linear.

The aforementioned scenarios are only two specific illustration cases in which
the proposed MID algorithm can be applied. Due to its change-point isolation
step prior to detection, our algorithm can be applied in more complicated
scenarios where each univariate signal could be, for example, piecewise-
polynomial or piecewise-exponential. In Sections 2.2.1 and 2.2.2, we provide
the main theorems for the consistency of our method in accurately estimating
the true number and the location of the change-points in Scenarios (S1) and
(S2), respectively. The theoretical results presented in this section are for
the $L_{\infty}$ norm being employed for the aggregation of the information
from the component data sequences; in practice, our algorithm exhibits very
good behaviour for other mean-dominant norms as well, something which is shown
in Sections 4 and 5.

####  2.2.1 Scenario (S1)

As already discussed, the first step of the detection process depends on an
appropriately chosen contrast function, which, for every component data
sequence, is applied to each change-point candidate. In Scenario (S1) , the
contrast function applied to the component data sequences $X_{t,j},\forall
j\in\left\\{1,\ldots,d\right\\}$ is the absolute value of the widely used
CUSUM statistic, with the latter being defined as

$\tilde{X}_{s,e}^{b,j}=\sqrt{\frac{e-b}{n(b-s+1)}}\sum_{t=s}^{b}X_{t,j}-\sqrt{\frac{b-s+1}{n(e-b)}}\sum_{t=b+1}^{e}X_{t,j},$
(2.4)

where $1\leq s\leq b<e\leq T$ and $n=e-s+1$. Before proceeding with the main
theoretical result for the consistency of our method, allow us to introduce
some more notation, as below.

$\displaystyle\delta_{T}:=\min_{j=1,\ldots,N+1}\mid r_{j}-r_{j-1}\mid,$
$\displaystyle\Delta_{j}^{q}:=\mid f_{r_{j+1},q}-f_{r_{j},q}\mid,\quad
j=1,\ldots,N,\quad q=1,\ldots,d$
$\displaystyle\underline{f}_{T}:=\inf_{j=1,\ldots,N}\left\\{\sup_{q=1,\ldots,d}\Delta_{j}^{q}\right\\}.$
(2.5)

For the consistency result with respect to the estimated number of change-
points and their estimated locations obtained by MID, we work under the
assumption (A1) as follows:

  * (A1)

The minimum distance, $\delta_{T}$, between two change-points and the minimum
magnitude of jumps, $\underline{f}_{T}$ as in (2.2.1), are connected by
$\sqrt{\delta_{T}}\underline{f}_{T}\geq\underline{C}\sqrt{\log\left(Td^{1/4}\right)}$,
for a large enough constant $\underline{C}$.

The number of change-points, $N$, is allowed to grow with the sample size $T$
and the dimensionality $d$. Theorem 2.1 provides the main theoretical result
for Scenario (S1). The proof is given in Appendix A.

######  Theorem 2.1.

Let $\left\\{\bm{X_{t}}\right\\}_{t=1,\ldots,T}$ follow model (1.1) with
$\bm{f_{t}}$ as in Scenario (S1) and $\bm{\epsilon_{t}}$ are i.i.d. random
vectors from the $d$-variate standard normal distribution. Let $N$ and
$r_{j},j=1,\ldots,N$ be the number and locations of the change-points, while
$\hat{N}$ and $\hat{r}_{j},j=1,\ldots,\hat{N}$ are their estimates sorted in
increasing order. In addition, for $[s_{j},e_{j}]$ being the interval where
$\hat{r}_{j}$ is obtained, we denote by $q_{j}:={\rm
argmax}_{k=1,\ldots,d}\mid\tilde{X}_{s_{j},e_{j}}^{\hat{r}_{j},k}\mid$.
Assuming that (A1) holds, then, there exist positive constants
$C_{1},C_{2},C_{3},$ and $C_{4}$, which do not depend on $T$ or $d$, such that
for
$C_{1}\sqrt{\log\left(Td^{1/4}\right)}\leq\zeta_{T,d}<C_{2}\sqrt{\delta_{T}}\underline{f}_{T}$,

${\rm pr}\left(V\leq C_{3}\log\left(Td^{\frac{1}{4}}\right)\right)\geq
1-\frac{C_{4}}{T},$

where
$V=\left\\{\hat{N}=N,\underset{j=1,\ldots,N}{\max}\left\\{\mid\hat{r}_{j}-r_{j}\mid\left(\Delta_{j}^{q_{j}}\right)^{2}\right\\}\right\\}$.

The lower bound for the probability in (2.1) does not depend on the
dimensionality $d$ and its order is $1-\mathcal{O}\left(T^{-1}\right)$.
Furthermore, the rate of convergence of the estimated change-point locations
does not depend on the minimum distance between two change-points,
$\delta_{T}$; the jump magnitude, $\Delta_{j}^{q_{j}}$, is the only quantity
that affects the rate. We notice, though, that to be able to match the
estimated change-point locations with the true ones, then $\delta_{T}$ should
be larger than the distance between the estimated and the true change-point
locations. Therefore, based on (2.1), we deduce that $\delta_{T}$ must be at
least $\mathcal{O}\left(\log\left(Td^{1/4}\right)\right)$. This, combined with
Assumption (A1) that requires $\delta_{T}(\underline{f}_{T})^{2}$ to be of
order at least $\mathcal{O}\left(\log\left(Td^{1/4}\right)\right)$, means that
$\underline{f}_{T}$ could decrease with $T$ in cases where $\delta_{T}$ is of
order higher than $\mathcal{O}\left(\log\left(Td^{1/4}\right)\right)$.

With respect to the threshold parameter, $\zeta_{T,d}$, the rate of its lower
bound is $\mathcal{O}\left(\sqrt{\log\left(Td^{1/4}\right)}\right)$; this will
also be used in practice as the default rate. Therefore, we have that

$\zeta_{T,d}=C\sqrt{\log\left(Td^{1/4}\right)},$ (2.6)

where $C$ is a positive constant. More details on the choice of $C$ are given
in Section 3.2.

####  2.2.2 Scenario (S2)

We are under the scenario where for any $j\in\left\\{1,\ldots,d\right\\}$, the
underlying signal $f_{t,j},t=1,\ldots,T$ has a continuous and piecewise-linear
structure as in Figure 4. In this case, the contrast function applied to the
component data sequences $X_{t,j},\forall j\in\left\\{1,\ldots,d\right\\}$ is

$C_{s,e}^{b}(\bm{X_{j}})=\mid\left\langle\bm{X_{j}},\bm{\phi_{s,e}^{b}}\right\rangle\mid,$
(2.7)

where for $n=e-s+1$ and

$\displaystyle\alpha_{s,e}^{b}=\sqrt{\frac{6}{n(n^{2}-1)(1+(e-b+1)(b-s+1)+(e-b)(b-s))}}$
$\displaystyle\beta_{s,e}^{b}=\sqrt{\frac{(e-b+1)(e-b)}{(b-s+1)(b-s)}},$

we have that the contrast vector, $\phi_{s,e}^{b}(t)$, is equal to

$\displaystyle\alpha_{s,e}^{b}\beta_{s,e}^{b}\left[(e+2b-3s+2)t-(be+bs-2s^{2}+2s)\right],\text{
for }t\in[s,b]$
$\displaystyle\frac{\alpha_{s,e}^{b}}{\beta_{s,e}^{b}}\left[(2e^{2}+2e-be-
bs)-(3e-2b-s+2)t\right],\text{ for }t\in[b+1,e]$ $\displaystyle
0,\text{otherwise}.$ (2.8)

For more details on how this vector is constructed for (S2), please see
section B.2 in the online supplementary material of Baranowski et al. (2019).
Proceeding now with the consistency result of MID applied to Scenario (S2), we
make the following assumption.

  * (A2)

The minimum distance, $\delta_{T}$, between two change-points and the minimum
magnitude of changes,
$\underline{f}_{T}:=\underset{j=1,\ldots,N}{\inf}\left\\{\underset{q=1,\ldots,d}{\sup}\mid
f_{r_{j}-1,q}+f_{r_{j}+1,q}-2f_{r_{j},q}\mid\right\\}$, are connected by
$\delta_{T}^{3/2}\underline{f}_{T}\geq\underline{C}^{*}\sqrt{\log\left(Td^{1/4}\right)}$,
where $\underline{C}^{*}$ is a large enough constant.

Theorem 2.2 provides the main theoretical result for Scenario (S2). The proof
can be found in Appendix B.

######  Theorem 2.2.

Let $\left\\{\bm{X_{t}}\right\\}_{t=1,\ldots,T}$ follow model (1.1) with
$\bm{f_{t}}$ as in Scenario (S2) and $\bm{\epsilon_{t}}$ are i.i.d. random
vectors from the $d$-variate standard normal distribution. Let $N$ and
$r_{j},j=1,\ldots,N$ be the number and locations of the change-points, while
$\hat{N}$ and $\hat{r}_{j},j=1,\ldots,\hat{N}$ are their estimates sorted in
increasing order. In addition, for $[s_{j},e_{j}]$ denoting the interval where
$\hat{r}_{j}$ is obtained during MID, we denote by $q_{j}:={\rm
argmax}_{k=1,\ldots,d}\left\\{C_{s_{j},e_{j}}^{\hat{r}_{j}}(\bm{X_{k}})\right\\}$.
With $\Delta_{j}^{q}:=\mid f_{r_{j}-1,q}+f_{r_{j}+1,q}-2f_{r_{j},q}\mid$ and
assuming that (A2) holds, then, there exist positive constants
$C_{1},C_{2},C_{3},$ and $C_{4}$, which do not depend on $T$ or $d$, such that
for
$C_{1}\sqrt{\log\left(Td^{1/4}\right)}\leq\zeta_{T,d}<C_{2}\delta_{T}^{3/2}\underline{f}_{T}$,
we obtain that

${\rm pr}\left(\tilde{V}\leq
C_{3}\left(\log\left(Td^{\frac{1}{4}}\right)\right)^{1/3}\right)\geq
1-\frac{C_{4}}{T},$ (2.9)

where
$\tilde{V}=\left\\{\hat{N}=N,\underset{j=1,\ldots,N}{\max}\left\\{\mid\hat{r}_{j}-r_{j}\mid\left(\Delta_{j}^{q_{j}}\right)^{2/3}\right\\}\right\\}$.

Similarly as in (S1), the upper bound for the probability in (2.9) does not
depend on the dimensionality of the data sequence $\bm{X_{t}}$; its rate of
convergence is $1-\mathcal{O}\left(T^{-1}\right)$. Furthermore, the rate of
convergence of the estimated change-point locations depends only on the change
magnitude, $\Delta_{j}^{q_{j}}$, as defined in the statement of Theorem 2.2.
The lower bound for the threshold, $\zeta_{T,d}$, is of order
$\mathcal{O}\left(\sqrt{\log\left(Td^{1/4}\right)}\right)$; this is the
default rate that is used in practise in Sections 4 and 5. Therefore,

$\zeta_{T,d}=C^{*}\sqrt{\log\left(Td^{1/4}\right)},$ (2.10)

where $C^{*}>0$. More details on the choice of the values for $C^{*}$ are
given in Section 3.2.

##  3 Practicalities and variants

###  3.1 Mean-dominant norms

The proposed MID methodology is based on an aggregation step of the contrast
function values obtained from each component data sequence using mean-dominant
norms in the way these are defined in Section 2 of Carlstein (1988). Even
though Theorems 2.1 and 2.2 cover the theoretical behaviour of MID under the
$L_{\infty}$ norm given in (1), similar results could be obtained when other
norms are employed. In the remaining sections, the discussion on the choice of
the parameter values as well as results on the practical performance of MID
will be focused when our method is combined with the $L_{\infty}$ or the
$L_{2}$ mean-dominant norm, as defined in (1), for the aggregation step of the
contrast function values.

###  3.2 Choice of parameter values

In order to choose the constants $C$ in (2.6) and $C^{*}$ in (2.10), we ran a
large scale simulation study involving data sequences
$\left\\{\bm{X_{t}}\right\\}_{t=1,\ldots,T}$, for $T=700,1400$ following model
(1.1), where $\bm{f_{t}}=\bm{0}$ while $\bm{\epsilon_{t}}$ are i.i.d. from the
$d$-variate standard normal distribution. Specifically, for each
$d\in\\{1,\ldots,50\\}$ we generated 500 replicates and applied MID under
scenarios (S1) and (S2) to each one of those replicates using various
threshold constant values $C$ and $C^{*}$, respectively, in order to estimate
the number of change-points. For each dimension, in order to control the Type
I error rate, $\alpha$, of falsely detecting change-points, we chose the
default constant to be the one that its number of times that successfully did
not detect any change-points was closer to $(1-\alpha)*500$. For $d>50$, we
keep the threshold that gave the best results for the simulated 50-dimensional
data sequence. Table 1 presents the results for
$\alpha\in\left\\{0.05,0.1\right\\}$ under Scenarios (S1) and (S2). From now
on, the obtained values for $C$ and $C^{*}$ will be referred to as the default
constants.

Table 1: The optimal values for the threshold constants, $C$ and $C^{*}$, which control the Type I error rate $\alpha$ under the Scenarios (S1) and (S2), respectively, for $d=1,\ldots,50$. Results are presented for the $L_{2}$ and the $L_{\infty}$ norms | Results for the (S1) scenario | Results for the (S2) scenario
---|---|---
Norm | $d$ |  $C$ ($\alpha=0.05$) |  $C$ ($\alpha=0.1$) | $d$ |  $C^{*}$ ($\alpha=0.05$) |  $C^{*}$ ($\alpha=0.1$)
| 1 | 1.7 | 1.55 | 1 | 1.65 | 1.55
| 2 | 1.25 | 1.25 | 2 | 1.25 | 1.2
| 3 | 1.1 | 1.05 | 3 | 1.05 | 1.05
| 4 | 1.05 | 0.95 | 4 | 0.95 | 0.95
| 5 | 0.95 | 0.9 | 5 | 0.9 | 0.9
| 6 | 0.9 | 0.9 | 6 | 0.9 | 0.85
$L_{2}$ | 7 | 0.9 | 0.8 | 7 | 0.8 | 0.8
| 8 | 0.8 | 0.8 | 8 | 0.8 | 0.75
| 9 | 0.8 | 0.75 | 9-11 | 0.75 | 0.75
| 10-13 | 0.75 | 0.75 | 12-16 | 0.7 | 0.7
| 14 | 0.75 | 0.65 | 17-19 | 0.65 | 0.6
| 15-20 | 0.7 | 0.65 | 20-22 | 0.6 | 0.6
| 21-23 | 0.65 | 0.6 | 24-42 | 0.6 | 0.55
| 24-39 | 0.6 | 0.6 | 43-50 | 0.55 | 0.55
| 40-50 | 0.6 | 0.55 |  |  |
| 1 | 1.7 | 1.55 | 1 | 1.65 | 1.55
| 2-3 | 1.75 | 1.7 | 2 | 1.7 | 1.6
| 4-6 | 1.8 | 1.7 | 3 | 1.75 | 1.6
$L_{\infty}$ | 7-13 | 1.85 | 1.75 | 4-5 | 1.75 | 1.65
| 14-25 | 1.9 | 1.8 | 6-13 | 1.75 | 1.7
| 26-28 | 1.9 | 1.85 | 14-25 | 1.8 | 1.75
| 29-50 | 1.95 | 1.85 | 26-38 | 1.85 | 1.8
|  |  |  | 39-50 | 1.9 | 1.85

###  3.3 Decision on the aggregation method

In Cho and Fryzlewicz (2015), Cho (2016), and Anastasiou et al. (2022), it has
been explained that the $L_{\infty}$-based aggregation of the contrast
functions for each component data sequence tends to exhibit a better behaviour
(compared to the $L_{1}$ and $L_{2}$-based aggregations) in scenarios where
the true change-points appear only in a small number of the component data
sequences. In contrast, due to spuriously large contrast function values, the
$L_{\infty}$ approach could suffer in situations where the change-points are
not sparse across the panel of the data sequences.

This difference in the behaviour between the $L_{2}$ and the $L_{\infty}$
norms examined in our paper is what has motivated us to introduce an
introductory step in MID, where the sparsity in a given $d$-dimensional data
sequence, $\bm{X_{t}}$, is first estimated. For $r_{1},\ldots,r_{N}$ being the
$N$ true change-points, allow us, for any $j\in\left\\{1,\ldots,N\right\\}$,
to define by $A_{j}\subseteq\left\\{1,\ldots,d\right\\}$ the set of indices
for the univariate component data sequences that contain the change-point
$r_{j}$. Then, for $\bm{X_{t}}$, the sparsity is given by

$sp=\max_{j=1,\ldots,N}\left\\{\mid A_{j}\mid\right\\}/d,$ (3.11)

where $\mid A_{j}\mid$ is the cardinality of the set $A_{j}$. It is
straightforward that $sp\in[0,1]$. Our proposed hybrid approach aims to data
adaptively decide on the aggregation method to be used and eliminates the
necessity of choosing between the $L_{2}$ and $L_{\infty}$ norms. The steps to
achieve our purpose are as follows.
Step 1: We apply MID paired with the $L_{\infty}$ aggregation rule in order to
obtain $\hat{r}_{1},\ldots,\hat{r}_{M}$. It has already been explained that
the $L_{\infty}$ norm is preferable and provides very good results in cases
with sparse change-points, while it tends to overestimate the number of
change-points (due to the contrast function taking spuriously large values)
under the scenario of dense change-points. Therefore, $M\geq N$.
Step 2: We will now estimate the sparsity in the given data as defined in
(3.11). With $\hat{r}_{0}=0$ and $\hat{r}_{M+1}=T$, we first collect the
triplets $(\hat{r}_{m-1}+1,\hat{r}_{m},\hat{r}_{m+1}),\forall
m\in\left\\{1,\ldots,M\right\\}$. After this, $\forall
i\in\left\\{1,\ldots,d\right\\}$, we calculate
$CS^{(i)}(\hat{r}_{m}):=C^{\hat{r}_{m}}_{\hat{r}_{m-1}+1,\hat{r}_{m+1}}(\bm{X_{i}})$,
with $C^{b}_{s,e}(\bm{X_{i}})$ being the relevant contrast function (based on
whether we are under Scenario (S1) of Section 2.2.1 or Scenario (S2) of
Section 2.2.2) value for the point $b$, when we are in the interval $[s,e]$,
for the univariate component data sequence $\bm{X_{i}}$. The $d$ contrast
function values for each $\hat{r}_{m}$ are collected in the sets

$S_{m}=\left\\{CS^{(1)}(\hat{r}_{m}),\ldots,CS^{(d)}(\hat{r}_{m})\right\\},\quad
j=1,\ldots,M.$ (3.12)

Step 3: For each $m=1,\ldots,M$, all the elements of $S_{m}$ are tested
against the relevant threshold value, $\zeta_{T}$, for univariate change-point
detection which is of the order $\mathcal{O}(\sqrt{\log T})$ in both scenarios
of piecewise-constant and continuous piecewise-linear signals (representing
scenarios (S1) and (S2), respectively, covered in this paper). In terms of the
threshold constants, we employ those of Anastasiou and Fryzlewicz (2022). We
denote by
$\hat{sp}_{m}:=\frac{\\#\left\\{i:CS^{(i)}(\hat{r}_{m})>\zeta_{T}\right\\}}{d}$.
The estimated sparsity is then
$\hat{sp}=\max_{m\in\left\\{1,\ldots,M\right\\}}\left\\{\hat{sp}_{m}\right\\}.$
Step 4: For cases where $\hat{sp}\leq 0.4$, we accept the result from Step 1,
where MID has been paired with the $L_{\infty}$ norm, whereas if $\hat{sp}\geq
0.6$, then MID is paired with the $L_{2}$ norm. Extensive simulations have
shown that there is no significant difference on the MID’s practical
performance with respect to accuracy (on both the estimated number of change-
points and the estimated locations) when $\hat{sp}\in(0.4,0.6)$. Therefore,
MID could be paired with any of the aforementioned two norms and give very
good results. For computational cost reasons, in such cases we accept the
result of Step 1. From now on, we denote by ${\rm MID}_{\rm opt}$ to be the
data-adaptive, sparsity-based MID version explained in this section, where an
aggregation method is chosen based on the estimated sparsity of the change-
points in the given data.

###  3.4 A permutation-based approach

MID is a threshold-based algorithm because at each step, the largest
aggregated contrast function value is tested against a predefined threshold in
order to decide whether there is a change-point at the corresponding location.
In Section 3.2 we have explained how the threshold constant is carefully
chosen to control the Type I error taking into account the dimensionality of
the data. However, misspecification of the threshold can possibly lead to the
misestimation of the number of change-points. To solve such issues, we propose
a variant of MID based on permutation.

The idea of a “data-adaptive threshold” through permutations or bootstrapping
is not new. In Cho (2016), a bootstrap procedure is proposed, which is
motivated by the representation theory developed for the Generalised Dynamic
Factor Model, in order to approximate the quantiles of the developed double
CUSUM test statistic under the null hypothesis of no change-points; the
obtained quantiles are used as a test criterion for detection. In Cabrieto et
al. (2018), a permutation-based approach is used to test the presence of
correlation changes in multivariate time series. Under a univariate framework,
in Antoch and Hušková (2001) a permutation scheme is proposed for deriving
critical values for test statistics based on functionals of the partial sums
$\sum_{i=1}^{k}(X_{i}-\bar{X}_{n}),k=1,\ldots,n$, where $X_{i}$ are the
observed data and $\bar{X}_{n}$ their mean value. The proposed scheme consists
of three steps: a) compute the test statistic using the original data, b)
construct the permutation distribution by computing the relevant test
statistic on permuted versions of the data, and c) reject the null hypothesis
of no change-points if the test statistic lies in the tails of the
distribution.

We propose a variant that combines the isolation technique of MID with an
extension of the permutation procedure used in Antoch and Hušková (2001) to
the multivariate framework. Although this permutation procedure tends to be
computationally expensive, it has a straightforward implementation. As
generally holds for permutation procedures, the test-statistic obtained from
the original data is compared to those obtained when applying the same steps
to several permuted versions of the data. For our proposed permutation scheme
all steps remain the same as MID with the only difference being the way the
algorithm chooses to accept or reject a change-point within each interval. To
be more precise, suppose that, for a given data sequence
$\left\\{\bm{X_{t}}\right\\}_{t=1,\ldots,T}$, the MID algorithm is at a step
where it looks for a change-point in the interval $I=[s^{*},e^{*}]$, where
$1\leq s^{*}<e^{*}\leq T$. As described in Section 2.1, MID returns a vector
$\bm{v}\in\mathbb{R}^{J}$, where $J$ is the amount of all change-point
candidates. The elements of $\bm{v}$ correspond to the aggregated contrast
function values for each candidate point in $I$. The next step is to store
$T_{I_{max}}=\underset{i\in\\{1,\ldots,J\\}}{\max}\left\\{v_{i}\right\\}$ and
repeat the following procedure a prespecified amount, denoted by $K$, of
times:

  1. 1.

Generate a random permutation from $(s^{*},\ldots,e^{*})$.

  2. 2.

Reorder the data according to the permutation.

  3. 3.

Calculate and store the maximum aggregated contrast function value for each
permutation.

The empirical distribution obtained from the maximum values is used to
construct our test. More precisely, we identify $\hat{b}_{I}={\rm
argmax}_{t}\\{v_{t}\\}$ as a change-point if, for given $\alpha\in(0,1)$,
$T_{I_{max}}>q_{1-\alpha}$, where $q_{1-\alpha}$ is the $100(1-\alpha)\%$
quantile of the empirical distribution.

The parameter $\alpha$ controls the probability of false detections. Small
values of $\alpha$ will make it harder for the algorithm to reject the null
hypothesis of no change-points, whereas large values can reduce the
probability of a Type II error. For the simulations in Section 4, we take
$\alpha=0.01$. Regarding the parameter $K$, the simulation results provided in
Antoch and Hušková (2001)–although for a univariate framework–suggest that the
empirical quantiles get stabilised quickly. Therefore, for our simulations, we
take $K=1000$.

##  4 Simulations

In this section, we investigate the performance of our method in Scenarios
(S1) and (S2) covered in Sections 2.2.1 and 2.2.2, respectively. Furthermore,
in (S1) MID is compared with state-of-the-art methods through a comprehensive
simulation study. The competitors are the Double Cusum (DC) method of Cho
(2016), the Sparsified Binary Segmentation (SBS) algorithm of Cho and
Fryzlewicz (2015), and the INSPECT algorithm of Wang and Samworth (2018). DC
and SBS are implemented in the hdbinseg R package, while for INSPECT we have
used the InspectChangepoint R package. For DC, the parameters used were the
ones that gave the best results in the simulation study carried out in the
relevant paper (Cho (2016)), while SBS and INSPECT were called with their
default arguments. Regarding our algorithm, we give results for the data-
adaptive, sparsity-based MID version of Section 3.3, denoted by ${\rm
MID}_{\rm opt}$, as well as for the permutation-based variants in Section 3.4;
these are denoted by ${\rm MIDPERL}_{2}$ and ${\rm MIDPERL}_{\infty}$ for the
$L_{2}$ and $L_{\infty}$ norms, respectively.

Table 2: Distribution of $\hat{N}-N$ over 100 simulated multivariate data sequences with 3 change-points under Scenario (S1) of Section 2.2.1. The average ARI, $d_{H}$, and computational times are also given Method | $sp$ | $d$ | $\hat{N}-N$ | ARI | $d_{h}$ | Time (s)
---|---|---|---|---|---|---
|  |  | $-2$ | $-1$ | 0 | 1 | 2 | 3 |  |  |
$\bm{{\rm MID}_{\rm opt}}$ | 0.2 | 30 | 0 | 0 | 95 | 5 | 0 | 0 | 1 | 0.01 | 0.43
DC | 0.2 | 30 | 0 | 0 | 90 | 10 | 0 | 0 | 0.98 | 0.03 | 5.91
SBS | 0.2 | 30 | 1 | 36 | 62 | 1 | 0 | 0 | 0.96 | 0.15 | 1.08
INSPECT | 0.2 | 30 | 0 | 0 | 96 | 4 | 0 | 0 | 1 | 0.01 | 0.31
${\rm MIDPERL}_{2}$ | 0.2 | 30 | 0 | 2 | 82 | 14 | 0 | 2 | 0.96 | 0.05 | 330.81
${\rm MIDPERL}_{\infty}$ | 0.2 | 30 | 0 | 0 | 79 | 18 | 2 | 1 | 0.96 | 0.05 | 172.72
$\bm{{\rm MID}_{\rm opt}}$ | 0.2 | 100 | 0 | 1 | 90 | 8 | 1 | 0 | 1 | 0.03 | 1.41
DC | 0.2 | 100 | 0 | 1 | 94 | 5 | 0 | 0 | 0.99 | 0.01 | 13.11
SBS | 0.2 | 100 | 5 | 25 | 70 | 0 | 0 | 0 | 0.97 | 0.14 | 2.31
INSPECT | 0.2 | 100 | 0 | 2 | 90 | 7 | 1 | 0 | 1 | 0.02 | 1.38
${\rm MIDPERL}_{2}$ | 0.2 | 100 | 0 | 4 | 81 | 13 | 2 | 0 | 0.97 | 0.04 | 1098.40
${\rm MIDPERL}_{\infty}$ | 0.2 | 100 | 0 | 1 | 79 | 16 | 4 | 0 | 0.95 | 0.07 | 615.85
$\bm{{\rm MID}_{\rm opt}}$ | 0.5 | 30 | 0 | 1 | 93 | 6 | 0 | 0 | 1 | 0.02 | 0.43
DC | 0.5 | 30 | 0 | 0 | 92 | 7 | 1 | 0 | 0.99 | 0.01 | 5.94
SBS | 0.5 | 30 | 2 | 23 | 75 | 0 | 0 | 0 | 0.98 | 0.10 | 1.08
INSPECT | 0.5 | 30 | 0 | 1 | 90 | 9 | 0 | 0 | 1 | 0.01 | 0.38
${\rm MIDPERL}_{2}$ | 0.5 | 30 | 0 | 2 | 80 | 15 | 2 | 1 | 0.97 | 0.06 | 331.91
${\rm MIDPERL}_{\infty}$ | 0.5 | 30 | 0 | 1 | 74 | 20 | 5 | 0 | 0.96 | 0.05 | 161.49
${\rm MID}_{\rm opt}$ | 0.5 | 100 | 0 | 0 | 83 | 16 | 1 | 0 | 0.98 | 0.02 | 1.38
DC | 0.5 | 100 | 0 | 0 | 90 | 9 | 1 | 0 | 0.99 | 0.02 | 13.12
SBS | 0.5 | 100 | 0 | 27 | 73 | 0 | 0 | 0 | 0.98 | 0.10 | 2.31
INSPECT | 0.5 | 100 | 0 | 0 | 93 | 7 | 0 | 0 | 0.99 | 0.02 | 1.31
$\bm{{\rm MIDPERL}_{2}}$ | 0.5 | 100 | 0 | 1 | 88 | 10 | 1 | 0 | 0.98 | 0.03 | 1069.54
$\bm{{\rm MIDPERL}_{\infty}}$ | 0.5 | 100 | 0 | 0 | 88 | 10 | 2 | 0 | 0.97 | 0.04 | 751.34
$\bm{{\rm MID}_{\rm opt}}$ | 0.8 | 30 | 0 | 1 | 97 | 2 | 0 | 0 | 1 | 0.01 | 0.44
DC | 0.8 | 30 | 0 | 1 | 90 | 9 | 0 | 0 | 0.99 | 0.02 | 5.94
SBS | 0.8 | 30 | 0 | 32 | 68 | 0 | 0 | 0 | 0.98 | 0.09 | 1.08
INSPECT | 0.8 | 30 | 0 | 1 | 94 | 4 | 1 | 0 | 1 | 0.01 | 0.31
${\rm MIDPERL}_{2}$ | 0.8 | 30 | 0 | 1 | 89 | 6 | 4 | 0 | 0.98 | 0.04 | 324.60
${\rm MIDPERL}_{\infty}$ | 0.8 | 30 | 0 | 0 | 80 | 17 | 2 | 1 | 0.97 | 0.07 | 161.49
$\bm{{\rm MID}_{\rm opt}}$ | 0.8 | 100 | 0 | 2 | 97 | 1 | 0 | 0 | 1 | 0.01 | 1.70
DC | 0.8 | 100 | 0 | 0 | 93 | 7 | 0 | 0 | 0.99 | 0.01 | 13.37
SBS | 0.8 | 100 | 4 | 24 | 72 | 0 | 0 | 0 | 0.97 | 0.14 | 2.35
INSPECT | 0.8 | 100 | 0 | 0 | 93 | 7 | 0 | 0 | 0.99 | 0.01 | 1.30
${\bm{{\rm MIDPERL}_{2}}}$ | 0.8 | 100 | 0 | 1 | 93 | 4 | 1 | 1 | 0.98 | 0.02 | 1085.59
${\rm MIDPERL}_{\infty}$ | 0.8 | 100 | 0 | 0 | 85 | 12 | 3 | 0 | 0.98 | 0.03 | 602.87

Simulation setup: We considered the settings where $T=1500$,
$d\in\left\\{30,100\right\\}$ and $N\in\left\\{3,20,50\right\\}$. With the
definition of sparsity as in (3.11), we took
$sp\in\left\\{0.2,0.5,0.8\right\\}$. In total, we tested the methods in 18
different setups covering a wide range of multivariate sequences. The change
magnitude at the component data sequences (either in the mean or in the slope
depending on whether we are under (S1) or (S2), respectively) follows the
uniform distribution ${\rm U}(1,2)$. Standard Gaussian noise was added to the
signals. We ran 100 replications for each setup. The frequency distribution of
$\hat{N}-N$ is provided, while the accuracy of the estimated locations is
evaluated through the Adjusted Rand Index (ARI) of the estimated segmentation
against the true one (Hubert and Arabie (1985)), and the scaled Hausdorff
distance,

$d_{H}=n_{s}^{-1}\max\bigg{\\{}\max_{j}\min_{k}\mid
r_{j}-\hat{r_{k}}\mid,\max_{k}\min_{k}\mid r_{j}-\hat{r_{k}}\mid,\bigg{\\}},$

where $n_{s}$ is the length of the largest segment. The average computational
times, in seconds, are also provided. The results, for MID and the competing
methods under Scenario (S1) are given in Tables 2-4. The method with the
highest empirical frequency of $\hat{N}-N$ being equal to zero (or close to
zero, depending on the example) and those within $5\%$ off the highest are
given in bold. For Scenario (S2), the results are presented in Table 5.

Table 3: Distribution of $\hat{N}-N$ over 100 simulated multivariate data sequences with 20 change-points under Scenario (S1) of Section 2.2.1. The average ARI, $d_{H}$, and computational times are also given Method | $sp$ | $d$ | $\hat{N}-N$ | ARI | $d_{h}$ | Time (s)
---|---|---|---|---|---|---
|  |  | $\leq-10$ | $(-10,-2)$ | $[-2,2]$ | $(2,10]$ |  |  |
$\bm{{\rm MID}_{\rm opt}}$ | 0.2 | 30 | 0 | 4 | 96 | 0 | 0.95 | 0.19 | 0.72
DC | 0.2 | 30 | 2 | 98 | 0 | 0 | 0.74 | 0.59 | 7.38
SBS | 0.2 | 30 | 100 | 0 | 0 | 0 | 0.50 | 1.19 | 1.36
INSPECT | 0.2 | 30 | 0 | 8 | 87 | 5 | 0.94 | 0.23 | 0.44
${\rm MIDPERL}_{2}$ | 0.2 | 30 | 0 | 16 | 84 | 0 | 0.94 | 0.25 | 476.36
$\bm{{\rm MIDPERL}_{\infty}}$ | 0.2 | 30 | 0 | 4 | 96 | 0 | 0.95 | 0.17 | 110.60
$\bm{{\rm MID}_{\rm opt}}$ | 0.2 | 100 | 0 | 1 | 98 | 1 | 0.97 | 0.08 | 2.20
DC | 0.2 | 100 | 16 | 84 | 0 | 0 | 0.76 | 0.58 | 18.36
SBS | 0.2 | 100 | 100 | 0 | 0 | 0 | 0.57 | 0.97 | 3.17
INSPECT | 0.2 | 100 | 0 | 3 | 91 | 6 | 0.97 | 0.20 | 1.86
$\bm{{\rm MIDPERL}_{2}}$ | 0.2 | 100 | 0 | 6 | 94 | 0 | 0.97 | 0.21 | 1488.64
$\bm{{\rm MIDPERL}_{\infty}}$ | 0.2 | 100 | 0 | 1 | 99 | 0 | 0.97 | 0.13 | 427.60
$\bm{{\rm MID}_{\rm opt}}$ | 0.5 | 30 | 0 | 0 | 100 | 0 | 0.97 | 0.11 | 0.70
DC | 0.5 | 30 | 17 | 83 | 0 | 0 | 0.75 | 0.62 | 10.11
SBS | 0.5 | 30 | 100 | 0 | 0 | 0 | 0.45 | 1.41 | 1.34
INSPECT | 0.5 | 30 | 0 | 1 | 94 | 5 | 0.97 | 0.16 | 0.44
$\bm{{\rm MIDPERL}_{2}}$ | 0.5 | 30 | 0 | 1 | 99 | 0 | 0.98 | 0.14 | 461.75
$\bm{{\rm MIDPERL}_{\infty}}$ | 0.5 | 30 | 0 | 0 | 99 | 1 | 0.97 | 0.12 | 103.31
$\bm{{\rm MID}_{\rm opt}}$ | 0.5 | 100 | 0 | 0 | 100 | 0 | 0.98 | 0.06 | 2.15
DC | 0.5 | 100 | 22 | 78 | 0 | 0 | 0.76 | 0.58 | 21.03
SBS | 0.5 | 100 | 100 | 0 | 0 | 0 | 0.49 | 1.20 | 3.21
INSPECT | 0.5 | 100 | 0 | 0 | 95 | 5 | 0.98 | 0.11 | 0.99
$\bm{{\rm MIDPERL}_{2}}$ | 0.5 | 100 | 0 | 0 | 100 | 0 | 0.98 | 0.12 | 1471.95
$\bm{{\rm MIDPERL}_{\infty}}$ | 0.5 | 100 | 0 | 0 | 100 | 0 | 0.97 | 0.09 | 450.34
$\bm{{\rm MID}_{\rm opt}}$ | 0.8 | 30 | 0 | 2 | 98 | 0 | 0.98 | 0.16 | 0.72
DC | 0.8 | 30 | 24 | 76 | 0 | 0 | 0.74 | 0.62 | 11.74
SBS | 0.8 | 30 | 100 | 0 | 0 | 0 | 0.41 | 1.53 | 1.34
INSPECT | 0.8 | 30 | 0 | 0 | 97 | 3 | 0.98 | 0.10 | 0.44
$\bm{{\rm MIDPERL}_{2}}$ | 0.8 | 30 | 0 | 1 | 99 | 0 | 0.98 | 0.14 | 417.58
$\bm{{\rm MIDPERL}_{\infty}}$ | 0.8 | 30 | 0 | 0 | 100 | 0 | 0.97 | 0.12 | 95.89
$\bm{{\rm MID}_{\rm opt}}$ | 0.8 | 100 | 0 | 2 | 98 | 0 | 0.98 | 0.13 | 3.66
DC | 0.8 | 100 | 16 | 84 | 0 | 0 | 0.77 | 0.57 | 21.20
SBS | 0.8 | 100 | 100 | 0 | 0 | 0 | 0.46 | 1.29 | 3.18
INSPECT | 0.8 | 100 | 0 | 0 | 96 | 4 | 0.99 | 0.08 | 1.86
$\bm{{\rm MIDPERL}_{2}}$ | 0.8 | 100 | 0 | 0 | 100 | 0 | 0.99 | 0.09 | 1467.22
$\bm{{\rm MIDPERL}_{\infty}}$ | 0.8 | 100 | 0 | 0 | 100 | 0 | 0.98 | 0.07 | 459.46
Table 4: Distribution of $\hat{N}-N$ over 100 simulated multivariate data sequences with 50 change-points under Scenario (S1) of Section 2.2.1. The average ARI, $d_{H}$, and computational times are also given Method | $sp$ | $d$ | $\hat{N}-N$ | ARI | $d_{h}$ | Time (s)
---|---|---|---|---|---|---
|  |  | $\leq-40$ | $(-40,-20)$ | $[-20,-10)$ | $[-10,10]$ |  |  |
$\bm{{\rm MID}_{\rm opt}}$ | 0.2 | 30 | 0 | 0 | 1 | 99 | 0.88 | 0.52 | 1.48
DC | 0.2 | 30 | 1 | 94 | 5 | 0 | 0.56 | 1.29 | 10.41
SBS | 0.2 | 30 | 100 | 0 | 0 | 0 | 0.19 | 4.34 | 1.91
INSPECT | 0.2 | 30 | 0 | 0 | 0 | 100 | 0.88 | 0.53 | 0.76
$\bm{{\rm MIDPERL}_{2}}$ | 0.2 | 30 | 0 | 0 | 3 | 97 | 0.88 | 0.50 | 650.05
$\bm{{\rm MIDPERL}_{\infty}}$ | 0.2 | 30 | 0 | 0 | 3 | 97 | 0.88 | 0.50 | 101.91
$\bm{{\rm MID}_{\rm opt}}$ | 0.2 | 100 | 0 | 0 | 0 | 100 | 0.93 | 0.40 | 4.49
DC | 0.2 | 100 | 0 | 72 | 28 | 0 | 0.59 | 1.13 | 22.60
SBS | 0.2 | 100 | 98 | 2 | 0 | 0 | 0.23 | 3.65 | 4.62
INSPECT | 0.2 | 100 | 0 | 0 | 0 | 100 | 0.93 | 0.51 | 2.87
$\bm{{\rm MIDPERL}_{2}}$ | 0.2 | 100 | 0 | 0 | 1 | 99 | 0.92 | 0.54 | 2188.56
$\bm{{\rm MIDPERL}_{\infty}}$ | 0.2 | 100 | 0 | 0 | 0 | 100 | 0.91 | 0.46 | 530.94
$\bm{{\rm MID}_{\rm opt}}$ | 0.5 | 30 | 0 | 0 | 0 | 100 | 0.92 | 0.44 | 1.44
DC | 0.5 | 30 | 1 | 79 | 20 | 0 | 0.56 | 1.26 | 11.62
SBS | 0.5 | 30 | 100 | 0 | 0 | 0 | 0.16 | 5.02 | 1.85
INSPECT | 0.5 | 30 | 0 | 0 | 0 | 100 | 0.94 | 0.43 | 0.76
$\bm{{\rm MIDPERL}_{2}}$ | 0.5 | 30 | 0 | 0 | 0 | 100 | 0.94 | 0.48 | 576.55
$\bm{{\rm MIDPERL}_{\infty}}$ | 0.5 | 30 | 0 | 0 | 0 | 100 | 0.92 | 0.43 | 85.13
$\bm{{\rm MID}_{\rm opt}}$ | 0.5 | 100 | 0 | 0 | 0 | 100 | 0.95 | 0.27 | 4.38
DC | 0.5 | 100 | 0 | 87 | 13 | 0 | 0.59 | 1.07 | 28.03
SBS | 0.5 | 100 | 100 | 0 | 0 | 0 | 0.18 | 4.49 | 4.46
INSPECT | 0.5 | 100 | 0 | 0 | 0 | 100 | 0.97 | 0.36 | 2.86
$\bm{{\rm MIDPERL}_{2}}$ | 0.5 | 100 | 0 | 0 | 0 | 100 | 0.96 | 0.42 | 2147.27
$\bm{{\rm MIDPERL}_{\infty}}$ | 0.5 | 100 | 0 | 0 | 0 | 100 | 0.93 | 0.39 | 524.08
$\bm{{\rm MID}_{\rm opt}}$ | 0.8 | 30 | 0 | 0 | 0 | 100 | 0.95 | 0.48 | 2.13
DC | 0.8 | 30 | 0 | 83 | 17 | 0 | 0.58 | 1.12 | 13.46
SBS | 0.8 | 30 | 100 | 0 | 0 | 0 | 0.15 | 5.44 | 1.85
INSPECT | 0.8 | 30 | 0 | 0 | 0 | 100 | 0.96 | 0.37 | 0.76
$\bm{{\rm MIDPERL}_{2}}$ | 0.8 | 100 | 0 | 0 | 0 | 100 | 0.96 | 0.42 | 513.94
$\bm{{\rm MIDPERL}_{\infty}}$ | 0.8 | 100 | 0 | 0 | 0 | 100 | 0.93 | 0.41 | 77.26
$\bm{{\rm MID}_{\rm opt}}$ | 0.8 | 100 | 0 | 0 | 0 | 100 | 0.96 | 0.44 | 6.92
DC | 0.8 | 100 | 0 | 90 | 10 | 0 | 0.60 | 1.07 | 29.73
SBS | 0.8 | 100 | 100 | 0 | 0 | 0 | 0.17 | 5.02 | 4.41
INSPECT | 0.8 | 100 | 0 | 0 | 0 | 100 | 0.98 | 0.29 | 2.86
$\bm{{\rm MIDPERL}_{2}}$ | 0.8 | 100 | 0 | 0 | 0 | 100 | 0.98 | 0.36 | 1723.12
$\bm{{\rm MIDPERL}_{\infty}}$ | 0.8 | 100 | 0 | 0 | 0 | 100 | 0.94 | 0.37 | 511.36

As the tables show, MID performs extremely well in all setups in both (S1) and
(S2). More specifically, for (S1) our method is 17 out of 18 times either the
best method overall or within $5\%$ off the best method with respect to the
estimated number of change-points. In all cases, MID attains very high values
for the ARI, and quite small (in most cases it actually attains the smallest
value) ones for the scaled Hausdorff distance; such results justify that apart
from being extremely accurate in estimating the correct number of change-
points, MID is also very accurate regarding the estimated change-point
locations. Regarding the permutation-based variants of MID, both ${\rm
MIDPERL}_{2}$ and ${\rm MIDPERL}_{\infty}$ show a very good behaviour in all
different scenarios in terms of accuracy with respect to both the estimated
number and the estimated change-point locations. Although these permutation-
based variants do not require specification of the threshold, their
computational cost is quite large (see Tables 2-4).

In regards to the competing methods, INSPECT has a very accurate behaviour in
all different scenarios. DC’s performance is very good only in those cases
with three change-points and it seems to heavily underestimate in cases with
more, regularly occurring change-points. SBS does not have a good performance
in the scenarios tested; it underestimates the number of change-points.

With respect to Scenario (S2) as in Section 2.2.2, the results in Table 5
exhibit MID’s very good performance. In conclusion, taking into account its
low computational time, we can deduce that our proposed method is reliable and
quick in accurately detecting change-points under various different (with
respect to the number of change-points, their magnitude, the sparsity, the
dimensionality of the data, and the structure of the changes) multivariate
settings. R code and instructions to replicate the results in this section are
available on https://github.com/apapan08/Simulations-MID.

Table 5: Distribution of $\hat{N}-N$ over 100 simulated multivariate data sequences under Scenario (S2) of Section 2.2.2. We present the results of MID for different levels of sparsity and dimension for the data sequence. The number of change-point is equal to 3, 20, or 50. The average ARI, $d_{H}$, and computational times are also given $\hat{N}$ | $sp$ | $d$ | $\hat{N}-N$ | ARI | $d_{h}$ | Time (s)
---|---|---|---|---|---|---
|  |  | $-2$ | $-1$ | 0 | 1 | 2 | 3 |  |  |
3 | 0.2 | 10 | 0 | 0 | 92 | 7 | 1 | 0 | 0.987 | 0.011 | 0.293
3 | 0.5 | 10 | 0 | 0 | 92 | 6 | 2 | 0 | 0.984 | 0.034 | 0.251
3 | 0.8 | 10 | 0 | 0 | 94 | 6 | 0 | 0 | 0.991 | 0.008 | 0.223
3 | 0.2 | 30 | 0 | 0 | 96 | 4 | 0 | 0 | 0.990 | 0.010 | 0.668
3 | 0.5 | 30 | 0 | 0 | 96 | 4 | 0 | 0 | 0.987 | 0.013 | 0.688
3 | 0.8 | 30 | 0 | 0 | 99 | 1 | 0 | 0 | 0.997 | 0.003 | 0.647
3 | 0.2 | 100 | 0 | 0 | 96 | 3 | 1 | 0 | 0.989 | 0.015 | 2.139
3 | 0.5 | 100 | 0 | 0 | 97 | 3 | 0 | 0 | 0.989 | 0.008 | 2.091
3 | 0.8 | 100 | 0 | 0 | 99 | 1 | 0 | 0 | 0.999 | 0.001 | 2.061
20 | 0.2 | 10 | 0 | 0 | 87 | 11 | 2 | 0 | 0.957 | 0.045 | 0.367
20 | 0.5 | 10 | 0 | 0 | 85 | 14 | 1 | 0 | 0.960 | 0.044 | 0.450
20 | 0.8 | 10 | 0 | 0 | 98 | 1 | 1 | 0 | 0.976 | 0.028 | 0.459
20 | 0.2 | 30 | 0 | 0 | 92 | 7 | 1 | 0 | 0.962 | 0.041 | 1.080
20 | 0.5 | 30 | 0 | 0 | 82 | 16 | 2 | 0 | 0.963 | 0.046 | 1.048
20 | 0.8 | 30 | 0 | 0 | 94 | 6 | 0 | 0 | 0.985 | 0.021 | 0.983
20 | 0.2 | 100 | 0 | 0 | 90 | 9 | 1 | 0 | 0.964 | 0.036 | 3.775
20 | 0.5 | 100 | 0 | 0 | 85 | 14 | 1 | 0 | 0.964 | 0.036 | 3.495
20 | 0.8 | 100 | 0 | 3 | 97 | 0 | 0 | 0 | 0.994 | 0.189 | 3.335
50 | 0.2 | 10 | 0 | 3 | 86 | 10 | 1 | 0 | 0.915 | 0.136 | 0.860
50 | 0.5 | 10 | 0 | 1 | 79 | 17 | 3 | 0 | 0.921 | 0.123 | 0.748
50 | 0.8 | 10 | 0 | 7 | 86 | 7 | 0 | 0 | 0.951 | 0.106 | 0.712
50 | 0.2 | 30 | 0 | 0 | 96 | 3 | 1 | 0 | 0.925 | 0.101 | 2.165
50 | 0.5 | 30 | 0 | 0 | 84 | 14 | 1 | 1 | 0.928 | 0.103 | 2.118
50 | 0.8 | 30 | 0 | 11 | 88 | 1 | 0 | 0 | 0.970 | 0.094 | 2.035
50 | 0.2 | 100 | 0 | 0 | 85 | 15 | 0 | 0 | 0.931 | 0.095 | 6.879
50 | 0.5 | 100 | 0 | 0 | 82 | 17 | 1 | 0 | 0.934 | 0.090 | 6.716
50 | 0.8 | 100 | 2 | 12 | 86 | 0 | 0 | 0 | 0.987 | 0.086 | 6.784

##  5 Real data examples

###  5.1 UK House Price Index

In this section, the performance of our method is studied on monthly
percentage changes in the UK house price index for all property types from
January 2000 to January 2022 in twenty London Boroughs. The data are available
from http://landregistry.data.gov.uk/app/ukhpi and they were last accessed in
May 2022. We have a multivariate data sequence
$\bm{X_{t}}=\left(X_{t,1},\ldots,X_{t,20}\right)^{\intercal}$, with
$t=1,\ldots,265$. The boroughs used are : Barnet (Ba), Bexley (Be), Bromley
(Br), Camden (Ca), Croydon (Cr), Ealing (Ea), Enfield (En), Greenwich (Gr),
Hackney (Ha), Hammersmith and Fulham (HF), Harrow (Har), Islington (Is),
Kensington and Chelsea (KC), Lambeth (La), Merton (Me), Newham (Ne), Richmond
upon Thames (RuT), Sutton (Su), Tower Hamlets (TH), and Wandsworth (Wa).
Similar data have been investigated in Baranowski et al. (2019), Fryzlewicz
(2014), Fryzlewicz (2020), and Anastasiou and Fryzlewicz (2022); however under
a univariate setting. Figure 5 indicates the results when ${\rm MID}_{{\rm
opt}}$ of Section 3.3 was employed for the detection of changes under Scenario
(S1) as explained in Section 2.2; with respect to the threshold constant, we
used the optimal values as in Table 1 for $\alpha=0.05$.

Figure 5: The monthly percentage changes in the UK house price index for the
twenty London boroughs under consideration. The estimated change-point
locations when ${\rm MID}_{{\rm opt}}$ was employed can be seen with red,
vertical lines.

Our method identifies 27 change-points in the mean structure of the
multivariate data sequence at hand. We have analysed the same data using the
competing methods of Section 4. INSPECT detects 26 change-points at locations
very similar to those detected by our method, DC detects one change-point at
location 41, while SBS does not detect any change-points. We need to
highlight, though, that in INSPECT, multiple change-points are estimated using
a wild binary segmentation scheme, which, due to the randomness involved, does
not necessarily detect the same change-points when it is employed more than
once over the same data set.

###  5.2 The COVID-19 outbreak

The performance of our method is also investigated on data from the COVID-19
pandemic. In this case, we focus though on the detection of changes under
Scenario (S2) as described in Section 2.2. The data under consideration
consist of the daily number of new lab-confirmed COVID-19 cases in the four
constituent countries of United Kingdom; England, Northern Ireland, Scotland,
and Wales. The period under investigation is from 01/04/2020 until 01/04/2022.
The data are available from https://coronavirus.data.gov.uk and they were last
accessed on the $30^{{\rm th}}$ of April 2022. Based on the description, in
this example we have a multivariate data sequence
$\bm{X_{t}}=\left(X_{t,1},\ldots,X_{t,4}\right)^{\intercal}$, with
$t=1,\ldots,731$. Due to the fact that the data are positive integer numbers,
we perform the Anscombe transform, $a:\mathbb{N}\rightarrow\mathbb{R}$, with
$a(x)=2\sqrt{x+3/8}$, to each $X_{t,j}$; this transform brings the
distribution of the component data sequences closer to the Gaussian one with
constant variance.

Our method has detected 29 change-points in the vector of the first partial
derivatives (changes in the slope for the component univariate data sequences)
which seem to capture the important movements in the data.

Figure 6: The daily number of COVID-19 cases for England, Northern Ireland,
Scotland, and Wales. The estimated change-point locations in the trend of each
component data sequence when our proposed ${\rm MID}_{{\rm opt}}$ method was
employed can be seen with red, vertical lines.

##  6 Conclusion

In this paper, the MID methodology has been proposed for multiple generalized
change-point detection in multivariate, possibly high-dimensional, data
sequences. The algorithm partly constitutes a generalization of the Isolate-
Detect algorithm of Anastasiou and Fryzlewicz (2022), which has been developed
for univariate change-point detection. Mean-dominant norms, as those in (1),
are employed for the aggregation of the information across the different
components of the multivariate data sequence.

The aggregated values for the relevant contrast function (depending on the
structure of the changes) are collected and compared to a threshold value,
$\zeta_{T,d}$. The rate of $\zeta_{T,d}$ with respect to both the length, $T$,
of the data sequence and its dimensionality, $d$, has been proven to be
$\mathcal{O}\left(\sqrt{\log(Td^{1/4})}\right)$. The optimal multiplicative
threshold constant, $C$, so that $\zeta_{T,d}=C\sqrt{\log(Td^{1/4})}$ has been
derived through a large scale simulation study, with special attention in
controlling the Type I error rate, $\alpha$, of falsely detecting change-
points for various values of the dimensionality of the data sequence.

Misspecification of the threshold is possible and could lead to misestimation
of the underlying signal. To solve such issues, in Section 3.4, a permutation-
based variant of our MID algorithm has been proposed with no threshold choice
requirement. The algorithmic steps of MID remain the same; the difference lies
in the way the method chooses to accept or reject a change-point within the
interval under consideration, suppose this is denoted by $I=[s^{*},e^{*}]$,
where $1\leq s^{*}<e^{*}\leq T$. In Section 2.1, it has been explained in
detail that MID will first return a vector $\bm{v}\in\mathbb{R}^{J}$, where
$J$ is the amount of all change-point candidates. The elements of $\bm{v}$
correspond to the aggregated contrast function values for each candidate point
in $I$. The next step is to store the maximum value of $\bm{v}$ as obtained
from the original data and compare it to the appropriate quantiles of the
empirical distribution of the values obtained when applying the same steps to
several permuted versions of the data; see the exact steps in Section 3.4. The
proposed permutation-based variant of MID, though computationally expensive,
performs very well in terms of accuracy with respect to the estimated number
and locations of the change-points; see the results in Section 4.

The choice of the mean-dominant norm to be employed in the aggregation step of
MID has already been discussed in the current paper, more specifically in
Section 3.3. Our aim has been to provide a method that in practice would
require minimal parameter choice to the user; towards this purpose, a data-
adaptive variant, named ${\rm MID}_{{\rm opt}}$, of the method has been
constructed. We first estimate the sparsity level in the given multivariate
data; the steps are explained in detail in Section 3.3. Depending on the value
of the estimated sparsity, ${\rm MID}_{\rm opt}$ estimates the change-points
employing either the $L_{\infty}$ or the $L_{2}$ mean-dominant norm as defined
in (1).

Through simulated and real-life data examples presented in Sections 4 and 5,
respectively, it has been shown that MID has very good performance in terms of
accuracy and speed. Specifically, MID lies in the top 5% (in terms of the
accurate estimation of the number and the location of the change-points) of
the best methods when compared to the state-of-the-art competitors. In
addition, ${\rm MID}_{\rm opt}$ is a very quick detection method which in a
few seconds can analyse signals with length in the range of thousands and
dimensionality in the range of hundreds; this is carried out automatically
with minimal decision making from the user on the aggregation method to be
employed.
Acknowledgments
This research occurred whilst Angelos Papanastasiou was employed as a Research
Special Scientist at the University of Cyprus, supported by the start-up
funding that Andreas Anastasiou has received from the University of Cyprus.

## References

  * Anastasiou et al. (2022) Anastasiou, A., I. Cribben, and P. Fryzlewicz (2022).  Cross-covariance isolate detect: A new change-point method for estimating dynamic functional connectivity.  Medical Image Analysis 75.
  * Anastasiou and Fryzlewicz (2022) Anastasiou, A. and P. Fryzlewicz (2022).  Detecting multiple generalized change-points by isolating single ones.  Metrika 85, 141–174.
  * Antoch and Hušková (2001) Antoch, J. and M. Hušková (2001).  Permutation tests in change point analysis.  Statistics & probability letters 53(1), 37–46.
  * Auger and Lawrence (1989) Auger, I. E. and C. E. Lawrence (1989).  Algorithms for the Optimal Identification of Segment Neighborhoods.  Bulletin of Mathematical Biology 51, 39–54.
  * Baranowski et al. (2019) Baranowski, R., Y. Chen, and P. Fryzlewicz (2019).  Narrowest-over-threshold detection of multiple change points and change-point-like features.  Journal of the Royal Statistical Society: Series B 81, 649–672.
  * Bücher et al. (2014) Bücher, A., I. Kojadinovic, T. Rohmer, and J. Segers (2014).  Detecting changes in cross-sectional dependence in multivariate time series.  Journal of Multivariate Analysis 132, 111–128.
  * Cabrieto et al. (2018) Cabrieto, J., F. Tuerlinckx, P. Kuppens, B. Hunyadi, and E. Ceulemans (2018).  Testing for the presence of correlation changes in a multivariate time series: A permutation based approach.  Scientific reports 8(1), 1–20.
  * Carlstein (1988) Carlstein, E. (1988).  Nonparametric change-point estimation.  Annals of Statistics 16, 188–197.
  * Cho (2016) Cho, H. (2016).  Change-point detection in panel data via double cusum statistic.  Electronic Journal of Statistics 10(2), 2000–2038.
  * Cho and Fryzlewicz (2015) Cho, H. and P. Fryzlewicz (2015).  Multiple-change-point detection for high dimensional time series via sparsified binary segmentation.  Journal of the Royal Statistical Society: Series B: Statistical Methodology, 475–507.
  * Cho and Kirch (2020) Cho, H. and C. Kirch (2020).  Data segmentation algorithms: Univariate mean change and beyond.  arXiv preprint arXiv:2012.12814.
  * Enikeeva and Harchaoui (2019) Enikeeva, F. and Z. Harchaoui (2019).  High-dimensional change-point detection under sparse alternatives.  Annals of Statistics 47(4), 2051–2079.
  * Fearnhead et al. (2019) Fearnhead, P., R. Maidstone, and A. Letchford (2019).  Detecting Changes in Slope With an ${L}_{0}$ Penalty.  Journal of Computational and Graphical Statistics 28, 265–275.
  * Frick et al. (2014) Frick, K., A. Munk, and H. Sieling (2014).  Multiscale change point inference.  Journal of the Royal Statistical Society: Series B: Statistical Methodology, 495–580.
  * Fryzlewicz (2014) Fryzlewicz, P. (2014).  Wild binary segmentation for multiple change-point detection.  Annals of Statistics 42, 2243–2281.
  * Fryzlewicz (2020) Fryzlewicz, P. (2020).  Detecting possibly frequent change-points: Wild binary segmentation 2 and steepest-drop model selection.  Journal of the Korean Statistical Society 49, 1027–1070.
  * Groen et al. (2013) Groen, J. J., G. Kapetanios, and S. Price (2013).  Multivariate methods for monitoring structural change.  Journal of Applied Econometrics 28(2), 250–274.
  * Hampel (1974) Hampel, F. R. (1974).  The influence curve and its role in robust estimation.  Journal of the American Statistical Association 69, 383–393.
  * Hocking et al. (2013) Hocking, T. D., G. Schleiermacher, I. Janoueix-Lerosey, V. Boeva, J. Cappo, O. Delattre, F. Bach, and J.-P. Vert (2013).  Learning smoothing models of copy number profiles using breakpoint annotations.  BMC bioinformatics 14(1), 1–15.
  * Horváth and Hušková (2012) Horváth, L. and M. Hušková (2012).  Change-point detection in panel data.  Journal of Time Series Analysis 33(4), 631–648.
  * Hubert and Arabie (1985) Hubert, L. and P. Arabie (1985).  Comparing partitions.  Journal of classification 2(1), 193–218.
  * Inclán and Tiao (1994) Inclán, C. and G. C. Tiao (1994).  Use of cumulative sums of squares for retrospective detection of changes of variance.  Journal of the American Statistical Association 89, 913–923.
  * Jackson et al. (2005) Jackson, B., J. D. Sargle, D. Barnes, S. Arabhi, A. Alt, P. Gioumousis, E. Gwin, P. Sangtrakulcharoen, L. Tan, and T. T. Tsai (2005).  An Algorithm for Optimal Partitioning of Data on an Interval.  IEEE Signal Processing Letters 12, 105–108.
  * Jirak (2015) Jirak, M. (2015).  Uniform change point tests in high dimension.  Annals of Statistics 43(6), 2451–2483.
  * Killick et al. (2012) Killick, R., P. Fearnhead, and I. A. Eckley (2012).  Optimal detection of changepoints with a linear computational cost.  Journal of the American Statistical Association 107(500), 1590–1598.
  * Kovács et al. (2020) Kovács, S., H. Li, P. Bühlmann, and A. Munk (2020).  Seeded Binary Segmentation: A general methodology for fast and optimal change point detection.  arXiv preprint arXiv:2002.06633.
  * Lavielle and Teyssiere (2006) Lavielle, M. and G. Teyssiere (2006).  Detection of multiple change-points in multivariate time series.  Lithuanian Mathematical Journal 46(3), 287–306.
  * Lavielle and Teyssiere (2007) Lavielle, M. and G. Teyssiere (2007).  Adaptive detection of multiple change-points in asset price volatility.  In Long memory in economics, pp. 129–156. Springer.
  * Maeng and Fryzlewicz (2019) Maeng, H. and P. Fryzlewicz (2019).  Detecting linear trend changes and point anomalies in data sequences.  arXiv preprint arXiv:1906.01939, 2019.
  * Matteson and James (2014) Matteson, D. S. and N. A. James (2014).  A Nonparametric Approach for Multiple Change Point Analysis of Multivariate Data.  Journal of the American Statistical Association 109, 334–345.
  * Ombao et al. (2005) Ombao, H., R. Von Sachs, and W. Guo (2005).  Slex analysis of multivariate nonstationary time series.  Journal of the American Statistical Association 100(470), 519–531.
  * Picard et al. (2011) Picard, F., E. Lebarbier, M. Hoebeke, G. Rigaill, B. Thiam, and S. Robin (2011).  Joint segmentation, calling, and normalization of multiple cgh profiles.  Biostatistics 12(3), 413–428.
  * Rigaill (2015) Rigaill, G. (2015).  A pruned dynamic programming algorithm to recover the best segmentations with 1 to Kmax change-points.  Journal de la Société Française de Statistique 156, 180–205.
  * Schröder and Fryzlewicz (2013) Schröder, A. L. and P. Fryzlewicz (2013).  Adaptive trend estimation in financial time series via multiscale change-point-induced basis recovery.  Statistics and its interface 6(4), 449–461.
  * Siris and Papagalou (2004) Siris, V. A. and F. Papagalou (2004).  Application of anomaly detection algorithms for detecting syn flooding attacks.  In IEEE Global Telecommunications Conference, 2004. GLOBECOM’04., Volume 4, pp. 2050–2054. IEEE.
  * Tibshirani (2014) Tibshirani, R. J. (2014).  Adaptive piecewise polynomial estimation via trend filtering.  Annals of Statistics 42, 285–323.
  * Truong et al. (2020) Truong, C., L. Oudre, and N. Vayatis (2020).  Selective review of offline change point detection methods.  Signal Processing 167, 107299.
  * Vert and Bleakley (2010) Vert, J.-P. and K. Bleakley (2010).  Fast detection of multiple change-points shared by many signals using group lars.  In NIPS, pp. 2343–2351.
  * Wang and Samworth (2018) Wang, T. and R. J. Samworth (2018).  High-dimensional changepoint estimation via sparse projection.  Journal of the Royal Statistical Society: Series B 80, 57–83.
  * Yu (2020) Yu, Y. (2020).  A review on minimax rates in change point detection and localisation.  arXiv preprint arXiv:2011.01857.
  * Yuan and Lin (2006) Yuan, M. and Y. Lin (2006).  Model selection and estimation in regression with grouped variables.  Journal of the Royal Statistical Society: Series B (Statistical Methodology) 68(1), 49–67.
  * Zou et al. (2014) Zou, C., G. Yin, L. Feng, and Z. Wang (2014).  Nonparametric maximum likelihood approach to multiple change-point problems.  The Annals of Statistics 42, 970–1002.

##  Appendix A Proof of Theorem 2.1

With $\tilde{X}_{s,e}^{b,j}$ as in (2.4), we denote, for $s\leq b<e$, by
$\tilde{f}_{s,e}^{b,j}=\sqrt{\frac{e-b}{n(b-s+1)}}\sum_{t=s}^{b}f_{t,j}-\sqrt{\frac{b-s+1}{n(e-b)}}\sum_{t=b+1}^{e}f_{t,j}$
the value of the CUSUM statistic for the $j^{th}$ underlying component signal.
Furthermore,

$\displaystyle C^{b}_{s,e}=\max\limits_{1\leq j\leq
d}\left|\tilde{X}_{s,e}^{b,j}\right|,\qquad\quad
D_{s,e}^{b}=\max\limits_{1\leq j\leq d}\left|\tilde{f}_{s,e}^{b,j}\right|,$
$\displaystyle A_{T}^{*}=\bigg{\\{}\max_{1\leq s\leq b<e\leq
T}\left\\{\max_{1\leq j\leq
d}\left|\tilde{X}_{s,e}^{t,j}-\tilde{f}_{s,e}^{t,j}\right|\right\\}\leq\sqrt{8\operatorname{log}\left(Td^{\frac{1}{4}}\right)}\bigg{\\}},$
$\displaystyle A_{T}=\bigg{\\{}\max_{1\leq s\leq b<e\leq
T}\left|C^{b}_{s,e}-D_{s,e}^{b}\right|\leq\sqrt{8\operatorname{log}\left(Td^{\frac{1}{4}}\right)}\bigg{\\}}.$
(A.13)

###### Proof.

Step 1: We will show that

${\rm pr}\left(A_{T}\right)\geq 1-1/(12\sqrt{\pi}T).$

Due to the fact that

$\displaystyle\left|C^{b}_{s,e}-D_{s,e}^{b}\right|$
$\displaystyle=\bigg{|}\max_{1\leq j\leq
d}\left|\tilde{X}_{s,e}^{b,j}\right|-\max_{1\leq j\leq
d}\left|\tilde{f}_{s,e}^{b,j}\right|\bigg{|}$ $\displaystyle\leq\max_{1\leq
j\leq d}\left|\tilde{X}_{s,e}^{b,j}-\tilde{f}_{s,e}^{b,j}\right|,$

it is straightforward to see that ${\rm pr}\left(A_{T}\right)\geq{\rm
pr}\left(A_{T}^{*}\right)$ and therefore, it suffices to show that ${\rm
pr}\left((A^{*}_{T})^{c}\right)\leq 1/(12\sqrt{\pi}T)$. From the definition of
our model in (1.1), we have that for any $j\in\left\\{1,\ldots,d\right\\}$,
$\tilde{X}_{s,e}^{b,j}-\tilde{f}^{b,j}_{s,e}=\tilde{\epsilon}_{s,e}^{b,j}$,
and simple calculations lead to
$\tilde{\epsilon}_{s,e}^{t,j}\sim\mathcal{N}(0,\,1)$. Therefore, for
$Z_{j}\sim\mathcal{N}(0,\,1)$, the Bonferroni inequality yields

$\displaystyle{\rm pr}\left((A^{*}_{T})^{c}\right)$ $\displaystyle={\rm
pr}\bigg{(}\max_{1\leq s\leq b<e\leq T}\left\\{\max_{1\leq j\leq
d}|Z_{j}|\right\\}>\sqrt{8\operatorname{log}\left(Td^{\frac{1}{4}}\right)}\bigg{)}$
$\displaystyle\leq\sum_{1\leq s\leq b<e\leq T}{\rm pr}\bigg{(}\max_{1\leq
j\leq
d}|Z_{j}|>\sqrt{8\operatorname{log}\left(Td^{\frac{1}{4}}\right)}\bigg{)}$
$\displaystyle\leq\frac{T^{3}}{6}{\rm pr}\bigg{(}\max_{1\leq j\leq
d}|Z_{j}|>\sqrt{8\operatorname{log}\left(Td^{\frac{1}{4}}\right)}\bigg{)}$
$\displaystyle\leq\frac{T^{3}}{6}\sum_{j=1}^{d}{\rm
pr}\bigg{(}|Z_{j}|>\sqrt{8\operatorname{log}\left(Td^{\frac{1}{4}}\right)}\bigg{)}$
$\displaystyle=\frac{T^{3}}{3}\sum_{j=1}^{d}{\rm
pr}\bigg{(}Z_{j}>\sqrt{8\operatorname{log}\left(Td^{\frac{1}{4}}\right)}\bigg{)}$
$\displaystyle=\frac{dT^{3}}{3}\frac{\phi\bigg{(}\sqrt{8\operatorname{log}\left(Td^{\frac{1}{4}}\right)}\bigg{)}}{\sqrt{8\operatorname{log}\left(Td^{\frac{1}{4}}\right)}}=\frac{1}{3T\sqrt{2\pi}\sqrt{8\operatorname{log}(Td^{\frac{1}{4}})}}$
$\displaystyle\leq\frac{1}{12\sqrt{\pi}T},$ (A.14)

where $\phi(x),x\in\mathbb{R}$ is the probability density function of the
standard normal distribution. We conclude that ${\rm pr}(A_{T}^{*})\geq
1-1/(12\sqrt{\pi}T)$, which leads to ${\rm pr}(A_{T})\geq
1-1/(12\sqrt{\pi}T)$.
Step 2 : Let us denote by
$\bm{\psi_{s,e}^{b}}=(\psi_{s,e}^{b}(1),\ldots,\psi_{s,e}^{b}(T))$ the
contrast vector, where

$\psi_{s,e}^{b}(t)=\begin{cases}\sqrt{\frac{e-b}{n(b-s+1)}},&t=s,\ldots,b,\\\
-\sqrt{\frac{b-s+1}{n(e-b)}},&t=b+1,\ldots,e,\\\
0,&\text{otherwise},\end{cases},$

for $n=e-s+1$. With $\bm{f_{j}}=\left(f_{1,j},\ldots,f_{T,j}\right)$ and for
$[s,e]$ being any interval that contains only one true change-point, namely
$r_{j}$, allow us for ease of presentation to introduce the notation

$A_{s,e}^{b}(k,r_{j}):=\frac{\left\langle\bm{\psi_{s,e}^{b}}\left\langle\bm{f_{k}},\bm{\psi_{s,e}^{b}}\right\rangle-\bm{\psi_{s,e}^{r_{j}}}\left\langle\bm{f_{k}},\bm{\psi_{s,e}^{r_{j}}}\right\rangle,\bm{\epsilon}\right\rangle}{\left\|\bm{\psi_{s,e}^{b}}\left\langle\bm{f_{k}},\bm{\psi_{s,e}^{b}}\right\rangle-\bm{\psi_{s,e}^{r_{j}}}\left\langle\bm{f_{k}},\bm{\psi_{s,e}^{r_{j}}}\right\rangle\right\|_{2}}.$
(A.15)

In (A.15) above, $k\in\left\\{1,\ldots,d\right\\}$, while $\|\cdot\|_{2}$ is
the Euclidean norm. Due to $\bm{\epsilon_{t}},t=1,\ldots,T$ following the
$d$-variate standard normal distribution, then straightforward calculations
lead to

$A_{s,e}^{b}(k,r_{j})\sim\mathcal{N}(0,1).$ (A.16)

For

$B_{T}=\left\\{\max_{\begin{subarray}{c}j=1,\ldots,N\\\
k=1,\ldots,d\end{subarray}}\max_{\begin{subarray}{c}r_{j-1}<s\leq r_{j}\\\
r_{j}<e\leq r_{j+1}\\\ s\leq
b<e\end{subarray}}\left|A_{s,e}^{b}(k,r_{j})\right|\leq\sqrt{8\operatorname{log}\left(Td^{\frac{1}{4}}\right)}\right\\},$
(A.17)

we will show that ${\rm pr}\left(B_{T}\right)\geq 1-1/(12\sqrt{\pi}T)$. Using
(A.16) and the Bonferroni inequality, then, for $Z\sim\mathcal{N}(0,1)$,
similar steps as in (A) yield

${\rm pr}\left(B_{T}^{c}\right)\leq\frac{1}{12\sqrt{\pi}T}.$ (A.18)

From (A) and (A.18) in Steps 1 and 2, respectively, we conclude that ${\rm
pr}\left(A_{T}\cap B_{T}\right)\geq 1-\frac{1}{6\sqrt{\pi}T}.$
Step 3: This is the main part of our proof, where we explain in detail how to
get the result in (2.1). For ease of understanding, we split this step into
two smaller parts. From now on, we assume that $A_{T}^{*}$ (and therefore,
also, $A_{T}$) and $B_{T}$ both hold. The constants we use are

$C_{1}=\sqrt{C_{3}}+\sqrt{8},\;\;C_{2}=\frac{1}{\sqrt{6}}-\frac{2\sqrt{2}}{\underline{C}},\;\;C_{3}=2(2\sqrt{2}+4)^{2},$

where $\underline{C}$ is as in condition (A1).
Step 3.1: For ease of presentation, we take $\lambda_{T}\leq\delta_{T}/3$; for
more information on the general case of $\lambda_{T}\leq\delta_{T}/m$, for an
$m>1$, see Remark 1 in the online supplementary material of Anastasiou and
Fryzlewicz (2022). Allow us now $\forall j\in\left\\{1,\ldots,N\right\\}$, to
define the intervals

$\displaystyle
I_{j}^{R}=\left[r_{j}+\frac{\delta_{T}}{3},r_{j}+2\frac{\delta_{T}}{3}\right)$
$\displaystyle
I_{j}^{L}=\left(r_{j}-2\frac{\delta_{T}}{3},r_{j}-\frac{\delta_{T}}{3}\right].$
(A.19)

It is apparent that in order for $I_{j}^{R}$ and $I_{j}^{L}$ to have at least
one point, then we actually implicitly require that $\delta_{T}>3$, which is
the case for sufficiently large $T$. Since the length of the intervals in (A)
is equal to $\delta_{T}/3$ and $\lambda_{T}\leq\delta_{T}/3$, then MID ensures
that for $K=\left\lceil T/\lambda_{T}\right\rceil$ and
$k,m\in\left\\{1,\ldots,K\right\\}$, there exists at least one
$c_{k}^{r}=k\lambda_{T}$ and at least one $c_{m}^{l}=T-m\lambda_{T}+1$ that
are in $I_{j}^{R}$ and $I_{j}^{L}$, respectively, $\forall j=1,\ldots,N$.

At the beginning of our algorithm, $s=1$, $e=T$ (see, for example, Figure 2)
and depending on whether $r_{1}\leq T-r_{N}$ then $r_{1}$ or $r_{N}$ will get
isolated in a right- or left-expanding interval, respectively. W.l.o.g.,
assume that $r_{1}\leq T-r_{N}$ and, as already mentioned, our algorithm
(which can be seen as an extension of the ID methodology in Anastasiou and
Fryzlewicz (2022)) naturally ensures that $\exists
k\in\left\\{1,\ldots,K\right\\}$ such that $c_{k}^{r}\in I_{1}^{R}$. We
highlight that there is no other change-point in $[1,c_{k}^{r}]$ apart from
$r_{1}$. With $C_{s,e}^{b}$ and $D_{s,e}^{b}$ as in (A), we will show that for
$\tilde{b}={\rm argmax}_{1\leq t<c_{k}^{r}}C_{1,c_{k}^{r}}^{t}$, then
$C_{1,c_{k}^{r}}^{\tilde{b}}>\zeta_{T,d}$. Using (A), we have that

$\displaystyle C_{1,c_{k}^{r}}^{\tilde{b}}\geq C_{1,c_{k}^{r}}^{r_{1}}\geq
D_{1,c_{k}^{r}}^{r_{1}}-\sqrt{8\log\left(Td^{\frac{1}{4}}\right)}$
$\displaystyle=\sqrt{\frac{r_{1}(c_{k}^{r}-r_{1})}{c_{k}^{r}}}\max\limits_{1\leq
j\leq
d}\left\\{\Delta_{1}^{j}\right\\}-\sqrt{8\log\left(Td^{\frac{1}{4}}\right)}$
$\displaystyle\geq\sqrt{\frac{\min\left\\{c_{k}^{r}-r_{1},r_{1}\right\\}}{2}}\max\limits_{1\leq
j\leq
d}\left\\{\Delta_{1}^{j}\right\\}-\sqrt{8\log\left(Td^{\frac{1}{4}}\right)},$
(A.20)

with $\Delta_{j}^{q}$ defined in (2.2.1). For $\delta_{T}$ as in (2.2.1), we
notice that $r_{1}\geq\delta_{T}$. Furthermore, since $c_{k}^{r}\in I_{1}^{R}$
as in (A), then $c_{k}^{r}-r_{1}\geq\delta_{T}/3$, meaning that
$\min\left\\{c_{k}^{r}-r_{1},r_{1}\right\\}\geq\frac{\delta_{T}}{3}$.
Therefore, using assumption (A1),

$\displaystyle C_{1,c_{k}^{r}}^{\tilde{b}}$
$\displaystyle\geq\sqrt{\frac{\delta_{T}}{6}}\max\limits_{1\leq j\leq
d}\left\\{\Delta_{1}^{j}\right\\}-\sqrt{8\log\left(Td^{\frac{1}{4}}\right)}$
$\displaystyle\geq\sqrt{\frac{\delta_{T}}{6}}\underline{f}_{T}-\sqrt{8\log\left(Td^{\frac{1}{4}}\right)}$
$\displaystyle=\left(\frac{1}{\sqrt{6}}-\frac{\sqrt{8\log\left(Td^{1/4}\right)}}{\sqrt{\delta_{T}}\underline{f}_{T}}\right)\sqrt{\delta_{T}}\underline{f}_{T}$
$\displaystyle\geq\left(\frac{1}{\sqrt{6}}-\frac{\sqrt{8\log\left(Td^{1/4}\right)}}{\underline{C}\sqrt{\log\left(Td^{1/4}\right)}}\right)\sqrt{\delta_{T}}\underline{f}_{T}$
$\displaystyle=\left(\frac{1}{\sqrt{6}}-\frac{2\sqrt{2}}{\underline{C}}\right)\sqrt{\delta_{T}}\underline{f}_{T}=C_{2}\sqrt{\delta_{T}}\underline{f}_{T}>\zeta_{T,d}.$
(A.21)

Therefore, there will be an interval of the form $[1,c_{k}^{r}]$, with
$c_{k}^{r}>r_{1}$, such that $[1,c_{k}^{r}]$ contains only $r_{1}$ and
$\max_{1\leq b<c_{k}^{r}}\left\\{C_{1,c_{k}^{r}}^{b}\right\\}>\zeta_{T,d}$.
Let us, for $k^{*}\in\left\\{1,\ldots,K\right\\}$, to denote by
$c_{k^{*}}^{r}\leq c_{k}^{r}$ the first right-expanding point where this
happens and let $b_{1}={\rm argmax}_{1\leq
t<c_{k^{*}}^{r}}\left\\{C_{1,c_{k^{*}}^{r}}^{t}\right\\}$ with
$\left\\{C_{1,c_{k^{*}}^{r}}^{b_{1}}\right\\}>\zeta_{T,d}$. In addition, for
ease of presentation, allow us to denote by $q_{1}:={\rm
argmax}_{j=1,\ldots,d}\left|\tilde{X}_{1,c_{k^{*}}^{r}}^{b_{1},j}\right|$;
this is the univariate component data sequence where the contrast function
value got maximised at location $b_{1}$. We will find $\gamma_{T}>0$ such that
for any $b^{*}\in\left\\{1,\ldots,c_{k^{*}}^{r}-1\right\\}$ with
$\left|b^{*}-r_{1}\right|\left(\Delta_{1}^{q_{1}}\right)^{2}>\gamma_{T}$, it
holds that

$\left(\tilde{X}_{1,c_{k^{*}}^{r}}^{r_{1},q_{1}}\right)^{2}>\left(\tilde{X}_{1,c_{k^{*}}^{r}}^{b^{*},q_{1}}\right)^{2}.$
(A.22)

Proving (A.22) and using the definition of $b_{1}$ we can conclude that
$|b_{1}-r_{1}|\left(\Delta_{1}^{q_{1}}\right)^{2}\leq\gamma_{T}$. Now, since
in our model, $X_{t,j}=f_{t,j}+\epsilon_{t,j}$, then (A.22) can be expressed
as

$\displaystyle\left(\tilde{f}_{1,c_{k^{*}}^{r}}^{r_{1},q_{1}}\right)^{2}-\left(\tilde{f}_{1,c_{k^{*}}^{r}}^{b^{*},q_{1}}\right)^{2}$
$\displaystyle>\left(\tilde{\epsilon}_{1,c_{k^{*}}^{r}}^{b^{*},q_{1}}\right)^{2}-\left(\tilde{\epsilon}_{1,c_{k^{*}}^{r}}^{r_{1},q_{1}}\right)^{2}$
$\displaystyle\;+2\left\langle\bm{\psi_{1,c_{k^{*}}^{r}}^{b^{*}}}\langle\bm{f_{q_{1}}},\bm{\psi_{1,c_{k^{*}}^{r}}^{b^{*}}}\rangle-\bm{\psi_{1,c_{k^{*}}^{r}}^{r_{1}}}\langle\bm{f_{q_{1}}},\bm{\psi_{1,c_{k^{*}}^{r}}^{r_{1}}}\rangle,\bm{\epsilon_{q_{1}}}\right\rangle.$
(A.23)

W.l.o.g. assume that $b^{*}\geq r_{1}$ and a similar approach as below holds
when $b^{*}<r_{1}$. Using Lemma 4 in the online supplementary material of
Baranowski et al. (2019), gives for the left-hand side of the inequality in
(A) that

$\left(\tilde{f}_{1,c_{k^{*}}^{r}}^{r_{1},q_{1}}\right)^{2}-\left(\tilde{f}_{1,c_{k^{*}}^{r}}^{b^{*},q_{1}}\right)^{2}=\frac{\left|b^{*}-r_{1}\right|r_{1}}{\left|b^{*}-r_{1}\right|+r_{1}}\left(\Delta_{1}^{q_{1}}\right)^{2}:=\Lambda.$
(A.24)

For the terms on the right-hand side of (A), since we are working under
$A_{T}^{*}$ as in (A), we obtain that

$\displaystyle\left(\tilde{\epsilon}_{1,c_{k^{*}}^{r}}^{b^{*},q_{1}}\right)^{2}-\left(\tilde{\epsilon}_{1,c_{k^{*}}^{r}}^{r_{1},q_{1}}\right)^{2}$
$\displaystyle\leq\left(\max_{s\leq b<e}\left\\{\max_{1\leq j\leq
d}\left|\tilde{\epsilon}_{s,e}^{b,j}\right|\right\\}\right)^{2}$
$\displaystyle\qquad-\left(\tilde{\epsilon}_{1,c_{k^{*}}^{r}}^{r_{1},q_{1}}\right)^{2}$
$\displaystyle\leq\left(\max_{s\leq b<e}\left\\{\max_{1\leq j\leq
d}\left|\tilde{\epsilon}_{s,e}^{b,j}\right|\right\\}\right)^{2}$
$\displaystyle\leq 8\log\left(Td^{\frac{1}{4}}\right),$

while from (A.17) and Lemma 4 from the online supplementary material of
Baranowski et al. (2019),

$\displaystyle
2\left\langle\bm{\psi_{1,c_{k^{*}}^{r}}^{b^{*}}}\langle\bm{f_{q_{1}}},\bm{\psi_{1,c_{k^{*}}^{r}}^{b^{*}}}\rangle-\bm{\psi_{1,c_{k^{*}}^{r}}^{r_{1}}}\langle\bm{f_{q_{1}}},\bm{\psi_{1,c_{k^{*}}^{r}}^{r_{1}}}\rangle,\bm{\epsilon_{q_{1}}}\right\rangle$
$\displaystyle\leq
2\|\bm{\psi_{1,c_{k^{*}}^{r}}^{b^{*}}}<\bm{f_{q_{1}}},\bm{\psi_{1,c_{k^{*}}^{r}}^{b^{*}}}>-\bm{\psi_{1,c_{k^{*}}^{r}}^{r_{1}}}<\bm{f_{q_{1}}},\bm{\psi_{1,c_{k^{*}}^{r}}^{r_{1}}}>\|_{2}$
$\displaystyle\qquad\times\sqrt{8\log\left(Td^{\frac{1}{4}}\right)}$
$\displaystyle=2\sqrt{8\Lambda\log\left(Td^{\frac{1}{4}}\right)}.$

Therefore, (A) is satisfied if the stronger inequality
$\Lambda>8\log\left(Td^{1/4}\right)+2\sqrt{8\Lambda\log\left(Td^{1/4}\right)}$
is satisfied, which has solution

$\Lambda>(2\sqrt{2}+4)^{2}\log\left(Td^{\frac{1}{4}}\right).$

From (A.24) and since

$\frac{\left|b^{*}-r_{1}\right|r_{1}}{\left|b^{*}-r_{1}\right|+r_{1}}\geq\frac{1}{2}\min\left\\{\left|b^{*}-r_{1}\right|,r_{1}\right\\},$

we deduce that (A.22) is implied by

$\displaystyle\min\left\\{\left|b^{*}-r_{1}\right|,r_{1}\right\\}$
$\displaystyle>\frac{2(2\sqrt{2}+4)^{2}\log\left(Td^{\frac{1}{4}}\right)}{\left(\Delta_{1}^{q_{1}}\right)^{2}}$
$\displaystyle=\frac{C_{3}\log\left(Td^{\frac{1}{4}}\right)}{\left(\Delta_{1}^{q_{1}}\right)^{2}}.$
(A.25)

However,

$\min\left\\{r_{1},c_{k^{*}}^{r}-r_{1}\right\\}>C_{3}\frac{\log\left(Td^{\frac{1}{4}}\right)}{\left(\Delta_{1}^{q_{1}}\right)^{2}}$
(A.26)

and this is because if we assume that
$\min\left\\{r_{1},c_{k^{*}}^{r}-r_{1}\right\\}\leq
C_{3}\log\left(Td^{\frac{1}{4}}\right)/\left(\Delta_{1}^{q_{1}}\right)^{2}$,
then

$\displaystyle
C_{1,c_{k^{*}}^{r}}^{b_{1}}=\left|\tilde{X}_{1,c_{k^{*}}^{r}}^{b_{1},q_{1}}\right|=\left|\tilde{X}_{1,c_{k^{*}}^{r}}^{b_{1},q_{1}}-\tilde{f}_{1,c_{k^{*}}^{r}}^{b_{1},q_{1}}+\tilde{f}_{1,c_{k^{*}}^{r}}^{b_{1},q_{1}}\right|$
$\displaystyle\leq\sqrt{8\log\left(Td^{\frac{1}{4}}\right)}+\left|\tilde{f}_{1,c_{k^{*}}^{r}}^{b_{1},q_{1}}\right|\leq\sqrt{8\log\left(Td^{\frac{1}{4}}\right)}+\left|\tilde{f}_{1,c_{k^{*}}^{r}}^{r_{1},q_{1}}\right|$
$\displaystyle=\sqrt{8\log\left(Td^{\frac{1}{4}}\right)}+\sqrt{\frac{r_{1}(c_{k^{*}}^{r}-r_{1})}{c_{k^{*}}^{r}}}\Delta_{1}^{q_{1}}$
$\displaystyle\leq\sqrt{8\log\left(Td^{\frac{1}{4}}\right)}+\sqrt{\min\left\\{c_{k^{*}}^{r}-r_{1},r_{1}\right\\}}\Delta_{1}^{q_{1}}$
$\displaystyle\leq\left(\sqrt{C_{3}}+\sqrt{8}\right)\sqrt{\log\left(Td^{\frac{1}{4}}\right)}=C_{1}\sqrt{\log\left(Td^{\frac{1}{4}}\right)}$
$\displaystyle\leq\zeta_{T,d}.$

This comes to a contradiction to $C_{1,c_{k^{*}}^{r}}^{b_{1}}>\zeta_{T,d}$,
which has already been proven in (A). Therefore, (A.26) holds and (A) is
restricted to
$\left|b^{*}-r_{1}\right|\left(\Delta_{1}^{q_{1}}\right)^{2}>C_{3}\log T$,
which implies (A.22). Therefore, necessarily

$\left|b_{1}-r_{1}\right|\left(\Delta_{1}^{q_{1}}\right)^{2}\leq
C_{3}\log\left(Td^{\frac{1}{4}}\right).$ (A.27)

So far, for $\lambda_{T}\leq\delta_{T}/3$ we have proven that working under
the assumption that $A_{T}^{*}$ (which implies that $A_{T}$ also holds) and
$B_{T}$ hold, there will be an interval $[1,c_{k^{*}}^{r}]$, with
$C_{1,c_{k^{*}}^{r}}^{b_{1}}>\zeta_{T,d}$, where $b_{1}=\underset{1\leq
t<c_{k^{*}}^{r}}{{\rm argmax}}C_{1,c_{k^{*}}^{r}}^{t}$ is an estimation of
$r_{1}$ that satisfies (A.27).
Step 3.2: After detecting the first change-point, MID follows the same process
as in Step 3.1 but now in the set $[c_{k^{*}}^{r},T]$, which contains the
change-points $r_{2},\ldots,r_{N}$. We do not check for possible change-points
the interval $[b_{1}+1,c_{k^{*}}^{r})$. However, this does not create any
issues because:

  * M.1

There is no change-point in $[b_{1}+1,c_{k^{*}}^{r})$, apart from maybe the
already detected $r_{1}$;

  * M.2

$c_{k^{*}}^{r}$ is at a location which allows for detection of $r_{2}$.

For M.1: We will split the explanation into two cases with respect to the
location of $b_{1}$.
Case 1: $b_{1}<r_{1}<c_{k^{*}}^{r}$. Using (A.27) and with

$\delta_{T}>3C_{3}\frac{\log\left(Td^{\frac{1}{4}}\right)}{\left(\Delta_{1}^{q_{1}}\right)^{2}},$
(A.28)

then since $c_{k}^{r}\in I_{1}^{R}$, we have that

$\displaystyle c_{k^{*}}^{r}-b_{1}$ $\displaystyle\leq
c_{k}^{r}-b_{1}=c_{k}^{r}-r_{1}+r_{1}-b_{1}$
$\displaystyle<2\frac{\delta_{T}}{3}+r_{1}-b_{1}$ $\displaystyle\leq
2\frac{\delta_{T}}{3}+\frac{C_{3}\log\left(Td^{\frac{1}{4}}\right)}{\left(\Delta_{1}^{q_{1}}\right)^{2}}<\delta_{T}.$

Since $r_{2}-r_{1}\geq\delta_{T}$ and $r_{1}$ is already in
$[b_{1}+1,c_{k^{*}}^{r})$, then there is no other change-point in
$[b_{1}+1,c_{k^{*}}^{r})$ apart from $r_{1}$.
Case 2: $r_{1}\leq b_{1}<c_{k^{*}}^{r}$. By construction
$c_{k^{*}}^{r}-r_{1}<2\delta_{T}/3$, which means that apart from $r_{1}$ there
is no other change-point in $[r_{1},c_{k^{*}}^{r})$. With $r_{1}\leq b_{1}$,
then $[b_{1}+1,c_{k^{*}}^{r})$ does not have any change-point.

Cases 1 and 2 above show that no matter the location of $b_{1}$, there is no
change-point in $[b_{1}+1,c_{k^{*}}^{r})$ other than possibly the previously
detected $r_{1}$. Similarly to the approach in Step 3.1, our method applied
now in $[s,e]=[c_{k^{*}}^{r},T]$, will first isolate, and then detect, $r_{2}$
or $r_{N}$ depending on whether $r_{2}-c_{k^{*}}^{r}$ is smaller or larger
than $T-r_{N}$. If $T-r_{N}<r_{2}-c_{k^{*}}^{r}$ then $r_{N}$ will get
isolated first in a left-expanding interval and the procedure to show its
detection is exactly the same as for the detection of $r_{1}$ in Step 3.1.
Therefore, for the sake of showing M.2 let us assume that
$r_{2}-c_{k^{*}}^{r}\leq T-r_{N}$.
For M.2: By construction, there exists a right expanding point
$c_{k_{2}}^{r}\in I_{2}^{R}$, with $I_{j}^{R}$ defined in (A). We will show
that $r_{2}$ gets detected in $[c_{k^{*}}^{r},c_{k^{*}_{2}}^{r}]$, for
$k_{2}^{*}\leq k_{2}$ and its detection is $b_{2}={\rm
argmax}_{c_{k^{*}}^{r}\leq
t<c_{k^{*}_{2}}^{r}}C_{c_{k^{*}}^{r},c_{k^{*}_{2}}^{r}}^{t}$, which satisfies
$\left|b_{2}-r_{2}\right|\left(\Delta_{2}^{q_{2}}\right)^{2}\leq
C_{3}\log\left(Td^{1/4}\right)$, where $q_{2}:={\rm
argmax}_{j=1,2,\ldots,d}\left|\tilde{X}_{c_{k^{*}}^{r},c_{k^{*}_{2}}^{r}}^{b_{2},j}\right|$.
Following similar steps as in (A), we have that

$C_{c_{k^{*}}^{r},c_{k_{2}^{*}}^{r}}^{b_{2}}>\zeta_{T,d}.$

We will now show that
$\left|b_{2}-r_{2}\right|\left(\Delta_{2}^{q_{2}}\right)^{2}\leq
C_{3}\log\left(Td^{1/4}\right)$. Following exactly the same process as in Step
3.1 and assuming now w.l.o.g. that $b_{2}<r_{2}$, we have that for
$b^{*}\in\left\\{c_{k^{*}}^{r},\ldots,c_{k^{*}_{2}}^{r}-1\right\\}$,

$\left(\tilde{X}_{c_{k^{*}}^{r},c_{k^{*}_{2}}^{r}}^{r_{2},q_{2}}\right)^{2}>\left(\tilde{X}_{c_{k^{*}}^{r},c_{k^{*}_{2}}^{r}}^{b^{*},q_{2}}\right)^{2}$
(A.29)

is implied by
$\min\left\\{\left|b^{*}-r_{2}\right|,c_{k^{*}_{2}}^{r}-r_{2}\right\\}>C_{3}\log
T/\left(\Delta_{2}^{f}\right)^{2}$. In the same way as in Step 3.1 and by
contradiction we can show that
$\min\left\\{c_{k^{*}_{2}}^{r}-r_{2},r_{2}-c_{k^{*}}^{r}+1\right\\}>C_{3}\log
T/\left(\Delta_{2}^{f}\right)^{2}$ and (A.29) is implied by
$\left|b^{*}-r_{2}\right|\left(\Delta_{2}^{q_{2}}\right)^{2}>C_{3}\log\left(Td^{1/4}\right)$.
Therefore
$\left|b_{2}-r_{2}\right|\left(\Delta_{2}^{q_{2}}\right)^{2}>C_{3}\log\left(Td^{1/4}\right)$
would mean that
$\left|\tilde{X}_{c_{k^{*}}^{r},c_{k^{*}_{2}}^{r}}^{r_{2},q_{2}}\right|>\left|\tilde{X}_{c_{k^{*}}^{r},c_{k^{*}_{2}}^{r}}^{b_{2},r_{2}}\right|$,
which is not true by the definition of $b_{2}$. Having said this, we conclude
that $\left|b_{2}-r_{2}\right|\left(\Delta_{2}^{q_{2}}\right)^{2}\leq
C_{3}\log\left(Td^{1/4}\right)$. Having detected $r_{2}$, then our algorithm
will proceed in the interval $[s,e]=[c_{k^{*}_{2}}^{r},T]$ and all the change-
points will get detected one by one since Step 3.2 will be applicable as long
as there are undetected change-points in $[s,e]$.

Denoting by $\hat{r}_{j}$ the estimation of $r_{j}$ as we did in the statement
of the theorem and for $[s^{*},e^{*}]$ being the interval where the isolation
and detection of $r_{j}$ occurs (in the way that we explained in Steps 3.1 and
3.2) allow us to denote by $q_{j}:={\rm
argmax}_{m=1,2,\ldots,d}\left|\tilde{X}_{s^{*},e^{*}}^{\hat{r}_{j},m}\right|$.
Then, Steps 3.1 and 3.2 have shown that MID will detect all the change-points
one by one and
$\left|\hat{r}_{j}-r_{j}\right|\left(\Delta_{j}^{q_{j}}\right)^{2}\leq
C_{3}\log\left(Td^{1/4}\right),\quad\forall
j\in\left\\{1,2,\ldots,N\right\\}$.
Step 4: The arguments given in Steps 1-3 hold in $A_{T}^{*}\cap B_{T}$. At the
beginning of the algorithm, $s=1,e=T$ and for $N\geq 1$, there exist
$k_{1}\in\left\\{1,\ldots,K\right\\}$ such that $s_{k_{1}}=s,e_{k_{1}}\in
I_{1}^{R}$ and $k_{2}\in\left\\{1,\ldots,K\right\\}$ such that $s_{k_{2}}\in
I_{N}^{L},e_{k_{2}}=e$. As in our previous steps, w.l.o.g. assume that
$r_{1}\leq T-r_{N}$ and $r_{1}$ gets isolated and detected first in an
interval $[s,c_{k^{*}}^{r}]$, where $c_{k^{*}}^{r}$ is less than or equal to
$e_{k_{1}}$. Then, $\hat{r}_{1}={\rm argmax}_{s\leq
t<c_{k^{*}}^{r}}C_{s,c_{k^{*}}^{r}}^{t}$ is the estimated location for $r_{1}$
and $\left|r_{1}-\hat{r}_{1}\right|\left(\Delta_{1}^{q_{1}}\right)^{2}\leq
C_{3}\log\left(Td^{1/4}\right)$. After this, the method continues in
$[c_{k^{*}}^{r},T]$ and keeps detecting all the change-points as explained in
Step 3. There will not be any double detection issues because naturally, at
each step of the algorithm, the new interval $[s,e]$ does not include any
previously detected change-points. Once all the change-points have been
detected one by one, then $[s,e]$ will contain no other change-points. Our
method will of course keep checking for possible change-points in right- and
left-expanding intervals denoted by $[s^{*},e^{*}]$. However, MID will not
detect anything in $[s^{*},e^{*}]$ because $\forall b\in[s^{*},e^{*})$,

$\displaystyle C_{s^{*},e^{*}}^{b}$ $\displaystyle\leq
D_{s^{*},e^{*}}^{b}+\sqrt{8\log\left(Td^{\frac{1}{4}}\right)}=\sqrt{8\log\left(Td^{\frac{1}{4}}\right)}$
$\displaystyle<C_{1}\sqrt{\log\left(Td^{\frac{1}{4}}\right)}\leq\zeta_{T,d}.$

After not detecting anything in all intervals of the above form, then the
algorithm concludes that there are not any more change-points to be detected
and stops. ∎

##  Appendix B Proof of Theorem 2.2

Before proving the result, we introduce some useful notation with

$\displaystyle G^{b}_{s,e}$ $\displaystyle:=\max\limits_{1\leq j\leq
d}\left\\{C_{s,e}^{b}(\bm{X_{j}})\right\\}$ $\displaystyle H_{s,e}^{b}$
$\displaystyle:=\max\limits_{1\leq j\leq
d}\left\\{C_{s,e}^{b}(\bm{f_{j}})\right\\},$ (B.30)

where $C_{s,e}^{b}(\cdot)$ is as in (2.7) of the main paper. For ease of
presentation, we define the sets

$\displaystyle\tilde{A}_{T}^{*}$ $\displaystyle=\bigg{\\{}\max_{1\leq s\leq
b<e\leq T}\left\\{\max_{1\leq j\leq
d}\left|C_{s,e}^{b}(\bm{X_{j}})-C_{s,e}^{b}(\bm{f_{j}})\right|\right\\}\leq\sqrt{8\operatorname{log}\left(Td^{\frac{1}{4}}\right)}\bigg{\\}},$
$\displaystyle\tilde{A}_{T}$ $\displaystyle=\bigg{\\{}\max_{1\leq s\leq
b<e\leq
T}\left|G^{b}_{s,e}-H_{s,e}^{b}\right|\leq\sqrt{8\operatorname{log}\left(Td^{\frac{1}{4}}\right)}\bigg{\\}},$
(B.31) $\displaystyle\tilde{B}_{T}$
$\displaystyle=\left\\{\max_{\begin{subarray}{c}j=1,\ldots,N\\\
k=1,\ldots,d\end{subarray}}\max_{\begin{subarray}{c}r_{j-1}<s\leq r_{j}\\\
r_{j}<e\leq r_{j+1}\\\ s\leq
b<e\end{subarray}}\left|\tilde{A}_{s,e}^{b}(k,r_{j})\right|\leq\sqrt{8\operatorname{log}\left(Td^{\frac{1}{4}}\right)}\right\\},$

where for $\bm{\phi_{s,e}^{b}}$ as in (2.2.2) of the main paper,

$\tilde{A}_{s,e}^{b}(k,r_{j}):=\frac{\left\langle\bm{\phi_{s,e}^{b}}\left\langle\bm{f_{k}},\bm{\phi_{s,e}^{b}}\right\rangle-\bm{\phi_{s,e}^{r_{j}}}\left\langle\bm{f_{k}},\bm{\phi_{s,e}^{r_{j}}}\right\rangle,\bm{\epsilon}\right\rangle|}{\left\|\bm{\phi_{s,e}^{b}}\left\langle\bm{f_{k}},\bm{\phi_{s,e}^{b}}\right\rangle-\bm{\phi_{s,e}^{r_{j}}}\left\langle\bm{f_{k}},\bm{\phi_{s,e}^{r_{j}}}\right\rangle\right\|_{2}}.$
(B.32)

###### Proof.

We will prove the more specific result

${\rm
pr}\left(\hat{N}=N,\max_{j=1,\ldots,N}\left(\left|\hat{r}_{j}-r_{j}\right|\left(\Delta_{j}^{q_{j}}\right)^{\frac{2}{3}}\right)\leq
C_{3}(\log\left(Td^{1/4}\right))^{\frac{1}{3}}\right)\leq
1-\frac{1}{6\sqrt{\pi}T},$ (B.33)

which implies the result in (2.9) of the main paper.
Steps 1 and 2: Due to the fact that ${\rm
pr}\left(\tilde{A}_{T}\right)\geq{\rm pr}\left(\tilde{A}_{T}^{*}\right)$, then
arguments as those used in Steps 1 and 2 of the proof of Theorem 2.1 can be
employed here as well in order to show that

$\displaystyle{\rm pr}\left(\tilde{A}_{T}\right)\geq{\rm
pr}\left(\tilde{A}_{T}^{*}\right)\geq 1-\frac{1}{12\sqrt{\pi}T}$
$\displaystyle{\rm pr}\left(\tilde{B}_{T}\right)\geq
1-\frac{1}{12\sqrt{\pi}T}.$ (B.34)

Therefore, ${\rm pr}\left(\tilde{A}_{T}\cap\tilde{B}_{T}\right)\geq
1-1/(6\sqrt{\pi}T)$.
Step 3: This is the main part of our proof, where we explain in detail how to
get the result in (B.33). From now on, we assume that $\tilde{A}_{T}^{*}$
(therefore, $\tilde{A}_{T}$ as well) and $\tilde{B}_{T}$ both hold. The
constants we use are

$C_{1}=\sqrt{\frac{2}{3}}C_{3}^{\frac{3}{2}}+\sqrt{8},\;\;C_{2}=\frac{1}{3\sqrt{72}}-\frac{2\sqrt{2}}{C^{*}},\;\;C_{3}=63^{\frac{1}{3}}(2\sqrt{2}+4)^{\frac{2}{3}},$

where $C^{*}$ is as in assumption (A2) of the main paper.
Step 3.1: First, $\forall j\in\left\\{1,\ldots,N\right\\}$, we define
$I_{j}^{R}$ and $I_{j}^{L}$ as in (A). At the beginning of our algorithm,
$s=1$, $e=T$ and depending on whether $r_{1}\leq T-r_{N}$, then $r_{1}$ or
$r_{N}$ will get isolated first in a right- or left-expanding interval,
respectively. W.l.o.g., assume that $r_{1}\leq T-r_{N}$. Because of the
structure of the ID methodology developed in Anastasiou and Fryzlewicz (2022)
and partly employed in our paper, it is ensured that
$\exists\tilde{k}\in\left\\{1,\ldots,K\right\\}$ such that
$c_{\tilde{k}}^{r}=\tilde{k}\lambda_{T}\in I_{1}^{R}$ and there is no other
change-point in $[1,c_{\tilde{k}}^{r}]$ apart from $r_{1}$. With, our aim now
is to show that for $\tilde{b}_{1}={\rm
argmax}_{1<t<c_{\tilde{k}}^{r}}G_{1,c_{\tilde{k}}^{r}}^{t}$, then
$G_{1,c_{\tilde{k}}^{r}}^{\tilde{b}_{1}}>\zeta_{T,d}$. Using (B), we obtain
that

$\displaystyle G_{1,c_{\tilde{k}}^{r}}^{\tilde{b}_{1}}$ $\displaystyle\geq
G_{1,c_{\tilde{k}}^{r}}^{r_{1}}\geq
H_{1,c_{\tilde{k}}^{r}}^{r_{1}}-\sqrt{8\log\left(Td^{\frac{1}{4}}\right)}=\max\limits_{1\leq
j\leq
d}\left\\{C_{1,c_{\tilde{k}}^{r}}^{r_{1}}(\bm{f_{j}})\right\\}-\sqrt{8\log\left(Td^{\frac{1}{4}}\right)}$
$\displaystyle\geq\max\limits_{1\leq j\leq
d}\left\\{\frac{1}{\sqrt{24}}\left(\min\left\\{r_{1}-1,c_{\tilde{k}}^{r}-r_{1}\right\\}\right)^{3/2}\Delta_{1}^{j}\right\\}-\sqrt{8\log\left(Td^{\frac{1}{4}}\right)}$
$\displaystyle=\frac{1}{\sqrt{24}}\left(\min\left\\{r_{1}-1,c_{\tilde{k}}^{r}-r_{1}\right\\}\right)^{3/2}\max\limits_{1\leq
j\leq
d}\left\\{\Delta_{1}^{j}\right\\}-\sqrt{8\log\left(Td^{\frac{1}{4}}\right)},$
(B.35)

where the third inequality is due to the result of Lemma 5 in the online
supplement of Baranowski et al. (2019). Now,
$r_{1}-1=r_{1}-r_{0}-1\geq\delta_{T}-1>\delta_{T}/3$, because in the case of
the continuous piecewise-linear signals of Scenario (S2) we necessarily have
that $\delta_{T}\geq 2$; otherwise change-point identifiability would not be
possible. In addition, since $c_{\tilde{k}}^{r}\in I_{1}^{R}$, then
$c_{\tilde{k}}^{r}-r_{1}\geq\delta_{T}/3$, meaning that

$\min\left\\{c_{\tilde{k}}^{r}-r_{1},r_{1}-1\right\\}\geq\frac{\delta_{T}}{3}.$
(B.36)

The result in (B), Assumption (A2) in Section 2.2.2 of the main paper, and
(B.36) yield

$\displaystyle G_{1,c_{\tilde{k}}^{r}}^{\tilde{b}_{1}}$
$\displaystyle\geq\frac{1}{\sqrt{24}}\left(\frac{\delta_{T}}{3}\right)^{3/2}\max\limits_{1\leq
j\leq
d}\left\\{\Delta_{1}^{j}\right\\}-\sqrt{8\log\left(Td^{\frac{1}{4}}\right)}$
$\displaystyle\geq\frac{1}{\sqrt{24}}\left(\frac{\delta_{T}}{3}\right)^{3/2}\underline{f}_{T}-\sqrt{8\log\left(Td^{\frac{1}{4}}\right)}$
$\displaystyle=\delta_{T}^{3/2}\underline{f}_{T}\left(\frac{1}{3\sqrt{72}}-\frac{\sqrt{8\log\left(Td^{\frac{1}{4}}\right)}}{\delta_{T}^{3/2}\underline{f}_{T}}\right)\geq\left(\frac{1}{3\sqrt{72}}-\frac{2\sqrt{2}}{C^{*}}\right)\delta_{T}^{3/2}\underline{f}_{T}$
$\displaystyle=C_{2}\delta_{T}^{3/2}\underline{f}_{T}>\zeta_{T,d}.$ (B.37)

The result in (B) above proves that there will be an interval of the form
$[1,c_{\tilde{k}}^{r}]$, with $c_{\tilde{k}}^{r}>r_{1}$, such that
$[1,c_{\tilde{k}}^{r}]$ contains only $r_{1}$ and $\max_{1\leq
b<c_{\tilde{k}}^{r}}G_{1,c_{\tilde{k}}^{r}}^{b}>\zeta_{T,d}$. Let us, for
$k^{*}\in\left\\{1,\ldots,K\right\\}$, denote by $c_{k^{*}}^{r}\leq
c_{\tilde{k}}^{r}$ the first right-expanding point where this happens and let
$b_{1}={\rm argmax}_{1\leq t<c_{k^{*}}^{r}}G_{1,c_{k^{*}}^{r}}^{t}$ with
$G_{1,c_{k^{*}}^{r}}^{b_{1}}>\zeta_{T,d}$. Note that $b_{1}$ can not be an
estimation of any other change-point as $[1,c_{k^{*}}^{r}]$ includes only
$r_{1}$.

Allow us, for ease of presentation, to denote by $q_{1}:={\rm
argmax}_{j=1,\ldots,d}\left\\{C_{1,c_{k^{*}}^{r}}^{b_{1}}(\bm{X_{j}})\right|$
and $q_{1}$ is basically the location index of the univariate component data
sequence where the contrast function value got maximised at point $b_{1}$.
Apparently, it holds that

$G_{1,c_{k^{*}}^{r}}^{b_{1}}=\max\limits_{1\leq j\leq
d}\left\\{C_{1,c_{k^{*}}^{r}}^{b_{1}}(\bm{X_{j}})\right\\}=C_{1,c_{k^{*}}^{r}}^{b_{1}}(\bm{X_{q_{1}}}).$

Our aim now is to find $\tilde{\gamma}_{T}>0$ such that for any
$b^{*}\in\left\\{1,\ldots,c_{k^{*}}^{r}-1\right\\}$ with
$\left|b^{*}-r_{1}\right|\left(\Delta_{1}^{q_{1}}\right)^{2/3}>\tilde{\gamma}_{T}$,
it holds that

$\left(C_{1,c_{k^{*}}^{r}}^{r_{1}}(\bm{X_{q_{1}}})\right)^{2}>\left(C_{1,c_{k^{*}}^{r}}^{b^{*}}(\bm{X_{q_{1}}})\right)^{2}.$
(B.38)

Proving (B.38) and using the definition of $b_{1}$, we can then conclude that
$|b_{1}-r_{1}|\left(\Delta_{1}^{q_{1}}\right)^{2/3}\leq\tilde{\gamma}_{T}$.
The inequality in (B.38) can be expressed as

$\displaystyle\left(C_{1,c_{k^{*}}^{r}}^{r_{1}}(\bm{f_{q_{1}}})\right)^{2}-\left(C_{1,c_{k^{*}}^{r}}^{b^{*}}(\bm{f_{q_{1}}})\right)^{2}$
$\displaystyle>\left(C_{1,c_{k^{*}}^{r}}^{b^{*}}(\bm{\epsilon_{q_{1}}})\right)^{2}-\left(C_{1,c_{k^{*}}^{r}}^{r_{1}}(\bm{\epsilon_{q_{1}}})\right)^{2}$
$\displaystyle+2\left\langle\bm{\phi_{1,c_{k^{*}}^{r}}^{b^{*}}}\langle\bm{f_{q_{1}}},\bm{\phi_{1,c_{k^{*}}^{r}}^{b^{*}}}\rangle-\bm{\phi_{1,c_{k^{*}}^{r}}^{r_{1}}}\langle\bm{f_{q_{1}}},\bm{\phi_{1,c_{k^{*}}^{r}}^{r_{1}}}\rangle,\bm{\epsilon_{q_{1}}}\right\rangle.$
(B.39)

W.l.o.g. assume that $b^{*}\geq r_{1}$ and a similar approach as below holds
when $b^{*}<r_{1}$. We denote by

$\Lambda:=\left(C_{1,c_{k^{*}}^{r}}^{r_{1}}(\bm{f_{q_{1}}})\right)^{2}-\left(C_{1,c_{k^{*}}^{r}}^{b^{*}}(\bm{f_{q_{1}}})\right)^{2}$

and for the terms in the right-hand side of (B), we get that

$\displaystyle\left(C_{1,c_{k^{*}}^{r}}^{b^{*}}(\bm{\epsilon_{q_{1}}})\right)^{2}-\left(C_{1,c_{k^{*}}^{r}}^{r_{1}}(\bm{\epsilon_{q_{1}}})\right)^{2}$
$\displaystyle\leq\left(\max_{s,e,b:s\leq b<e}\left\\{\max_{1\leq j\leq
d}\left\\{C_{s,e}^{b}(\bm{\epsilon_{j}})\right\\}\right\\}\right)^{2}$
$\displaystyle\;\;-\left(C_{1,c_{k}^{r}}^{r_{1}}(\bm{\epsilon_{q_{1}}})\right)^{2}$
$\displaystyle\leq\left(\max_{s,e,b:s\leq b<e}\left\\{\max_{1\leq j\leq
d}\left\\{C_{s,e}^{b}(\bm{\epsilon_{j}})\right\\}\right\\}\right)^{2}$
$\displaystyle\leq 8\log\left(Td^{\frac{1}{4}}\right),$

while from (B) and Lemma 7 in the online supplement of Baranowski et al.
(2019),

$\displaystyle
2\left\langle\bm{\phi_{1,c_{k^{*}}^{r}}^{b^{*}}}\langle\bm{f_{q_{1}}},\bm{\phi_{1,c_{k^{*}}^{r}}^{b^{*}}}\rangle-\bm{\phi_{1,c_{k^{*}}^{r}}^{r_{1}}}\langle\bm{f_{q_{1}}},\bm{\phi_{1,c_{k^{*}}^{r}}^{r_{1}}}\rangle,\bm{\epsilon_{q_{1}}}\right\rangle$
$\displaystyle\leq
2\|\bm{\phi_{1,c_{k^{*}}^{r}}^{b^{*}}}<\bm{f_{q_{1}}},\bm{\phi_{1,c_{k^{*}}^{r}}^{b^{*}}}>-\bm{\phi_{1,c_{k^{*}}^{r}}^{r_{1}}}<\bm{f_{q_{1}}},\bm{\phi_{1,c_{k^{*}}^{r}}^{r_{1}}}>\|_{2}\sqrt{8\log\left(Td^{\frac{1}{4}}\right)}$
$\displaystyle=2\sqrt{8\Lambda\log\left(Td^{\frac{1}{4}}\right)}.$

Therefore (B) is satisfied if the stronger inequality
$\Lambda>8\log\left(Td^{1/4}\right)+2\sqrt{8\Lambda\log\left(Td^{1/4}\right)}$
is satisfied, which has solution

$\Lambda>(2\sqrt{2}+4)^{2}\log\left(Td^{\frac{1}{4}}\right).$ (B.40)

Using now Lemma 7 in the online supplement of Baranowski et al. (2019), we
have that (B.40) is implied by

$\displaystyle\frac{1}{63}\left(\min\left\\{\left|r_{1}-b^{*}\right|,r_{1}-1\right\\}\right)^{3}\left(\Delta_{1}^{q_{1}}\right)^{2}>\left(2\sqrt{2}+4\right)^{2}\log\left(Td^{\frac{1}{4}}\right)$
$\displaystyle\Leftrightarrow\min\left\\{\left|r_{1}-b^{*}\right|,r_{1}-1\right\\}>\frac{\left(63\log\left(Td^{\frac{1}{4}}\right)\right)^{1/3}(2\sqrt{2}+4)^{2/3}}{\left(\Delta_{1}^{q_{1}}\right)^{2/3}}$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\quad\;\;=\frac{C_{3}\left(\log\left(Td^{\frac{1}{4}}\right)\right)^{1/3}}{\left(\Delta_{1}^{q_{1}}\right)^{2/3}}.$
(B.41)

However,

$\displaystyle\min\left\\{r_{1}-1,c_{k^{*}}^{r}-r_{1}\right\\}>2^{1/3}C_{3}\frac{\left(\log\left(Td^{\frac{1}{4}}\right)\right)^{1/3}}{\left(\Delta_{1}^{q_{1}}\right)^{2/3}}-1$
(B.42)

and this is because if we assume that

$\min\left\\{r_{1}-1,c_{k^{*}}^{r}-r_{1}\right\\}\leq
2^{1/3}C_{3}\left(\log\left(Td^{\frac{1}{4}}\right)\right)^{1/3}/\left(\Delta_{1}^{q_{1}}\right)^{2/3}-1$

then we obtain that

$\displaystyle G_{1,c_{k^{*}}^{r}}^{b_{1}}$
$\displaystyle=C_{1,c_{k^{*}}^{r}}^{b_{1}}(\bm{X_{q_{1}}})\leq\left|C_{1,c_{k^{*}}^{r}}^{b_{1}}(\bm{X_{q_{1}}})-C_{1,c_{k^{*}}^{r}}^{b_{1}}(\bm{f_{q_{1}}})\right|+C_{1,c_{k^{*}}^{r}}^{b_{1}}(\bm{f_{q_{1}}})$
$\displaystyle\leq\sqrt{8\log\left(Td^{\frac{1}{4}}\right)}+C_{1,c_{k^{*}}^{r}}^{b_{1}}(\bm{f_{q_{1}}})$
$\displaystyle\leq\frac{1}{\sqrt{3}}\left(\min\left\\{r_{1}-1,c_{k^{*}}^{r}-r_{1}\right\\}+1\right)^{3/2}\Delta_{1}^{q_{1}}+\sqrt{8\log\left(Td^{\frac{1}{4}}\right)}$
$\displaystyle\leq\frac{1}{\sqrt{3}}\left(2^{1/3}C_{3}\frac{\left(\log\left(Td^{\frac{1}{4}}\right)\right)^{1/3}}{\left(\Delta_{1}^{q_{1}}\right)^{2/3}}\right)^{3/2}\Delta_{1}^{q_{1}}+\sqrt{8\log\left(Td^{\frac{1}{4}}\right)}$
$\displaystyle=\sqrt{\frac{2}{3}}C_{3}^{3/2}\sqrt{\log\left(Td^{\frac{1}{4}}\right)}+\sqrt{8\log\left(Td^{\frac{1}{4}}\right)}$
$\displaystyle=\left(\sqrt{\frac{2}{3}}C_{3}^{3/2}+\sqrt{8}\right)\sqrt{\log\left(Td^{\frac{1}{4}}\right)}=C_{1}\sqrt{\log\left(Td^{\frac{1}{4}}\right)}\leq\zeta_{T,d},$
(B.43)

where the third inequality is a direct application of Lemma 5 in the online
supplement of Baranowski et al. (2019). The result in (B) comes to a
contradiction to $G_{1,c_{k^{*}}^{r}}^{b_{1}}>\zeta_{T,d}$ which has already
been proven in (B). Therefore, (B.42) holds and for $Td^{1/4}$ being
sufficiently large,

$\min\left\\{r_{1}-1,c_{k^{*}}^{r}-r_{1}\right\\}>2^{1/3}C_{3}\frac{\left(\log\left(Td^{\frac{1}{4}}\right)\right)^{1/3}}{\left(\Delta_{1}^{q_{1}}\right)^{2/3}}-1>C_{3}\frac{\left(\log\left(Td^{\frac{1}{4}}\right)\right)^{1/3}}{\left(\Delta_{1}^{q_{1}}\right)^{2/3}}.$
(B.44)

From (B.44) we deduce that (B) is restricted to

$\left|r_{1}-b^{*}\right|>C_{3}\frac{\left(\log\left(Td^{\frac{1}{4}}\right)\right)^{1/3}}{\left(\Delta_{1}^{q_{1}}\right)^{2/3}},$

which implies (B.38). Therefore, necessarily,

$\left|b_{1}-r_{1}\right|\left(\Delta_{1}^{q_{1}}\right)^{2/3}\leq
C_{3}\left(\log\left(Td^{\frac{1}{4}}\right)\right)^{1/3}.$ (B.45)

So far, we have proven that working under the sets $\tilde{A}_{T}^{*}$ and
$\tilde{B}_{T}$, there will be an interval of the form $[1,c_{k^{*}}^{r}]$,
with $G_{1,c_{k^{*}}^{r}}^{b_{1}}>\zeta_{T,d}$, where $b_{1}=\underset{1\leq
t<c_{k^{*}}^{r}}{\rm argmax}\;G_{1,c_{k^{*}}^{r}}^{t}$ is an estimation of
$r_{1}$ that satisfies (B.45).
Step 3.2: After detecting the first change-point, MID follows the same process
as in Step 3.1 in the set $[c_{k^{*}}^{r},T]$, which contains
$r_{2},\ldots,r_{N}$. This means that we do not check for possible change-
points in the interval $[b_{1}+1,c_{k^{*}}^{r})$ and we need to prove that:

  * (M.1)

There is no other change-point in $[b_{1}+1,c_{k^{*}}^{r})$, apart from
possibly the already detected $r_{1}$;

  * (M.2)

$c_{k^{*}}^{r}$ is at a location which allows for detection of $r_{2}$.

For (M.1): The approach is the same as the one in Step 3.2 in the proof of
Theorem 2.1 and will not be repeated here.

Similarly to the approach in Step 3.1, our method applied now to
$[c_{k^{*}}^{r},T]$, will first detect $r_{2}$ or $r_{N}$ depending on whether
$r_{2}-c_{k^{*}}^{r}$ is smaller or larger than $T-r_{N}$. If
$T-r_{N}<r_{2}-c_{k^{*}}^{r}$ then $r_{N}$ will get isolated first and the
procedure to show its detection is exactly the same as in Step 3.1 where we
explained the detection of $r_{1}$. Therefore, w.l.o.g. and also for the sake
of showing (M.2) let us assume that $r_{2}-c_{k^{*}}^{r}\leq T-r_{N}$.
For (M.2): Due to the isolation aspect of MID, it is certain that there exists
a right expanding point $c_{k_{2}}^{r}$, such that $c_{k_{2}}^{r}\in
I_{2}^{R}$, with $I_{j}^{R}$ being defined in (A) of the main paper. We will
show that $r_{2}$ gets detected in $[c_{k^{*}}^{r},c_{k^{*}_{2}}^{r}]$, for
$k_{2}^{*}\leq k_{2}$ and its detection is $b_{2}={\rm
argmax}_{c_{k^{*}}^{r}\leq
t<c_{k^{*}_{2}}^{r}}G_{c_{k^{*}}^{r},c_{k^{*}_{2}}^{r}}^{t}$, which satisfies
$\left|b_{2}-r_{2}\right|\left(\Delta_{2}^{q_{2}}\right)^{2/3}\leq
C_{3}\left(\log\left(Td^{1/4}\right)\right)^{1/3}$, where $q_{2}:={\rm
argmax}_{j=1,\ldots,d}\left\\{C_{c_{k^{*}}^{r},c_{k^{*}_{2}}^{r}}^{b_{2}}(\bm{X_{j}})\right\\}$.
Using again Lemma 5 from the online supplement of Baranowski et al. (2019) and
for $\tilde{b}_{2}={\rm argmax}_{c_{k^{*}}^{r}\leq
t<c_{k_{2}}^{r}}G_{c_{k^{*}}^{r},c_{k_{2}}^{r}}^{t}$, we have that

$\displaystyle G_{c_{k^{*}}^{r},c_{k_{2}}^{r}}^{\tilde{b}_{2}}$
$\displaystyle\geq G_{c_{k^{*}}^{r},c_{k_{2}}^{r}}^{r_{2}}\geq
H_{c_{k^{*}}^{r},c_{k_{2}}^{r}}^{r_{2}}-\sqrt{8\log\left(Td^{\frac{1}{4}}\right)}$
$\displaystyle\geq\frac{1}{\sqrt{24}}\left(\min\left\\{r_{2}-c_{k^{*}}^{r},c_{k_{2}}^{r}-r_{2}\right\\}\right)^{3/2}\underline{f}_{T}-\sqrt{8\log\left(Td^{\frac{1}{4}}\right)}.$
(B.46)

By construction,

$\displaystyle c_{k_{2}}^{r}-r_{2}\geq\frac{\delta_{T}}{3}$ $\displaystyle
r_{2}-c_{k^{*}}^{r}\geq
r_{2}-c_{k}^{r}=r_{2}-r_{1}-(c_{k}^{r}-r_{1})\geq\delta_{T}-(c_{k}^{r}-r_{1})>\delta_{T}-2\frac{\delta_{T}}{3}=\frac{\delta_{T}}{3},$

which means that
$\min\left\\{c_{k_{2}}^{r}-r_{2},r_{2}-c_{k^{*}}^{r}\right\\}\geq\delta_{T}/3$.
Therefore, continuing from (B) and using the exact same calculations as in
(B), we have that

$\displaystyle G_{c_{k^{*}}^{r},c_{k_{2}^{*}}^{r}}^{b_{2}}\geq
C_{2}\delta_{T}^{3/2}\underline{f}_{T}>\zeta_{T,d}.$

We will now show that
$\left|b_{2}-r_{2}\right|\left(\Delta_{2}^{q_{2}}\right)^{2/3}\leq
C_{3}\left(\log T\right)^{1/3}$. Following the same process as in Step 3.1 and
assuming now that $b_{2}<r_{2}$, we have that for
$b^{*}\in\left\\{c_{k^{*}}^{r},\ldots,c_{k^{*}_{2}}^{r}-1\right\\}$,

$\left(C_{c_{k^{*}}^{r},c_{k^{*}_{2}}^{r}}^{r_{2}}(\bm{X_{q_{2}}})\right)^{2}>\left(C_{c_{k^{*}}^{r},c_{k^{*}_{2}}^{r}}^{b^{*}}(\bm{X_{q_{2}}})\right)^{2}$
(B.47)

is implied by
$\min\left\\{\left|b^{*}-r_{2}\right|,c_{k_{2}^{*}}^{r}-r_{2}\right\\}>C_{3}\left(\log\left(Td^{1/4}\right)\right)^{1/3}/\left(\Delta_{2}^{q_{2}}\right)^{2/3}$.
However, following the same procedure as in Step 3.1 we can show that for
sufficiently large $Td^{1/4}$,

$\min\left\\{c_{k_{2}^{*}}^{r}-r_{2},r_{2}-c_{k^{*}}^{r}\right\\}>C_{3}\frac{\left(\log\left(Td^{1/4}\right)\right)^{1/3}}{\left(\Delta_{2}^{q_{2}}\right)^{2/3}}.$

Thus, (B.47) is implied by
$\left|b^{*}-r_{2}\right|\left(\Delta_{2}^{q_{2}}\right)^{2/3}>C_{3}(\log\left(Td^{1/4}\right))^{1/3}$.
Therefore,
$\left|b_{2}-r_{2}\right|\left(\Delta_{2}^{q_{2}}\right)^{2/3}>C_{3}\left(\log\left(Td^{1/4}\right)\right)^{1/3}$
would necessarily mean that
$G_{c_{k^{*}}^{r},c_{k^{*}_{2}}^{r}}^{r_{2}}>G_{c_{k^{*}}^{r},c_{k^{*}_{2}}^{r}}^{b_{2}}$,
which is not true by the definition of $b_{2}$. Having said this, we conclude
that $\left|b_{2}-r_{2}\right|\left(\Delta_{2}^{q_{2}}\right)^{2/3}\leq
C_{3}\left(\log\left(Td^{1/4}\right)\right)^{1/3}$.

Having detected $r_{2}$, then our algorithm will proceed in the interval
$[s,e]=[c_{k^{*}_{2}}^{r},T]$ and all the change-points will get detected one
by one since Step 3.2 will be applicable as long as there are previously
undetected change-points in $[s,e]$. Denoting by $\hat{r}_{j}$ the estimation
of $r_{j}$ as we did in the statement of the theorem and for $[s^{*},e^{*}]$
being the interval where the isolation and detection of $r_{j}$ occurs (in the
way that we explained in Steps 3.1 and 3.2) allow us to denote by $q_{j}:={\rm
argmax}_{m=1,\ldots,d}\left\\{C_{s^{*},e^{*}}^{\hat{r}_{j}}(\bm{X_{m}})\right\\}$.
Then, Steps 3.1 and 3.2 have shown that MID will detect all change-points one
by one and
$\left|\hat{r}_{j}-r_{j}\right|\left(\Delta_{j}^{q_{j}}\right)^{2/3}\leq
C_{3}\left(\log\left(Td^{\frac{1}{4}}\right)\right)^{1/3},\quad\forall
j\in\left\\{1,\ldots,N\right\\}.$
Step 4: The arguments given in Steps 1-3 hold in
$\tilde{A}_{T}^{*}\cap\tilde{B}_{T}$. At the beginning of the algorithm,
$s=1,e=T$ and for $N\geq 1$, there exist $k_{1}\in\left\\{1,\ldots,K\right\\}$
such that $s_{k_{1}}=s,e_{k_{1}}\in I_{1}^{R}$ and
$k_{2}\in\left\\{1,\ldots,K\right\\}$ such that $s_{k_{2}}\in
I_{N}^{L},e_{k_{2}}=e$, where $I_{j}^{R}$ and $I_{j}^{L}$ are defined in (A)
of the main paper. As in our previous steps, w.l.o.g. assume that $r_{1}\leq
T-r_{N}+1$, meaning that $r_{1}$ gets isolated and detected first in an
interval $[s,c_{k^{*}}^{r}]$, where $c_{k^{*}}^{r}\leq e_{k_{1}}$. Then,
$\hat{r}_{1}={\rm argmax}_{s\leq t<c_{k^{*}}^{r}}G_{s,c_{k^{*}}^{r}}^{t}$ is
the estimated location for $r_{1}$ and
$\left|r_{1}-\hat{r}_{1}\right|\left(\Delta_{1}^{q_{1}}\right)^{2/3}\leq
C_{3}\left(\log\left(Td^{1/4}\right)\right)^{1/3}$. After this, the algorithm
continues in $[c_{k^{*}}^{r},T]$ and keeps detecting all the change-points as
explained in Step 3. It is important to note that there will not be any double
detection issues because naturally, at each step of the algorithm, the new
interval $[s,e]$ does not include any previously detected change-points.

Once all the change-points have been detected one by one, then $[s,e]$ will
have no other change-points in it. Our method will keep interchangeably
checking for possible change-points in right- and left-expanding intervals
denoted by $[s^{*},e^{*}]$. MID will not detect anything in $[s^{*},e^{*}]$
since $\forall b\in[s^{*},e^{*})$, using (B) we have that

$G_{s^{*},e^{*}}^{b}\leq
H_{s^{*},e^{*}}^{b}+\sqrt{8\log\left(Td^{\frac{1}{4}}\right)}=\sqrt{8\log\left(Td^{1/4}\right)}<C_{1}\sqrt{\log\left(Td^{1/4}\right)}\leq\zeta_{T,d}.$

After not detecting anything in all intervals of the above form, then the
algorithm concludes that there are not any change-points in $[s,e]$ and stops.
∎