# Random Forests for Change Point Detection

Malte Londschien1,2, Peter Bühlmann1, Solt Kovács1

1Seminar for Statistics, ETH Zürich, Switzerland
2AI Center, ETH Zürich, Switzerland

(May 2022)

###### Abstract

We propose a novel multivariate nonparametric multiple change point detection
method using classifiers. We construct a classifier log-likelihood ratio that
uses class probability predictions to compare different change point
configurations. We propose a computationally feasible search method that is
particularly well suited for random forests, denoted by changeforest. However,
the method can be paired with any classifier that yields class probability
predictions, which we illustrate by also using a $k$-nearest neighbor
classifier. We provide theoretical results motivating our choices. In a large
simulation study, our proposed changeforest method achieves improved empirical
performance compared to existing multivariate nonparametric change point
detection methods. An efficient implementation of our method is made available
for R, Python, and Rust users in the changeforest software package.

_Keywords_ break point detection, classifiers, multivariate time series,
nonparametric

##  1 Introduction

Change point detection considers the localization of abrupt distributional
changes in ordered observations, often time series. We focus on offline
problems, retrospectively detecting changes after all samples have been
observed. Inferring abrupt structural changes has a wide range of
applications, including bioinformatics (Olshen et al., 2004; Picard et al.,
2005), neuroscience (Kaplan et al., 2001), biochemistry (Hotz et al., 2013),
climatology (Reeves et al., 2007), and finance (Kim et al., 2005). A rich
literature has developed around statistical methods that recover change points
in different scenarios, see Truong et al. (2020) for a recent review.

Parametric change point detection methods typically assume that observations
between change points stem from a finite-dimensional family of distributions.
Change points can then be estimated by maximizing a regularized log-likelihood
over various segmentations. The classical scenario of independent univariate
Gaussian variables with constant variance and changes in the mean goes back to
the 1950s (Page, 1954, 1955). It has recently been studied by, among others,
Frick et al. (2014), Fryzlewicz (2014), and references therein. Pein et al.
(2017) consider a relaxed version that allows certain changes in the variance
in addition to the mean shifts. Generalizations exist for multivariate
scenarios. Recently, even high-dimensional scenarios have been studied: Wang
and Samworth (2018) consider multivariate Gaussian observations with sparse
mean shifts, Wang et al. (2021) study the problem of changing covariance
matrices of sub-Gaussian random vectors, Roy et al. (2017) study estimation in
Markov random field models, and Londschien et al. (2021) consider time-varying
Gaussian graphical models.

Nonparametric change point detection methods try to use measures that do not
rely on parametric forms of the distribution or the nature of change.
Proposals for univariate nonparametric change point detection methods include
Pettitt (1979), Carlstein (1988), Dümbgen (1991), and, more recently, Zou et
al. (2014) and Madrid–Padilla et al. (2021a). Multivariate setups are
challenging even in parametric scenarios. The few multivariate nonparametric
change point detection methods we are aware of are based on ranks (Lung-Yut-
Fong et al., 2015), distances (Matteson and James, 2014; Chen and Zhang, 2015;
Zhang and Chen, 2021), kernel distances (Arlot et al., 2019; Garreau and
Arlot, 2018), kernel densities (Madrid–Padilla et al., 2021b), and relative
density ratio estimation (Liu et al., 2013).

Many of these nonparametric proposals are related to hypothesis testing. For a
single change point, they maximize a test statistic that measures the
dissimilarity of distributions. In the modern era of statistics and machine
learning, many nonparametric methods are available to learn complex
conditional class probability distributions, such as random forests (Breiman,
2001) and neural networks (McCulloch and Pitts, 1943). These have proven to
often outperform simple rank and distance-based methods. Friedman (2004)
proposed to use binary classifiers for two-sample testing. Recently, Lopez-Paz
and Oquab (2017) applied this framework in combination with neural networks
and Hediger et al. (2022) with random forests. Similar to such two-sample
testing approaches, we use classifiers to construct a novel class of
multivariate nonparametric multiple change point detection methods.

###  1.1 Our contribution

We propose a novel classifier-based nonparametric change point detection
framework. Motivated by parametric change point detection, we construct a
classifier log-likelihood ratio that uses class probability predictions to
compare different change point configurations. Theoretical results for the
population case motivate the development of our algorithm.

We present a novel search method based on binary segmentation that requires a
constant number of classifier fits to find a single change point. For multiple
change points, the number of classifier fits required scales linearly with the
number of change points, making the algorithm highly efficient. Our method is
implemented for random forests and $k$-nearest neighbor classifiers in the
changeforest software package, available for R, Python, and Rust users. The
algorithm achieves improved empirical performance compared to existing
multivariate nonparametric change point detection methods.

##  2 A nonparametric classifier log-likelihood ratio

We are interested in detecting structural breaks in the probability
distribution of a time series. More formally, consider a sequence of
independent random variables $(X_{i})_{i=1}^{n}\subset{\mathbb{R}}^{p}$ with
distributions $\tilde{\mathbb{P}}_{1},\dotsc,\tilde{\mathbb{P}}_{n}$ such that
the map $i\mapsto\tilde{\mathbb{P}}_{i}$ is piecewise constant. Let

$\alpha^{0}:=\\{0,n\\}\cup\\{i:\tilde{\mathbb{P}}_{i}\neq\tilde{\mathbb{P}}_{i+1}\\}$

be the set of _segment boundaries_ and denote with $K^{0}:=|\alpha^{0}|-1$ the
total number of segments. We label the elements of $\alpha^{0}$ by their
natural order starting with zero. Then consecutive elements in $\alpha^{0}$
define segments ${{(\alpha_{k\shortminus 1}^{0},\alpha_{k}^{0}]}}$ for
$k=1,\dotsc,K^{0}$ within which the
$X_{i}\sim{\mathbb{P}}_{k}\coloneqq\tilde{\mathbb{P}}_{\alpha_{k}}$ are i.i.d.
We call $\alpha^{0}_{1},\dotsc,\alpha^{0}_{K^{0}\shortminus 1}$ the _change
points_ of the sequence $X_{1},\dotsc,X_{n}$. We aim to estimate the change
points (or equivalently $\alpha^{0}$) of the time series process
$X_{1},\dotsc,X_{n}$ upon observing a realization $x_{1},\dotsc,x_{n}$. We
construct a classifier log-likelihood ratio (1) and use it for change point
detection. We motivate this ratio later in Section 2.2, drawing parallels to
parametric methods that we introduce in Section 2.1.

For any segmentation $\alpha$, let ${\bf y}_{\alpha}=(y_{i})_{i=1}^{n}$ be the
sequence such that $y_{i}=k$ whenever $i\in{{(\alpha_{k\shortminus
1},\alpha_{k}]}}$. Consider a classification algorithm $\hat{p}$ that can
produce class probability predictions and denote with $\hat{p}_{\alpha}$ the
classifier trained on covariates ${\bf X}=(x_{i})_{i=1}^{n}$ and labels ${\bf
y}_{\alpha}$. Write $\hat{p}_{\alpha}(x_{i})_{k}$ for the trained classifier’s
class $k$ probability prediction for observation $x_{i}$. The _classifier log-
likelihood ratio_ for change point detection is

$G\left((x_{i})_{i=1}^{n}\mid\alpha,\hat{p}\right)\coloneqq\sum_{k=1}^{K}\sum_{i=\alpha_{k\shortminus
1}+1}^{\alpha_{k}}\log(\frac{n}{\alpha_{k}-\alpha_{k\shortminus
1}}\hat{p}_{\alpha}(x_{i})_{k}).$ (1)

Consider a setup with a single change point. If a split candidate $s$ is close
to that change point, we expect the trained classifier $\hat{p}_{\\{0,s,n\\}}$
to perform better than random at separating $x_{1},\ldots x_{s}$ from
$x_{s+1},\ldots x_{n}$ and
$G\left((x_{i})_{i=1}^{n}\mid\\{0,s,n\\},\hat{p}\right)$ will take a positive
value. We will later relate
$G\left((x_{i})_{i=1}^{n}\mid\\{0,s,n\\},\hat{p}\right)$ to the gain
$G_{(0,n]}(s)$ from Equation (5) used in binary segmentation, where
$G_{(0,n]}(s)$ is the reduction in loss when splitting the segment ${(0,n]}$
at $s$, used in binary segmentation style algorithms. For a segmentation
$\alpha$ with multiple change points, the classifier log-likelihood ratio (1)
measures the separability of the segments
${(\alpha_{0},\alpha_{1}]},\ldots,{(\alpha_{K\shortminus 1},\alpha_{K}]}$. We
expect this separability to be highest for $\alpha=\alpha^{0}$, which we
formalize in Proposition 2.

###  2.1 The parametric setting

Since first proposals from Page (1954), parametric approaches are typically
based on maximizing a parametric log-likelihood, see also Frick et al. (2014)
and references therein. For this, one assumes that the
$X_{i}\sim{\mathbb{P}}_{k}$ are i.i.d. within segments
$i\in{{(\alpha_{k\shortminus 1}^{0},\alpha_{k}^{0}]}}$ and that the
distributions
${\mathbb{P}}_{k}\in\\{{\mathbb{P}}_{\vartheta}\mid\vartheta\in\Theta\\}$ for
some finite-dimensional parameter space $\Theta$. The distribution of
$(X_{i})_{i=1}^{n}$ can then be parametrized with the tuple
$(\alpha^{0},(\vartheta^{0}_{k})_{k=1}^{K^{0}})$ such that for all
$k=1,\dotsc,K^{0}$ and $i=\alpha^{0}_{k\shortminus 1}+1,\dotsc,\alpha^{0}_{k}$
the $X_{i}$ are ${\mathbb{P}}_{\vartheta^{0}_{k}}$-distributed.

Assuming the ${\mathbb{P}}_{\vartheta}$ have densities $p_{\vartheta}$, we can
express the log-likelihood of observing some sequence $(x_{i})_{i=1}^{n}$ in a
setup parametrized by $(\alpha,(\vartheta_{k})_{k=1}^{K})$ as

$\ell\left((x_{i})_{i=1}^{n}\mid\alpha,(\vartheta_{k})_{k=1}^{K}\right)=\sum_{k=1}^{K}\sum_{i=\alpha_{k\shortminus
1}+1}^{\alpha_{k}}\log(p_{\vartheta_{k}}(x_{i})).$ (2)

An estimate of $\alpha_{0}$ can be obtained by maximizing (2). Since the log-
likelihood is increasing in $K=|\alpha|-1$, an additional penalty term
$\gamma>0$ has to be subtracted whenever the number of change points is
unknown. Let

$\hat{\alpha}_{\gamma}:=\underset{\alpha\in{\cal A}}{\operatorname*{argmax}}\
\underset{\vartheta_{1},\ldots,\vartheta_{K}}{\max}\
\ell\left((x_{i})_{i=1}^{n}\mid\alpha,(\vartheta_{k})_{k=1}^{K}\right)-|\alpha|\gamma\
,$ (3)

where ${\cal A}$ is a set of possible segmentations, typically restricted such
that each segment has a minimum length of $\delta n$ for some $\delta>0$.

The most popular parametrization of in-segment distributions are Gaussian
distributions with known variance and a changing mean, see for example Yao
(1988) and Fryzlewicz (2014). Next to changes in the mean, Killick et al.
(2012) also study changes in the variance, where the mean of the Gaussian is
unknown but constant over all segments. Kim and Siegmund (1989) study changes
in the parameters of a linear regression model. Shen and Zhang (2012) and
Cleynen and Lebarbier (2014) consider changes in the rate of a Poisson and
changes in the success rate in a negative binomial, with applications
detecting copy number alterations in DNA. Frick et al. (2014) consider the
general case where in-segment distributions belong to an exponential family.

###  2.2 The nonparametric setting

Often, a good parametric model is hard to justify and nonparametric methods
can be applied instead. Again, consider a classification algorithm $\hat{p}$
that yields class probability predictions after training and recall that ${\bf
y}_{\alpha}=(y_{i})_{i=1}^{n}$ with $y_{i}=k$ whenever
$i\in{{(\alpha_{k\shortminus 1},\alpha_{k}]}}$ and that ${{\bf
X}=(x_{i})_{i=1}^{n}}$. For any $\alpha$, training $\hat{p}$ on $({\bf X},{\bf
y}_{\alpha})$ yields the function $\hat{p}_{\alpha}\coloneqq\hat{p}({\bf
X},{\bf y}_{\alpha})$ that assigns class $k$ probability estimates
$\hat{p}_{\alpha}(x)_{k}$ corresponding to segments ${{(\alpha_{k\shortminus
1},\alpha_{k}]}}$ to observations $x$. Write
${\mathbb{P}}_{(u,v]}:=\frac{1}{v-u}\sum_{i=u+1}^{v}\tilde{\mathbb{P}}_{i}$
for the mixture distribution of $X_{u+1},\ldots,X_{v}$ and assume that
densities $p_{(u,v]}$ exist. Assuming a uniform prior on class probabilities,
the trained classifier $\hat{p}_{\alpha}$ estimates

$\hat{p}_{\alpha}(x)_{k}\approx{\mathbb{P}}(Y=k\mid
X=x)=\frac{\differential{\mathbb{P}}(X=x\mid
Y=k)}{\differential{\mathbb{P}}(X=x)}{\mathbb{P}}(Y=k)=\frac{p_{(\alpha_{k\shortminus
1},\alpha_{k}]}(x_{i})}{p_{(0,n]}(x_{i})}\frac{\alpha_{k}-\alpha_{k\shortminus
1}}{n}$

for $k=1,\ldots K$. This yields an approximation

$\frac{p_{(\alpha_{k\shortminus
1},\alpha_{k}]}(x)}{p_{(0,n]}(x)}\approx\frac{n}{\alpha_{k}-\alpha_{k\shortminus
1}}\hat{p}_{\alpha}(x)_{k}\,,$

which we use in the definition of the classifier log-likelihood ratio (1):

$\displaystyle G\left((x_{i})_{i=1}^{n}\mid\alpha,\hat{p}\right)$
$\displaystyle=\sum_{k=1}^{K}\sum_{i=\alpha_{k\shortminus
1}+1}^{\alpha_{k}}\log(\frac{n}{\alpha_{k}-\alpha_{k\shortminus
1}}\hat{p}_{\alpha}(x_{i})_{k})\approx\sum_{k=1}^{K}\sum_{i=\alpha_{k\shortminus
1}+1}^{\alpha_{k}}\log(\frac{p_{{(\alpha_{k\shortminus
1},\alpha_{k}]}}(x_{i})}{p_{(0,n]}(x_{i})})\,.$

This resembles the parametric log-likelihood (2), where
$p_{(\alpha^{0}_{k\shortminus 1},\alpha^{0}_{k}]}=p_{\vartheta^{0}_{k}}$, up
to $-\sum_{i=1}^{n}\linebreak[1]\log(p_{(0,n]}(x_{i}))$, which is independent
of $\alpha$. In Equation (1), the classification algorithm $\hat{p}$ takes the
role of the parametrization $\vartheta\in\Theta$ in Equation (2). Analogously
to (3), a nonparametric estimator for $\alpha^{0}$ could be given by

$\hat{\alpha}_{\gamma}=\operatorname*{argmax}_{\alpha\in{\cal
A}}G\left((x_{i})_{i=1}^{n}\mid\alpha,\hat{p}\right)-\gamma|\alpha|\ ,$ (4)

where ${\cal A}$ is a suitable set of possible segmentations.

######  Remark 1.

To further motivate (1), assume that $\hat{p}$ is a generative classifier (as
in linear discriminant analysis). Write
$\hat{\vartheta}({(u,v]})\coloneqq\operatorname*{argmax}_{\vartheta}\sum_{i=u+1}^{v}\log(p_{\vartheta}(x_{i}))$
and for any segmentation $\alpha$ let
$\pi_{k}\coloneqq\frac{\alpha_{k\shortminus 1}-\alpha_{k}}{n}$, such that
$\hat{p}_{\alpha}(x_{i})_{k}=\frac{\pi_{k}\
p_{\hat{\vartheta}({{(\alpha_{k\shortminus
1},\alpha_{k}]}})}(x_{i})}{\sum_{j=1}^{K}\pi_{j}\
p_{\hat{\vartheta}({{(\alpha_{j\shortminus 1},\alpha_{j}]}})}(x_{i})}$. Then
the classifier log-likelihood ratio is equal to the parametric log-likelihood
ratio between the segment distributions
$p_{\hat{\vartheta}({{(\alpha_{k\shortminus 1},\alpha_{k}]}})}$ and the
mixture model $\sum_{j=1}^{K}\pi_{j}\
p_{\hat{\vartheta}({{(\alpha_{j\shortminus 1},\alpha_{j}]}})}(\cdot)_{j}$:

$\displaystyle G\left((x_{i})_{i=1}^{n}\mid\alpha,\hat{p}\right)$
$\displaystyle=\sum_{k=1}^{K}\sum_{i=\alpha_{k\shortminus
1}+1}^{\alpha_{k}}\log(\frac{p_{\hat{\vartheta}({{(\alpha_{k\shortminus
1},\alpha_{k}]}})}(x_{i})}{\pi_{j}\,\sum_{j=1}^{K}p_{\hat{\vartheta}({{(\alpha_{j\shortminus
1},\alpha_{j}]}})}(x_{i})})$
$\displaystyle=\ell\left((x_{i})_{i=1}^{n}\mid(\alpha,(\hat{\vartheta}({{(\alpha_{k\shortminus
1},\alpha_{k}]}}))_{k=1}^{K})\right)-\sum_{n=1}^{n}\log(\sum_{k=1}^{K}\pi_{k}\
p_{\hat{\vartheta}({{(\alpha_{k\shortminus 1},\alpha_{k}]}})}(x_{i}))\ .$

The following proposition suggests that for any estimator approximating the
Bayes classifier, minimizing $G$ over segmentations $\alpha$ yields a
reasonable estimator for $\alpha^{0}$.

######  Proposition 2.

Let

$p^{\infty}_{\alpha}(x)\coloneqq{\mathbb{P}}\left(Y=k\mid
X=x\right)=\left(\frac{\alpha_{k\shortminus
1}-\alpha_{k}}{n}\frac{p_{(\alpha_{k\shortminus
1},\alpha_{k}]}}{p_{(0,n]}}\right)_{k=1}^{|\alpha|-1}$

be the infinite sample Bayes classifier corresponding to the segmentation
$\alpha$. If ${\cal A}$ is a set of segmentations containing $\alpha^{0}$ and
$D_{KL}(P\;\|\;Q)$ is the Kullback-Leibler divergence of $Q$ with respect to
$P$, then

$\displaystyle\max_{\alpha\in{\cal
A}}\,{\mathbb{E}}\left[G\left((X_{i})_{i=1}^{n}\mid\alpha,p^{\infty}\right)\right]$
$\displaystyle={\mathbb{E}}\left[G\left((X_{i})_{i=1}^{n}\mid\alpha^{0},p^{\infty}\right)\right]$
$\displaystyle=\sum_{k=1}^{K}(\alpha^{0}_{k}-\alpha^{0}_{k\shortminus
1})D_{KL}({\mathbb{P}}_{k}\;\|\;{\mathbb{P}}_{(0,n]})$

and any maximizer of
${\mathbb{E}}\left[G\left((X_{i})_{i=1}^{n}\mid\alpha,p^{\infty}\right)\right]$
is equal to $\alpha^{0}$ or is an oversegmentation of $\alpha^{0}$. In
particular, the maximizer $\alpha$ with the smallest cardinality $|\alpha|$ is
equal to $\alpha^{0}$.

##  3 The changeforest algorithm

Finding a solution to the nonparametric estimator (4) by evaluating the
classifier log-likelihood ratio (1) for all possible candidate segmentations
$\alpha\in{\cal A}$ is infeasible. We present an efficient search procedure
based on binary segmentation that finds an approximate solution to (4). We
furthermore motivate the use of random forests as classifiers and present a
model selection procedure adapted to our search procedure. These building
blocks combine to form the changeforest algorithm, which we summarize in
Section 3.5.

###  3.1 Binary segmentation

Binary segmentation (Vostrikova, 1981) is a popular greedy algorithm to obtain
an approximate solution to the parametric maximum log-likelihood estimator
(3). For
$\hat{\vartheta}({(u,v]})\coloneqq\operatorname*{argmax}_{\vartheta}\linebreak[1]\sum_{i=u+1}^{v}\log(p_{\vartheta}(x_{i}))$,
binary segmentation recursively splits segments at the split $s$ maximizing
the gain

$G_{(u,v]}(s):=\sum_{i=u+1}^{s}\log(\frac{p_{\hat{\vartheta}({(u,s]})}(x_{i})}{p_{\hat{\vartheta}({(u,v]})}(x_{i})})+\sum_{i=s+1}^{v}\log(\frac{p_{\hat{\vartheta}({(s,v]})}(x_{i})}{p_{\hat{\vartheta}({(u,v]})}(x_{i})})\,,$
(5)

the increase in log-likelihood, a parametric log-likelihood ratio, until a
stopping criterion is met. For changes in the mean, where
$\hat{\vartheta}({(u,v]})=\frac{1}{v-u}\sum_{i=u+1}^{v}x_{i}$ and
${\log(p_{\vartheta}(x))=\frac{1}{\sqrt{2\pi}}-\frac{1}{2}(x-\vartheta)^{2}}$,
the normalized gain $\frac{2}{v-u}G_{(u,v]}(s)$ is equal to
$\left(\sqrt{\frac{v-s}{s-u}}\sum_{i=u+1}^{s}x_{i}-\sqrt{\frac{s-u}{v-s}}\sum_{i=s+1}^{v}x_{i}\right)^{2}$,
the square of the CUSUM statistic, first presented by Page (1954). Binary
segmentation typically requires ${\cal O}(n\log(n))$ evaluations of the gain
and is typically faster than search methods based on dynamic programming such
as PELT (Killick et al., 2012).

We replace the parametric gain $G_{(u,v]}(s)$ in binary segmentation with the
nonparametric classifier log-likelihood ratio from (1)

$\displaystyle G\left((x_{i})_{i=u+1}^{s}\mid\\{u,s,v\\},\hat{p}\right)$
$\displaystyle=\sum_{i=u+1}^{s}\log(\frac{v-u}{s-u}\hat{p}_{\\{u,s,v\\}}(x_{i})_{1})+\sum_{i=s+1}^{v}\log(\frac{v-u}{v-s}\hat{p}_{\\{u,s,v\\}}(x_{i})_{2})$
(6)
$\displaystyle\approx\sum_{i=u+1}^{s}\log(\frac{p_{{(u,s]}}(x_{i})}{p_{{(u,v]}}(x_{i})})+\sum_{i=s+1}^{v}\log(\frac{p_{{(s,v]}}(x_{i})}{p_{{(u,v]}}(x_{i})})\,.$

In many parametric settings, the expected gain curve $G_{(u,v]}$ will be
piecewise convex between the underlying change points (Kovács et al., 2020c).
Proposition 3 shows that the same holds for the nonparametric variant (6) in
the population case. Figure 1 shows the nonparametric gains for a random
forest and a $k$-nearest neighbor classifier applied to simulated data with
two change points.

######  Proposition 3.

For the Bayes classifier $p^{\infty}$, the expected classifier log-likelihood
ratio

$s\mapsto{\mathbb{E}}\left[G\left((X_{i})_{i=u+1}^{v}\mid\\{u,s,v\\},p^{\infty}\right)\right]$

is piecewise convex between the underlying change points $\alpha^{0}$, with
strict convexity in the segment ${(\alpha^{0}_{k\shortminus
1},\alpha^{0}_{k}]}$ if ${\mathbb{P}}_{(u,\alpha^{0}_{k\shortminus
1}]}\neq{\mathbb{P}}_{(\alpha^{0}_{k\shortminus 1},\alpha^{0}_{k}]}$ or
${\mathbb{P}}_{(\alpha^{0}_{k\shortminus
1},\alpha^{0}_{k}]}\neq{\mathbb{P}}_{(\alpha^{0}_{k},v]}$.

Figure 1:  Classifier-based gain curves, using a random forest (left) and a
$k$-nearest neighbor classifier (right). The time series with 600 observations
of dimension 5 and change points at $t=200,400$ (green crosses) was simulated
with the challenging change in covariance setup described in Section 4.2. The
maximizer is marked in red.

###  3.2 The two-step search algorithm

Binary segmentation relies on the full grid search, where $G_{(u,v]}(s)$ is
evaluated for all ${s=u+1,\ldots,v}$, to find the maximizer of the gain
$G_{(u,v]}(s)$. In many traditional parametric settings, such as change in
mean, log-likelihoods of neighboring segments can be recovered using cheap
${\cal O}(1)$ updates, enabling change point detection with binary
segmentation in typically ${\cal O}(n\log(n))$ time.

For many classifiers, random forests in particular, such updates are
unavailable and the classifiers need to be recomputed from scratch, making the
computational cost of the grid search in binary segmentation prohibitive.
Similar computational issues arise in high-dimensional regression models.
Here, Kaul et al. (2019) propose a two-step approach. Starting with an initial
guess $s^{(0)}$, they fit a single high-dimensional linear regression for each
of the segments ${(u,s^{(0)}]}$ and ${(s^{(0)},v]}$ and generate an improved
estimate of the optimal split using the resulting residuals. They repeat the
procedure two times and show a result about consistency in the high-
dimensional regression setting with a single change point if the change point
is sufficiently close to the first guess $s^{(0)}$.

We apply a variant of the two-step search paired with the classifier log-
likelihood ratio for multiple change point scenarios. Instead of residuals, we
recycle the class probability predictions and the resulting classifier log-
likelihood ratios for individual observations from a single classifier fit
$\hat{p}_{\\{u,s^{(0)},v\\}}$, approximating

$G\left((x_{i})_{i=u+1}^{v}\mid\\{u,s,v\\},\hat{p}\right)\approx\sum_{i=u+1}^{s}\log(\frac{v-u}{s^{(0)}-u}\hat{p}_{\\{u,s^{(0)},v\\}}(x_{i})_{1})+\sum_{i=s+1}^{v}\log(\frac{v-u}{v-s^{(0)}}\hat{p}_{\\{u,s^{(0)},v\\}}(x_{i})_{2}).$

We call this the _approximate gain_ in the following. Note that the classifier
and normalization factors are fixed, but the summation varies with the split
$s$. Like Kaul et al. (2019), we compute this twice, using the first maximizer
of the approximate gain as a second guess $s^{(1)}$. This allows us to find
local maxima of the nonparametric gain (6) in a constant number of classifier
fits.

Kaul et al. (2019) only apply the two-step search in scenarios with a single
change point. In scenarios with multiple change points, as encountered in
binary segmentation, the class probability predictions from the initial
classifier fit might contain little information if
${\mathbb{P}}_{(u,s^{(0)}]}\approx{\mathbb{P}}_{(s^{(0)},v]}$. To avoid this,
we start with multiple initial guesses and select the split point
corresponding to the overall highest approximate gain as the second guess. The
resulting two-step algorithm, as implemented in changeforest, using three
initial guesses at the segment’s $\frac{1}{4},\frac{1}{2}$ and
$\frac{3}{4}$-quantiles, is presented in Algorithm 2. Proposition 4 suggests
good estimation performance when coupling the two-step search with our
classifier-based approach in a single change point scenario.

######  Proposition 4.

For the Bayes classifier $p^{\infty}$ and any initial guess $u<s^{(0)}<v$, the
expected approximate gain

$s\mapsto{\mathbb{E}}\left[\sum_{i=u+1}^{s}\log(\frac{v-u}{s^{(0)}-u}p_{\\{u,s^{(0)},v\\}}^{\infty}(X_{i})_{1})+\sum_{i=s+1}^{v}\log(\frac{v-u}{v-s^{(0)}}p_{\\{u,s^{(0)},v\\}}^{\infty}(X_{i})_{2})\right]$

is piecewise linear between the underlying change points $\alpha^{0}$. If
there is a single change point $a^{0}\in{(u,v]}$, the expected approximate
gain has a unique maximum at $s=a^{0}$.

Figure 2 shows approximate gain curves and the underlying class probability
predictions for a simulated dataset with two change points. The approximate
gain curves have a piecewise linear shape with kinks at the underlying change
points, as predicted by Proposition 4. One can also observe how the
distributions of the probability predictions change at the true change points.
For the second row, the choice of the initial guess $s^{(0)}=300$ is precisely
such that ${\mathbb{P}}_{(0,300]}={\mathbb{P}}_{(300,600]}$, and the
approximate gain curve is flat. Figure 4 shows approximate gain curves in
different stages of binary segmentation for a different simulated time series.
Again, the approximate gain curves are roughly piecewise linear with kinks at
the underlying change points.

Figure 2:  Approximate gain curves (left) from single random forest fits for
initial guesses $s^{(0)}=150$ (top), $300$ (second row), and $450$ (third row)
and the underlying out-of-bag class probability predictions (right). The last
row shows the approximate gain curve for the second guess $\hat{s}^{(1)}=400$.
The underlying data is the same as in Figure 1, with true change points at
$t=200,400$ (green crosses). The initial guess and the curve’s maximizer are
marked in blue and red vertical dashed lines.

###  3.3 Choice of classifier

We recommend pairing our classifier-based method with random forests (Breiman,
2001). We motivate this choice.

A classifier should satisfy three properties when paired with our change point
detection methodology to enjoy good estimation performance and computational
feasibility: (i) It needs to have a low dependence on hyperparameter tuning,
(ii) it needs to be able to generate unbiased class probability predictions,
and (iii) the training cost should be reasonable even for a large number of
samples. The hyperparameters optimal for classification can vary between
segments encountered in binary segmentation. Retuning them, as is done by
Londschien et al. (2021), is infeasible. Thus, the classifier needs to provide
suitable class probability predictions even for suboptimal hyperparameters. If
the class probability predictions have a strong bias due to overfitting, the
approximate gain curve will have a kink at the initial guess, see Figure 3.
For a weak underlying signal, the initial guess might even be the maximizer of
the curve, causing the two-step algorithm to fail in finding a true change
point. As for computational cost: Our optimization algorithm based on binary
segmentation and the two-step search massively reduces the number of
classifier fits required for change point estimation. Still, hundreds of fits
might be necessary. For this not to become restrictive, the classifier’s
training cost should not explode when the number of samples gets large.

Random forests are a perfect fit: They are known to be an off-the-shelf
classifier that fits most data without much hyperparameter tuning, they can
generate unbiased probability estimates through out-of-bag prediction, and
their training complexity is nearly linear in the sample size for a fixed
maximal depth.

We highlight another attractive attribute of random forests: They are
applicable to a wide range of types of data. The need for nonparametric change
point detection arises primarily when little is known about the data to be
analyzed. Else, a smartly chosen parametric method might be better suited. In
such situations, classifiers with broad applicability are optimal. As a tree-
based method, random forests are invariant to feature rescaling and robust to
the inclusion of irrelevant features. Due to these properties, random forests
are popular as black box or off-the-shelf learning algorithms and thus well
suited for many datasets encountered in nonparametric change point detection.

While random forests are our classifier of choice, we believe that our method,
including the two-step search, can be paired with any classification algorithm
that can produce unbiased class probability predictions. As an illustration,
we pair our approach with the Euclidean norm $k$-nearest neighbor classifier.
It has a single hyperparameter $k$, which we set to the square root of the
segment length. Furthermore, the algorithm can be adapted to recover leave-
one-out cross-validated probability estimates. Our implementation has a
computational complexity of ${\cal O}(n^{2})$ to compute all pairwise
distances once, which we can recycle for all segments. More efficient
implementations exist (Arya et al., 1998). However, the quadratic runtime and
memory requirements are not an issue for the reasonably sized datasets we use
for benchmarking.

Figure 3:  In-sample (not out-of-bag) class probability predictions of a
random forest (right) and the corresponding approximate gain curve (left).
Predictions are from a single random forest fit for the initial guess
$s^{(0)}=450$ (blue line). The underlying data is the same as in Figures 1 and
2, with true change points at $t=200,400$ (green crosses). The maximizer is
marked in red. The predictions show a clear bias, dominating the underlying
signal (kink around the true change point at $t=400$), and the approximate
gain curve has a peak at the initial guess. Compare this to the third row of
Figure 2, where out-of-bag class probability predictions were used, resulting
in a maximum around a true change point at $t=400$.

###  3.4 Model selection

So far, we have discussed how to accurately and efficiently estimate the
location of change points in a signal. We now discuss how to decide whether to
keep or discard a candidate split.

Separating true change points from false ones is crucial for good estimation
performance and is a task as complex as finding candidate splits. Two
approaches for model selection are particularly popular in change point
detection: thresholding and permutation tests. For thresholding, a candidate
split is kept if its inclusion results in an increase of the log-likelihood at
least as large as the threshold $\gamma$. Thresholding is often applied in
parametric settings, where BIC-type criteria for the minimal gain to split can
be derived with asymptotic theory, see for example Yao (1988). Permutation
tests are an alternative to thresholding when such asymptotic theory is
unavailable. Here, a statistic gathered for a candidate split is compared to
values of the statistic gathered after the underlying data was permuted, see
for example Matteson and James (2014).

There is no asymptotic theory available to suggest a threshold for the
(approximate) nonparametric classifier log-likelihood ratio presented in
Section 2. Meanwhile, permutation tests require a high number of classifier
fits in each splitting step, with a severe effect on the overall computational
cost. We apply a leaner pseudo-permutation test. Instead of permuting the
observations and refitting the classifier, we only permute the classifier’s
predictions gathered in the first step of the two-step search, and thus the
classifier log-likelihood ratios. We then compute the approximate gains after
shuffling and take their maxima. In the absence of change points, there is no
underlying signal for the classifier to learn other than the relative class
sizes. We thus expect the distribution of the classifier log-likelihood ratios
to be approximately invariant under permutations. This allows us to compare
the maximal gain obtained in the first step of the two-step search to the
maximal gains computed after permuting. Algorithm 3 summarizes the procedure.

The significance threshold for the pseudo-permutation test is a tuning
parameter rather than a valid significance level. We choose 0.02 by default.
This will typically lead to many false positives due to multiple testing, as
the permutation test is applied to each segment encountered in binary
segmentation. We use 199 permutations. We analyze the empirical false positive
rate of the pseudo-permutation test in Section 4.7.

###  3.5 Our proposal

We present changeforest in Algorithm 1.

We comment on two details: We use $\log_{\eta}(x)$ instead of the natural
logarithm and rescale probabilities with
$1\,/\,\pi_{i}^{\\{u,s,v\\}}=\frac{v-u-1}{s-u-\mathds{1}_{\\{i\leqslant
s\\}}}$ instead of $\frac{v-u}{s-u}$. We introduce
$\log_{\eta}(x)\coloneqq\log((1-\eta)x+\eta)$ to use instead of $\log(x)$ in
the definition of the classifier log-likelihood ratio. Class probability
predictions can take values 0 and 1, prohibiting us from using the logarithm
directly to compute individual classifier log-likelihood ratios. With
$\eta=\exp(-6)$, we effectively cap individual classifier log-likelihood
ratios from below by $-6$. Concerning the rescaling: the out-of-bag
probability predictions of a random forest behave like those from leave-one-
out cross-validation. The prediction for observation $x_{i}$ was generated by
a forest that was trained on the $v-u-1$ observations
$(x_{j})_{j\in\\{u+1,\ldots,v\\}\backslash\\{i\\}}$. If $i\leqslant s$, then
$s-u-1$ of these observations belong to the first class, and if $i>s$, $s-u$
do. Thus, in the absence of any change points, we expect
$\hat{p}_{\\{u,s,v\\}}(x_{i})\approx\frac{s-u-\mathds{1}_{\\{i\leqslant
s\\}}}{v-u-1}=\pi_{i}^{\\{u,s,v\\}}$ and thus scale predictions by the inverse
of $\pi_{i}^{\\{u,s,v\\}}$ before taking the logarithm to recover the
classifier log-likelihood ratio estimates.

Our algorithm is visualized in Figure 4, where we display the approximate gain
curves of the two-step search for segments encountered in binary segmentation.

Algorithm 1 changeforest

Input:  Observations $(x_{i})_{i=1}^{n}$, a classifier $\hat{p}$, and a
minimum relative segment length $\delta>0$.

Output:  Change point estimate
$\hat{\alpha}\leftarrow\textrm{BinarySegmentation}((x_{i})_{i=1}^{n},\hat{p},\delta)$.

function $\textrm{BinarySegmentation}((x_{i})_{i=u+1}^{v},\hat{p},\delta)$

if $v-u<2\delta n$ then

return $\varnothing$

end if

$\hat{s},((\ell_{i,k,j})_{i=u+1,\ldots,v}^{k=1,2})_{j=1,2,3}\leftarrow\textrm{TwoStepSearch}((x_{i})_{i=u+1}^{v},\hat{p},\delta)$
$\triangleright$ c.f. Algorithm 2

$q=\textrm{ModelSelection}(((\ell_{i,k,j})_{i=u+1,\ldots,v}^{k=1,2})_{j=1,2,3},\delta)$
$\triangleright$ c.f. Algorithm 3

if $q\leqslant 0.02$ then

$\hat{\alpha}_{\textrm{left}}\leftarrow\textrm{BinarySegmentation}((x_{i})_{i=u+1}^{\hat{s}},\hat{p},\delta)$

$\hat{\alpha}_{\textrm{right}}\leftarrow\textrm{BinarySegmentation}((x_{i})_{i=\hat{s}+1}^{v},\hat{p},\delta)$

return
$\hat{\alpha}_{\textrm{left}}\cup\\{\hat{s}\\}\cup\hat{\alpha}_{\textrm{right}}$

else

return $\varnothing$

end if

end function

Algorithm 2 TwoStepSearch

Input: Observations $(x_{i})_{i=u+1}^{v}$, a classifier $\hat{p}$, and a
minimum relative segment length $\delta>0$.

Output: A single change point estimate $\hat{s}$ and classifier log-likelihood
ratios $((\ell_{i,k,j})_{i=u+1,\ldots,v}^{k=1,2})_{j=1,2,3}$ for model
selection.

Step 1:  Define
$(s^{(0)}_{1},s^{(0)}_{2},s^{(0)}_{3})=(\lfloor\frac{3u+v}{4}\rfloor,\lfloor\frac{v+u}{2}\rfloor,\lfloor\frac{u+3v}{4}\rfloor)$
and $\log_{\eta}(x)=\log((1-\eta)x+\eta)$ for $\eta=\exp(-6)$ and set
$\pi^{\\{u,s,v\\}}_{i}\leftarrow\frac{s-u-\mathds{1}_{\\{i\leqslant
s\\}}}{v-u-1}$. For $j=1,2,3$, train binary classifiers
$\hat{p}_{\\{u,s^{(0)}_{j},v\\}}$ on $({\bf X},{\bf
y}_{\\{u,s^{(0)}_{j},v\\}})=((x_{i})_{i=u+1}^{v},(1+\mathds{1}_{\\{i>s^{(0)}_{j}\\}})_{i=u+1}^{v})$
and let

$\ell_{i,1,j}=\log_{\eta}\left(\hat{p}_{\\{u,s^{(0)}_{j},v\\}}(x_{i})_{1}/\pi^{\\{u,s^{(0)}_{j},v\\}}_{i}\right),\
\ell_{i,2,j}=\log_{\eta}\left(\hat{p}_{\\{u,s^{(0)}_{j},v\\}}(x_{i})_{2}/(1-\pi^{\\{u,s^{(0)}_{j},v\\}}_{i})\right).$

Set  $\hat{s}^{(1)}=\operatorname*{argmax}_{s=u+1+\lceil\delta
n\rceil,\dotsc,v-\lceil\delta
n\rceil}\max_{j=1,2,3}\sum_{i=u+1}^{s}\ell_{i,1,j}+\sum_{i=s+1}^{v}\ell_{i,2,j}$

Step 2:  Train a binary classifier $\hat{p}_{\\{u,\hat{s}^{(1)},v\\}}$ and let

$\ell_{i,1}=\log_{\eta}\left(\hat{p}_{\\{u,s^{(1)},v\\}}(x_{i})_{1}/\pi^{\\{u,s^{(1)},v\\}}_{i}\right),\
\ell_{i,2}=\log_{\eta}\left(\hat{p}_{\\{u,s^{(1)},v\\}}(x_{i})_{2}/(1-\pi^{\\{u,s^{(1)},v\\}}_{i})\right).$

Set  $\hat{s}=\hat{s}^{(2)}=\operatorname*{argmax}_{s=u+1+\lceil\delta
n\rceil,\dotsc,v-\lceil\delta
n\rceil}\sum_{i=u+1}^{s}\ell_{i,1}+\sum_{i=s+1}^{v}\ell_{i,2}.$

Return:
$\left(\hat{s},((\ell_{i,k,g})_{i=u+1,\ldots,v}^{k=1,2})_{j=1,2,3}\right)$

Algorithm 3 ModelSelection

Input: Classifier log-likelihood ratios
$((\ell_{i,k,j})_{i=u+1,\ldots,v}^{k=1,2})_{j=1,2,3}$ from the two-step search
and a minimum relative segment length $\delta>0$.

Output: An approximate $p$-value from the pseudo-permutation test.

Step 1:  Recover the maximal gain encountered in the first step of the two-
step search.  $G_{0}=\max_{s=u+\lceil\delta n\rceil+1,\ldots,v-\lceil\delta
n\rceil}\ \max_{j=1,2,3}\
\sum_{i=u+1}^{s}l_{i,1,j}+\sum_{i=s+1}^{v}l_{i,2,j}.$

Step 2:  For random permutations
$\sigma_{l}:\\{u+1,\ldots,v\\}\rightarrow\\{u+1,\ldots,v\\}$ for
$l=1\ldots,199$, compute maximal gains after permuting likelihoods
$G_{l}=\max_{j=1,2,3}\ \max_{s=u+\lceil\delta n\rceil+1,\ldots,v-\lceil\delta
n\rceil}\
\sum_{i=u+1}^{s}l_{\sigma_{l}(i),1,j}+\sum_{i=s+1}^{v}l_{\sigma_{l}(i),2,j}.$

Return:  $\\#\\{l=0,\ldots,199\mid G_{l}\geqslant G_{0}\\}\,/\,200$

Figure 4:  Approximate gain curves from the last step of the two-step search
in the different steps of binary segmentation as encountered in changeforest.
The underlying dataset was simulated with the abalone setup, see Section 4.2
for details. The true change points at
$t=568,635,1122,1225,1914,2305,2508,2775,2817,2875,3134,3249,3306,3432$ are
marked with green crosses. Change point estimates, corresponding to maxima of
the approximate gain curves that are accepted in our model selection
procedure, are marked in red. The second guess $\hat{s}^{(1)}$ for the two-
step search is marked in blue. Current segment boundaries are marked with
black vertical lines. Here, changeforest estimates change points at
$t=568,636,1122,1229,1917,2305,2444,2798,2875,3132,3245,3306,3425$,
corresponding to an adjusted Rand index of 0.97.

###  3.6 The changeforest software package

A substantial part of this work is the changeforest software package,
implementing Algorithm 1 for random forests and $k$-nearest neighbor
classifiers. It also implements the parametric change in mean setup, paired
with binary segmentation, which we use as a baseline in the simulations.

The changeforest package is available for Python users on PyPI and conda-forge
(conda-forge community, 2015), for R users on conda-forge, and for Rust users
on crates.io. Its backend is implemented in the system programming language
Rust (Matsakis and Klock, 2014). For inquiries, installation instructions, a
tutorial, and more information, please visit
github.com/mlondschien/changeforest.

##  4 Empirical results

We present results from a simulation study comparing the empirical performance
of our classifier-based method changeforest to available multivariate
nonparametric competitors. The source code for the simulations, together with
guidance on reproducing tables and figures, is available at
github.com/mlondschien/changeforest-simulations. The repository can easily be
expanded to benchmark new methods provided by users.

###  4.1 Competing methods

We are aware of the following nonparametric methods for multivariate multiple
change point detection with existing reasonably efficient implementations:
Matteson and James’ (2014) _e-divisive_ (ECP), Lung-Yut-Fong et al.’s (2015)
_MultiRank_ , and kernel change point methods such as Arlot et al.’s (2019)
_KCP_.

ECP searches for a single change point minimizing energy distances within
segments. The method generalizes to multiple change point scenarios through
binary segmentation, and the significance for change points is assessed with a
permutation test. An implementation of ECP is available in the R-package ecp
(James and Matteson, 2015), which we use in the simulations.

MultiRank uses a rank-based multivariate homogeneity statistic combined with
dynamic programming. The significance of change points is assessed based on
asymptotic theory. We used code made available to us by the authors to run
simulations with MultiRank, which is included in the simulations repository.

Kernel change point methods minimize a within-segment average kernelized
dissimilarity measure. Arlot et al. (2019) propose KCP, an algorithm based on
dynamic programming. An efficient implementation can be found in the Python
package ruptures (Truong et al., 2020). We used the slope heuristic as
proposed by Arlot et al. (2019) for model selection.

For all packages, we used default parameters, if available. For ECP, this
results in $\alpha=1$ and using 199 permutations at significance level 0.05.
For KCP, we use a Gaussian kernel with a bandwidth of 0.1, optimal for the
simulated scenarios. See Section 4.6 and Table 9 for details. We supplied
information about the minimum relative segment length $\delta$, equal to 0.01
in the main simulations, to each method, either by informing about the minimum
number of observations per segment $\lceil\delta n\rceil$ or by informing
about the maximal number of change points $\lfloor 1/\delta\rfloor$.

Other multivariate nonparametric multiple change point detection methods that
we did not include in our simulations include Liu et al.’s (2013) _RuLSIF_ ,
Cabrieto et al.’s (2017) _Decon_ , Madrid–Padilla et al.’s (2021b) _MNWBS_ ,
and Zhang and Chen’s (2021) _gMulti_. For gMulti, no implementation has been
made available by the time of writing. Implementations for MNWBS and Decon
exist. However, these ran for over one hour for only medium-sized datasets
with 5000 observations. The implementation for RuLSIF is only available in
Matlab and does not include an automatic model selection procedure.

###  4.2 Simulation setups

We evaluate all methods on three parametric and six nonparametric scenarios.
We display their key characteristics below the header in Table 2.

As for parametric setups, we use the _change in mean_ (CIM) and _change in
covariance_ (CIC) setups from (Matteson and James, 2014, Section 4.3) and the
_Dirichlet_ setup from (Arlot et al., 2019, Section 6.1). For both works, the
authors present empirical results of their methods on datasets for three sets
of parameters. We use the set of parameters corresponding to medium
difficulty.

For the change in mean and change in covariance setups, we generate time
series of dimension $d=5$ with $n=600$ observations and change points at
$t=200,400$. Observations in the first and last segment are independently
drawn from a standard multivariate Gaussian distribution. In the change in
mean setup, entries in the second segment are i.i.d. Gaussian with a mean
shift of $\mu=2$ in each coordinate. In the change in covariance setup,
entries in the second segment are i.i.d. normal with mean zero and unit
variance, but with a covariance of $\rho=0.7$ between coordinates. Lastly, the
Dirichlet dataset consists of $n=1000$ observations and has dimension $d=20$.
Ten change points are located at $t=100,130,220,320,370,520,620,740,790,870$.
In each segment, the observations are distributed according to the Dirichlet
distribution, with the $d$ parameters drawn i.i.d. uniformly from $[0,0.2]$.

We also generate time series with nonparametric in-segment distributions. For
this, we use popular multiclass classification datasets and treat observations
corresponding to a single class as i.i.d. draws from a class-specific
distribution. We discard classes with fewer than $\delta n=\frac{n}{100}$
observations and shuffle and randomly concatenate the remaining classes.
Categorical variables are dummy encoded. Finally, each covariate gets
normalized by the median absolute deviation of absolute consecutive
differences, as is also used by Fryzlewicz (2014) and further discussed by
Kovács et al. (2020a). This normalization does not affect our random forest-
based method changeforest and rank-based MultiRank but increases performance
for distance-based competitors ECP and KCP.

We use the following datasets: The _iris_ flower dataset (Fisher, 1936)
contains 150 samples of four measurements from three species of the iris
flower. The _glass_ identification dataset (Evett and Spiehler, 1989),
initially motivated by criminological investigation, contains measurements of
chemical compounds for different types of glass. The _wine_ dataset (Cortez et
al., 2009) is the result of physicochemical analysis of both red and white
vinho verde from north-western Portugal. It includes quality scores from blind
testing, which we use as class labels. We encode the wine color as binary. The
Wisconsin _breast cancer_ dataset (Street et al., 1993) contains nominal
characteristics of cell nuclei computed from a fine needle aspirate of breast
tissue labeled as malignant or benign. The resulting simulated time series has
a single change point. The _abalone_ dataset (Waugh, 1995) contains easily
obtainable physical measurements of snails and their age, determined by
cutting through the shell cone and counting the number of rings. Two
categorical variables are dummy encoded. We use the number of rings as class
labels. The _dry beans_ dataset (Koklu and Ozkan, 2020) contains shape
measurements derived from images of seven different types of dry beans.

The change in mean and change in covariance and possibly also the abalone and
wine datasets feature similar or equal distributions in non-neighboring
segments. We expect these scenarios to be challenging for our methodology, see
Section 3.2 and Figure 2.

###  4.3 Performance measures

We mainly use the adjusted Rand index (Hubert and Arabie, 1985), a standard
measure to compare clusterings, to evaluate the goodness of fit of estimated
segmentations. Given two partitionings of $n$ observations, the Rand index
(Rand, 1971) is the number of agreements, pairs of observations either in the
same or in different subsets for both partitionings, divided by the total
number of pairs $n\choose 2$. The adjusted Rand index is the normalized
difference between the Rand index and its expectation when choosing partitions
randomly. Consequently, the adjusted Rand index can take negative values but
no values greater than one. Its expected value is zero for random class
assignments, but the expected value for random segmentations, namely class
assignments consistent with the time structure, is significantly higher. We
present a set of segmentations together with their adjusted Rand indices in
Table 1.

true segmentation estimated segmentation ARI $d_{H}$ comment 0, 17, 46, 55,
68, 144, 214 0, 17, 46, 55, 68, 144, 214 1.00 0.000 perfect fit 0, 15, 45, 55,
68, 142, 214 0.95 0.009 almost perfect fit 0, 46, 55, 68, 144, 214 0.95 0.079
undersegmentation 0, 17, 46, 55, 144, 214 0.89 0.061 undersegmentation 0, 17,
46, 55, 68, 80, 144, 214 0.91 0.056 oversegmentation 0, 17, 46, 55, 68, 100,
144, 214 0.83 0.150 oversegmentation 0, 50, 100, 150, 214 0.61 0.150 random
segmentation 0, 214 0.00 0.327 no segmentation

Table 1: Adjusted Rand indices (ARI) and Hausdorff distances ($d_{H}$) for
different segmentations. The true segmentation comes from the nonparametric
glass setup.

We believe that an adjusted Rand index above 0.95 corresponds to a
segmentation that is close to perfect. An index of 0.9 corresponds to a good
segmentation, possibly oversegmenting with a false positive close to a true
change point or missing one of two close change points. On the other hand, we
expect a segmentation with an adjusted Rand index significantly below 0.8 not
to provide much scientific value.

We additionally use the Hausdorff distance to evaluate change point estimates.
For two segmentations $\alpha,\alpha^{\prime}$, set
$d(\alpha,\alpha^{\prime})\coloneqq\frac{1}{n}\max_{a\in\alpha}\min_{a^{\prime}\in\alpha^{\prime}}|a-a^{\prime}|$.
Then $d(\alpha^{0},\hat{\alpha})$ is the largest relative distance of a true
change point in $\alpha^{0}$ to its closest counterpart in $\hat{\alpha}$, and
$d(\hat{\alpha},\alpha^{0})$ is the largest relative distance of a found
change point in $\hat{\alpha}$ to the closest entry in $\alpha^{0}$. Thus
$d(\alpha^{0},\hat{\alpha})$ especially penalizes undersegmentation while
$d(\hat{\alpha},\alpha)$ penalizes oversegmentation. Define the Hausdorff
distance as
$d_{H}(\alpha,\alpha^{\prime})\coloneqq\max(d(\alpha,\alpha^{\prime}),d(\alpha^{\prime},\alpha))$.
Table 1 also includes the Hausdorff distance for the selected segmentations.

###  4.4 Main simulation results

We present the results of our main simulation study in Table 2, where we
display average adjusted Rand indices over 500 simulations.

CIM CIC Dirichlet iris glass $n,d,K$ $600,5,3$ $600,5,3$ $1000,20,11$
$150,4,3$ $214,8,6$ change in mean 0.999 (0.003) 0.000 (0.000) 0.000 (0.000)
0.935 (0.082) 0.490 (0.171) changeforest 0.987 (0.035) 0.925 (0.122) 0.986
(0.019) 0.983 (0.040) 0.923 (0.072) changekNN 0.993 (0.023) 0.011 (0.077)
0.688 (0.245) 0.990 (0.031) 0.879 (0.089) ECP 0.994 (0.025) 0.430 (0.317)
0.836 (0.080) 0.991 (0.031) 0.608 (0.243) KCP 0.998 (0.012) 0.955 (0.057)
0.870 (0.078) 0.994 (0.023) 0.744 (0.234) MultiRank 1.000 (0.001) 0.003
(0.036) 0.966 (0.025) 0.569 (0.256) 0.598 (0.093)

breast cancer abalone wine dry beans average worst $n,d,K$ $699,9,2$
$4066,9,15$ $6462,12,5$ $13611,16,7$ change in mean 0.791 (0.027) 0.885
(0.055) 0.983 (0.018) 0.991 (0.014) 0.675 (0.067) 0.000 changeforest 0.983
(0.074) 0.932 (0.049) 0.993 (0.023) 0.997 (0.010) 0.968 (0.061) 0.923
changekNN 0.905 (0.166) 0.880 (0.097) 0.972 (0.040) 0.996 (0.011) 0.813
(0.113) 0.011 ECP 0.995 (0.021) 0.865 (0.064) 0.982 (0.032) 1.000 (0.002)
0.854 (0.136) 0.430 KCP 0.996 (0.033) 0.916 (0.043) 0.988 (0.011) 1.000
(0.003) 0.940 (0.087) 0.744 MultiRank 0.923 (0.119) 0.813 (0.115) 0.957
(0.032) 0.729 (0.114) 0.003

Table 2:  Average adjusted Rand indices over 500 simulations, with standard
deviations in parentheses. Optimal scores are marked in bold. The MultiRank
method was unable to produce results for the dry beans simulation setup.
Average and worst denote the average and worst performances across all setups.

Our method, changeforest, scores above 0.9 on average for each simulation
scenario studied. All other methods have an average score lower than 0.75 for
at least one scenario. changeforest is the best performing method for four
setups: Dirichlet, glass, abalone, and wine. Furthermore, in four further
setups, change in mean, iris, breast cancer, and dry beans, our method
changeforest scores above 0.975 on average, corresponding to an almost perfect
segmentation. The only setup for which changeforest does not perform best or
above 0.975 is change in covariance. Here, changeforest scores second-best
behind KCP with an average adjusted Rand index of 0.925, and no other method
scores above 0.5. We compute the average score for each method to combine the
results for the nine simulation scenarios. changeforest scores highest.
Overall, changeforest is widely applicable, with robust performance on all
parametric and nonparametric setups.

We include Table 5 with median Hausdorff distances in the Appendix. The
results are similar, with changeforest being among the methods with the best
score, except for change in covariance, where it achieves a median Hausdorff
distance of 0.017, behind KCP.

Table 6 in the Appendix shows the average number of estimated change points
for each method and simulation setup. changeforest slightly oversegments for
most setups, but not for the glass and wine setups.

Figure 7 in the Appendix shows histograms of change point estimates for the
Dirichlet setup, where the underlying change points are constant.
Corresponding to high adjusted Rand indices of 0.986 and 0.966, change point
estimates by changeforest and MultiRank are sharply concentrated around the
true change point locations. Meanwhile, estimates by changekNN, ECP, and KCP
are visually more spread out, corresponding to worse adjusted Rand indices
smaller than 0.9. Lastly, the change in mean method fails to produce any
change point estimates.

We present average computation times for each method in Table 3. Simulations
were run on eight Intel Xeon 2.3 GHz cores with 4 GB of RAM available per core
(32 GB in total). The simple change in mean method is blazingly fast. On the
largest simulation setup dry beans with $n=13611$ and $d=16$, the computation
time was around 0.002 seconds on average. changeforest is the fastest
nonparametric method on the dry beans dataset, requiring around 2.6 seconds on
average. On the other hand, ECP is prohibitively slow, requiring more than one
hour to estimate change points on the dry beans dataset on average. The code
for MultiRank is experimental and broke for the dry beans dataset due to the
multiple dummy encoded columns. We provide further analysis of the
computational efficiency of available implementations in Section 4.5.

CIM CIC Dirichlet iris glass breast abalone wine dry cancer beans $n,d$ 600, 5
600, 5 1000, 20 150, 4 214, 8 699, 9 4066, 9 6462, 12 13611, 16 change in mean
<0.01 <0.01 <0.01 <0.01 <0.01 <0.01 <0.01 <0.01 <0.01 changeforest 0.08 0.09
0.43 0.03 0.06 0.07 0.86 0.89 2.6 changekNN 0.04 0.04 0.15 <0.01 0.01 0.05 1.8
4.7 24 ECP 4.2 2.9 18 0.37 0.72 4.3 386 741 5356 KCP 0.08 0.08 0.18 0.02 0.04
0.10 2.5 6.6 31 MultiRank 0.77 0.76 2.1 0.10 0.17 1.0 34 90

Table 3:  Average computation times in seconds on 8 Intel Xeon 2.3 GHz cores.

###  4.5 Performance for a varying number of segments and sample size

The performance of change point estimation methods depends on an interplay
between the number of observations, the size of distributional shifts between
segments, and the number of change points. In Table 2 we see that for some
datasets, such as CIM, iris, and dry beans, all methods achieve average
adjusted Rand indices above 0.975. Here, we believe that the signal-to-noise
ratio is so high that change point estimation is _easy_ , and the difference
between the average adjusted Rand index and the perfect score 1.0 is mostly a
result of random false positives. On the other hand, the change in covariance
and glass setups appear truly difficult, clearly separating methods that
perform well (above 0.9) and those performing not much better than random
guessing (below 0.8).

We introduce two simulation settings that let us generate time series of any
length $n$ with any number of change points $K$, allowing us to control the
underlying signal-to-noise ratio. We draw
$\tilde{N}_{k}\sim\operatorname{Exp}(1)$ i.i.d. for $k=1,\ldots,K$ and define
$N_{i}=\frac{1}{10K}+\frac{0.9}{\sum_{j=1}^{K}\tilde{N}_{j}}\tilde{N}_{k}$,
resulting in a minimum relative segment length of $\delta=\frac{1}{10K}$. We
use $\textrm{round}(nN_{k})$ as segment lengths, where we round such that
$\sum_{k=1}^{K}\textrm{round}(nN_{k})=n$.

We present two simulation setups: One based on the Dirichlet distribution and
one on the dry beans dataset. We construct the former as in Section 4.2. For
each segment, we draw i.i.d. observations from a Dirichlet distribution of
dimension $d=20$, after sampling its parameters uniformly from $[0,0.2]$. To
generate time series of arbitrary length based on the dry beans dataset, we
first normalize each covariate to have an average within-class variance of
one. Then, for each segment, we draw a class label unequal to that of the
previous segment and draw $\textrm{round}(nN_{k})$ observations with
replacement from the corresponding class. Finally, we add i.i.d. standard
Gaussian noise to each observation. We simulate 500 datasets with $K=20,80$
and vary $n=250,354,500,707,1000,\dotsc,64000$ for each setup. The mean
adjusted Rand indices of change point estimates from our and competing methods
are displayed in Figure 5. Each method was supplied with the minimum relative
segment length $\delta=\frac{1}{10K}$.

Figure 5:  Evolution of the average adjusted Rand index for methods by the
number of observations for 20 and 80 segments. The $y$-axis was scaled to
better pronounce values in $[0.6,1]$. Two times the standard deviation of the
mean score is marked with vertical bars. All methods were run until
computational time got prohibitive (ECP) or 32 GB of memory was not sufficient
for computation (KCP, MultiRank, changekNN).

The change in mean method does not perform well for either setup. Our method
changeforest performs best, except for very small sample sizes. For the dry
beans simulation setup when $n/K<20$, ECP performs best, but not much better
than 0.8. For the Dirichlet simulation setup for $K=20,n=354$ and
$K=80,n=1000$, MultiRank performs best. Surprisingly, KCP performs very badly
for the dry beans-based simulation, not selecting any change points at all.
This might be due to the higher number of change points compared to the main
simulations, or the additional Gaussian noise added.

We display the average computational time of the different methods in Figure
6, together with dashed grey lines corresponding to linear and quadratic time
complexity. Again, the parametric change in mean method is very fast. For both
change in mean and changeforest, time complexity appears to scale
approximately linear in the number of observations. Meanwhile, changekNN,
MultiRank, and KCP scale approximately quadratic, while ECP appears to have
worse than quadratic time complexity. changeforest is the quickest
nonparametric method for $n\gtrsim 5000$.

Figure 6:  Average time for change point detection on 8 Intel Xeon 2.3 GHz
cores with 4 GB of RAM available per core for the Dirichlet simulation setup.
All methods were run until computational time got prohibitive (ECP) or the in
total 32 GB of RAM we not sufficient (KCP, MultiRank, changekNN).

###  4.6 Dependence on hyperparameters

We investigate how the hyperparameters of random forests affect changeforest’s
performance. For this, we vary the three most important parameters for random
forests: the number of trees from $20$ to $100$ (default) to $500$, the
maximal depth of individual trees from $2$ to $8$ (default) to infinity, and
the number of features considered at each split (mtry) from $1$ to $\sqrt{d}$
(default, rounded down) to $d$. The average adjusted Rand indices for each set
of parameters for the simulation setups introduced in Section 4.2 are in Table
7 in the Appendix. Average computational costs are in Table 8. We summarize
the results here.

We observe that choosing even very extreme parameters, such as $\textrm{number
of trees}=20$, $\textrm{maximal depth}=2$, and $\textrm{mtry}=1,$ yields
acceptable change point estimation performance for simple simulation setups.
With these parameters, changeforest performs on average better than ECP,
MultiRank, and changekNN, our method paired with a $k$-nearest neighbor
classifier. It even outperforms the default configuration (100, 8, $\sqrt{d}$)
for the iris, breast cancer, and wine simulation setups, possibly due to
selecting fewer false positives because of higher noise in the predictions. On
the other hand, increasing the number of trees from 100 to 500 and splitting
individual trees to purity does not significantly increase the performance
compared to the default values while increasing computational cost fivefold.

This low dependence on the choice of tuning parameters is in contrast to KCP,
where the bandwidth of the Gaussian kernel has to be fine-tuned for optimal
performance. Table 9 in the Appendix shows the estimation performance of KCP
for different kernels and bandwidths.

###  4.7 False positive rates

We empirically analyze our model selection procedure. For this, we select the
largest class of each dataset from Section 4.2, shuffle the observations
before applying the change point detection methods and collect the percentage
out of 2500 simulations where at least one change point was detected. See
Table 4 for the results.

dataset CIM CIC Dirichlet iris glass breast abalone wine dry worst cancer
beans change in mean 0.12 0.12 100.0 27.04 91.44 3.44 7.68 1.20 1.40 100.00
changeforest 3.36 3.36 3.76 5.00 3.80 3.80 3.00 3.52 2.48 5.00 changekNN 2.44
2.44 1.80 3.64 2.40 29.80 3.28 6.80 2.04 29.80 ECP 5.76 6.20 5.08 6.80 5.72
4.76 4.76 5.16 5.16 6.80 KCP 0.00 0.00 0.00 3.04 0.32 5.24 11.32 0.08 0.04
11.32 MultiRank 0.20 0.20 12.08 93.64 90.44 100.0 100.0 100.0 100.00

Table 4:  Percentage of simulations where at least one change point was
detected in a homogeneous dataset. For each dataset and method, 2500
simulations were performed.

As our model selection procedure is not a true permutation test, we observe
more false positives than expected, given the approximate significance
threshold of $0.02$ for both changeforest and the $k$-nearest neighbor method
changekNN. However, changeforest has a similar false positive rate across all
simulation setups. The methods based on heuristics (KCP) or asymptotic theory
(MultiRank, change in mean) have varying degrees of false positive rates by
dataset. ChangekNN, which uses the same model selection procedure as
changeforest, has a false positive rate of almost 30% for the breast cancer
dataset.

##  5 Conclusion

Our proposed changeforest algorithm uses random forests to efficiently and
accurately perform multiple change point detection across a wide range of
dataset types. Motivated by parametric change point detection methodology, the
algorithm optimizes a nonparametric classifier-based log-likelihood ratio.
When coupled with random forests, it performs very well or better than
competitors for all simulation setups studied. The computation time of
changeforest scales almost linearly with the number of observations, enabling
nonparametric change point detection even for very large datasets. The
changeforest software package is available to R, Python, and Rust users. See
github.com/mlondschien/changeforest for installation instructions.

As argued in Section 3.3, random forests are a natural choice as a classifier
to pair with our methodology due to the possibility of retrieving unbiased
out-of-bag class probability predictions, the low dependence on
hyperparameters, and almost linear time complexity. However, in principle, any
classifier yielding class probability predictions can be used to generate
classifier log-likelihood ratios, which we demonstrated by pairing our
methodology with a $k$-nearest neighbor classifier. Other choices, such as
neural networks, are also possible. Here, using cross-validation to generate
unbiased class probability predictions would be necessary.

The performance of the changeforest algorithm could be improved by pairing it
with wild binary segmentation of Fryzlewicz (2014) or seeded binary
segmentation of Kovács et al. (2020b). This would help particularly in
scenarios with very short segments and equal or similar distributions in
surrounding segments. However, this would require changes to the model
selection procedure presented to avoid overfitting. Since we already observe
very good performance with binary segmentation, we do not include results from
wild or seeded binary segmentation here.

Acknowledgments

Malte Londschien is supported by the ETH Foundations of Data Science and the
ETH AI Center. Solt Kovács and Peter Bühlmann are supported by the European
Research Council (ERC) under the European Union’s Horizon 2020 research and
innovation programme (Grant agreement No. 786461 CausalStats - ERC-2017-ADG).

Authors’ Contributions

S.K. conceived and co-supervised the project. M.L. and S.K. developed the
methodology. M.L. designed and implemented the changeforest software package,
designed, implemented, and ran the simulations, developed the theory, and
wrote the manuscript. P.B. co-supervised the project and provided feedback and
ideas during the development. S.K. and P.B. provided feedback on the
manuscript.

## References

  * Arlot et al. (2019) Sylvain Arlot, Alain Celisse, and Zaid Harchaoui.  A kernel multiple change-point algorithm via model selection.  _Journal of Machine Learning Research_ , 20(162), 2019.
  * Arya et al. (1998) Sunil Arya, David M. Mount, Nathan S. Netanyahu, Ruth Silverman, and Angela Y. Wu.  An optimal algorithm for approximate nearest neighbor searching fixed dimensions.  _Journal of the ACM_ , 45(6):891–923, 1998.
  * Breiman (2001) Leo Breiman.  Random forests.  _Machine Learning_ , 45(1):5–32, 2001.
  * Cabrieto et al. (2017) Jedelyn Cabrieto, Francis Tuerlinckx, Peter Kuppens, Mariel Grassmann, and Eva Ceulemans.  Detecting correlation changes in multivariate time series: A comparison of four non-parametric change point detection methods.  _Behavior Research Methods_ , 49(3):988–1005, 2017.
  * Carlstein (1988) Edward Carlstein.  Nonparametric change-point estimation.  _The Annals of Statistics_ , 16(1):188–197, 1988\.
  * Chen and Zhang (2015) Hao Chen and Nancy Zhang.  Graph-based change-point detection.  _The Annals of Statistics_ , 43(1):139–176, 2015\.
  * Cleynen and Lebarbier (2014) Alice Cleynen and Émilie Lebarbier.  Segmentation of the Poisson and negative binomial rate models: a penalized estimator.  _ESAIM: Probability and Statistics_ , 18:750–769, 2014\.
  * conda-forge community (2015) conda-forge community.  The conda-forge Project: Community-based Software Distribution Built on the conda Package Format and Ecosystem, July 2015.
  * Cortez et al. (2009) Paulo Cortez, António Cerdeira, Fernando Almeida, Telmo Matos, and José Reis.  Modeling wine preferences by data mining from physicochemical properties.  _Decision support systems_ , 47(4):547–553, 2009\.
  * Dümbgen (1991) Lutz Dümbgen.  The asymptotic behavior of some nonparametric change-point estimators.  _The Annals of Statistics_ , 19(3):1471–1495, 1991.
  * Evett and Spiehler (1989) Ian W. Evett and Ernest J. Spiehler.  Rule induction in forensic science.  In _Knowledge Based Systems_ , pages 152–160. Online Publications, 1989.
  * Fisher (1936) Ronald A. Fisher.  The use of multiple measurements in taxonomic problems.  _Annals of Eugenics_ , 7(2):179–188, 1936.
  * Frick et al. (2014) Klaus Frick, Axel Munk, and Hannes Sieling.  Multiscale change point inference.  _Journal of the Royal Statistical Society: Series B_ , 76(3):495–580, 2014.
  * Friedman (2004) Jerome Friedman.  On multivariate goodness-of-fit and two-sample testing.  Technical report, Stanford University, 2004.
  * Fryzlewicz (2014) Piotr Fryzlewicz.  Wild binary segmentation for multiple change-point detection.  _The Annals of Statistics_ , 42(6):2243–2281, 2014.
  * Garreau and Arlot (2018) Damien Garreau and Sylvain Arlot.  Consistent change-point detection with kernels.  _Electronic Journal of Statistics_ , 12(2):4440–4486, 2018.
  * Hediger et al. (2022) Simon Hediger, Loris Michel, and Jeffrey Näf.  On the use of random forest for two-sample testing.  _Computational Statistics & Data Analysis_, 170:170435, 2022.
  * Hotz et al. (2013) Thomas Hotz, Ole Mathis Schütte, Hannes Sieling, Tatjana Polupanow, Ulf Diederichsen, Claudia Steinem, and Axel Munk.  Idealizing ion channel recordings by a jump segmentation multiresolution filter.  _IEEE Transactions on NanoBioscience_ , 12(4):376–386, 2013.
  * Hubert and Arabie (1985) Lawrence Hubert and Phipps Arabie.  Comparing partitions.  _Journal of Classification_ , 2(1):193–218, 1985\.
  * James and Matteson (2015) Nicholas James and David Matteson.  ecp: An R package for nonparametric multiple change point analysis of multivariate data.  _Journal of Statistical Software_ , 62(7):1–25, 2015.
  * Kaplan et al. (2001) Alexander Kaplan, Joachim Röschke, Boris Darkhovsky, and Juergen Fell.  Macrostructural EEG characterization based on nonparametric change point segmentation: application to sleep analysis.  _Journal of Neuroscience Methods_ , 106(1):81–90, 2001.
  * Kaul et al. (2019) Abhishek Kaul, Venkata K. Jandhyala, and Stergios B. Fotopoulos.  An efficient two step algorithm for high dimensional change point regression models without grid search.  _Journal of Machine Learning Research_ , 20:1–40, 2019\.
  * Killick et al. (2012) Rebecca Killick, Paul Fearnhead, and Idris A. Eckley.  Optimal detection of changepoints with a linear computational cost.  _Journal of the American Statistical Association_ , 107(500):1590–1598, 2012.
  * Kim et al. (2005) Chang-Jin Kim, James C. Morley, and Charles R. Nelson.  The structural break in the equity premium.  _Journal of Business & Economic Statistics_, 23(2):181–191, 2005.
  * Kim and Siegmund (1989) Hyune-Ju Kim and David Siegmund.  The likelihood ratio test for a change-point in simple linear regression.  _Biometrika_ , 76(3):409–423, 1989.
  * Koklu and Ozkan (2020) Murat Koklu and Ilker Ali Ozkan.  Multiclass classification of dry beans using computer vision and machine learning techniques.  _Computers and Electronics in Agriculture_ , 174:105507, 2020.
  * Kovács et al. (2020a) Solt Kovács, Housen Li, and Peter Bühlmann.  Seeded intervals and noise level estimation in change point detection: a discussion of Fryzlewicz (2020).  _Journal of the Korean Statistical Society_ , 49(4):1081–1089, 2020a.
  * Kovács et al. (2020b) Solt Kovács, Housen Li, Peter Bühlmann, and Axel Munk.  Seeded binary segmentation: A general methodology for fast and optimal change point detection.  _arXiv:2002.06633_ , 2020b.  To appear in Biometrika.
  * Kovács et al. (2020c) Solt Kovács, Housen Li, Lorenz Haubner, Axel Munk, and Peter Bühlmann.  Optimistic search strategy: Change point detection for large-scale data via adaptive logarithmic queries.  _arXiv:2010.10194_ , 2020c.
  * Liu et al. (2013) Song Liu, Makoto Yamada, Nigel Collier, and Masashi Sugiyama.  Change-point detection in time-series data by relative density-ratio estimation.  _Neural Networks_ , 43:72–83, 2013.
  * Londschien et al. (2021) Malte Londschien, Solt Kovács, and Peter Bühlmann.  Change-point detection for graphical models in the presence of missing values.  _Journal of Computational and Graphical Statistics_ , 30(3):768–779, 2021.
  * Lopez-Paz and Oquab (2017) David Lopez-Paz and Maxime Oquab.  Revisiting classifier two-sample tests.  In _International Conference on Learning Representations_ , 2017.
  * Lung-Yut-Fong et al. (2015) Alexandre Lung-Yut-Fong, Céline Lévy-Leduc, and Olivier Cappé.  Homogeneity and change-point detection tests for multivariate data using rank statistics.  _Journal de la Société Française de Statistique_ , 156(4):133–162, 2015.
  * Madrid–Padilla et al. (2021a) Oscar–Hernan Madrid–Padilla, Yi Yu, Daren Wang, and Alessandro Rinaldo.  Optimal nonparametric change point analysis.  _Electronic Journal of Statistics_ , 15(1):1154 – 1201, 2021a.
  * Madrid–Padilla et al. (2021b) Oscar–Hernan Madrid–Padilla, Yi Yu, Daren Wang, and Alessandro Rinaldo.  Optimal nonparametric multivariate change point detection and localization.  _IEEE Transactions on Information Theory_ , 68(3):1922–1944, 2021b.
  * Matsakis and Klock (2014) Nicholas D. Matsakis and Felix S. Klock.  The rust language.  In _ACM SIGAda Ada Letters_ , volume 34, pages 103–104. ACM, 2014\.
  * Matteson and James (2014) David S. Matteson and Nicholas A. James.  A nonparametric approach for multiple change point analysis of multivariate data.  _Journal of the American Statistical Association_ , 109(505):334–345, 2014.
  * McCulloch and Pitts (1943) Warren S. McCulloch and Walter Pitts.  A logical calculus of the ideas immanent in nervous activity.  _The Bulletin of Mathematical Biophysics_ , 5(4):115–133, 1943.
  * Olshen et al. (2004) Adam B. Olshen, E. Seshan Venkatraman, Robert Lucito, and Michael Wigler.  Circular binary segmentation for the analysis of array-based DNA copy number data.  _Biostatistics_ , 5(4):557–572, 2004.
  * Page (1954) Ewan S. Page.  Continuous inspection schemes.  _Biometrika_ , 41(1/2):100–115, 1954.
  * Page (1955) Ewan S. Page.  A test for a change in a parameter occurring at an unknown point.  _Biometrika_ , 42(3/4):523–527, 1955.
  * Pein et al. (2017) Florian Pein, Hannes Sieling, and Axel Munk.  Heterogeneous change point inference.  _Journal of the Royal Statistical Society: Series B_ , 79(4):1207–1227, 2017.
  * Pettitt (1979) Anthony N. Pettitt.  A non-parametric approach to the change-point problem.  _Journal of the Royal Statistical Society: Series C_ , 28(2):126–135, 1979.
  * Picard et al. (2005) Franck Picard, Stephane Robin, Marc Lavielle, Christian Vaisse, and Jean-Jacques Daudin.  A statistical approach for array CGH data analysis.  _BMC Bioinformatics_ , 6(1):1–14, 2005.
  * Rand (1971) William M. Rand.  Objective criteria for the evaluation of clustering methods.  _Journal of the American Statistical association_ , 66(336):846–850, 1971.
  * Reeves et al. (2007) Jaxk Reeves, Jien Chen, Xiaolan L. Wang, Robert Lund, and Qi Qi Lu.  A review and comparison of changepoint detection techniques for climate data.  _Journal of Applied Meteorology and Climatology_ , 46(6):900–915, 2007.
  * Roy et al. (2017) Sandipan Roy, Yves Atchadé, and George Michailidis.  Change point estimation in high dimensional Markov random-field models.  _Journal of the Royal Statistical Society: Series B_ , 79(4):1187–1206, 2017.
  * Shen and Zhang (2012) Jeremy J. Shen and Nancy R. Zhang.  Change-point model on nonhomogeneous Poisson processes with application in copy number profiling by next-generation DNA sequencing.  _The Annals of Applied Statistics_ , 6(2):476–496, 2012.
  * Street et al. (1993) Nick Street, William H. Wolberg, and Olvi L. Mangasarian.  Nuclear feature extraction for breast tumor diagnosis.  In _Biomedical Image Processing and Biomedical Visualization_ , volume 1905. SPIE, 1993.
  * Truong et al. (2020) Charles Truong, Laurent Oudre, and Nicolas Vayatis.  Selective review of offline change point detection methods.  _Signal Processing_ , 167:107299, 2020.
  * Vostrikova (1981) Lyudmila Y. Vostrikova.  Detecting ’disorder’ in multidimensional random processes.  _Soviet Mathematics Doklady_ , 24:55–59, 1981.
  * Wang et al. (2021) Daren Wang, Yi Yu, and Alessandro Rinaldo.  Optimal covariance change point localization in high dimensions.  _Bernoulli_ , 27(1):554–575, 2021.
  * Wang and Samworth (2018) Tengyao Wang and Richard J. Samworth.  High dimensional change point estimation via sparse projection.  _Journal of the Royal Statistical Society: Series B_ , 80(1):57–83, 2018.
  * Waugh (1995) Samuel G. Waugh.  _Extending and benchmarking Cascade-Correlation: extensions to the Cascade-Correlation architecture and benchmarking of feed-forward supervised artificial neural networks_.  PhD thesis, University of Tasmania, 1995.
  * Yao (1988) Yi-Ching Yao.  Estimating the number of change-points via Schwarz’ criterion.  _Statistics & Probability Letters_, 6(3):181–189, 1988.
  * Zhang and Chen (2021) Yuxuan Zhang and Hao Chen.  Graph-based multiple change-point detection.  _arXiv:2110.01170_ , 2021.
  * Zou et al. (2014) Changliang Zou, Guosheng Yin, Long Feng, and Zhaojun Wang.  Nonparametric maximum likelihood approach to multiple change-point problems.  _The Annals of Statistics_ , 42(3):970–1002, 2014\.

##  Appendix A Proofs

See 2

###### Proof.

Let $u<v$ and define
$\kappa({(u,v]})\coloneqq\\{k\mid{(\alpha^{0}_{k\shortminus
1},\alpha^{0}_{k}]}\cap{(u,v]}\neq\varnothing\\}$. If there exists some $k$
such that $\alpha_{k\shortminus 1}^{0}\leqslant u<v\leqslant\alpha_{k}^{0}$,
then $\kappa({(u,v]})=\\{k\\}$ and $p_{(u,v]}=p_{(\alpha_{k\shortminus
1}^{0},\alpha_{k}^{0}]}$. Consequently,

${\mathbb{E}}\left[\sum_{i=u+1}^{v}\log(\frac{p_{(u,v]}(X_{i})}{p_{(0,n]}(X_{i})})\right]=(v-u)D_{KL}({\mathbb{P}}_{k}\;\|\;{\mathbb{P}}_{(0,n]}).$

This implies, that

${\mathbb{E}}\left[G_{n}\left((X_{i})_{i=1}^{n}\mid\alpha^{0},p^{\infty}\right)\right]=\sum_{k=1}^{K}\frac{\alpha^{0}_{k}-\alpha^{0}_{k\shortminus
1}}{n}D_{KL}({\mathbb{P}}_{k}\;\|\;{\mathbb{P}}_{(0,n]}).$

Now assume that $\kappa({(u,v]})$ contains more than a single element. Set
$\pi_{k}({(u,v]})\coloneqq\frac{|{(\alpha^{0}_{k\shortminus
1},\alpha^{0}_{k}]}\cap{(u,v]}|}{v-u}$. Then

$\displaystyle{\mathbb{E}}\left[\sum_{i=u+1}^{v}\log(\frac{p_{(u,v]}(X_{i})}{p_{(0,n]}(X_{i})})\right]$
$\displaystyle=(v-u){\mathbb{E}}_{X\sim{\mathbb{P}}_{(u,v]}}\left[\log(\frac{p_{(u,v]}(X)}{p_{(0,n]}(X)})\right]$
$\displaystyle=(v-u){\mathbb{E}}_{X\sim{\mathbb{P}}_{(u,v]}}\left[\log(\frac{\sum_{k\in\kappa({(u,v]})}\pi_{k}({(u,v]})p_{k}(X)}{p_{(0,n]}(X)})\right]$
$\displaystyle<(v-u)\sum_{k\in\kappa({(u,v]})}\pi_{k}({(u,v]}){\mathbb{E}}_{X\sim{\mathbb{P}}_{k}}\left[\log(\frac{p_{k}(X)}{p_{(0,n]}(X)})\right]$
$\displaystyle=(v-u)\sum_{k\in\kappa({(u,v]})}\pi_{k}({(u,v]})D_{KL}({\mathbb{P}}_{k}\;\|\;{\mathbb{P}}_{(0,n]})\,,$

where we used the fact that the continuous extension of $x\mapsto x\log(x)$ on
$[0,\infty)$ is strictly convex in the third line. Consider any segmentation
$\alpha$ that is not an oversegmentation of $\alpha^{0}$. Then there exists at
least one segment ${(\alpha_{k^{\prime}\shortminus 1},\alpha_{k^{\prime}}]}$
such that $\kappa({(\alpha_{k^{\prime}\shortminus 1},\alpha_{k^{\prime}}]})$
contains more than one element and thus

$\displaystyle{\mathbb{E}}\left[G_{n}\left((X_{i})_{i=1}^{n}\mid\alpha,p^{\infty}\right)\right]$
$\displaystyle=\sum_{k=1}^{K}{\mathbb{E}}\left[\sum_{i=\alpha_{k\shortminus
1}+1}^{\alpha_{k}}\log(\frac{p_{{(\alpha_{k\shortminus
1},\alpha_{k}]}}(X_{i})}{p_{(0,n]}(X_{i})})\right]$
$\displaystyle<\sum_{k=1}^{K}(\alpha_{k}-\alpha_{k\shortminus
1})\sum_{l\in\kappa({{(\alpha_{k\shortminus
1},\alpha_{k}]}})}\pi_{l}({{(\alpha_{k\shortminus
1},\alpha_{k}]}})D_{KL}({\mathbb{P}}_{l}\;\|\;{\mathbb{P}}_{(0,n]})$
$\displaystyle=\sum_{k=1}^{K}\frac{\alpha^{0}_{k}-\alpha^{0}_{k\shortminus
1}}{n}D_{KL}({\mathbb{P}}_{k}\;\|\;{\mathbb{P}}_{(0,n]})\,,$

applying the above inequality with $u=\alpha_{k^{\prime}\shortminus 1}$ and
$v=\alpha_{k^{\prime}}$. ∎

See 3

###### Proof.

Let $s\in{(u,v]}\backslash\alpha^{0}$. We show that the expected classifier
log-likelihood ratio

$G(s):={\mathbb{E}}\left[G\left((X_{i})_{i=u+1}^{v}\mid\\{u,s,v\\},p^{\infty}\right)\right]$

is convex at $s$, that is, $G(s+1)-2G(s)+G(s+1)>0$. By definition

$G(s)={\mathbb{E}}\left[\sum_{i=u+1}^{s}\log(\frac{p_{(u,s]}(X_{i})}{p_{(u,v]}(X_{i})})+\sum_{i=s+1}^{v}\log(\frac{p_{(s,v]}(X_{i})}{p_{(u,v]}(X_{i})})\right]$

and

$\displaystyle G(s+1)-G(s)$
$\displaystyle={\mathbb{E}}\left[\sum_{i=u+1}^{s+1}\log(\frac{p_{(u,s+1]}(X_{i})}{p_{(u,s]}(X_{i})})+\sum_{i=s+2}^{v}\log(\frac{p_{(s+1,v]}(X_{i})}{p_{(s,v]}(X_{i})})+\log(\frac{p_{(u,s]}(X_{s+1})}{p_{(s,v]}(X_{s+1})})\right]$
$\displaystyle=(s-u+1)D_{KL}({\mathbb{P}}_{(u,s+1]}\;\|\;{\mathbb{P}}_{(u,s]})+(v-s-1)D_{KL}({\mathbb{P}}_{(s+1,v]}\;\|\;{\mathbb{P}}_{(s,v]})+$
$\displaystyle\hphantom{\,\,\,\,\,\,\,\,}{\mathbb{E}}\left[\log(\frac{p_{(u,v]}(X_{s+1})}{p_{(s,v]}(X_{s+1})})\right]\geqslant{\mathbb{E}}\left[\log(\frac{p_{(u,s]}(X_{s+1})}{p_{(s,v]}(X_{s+1})})\right]\,,$

with a strict inequality if ${\mathbb{P}}_{(u,s+1]}\neq{\mathbb{P}}_{(u,s]}$
or if $s<v-1$ and ${\mathbb{P}}_{(s+1,v]}\neq{\mathbb{P}}_{(s,v]}$. Similarly,

$\displaystyle G(s-1)-G(s)$
$\displaystyle={\mathbb{E}}\left[\sum_{i=u+1}^{s-1}\log(\frac{p_{(u,s-1]}(X_{i})}{p_{(u,s]}(X_{i})})+\sum_{i=s+1}^{v}\log(\frac{p_{(s-1,v]}(X_{i})}{p_{(s,v]}(X_{i})})+\log(\frac{p_{(s,v]}(X_{s-1})}{p_{(u,s]}(X_{s-1})})\right]$
$\displaystyle=(s-u-1)D_{KL}({\mathbb{P}}_{(u,s-1]}\;\|\;{\mathbb{P}}_{(u,s]})+(v-s+1)D_{KL}({\mathbb{P}}_{(s-1,v]}\;\|\;{\mathbb{P}}_{(s,v]})+$
$\displaystyle\hphantom{\,\,\,\,\,\,\,\,}{\mathbb{E}}\left[\log(\frac{p_{(s,v]}(X_{s})}{p_{(u,v]}(X_{s})})\right]\geqslant{\mathbb{E}}\left[\log(\frac{p_{(s,v]}(X_{s})}{p_{(u,s]}(X_{s})})\right]\,,$

with a strict inequality if ${\mathbb{P}}_{(s-1,v]}\neq{\mathbb{P}}_{(s,v]}$
or if $s>u+1$ and ${\mathbb{P}}_{(u,s-1]}\neq{\mathbb{P}}_{(s,v]}$. Note that
$X_{s},X_{s+1}$ are i.i.d. and thus

${\mathbb{E}}\left[\log(\frac{p_{(u,s]}(X_{s+1})}{p_{(s,v]}(X_{s+1})})+\log(\frac{p_{(s,v]}(X_{s})}{p_{(u,s]}(X_{s})})\right]={\mathbb{E}}\left[\log(\frac{p_{(u,s]}(X_{s+1})}{p_{(u,s]}(X_{s})})+\log(\frac{p_{(s,v]}(X_{s})}{p_{(s,v]}(X_{s+1})})\right]=0.$

Combining everything together,

$G(s+1)-2G(s)+G(s-1)=(G(s+1)-G(s))+(G(s-1)-G(s))\geqslant 0,$

with strict inequality if ${\mathbb{P}}_{(u,s+1]}\neq{\mathbb{P}}_{(u,s]}$ or
${\mathbb{P}}_{(s-1,v]}\neq{\mathbb{P}}_{(s,v]}$. Choose $k$ such that
$s\in{(\alpha_{k\shortminus 1}^{0},\alpha_{k}^{0}]}$. Then also
$s+1\in{(\alpha_{k\shortminus 1}^{0},\alpha_{k}^{0}]}$. Then this is
equivalent to ${\mathbb{P}}_{(u,\alpha_{k\shortminus
1}^{0}]}\neq{\mathbb{P}}_{(\alpha^{0}_{k\shortminus 1},\alpha^{0}_{k}]}$ or
${\mathbb{P}}_{(\alpha^{0}_{k\shortminus
1},\alpha^{0}_{k}]}\neq{\mathbb{P}}_{(\alpha_{k}^{0},v]}$. ∎

See 4

###### Proof.

Write

$\displaystyle G(s)$
$\displaystyle\coloneqq\sum_{i=u+1}^{s}\log(\frac{v-u}{s^{(0)}-u}p_{\\{u,s^{(0)},v\\}}^{\infty}(X_{i}))+\sum_{i=s+1}^{v}\log(\frac{v-u}{v-s^{(0)}}(1-p_{\\{u,s^{(0)},v\\}}^{\infty}(X_{i})))$
(7)
$\displaystyle=\sum_{i=u+1}^{s}\log(\frac{p_{(u,s^{(0)}]}(X_{i})}{p_{(u,v]}(X_{i})})+\sum_{i=s+1}^{v}\log(\frac{p_{(s^{(0)},v]}(X_{i})}{p_{(u,v]}(X_{i})})$
(8)

and define

$U_{j}\coloneqq
G(j)-G(j-1)=\log(\frac{p_{(u,s^{(0)}]}(X_{j})}{p_{(s^{(0)},v]}(X_{j})}),$

such that $G(s)=G(u)+\sum_{j=u+1}^{s}U_{j}$. In particular, for any segment
${(\alpha_{k\shortminus 1}^{0},\alpha_{k}]}$ intersecting the segment
${(u,v]}$ and any $s\in{(\alpha_{k\shortminus 1}^{0},\alpha_{k}]}$, we have
that ${\mathbb{E}}[G(s)]={\mathbb{E}}[G(\alpha_{k\shortminus
1}^{0})]+(s-\alpha_{k\shortminus 1}^{0}){\mathbb{E}}[U_{\alpha_{k\shortminus
1}^{0}+1}]$ as the $U_{\alpha_{k\shortminus 1}^{0}+1},\ldots,U_{\alpha_{k}}$
are i.i.d., showing that ${\mathbb{E}}[G(s)]$ is piecewise linear between the
change points in $\alpha^{0}$.

Now assume that there is a single change point $a^{0}$ in ${(u,v]}$. Assume
without loss of generality that $s^{(0)}\geqslant a^{0}$ such that
$p_{(u,s^{(0)}]}=\frac{a^{0}-u}{s^{(0)}-u}p_{(u,a^{0}]}+\frac{s^{(0)}-a^{0}}{s^{(0)}-u}p_{(a^{0},v]}$
and $p_{(s^{(0)},v]}=p_{(a^{0},v]}$. Then, for $j\leqslant a^{0}$,

$\displaystyle{\mathbb{E}}[U_{j}]$
$\displaystyle={\mathbb{E}}_{{\mathbb{P}}_{(u,a^{0}]}}\left[\log(\frac{p_{(u,s^{(0)}]}(X_{j})}{p_{(s^{(0)},v]}(X_{j})})\right]$
$\displaystyle=\frac{s^{(0)}-u}{a^{0}-u}{\mathbb{E}}_{{\mathbb{P}}_{(u,s^{(0)}]}}\left[\log(\frac{p_{(u,s^{(0)}]}(X_{j})}{p_{(s^{(0)},v]}(X_{j})})\right]-\frac{s^{(0)}-a^{0}}{a^{0}-u}{\mathbb{E}}_{{\mathbb{P}}_{(a^{0},v]}}\left[\log(\frac{p_{(u,s^{(0)}]}(X_{j})}{p_{(s^{(0)},v]}(X_{j})})\right]$
$\displaystyle=\frac{s^{(0)}-u}{a^{0}-u}D_{KL}({\mathbb{P}}_{(u,s^{(0)}]}\;\|\;{\mathbb{P}}_{(a^{0},v]})+\frac{s^{(0)}-a^{0}}{a^{0}-u}D_{KL}({\mathbb{P}}_{(s^{(0)},v]}\;\|\;{\mathbb{P}}_{(u,s^{(0)}]})>0.$

Similarly, for $j>a^{0}$

${\mathbb{E}}[U_{s}]={\mathbb{E}}_{{\mathbb{P}}_{(a^{0},v]}}\left[\log(\frac{p_{(u,s]}(X_{j})}{p_{(s,v]}(X_{j})})\right]=-D_{KL}({\mathbb{P}}_{2}\;\|\;{\mathbb{P}}_{(u,s]})<0,$

which shows that
$s\mapsto{\mathbb{E}}[G(s)]={\mathbb{E}}[G(u)]+\sum_{i=u+1}^{s}{\mathbb{E}}[U_{i}]$
has a unique maximum at $s=a^{0}$. ∎

##  Appendix B Figures

Figure 7:  Histograms of cumulative change point estimates for 500 simulation
runs on the Dirichlet simulation setup. The true change point locations are
marked with green crosses. Note the different scales on the $y$-axis.

##  Appendix C Tables

dataset CIM CIC Dirichlet iris glass breast abalone wine dry average cancer
beans change in mean 0.0 33.3 48.0 1.3 16.8 16.7 5.0 0.4 1.3 13.7 changeforest
0.0 1.7 0.2 0.0 2.8 0.0 3.1 0.0 0.0 0.9 changekNN 0.0 33.3 6.6 0.0 7.9 0.0 5.0
1.5 0.0 6.0 ECP 0.0 30.7 5.1 0.0 31.1 0.0 6.2 0.3 0.0 8.1 KCP 0.0 0.7 5.0 0.0
7.9 0.0 4.2 0.2 0.0 2.0 MultiRank 0.0 33.3 3.0 15.3 15.9 0.0 9.1 3.0 10.0

Table 5:  Median relative Hausdorff distance over 500 simulations (in
percent). Optimal scores are marked in bold.

dataset CIM CIC Dirichlet iris glass breast abalone wine dry cancer beans
ground truth 2 2 10 2 5 1 14 4 6 change in mean 2.00 0.00 0.00 2.99 22.06
16.62 11.07 3.65 7.04 changeforest 2.13 2.36 10.54 2.15 4.93 1.08 12.09 4.21
6.23 changekNN 2.12 0.59 8.01 2.11 4.16 1.49 11.74 4.52 6.21 ECP 2.07 1.26
7.70 2.09 2.58 1.07 9.73 3.71 6.05 KCP 2.03 2.38 8.39 2.05 3.44 1.01 10.62
3.60 6.00 MultiRank 2.00 0.45 9.18 13.31 21.74 8.28 11.77 3.00

Table 6:  Average number of change points estimated over 500 simulations.

trees max. depth mtry CIM CIC Dirichlet iris glass 20 2 1 0.987 (0.033) 0.133
(0.290) 0.972 (0.020) 0.984 (0.039) 0.889 (0.084) $\sqrt{d}$ 0.988 (0.032)
0.185 (0.339) 0.982 (0.017) 0.985 (0.038) 0.903 (0.081) $d$ 0.988 (0.033)
0.189 (0.339) 0.982 (0.016) 0.989 (0.030) 0.913 (0.080) 8 1 0.992 (0.017)
0.286 (0.411) 0.982 (0.014) 0.991 (0.025) 0.892 (0.081) $\sqrt{d}$ 0.991
(0.016) 0.266 (0.403) 0.988 (0.012) 0.993 (0.021) 0.889 (0.098) $d$ 0.987
(0.048) 0.206 (0.370) 0.988 (0.011) 0.992 (0.024) 0.869 (0.131) $\infty$ 1
0.871 (0.268) 0.000 (0.000) 0.979 (0.016) 0.993 (0.022) 0.855 (0.148)
$\sqrt{d}$ 0.849 (0.296) 0.000 (0.000) 0.988 (0.012) 0.993 (0.022) 0.872
(0.121) $d$ 0.777 (0.361) 0.000 (0.000) 0.988 (0.012) 0.992 (0.023) 0.848
(0.157) 100 2 1 0.983 (0.040) 0.498 (0.438) 0.983 (0.019) 0.980 (0.042) 0.899
(0.079) $\sqrt{d}$ 0.982 (0.042) 0.521 (0.439) 0.981 (0.021) 0.982 (0.040)
0.912 (0.071) $d$ 0.983 (0.039) 0.499 (0.434) 0.976 (0.021) 0.982 (0.040)
0.919 (0.068) 8 1 0.987 (0.034) 0.929 (0.130) 0.985 (0.017) 0.981 (0.041)
0.923 (0.069) $\sqrt{d}$ 0.987 (0.035) 0.925 (0.122) 0.986 (0.019) 0.983
(0.040) 0.923 (0.072) $d$ 0.990 (0.031) 0.891 (0.204) 0.984 (0.019) 0.986
(0.034) 0.917 (0.079) $\infty$ 1 0.984 (0.038) 0.880 (0.217) 0.986 (0.018)
0.983 (0.039) 0.917 (0.075) $\sqrt{d}$ 0.985 (0.035) 0.882 (0.222) 0.987
(0.019) 0.984 (0.039) 0.920 (0.073) $d$ 0.988 (0.031) 0.815 (0.312) 0.985
(0.018) 0.986 (0.035) 0.918 (0.076) 500 2 1 0.977 (0.048) 0.529 (0.435) 0.979
(0.024) 0.977 (0.047) 0.910 (0.074) $\sqrt{d}$ 0.978 (0.045) 0.565 (0.432)
0.977 (0.025) 0.981 (0.042) 0.916 (0.072) $d$ 0.980 (0.045) 0.519 (0.433)
0.974 (0.023) 0.982 (0.040) 0.918 (0.070) 8 1 0.982 (0.042) 0.948 (0.051)
0.986 (0.019) 0.980 (0.046) 0.925 (0.071) $\sqrt{d}$ 0.983 (0.040) 0.942
(0.068) 0.985 (0.019) 0.981 (0.044) 0.924 (0.072) $d$ 0.986 (0.037) 0.929
(0.110) 0.983 (0.020) 0.986 (0.035) 0.921 (0.076) $\infty$ 1 0.983 (0.040)
0.941 (0.080) 0.988 (0.017) 0.981 (0.045) 0.921 (0.074) $\sqrt{d}$ 0.981
(0.042) 0.944 (0.068) 0.986 (0.018) 0.981 (0.043) 0.925 (0.072) $d$ 0.984
(0.039) 0.935 (0.110) 0.984 (0.019) 0.986 (0.037) 0.918 (0.078)

trees max. depth mtry breast cancer abalone wine dry beans average 20 2 1
0.992 (0.050) 0.868 (0.083) 0.984 (0.022) 0.997 (0.008) 0.867 (0.107)
$\sqrt{d}$ 0.991 (0.044) 0.907 (0.062) 0.991 (0.014) 0.997 (0.009) 0.881
(0.120) $d$ 0.988 (0.053) 0.923 (0.054) 0.991 (0.021) 0.997 (0.010) 0.884
(0.120) 8 1 0.996 (0.026) 0.843 (0.106) 0.996 (0.009) 0.998 (0.004) 0.886
(0.145) $\sqrt{d}$ 0.994 (0.039) 0.852 (0.108) 0.996 (0.009) 0.998 (0.004)
0.885 (0.144) $d$ 0.993 (0.042) 0.809 (0.131) 0.997 (0.004) 0.998 (0.004)
0.871 (0.140) $\infty$ 1 0.998 (0.006) 0.081 (0.139) 0.991 (0.010) 0.998
(0.004) 0.752 (0.113) $\sqrt{d}$ 0.998 (0.010) 0.089 (0.143) 0.992 (0.009)
0.998 (0.004) 0.753 (0.117) $d$ 0.998 (0.007) 0.070 (0.118) 0.992 (0.009)
0.998 (0.004) 0.740 (0.137) 100 2 1 0.975 (0.089) 0.914 (0.052) 0.976 (0.037)
0.993 (0.017) 0.911 (0.154) $\sqrt{d}$ 0.975 (0.084) 0.931 (0.040) 0.984
(0.037) 0.990 (0.023) 0.917 (0.154) $d$ 0.975 (0.084) 0.937 (0.038) 0.984
(0.040) 0.987 (0.026) 0.916 (0.152) 8 1 0.979 (0.075) 0.931 (0.047) 0.991
(0.031) 0.995 (0.015) 0.967 (0.061) $\sqrt{d}$ 0.983 (0.074) 0.932 (0.049)
0.993 (0.023) 0.997 (0.010) 0.968 (0.059) $d$ 0.984 (0.069) 0.927 (0.055)
0.995 (0.019) 0.996 (0.013) 0.963 (0.081) $\infty$ 1 0.996 (0.017) 0.879
(0.096) 0.997 (0.005) 0.998 (0.006) 0.958 (0.085) $\sqrt{d}$ 0.994 (0.026)
0.902 (0.075) 0.997 (0.005) 0.996 (0.014) 0.961 (0.085) $d$ 0.996 (0.019)
0.892 (0.077) 0.997 (0.005) 0.992 (0.022) 0.952 (0.112) 500 2 1 0.976 (0.085)
0.916 (0.046) 0.971 (0.056) 0.984 (0.033) 0.913 (0.154) $\sqrt{d}$ 0.972
(0.091) 0.932 (0.040) 0.972 (0.057) 0.984 (0.032) 0.920 (0.153) $d$ 0.973
(0.088) 0.938 (0.039) 0.972 (0.061) 0.982 (0.034) 0.915 (0.153) 8 1 0.969
(0.098) 0.934 (0.045) 0.984 (0.043) 0.992 (0.020) 0.967 (0.053) $\sqrt{d}$
0.978 (0.082) 0.938 (0.044) 0.987 (0.037) 0.993 (0.019) 0.968 (0.052) $d$
0.982 (0.075) 0.933 (0.050) 0.994 (0.021) 0.995 (0.015) 0.968 (0.057) $\infty$
1 0.991 (0.030) 0.914 (0.062) 0.997 (0.006) 0.996 (0.012) 0.968 (0.048)
$\sqrt{d}$ 0.994 (0.026) 0.927 (0.052) 0.997 (0.007) 0.991 (0.023) 0.970
(0.045) $d$ 0.994 (0.032) 0.919 (0.061) 0.996 (0.009) 0.994 (0.019) 0.968
(0.054)

Table 7:  Average adjusted Rand indices of changeforest for different tuning
parameters over 500 simulations, with standard deviations in parentheses.
Optimal scores and those within two standard deviations are marked in bold.
When necessary, mtry, the number of features evaluated at each split, was
rounded down to the next integer.

trees max. depth mtry CIM CIC Dirichlet iris glass wine breast cancer abalone
dry beans 20 2 1 0.03 0.01 0.10 0.01 0.02 0.02 0.23 0.24 0.64 $\sqrt{d}$ 0.02
0.01 0.11 0.01 0.02 0.02 0.24 0.25 0.68 $d$ 0.03 0.01 0.14 0.01 0.02 0.02 0.27
0.27 0.86 8 1 0.03 0.02 0.16 0.01 0.02 0.03 0.32 0.34 0.84 $\sqrt{d}$ 0.03
0.02 0.16 0.01 0.03 0.03 0.35 0.36 0.97 $d$ 0.04 0.03 0.21 0.01 0.03 0.03 0.42
0.45 1.5 $\infty$ 1 0.04 0.02 0.18 0.01 0.03 0.03 0.23 0.57 1.2 $\sqrt{d}$
0.04 0.02 0.17 0.01 0.03 0.03 0.27 0.57 1.4 $d$ 0.05 0.02 0.22 0.01 0.03 0.03
0.31 0.72 2.1 100 2 1 0.04 0.03 0.22 0.02 0.03 0.04 0.43 0.45 1.3 $\sqrt{d}$
0.04 0.03 0.24 0.02 0.04 0.04 0.46 0.48 1.4 $d$ 0.05 0.04 0.33 0.02 0.04 0.04
0.56 0.56 2.1 8 1 0.07 0.08 0.43 0.02 0.06 0.07 0.77 0.79 2.1 $\sqrt{d}$ 0.08
0.09 0.43 0.02 0.06 0.07 0.85 0.87 2.6 $d$ 0.10 0.10 0.57 0.02 0.07 0.07 1.1
1.2 4.7 $\infty$ 1 0.09 0.10 0.50 0.03 0.06 0.09 1.3 1.6 3.6 $\sqrt{d}$ 0.10
0.11 0.46 0.02 0.06 0.09 1.4 1.6 4.2 $d$ 0.11 0.12 0.60 0.03 0.07 0.09 1.7 2.2
7.4 500 2 1 0.11 0.08 0.67 0.04 0.10 0.10 1.2 1.4 4.3 $\sqrt{d}$ 0.12 0.09
0.75 0.04 0.10 0.11 1.3 1.5 5.0 $d$ 0.14 0.11 1.1 0.05 0.12 0.12 1.7 1.9 8.2 8
1 0.25 0.27 1.6 0.08 0.20 0.22 2.7 3.0 8.4 $\sqrt{d}$ 0.27 0.30 1.6 0.08 0.20
0.22 3.1 3.4 11 $d$ 0.33 0.37 2.2 0.08 0.22 0.25 4.2 5.0 21 $\infty$ 1 0.32
0.36 1.9 0.08 0.22 0.30 4.9 6.8 15 $\sqrt{d}$ 0.33 0.38 1.7 0.08 0.21 0.29 5.3
7.0 18 $d$ 0.38 0.44 2.3 0.08 0.22 0.31 6.7 9.7 33

Table 8:  Average computational times in seconds on 8 Intel Xeon 2.3 GHz
cores. When necessary, mtry, the number of features evaluated at each split,
was rounded down to the next integer.

kernel bandwidth CIM CIC Dirichlet iris glass cosine 0.801 (0.166) 0.360
(0.250) 0.858 (0.082) 0.028 (0.000) 0.031 (0.000) linear 0.988 (0.033) 0.343
(0.280) 0.869 (0.078) 0.849 (0.163) 0.367 (0.147) Gaussian 0.025 0.993 (0.025)
0.376 (0.357) 0.799 (0.112) 0.943 (0.105) 0.828 (0.157) Gaussian 0.05 0.996
(0.018) 0.731 (0.333) 0.866 (0.080) 0.977 (0.063) 0.830 (0.171) Gaussian 0.1
0.998 (0.012) 0.955 (0.057) 0.870 (0.078) 0.994 (0.023) 0.744 (0.234) Gaussian
0.2 0.999 (0.005) 0.978 (0.032) 0.871 (0.079) 0.996 (0.013) 0.588 (0.258)
Gaussian 0.4 1.000 (0.002) 0.983 (0.022) 0.870 (0.079) 0.994 (0.023) 0.465
(0.278) Gaussian 0.8 0.999 (0.002) 0.973 (0.090) 0.868 (0.079) 0.831 (0.218)
0.153 (0.246) oracle 1.000 (0.002) 0.983 (0.022) 0.871 (0.079) 0.996 (0.013)
0.830 (0.171)

kernel bandwidth breast cancer abalone wine dry beans average cosine 0.036
(0.005) 0.390 (0.037) 0.157 (0.024) 0.794 (0.078) 0.384 (0.108) linear 0.968
(0.054) 0.840 (0.062) 0.630 (0.158) 0.965 (0.021) 0.758 (0.136) Gaussian 0.025
1.000 (0.002) 0.892 (0.053) 0.986 (0.016) 1.000 (0.003) 0.868 (0.141) Gaussian
0.05 0.999 (0.003) 0.905 (0.045) 0.987 (0.014) 1.000 (0.003) 0.921 (0.130)
Gaussian 0.1 0.996 (0.033) 0.916 (0.043) 0.988 (0.011) 1.000 (0.003) 0.940
(0.087) Gaussian 0.2 0.989 (0.055) 0.918 (0.045) 0.988 (0.011) 1.000 (0.000)
0.925 (0.094) Gaussian 0.4 0.987 (0.055) 0.878 (0.076) 0.970 (0.027) 1.000
(0.000) 0.905 (0.102) Gaussian 0.8 0.988 (0.046) 0.819 (0.101) 0.916 (0.026)
0.967 (0.039) 0.835 (0.123) oracle 1.000 (0.002) 0.918 (0.045) 0.988 (0.011)
1.000 (0.003) 0.954 (0.065)

Table 9:  Average adjusted Rand indices of KCP for different tuning parameters
over 500 simulations, with standard deviations in parentheses. Optimal scores
and those within two standard deviations are marked in bold. The oracle is the
best performing hyperparameter combination. We used a Gaussian kernel with
bandwidth 0.1 in the main simulations.