#  Change-Point Detection and Data Segmentation Chapter: The Multiple Change-
in-Gaussian-Mean Problem

Paul Fearnhead Department of Mathematics and Statistics, Lancaster University,
Lancaster LA1 4YF, UK. Email<EMAIL_ADDRESS>Piotr Fryzlewicz
Department of Statistics, London School of Economics, London WC2A 2AE, UK.
Email<EMAIL_ADDRESS>
###### Abstract

A manuscript version of the chapter “The Multiple Change-in-Gaussian-Mean
Problem” from the book “Change-Point Detection and Data Segmentation” by
Fearnhead and Fryzlewicz, currently in preparation. All R code and data to
accompany this chapter and the book are gradually being made available through
https://github.com/pfryz/cpdds.

###### Contents

  1. 1 Impact of Loss Function and Typical Forms of Assumptions and Results
  2. 2 Relationship to Model Choice in Linear Modelling: Similarities and Differences
  3. 3 Estimators Defined as Optima
    1. 3.1 Estimating Change-point Locations When $N_{0}$ Is Known
    2. 3.2 (*) Estimating Change-point Locations When $N_{0}$ Is Known – Theoretical Details
    3. 3.3 When $N_{0}$ Is Unknown: General Principles of Penalised Estimation
    4. 3.4 When $N_{0}$ Is Unknown: Schwarz Information Criterion
    5. 3.5 When $N_{0}$ Is Unknown: Penalised Approaches Derived from Other Information Criteria
    6. 3.6 Computation for Estimators Defined as Optima
  4. 4 Binary Segmentation and Other Hierarchical Methods
    1. 4.1 Binary Segmentation
    2. 4.2 Wild Binary Segmentation
    3. 4.3 From the Bottom Up: Tail-Greedy Unbalanced Haar
  5. 5 Moving Sum and Isolation Approaches
    1. 5.1 MOSUM
    2. 5.2 Isolate-Detect
  6. 6 $L_{1}$-penalised Methods
  7. 7 Estimating the Noise Variance
  8. 8 Adaptive Estimation of $N_{0}$ Based on Output of Solution Path
    1. 8.1 Dimension Jump, Slope Estimation and Lavielle’s Formalisation
    2. 8.2 Steepest Drop to Low Levels
  9. 9 Software: Speed, Design and Accuracy
    1. 9.1 Speed Considerations
    2. 9.2 Design and Functionality
    3. 9.3 Accuracy and Performance Considerations
  10. 10 Example: Well-log Data
    1. 10.1 Exploratory Data Analysis
    2. 10.2 Penalised Cost Approach
    3. 10.3 Wild Binary Segmentation
    4. 10.4 MOSUM
    5. 10.5 Robust Implementation
  11. 11 Bibliographical Notes
  12. 12 Summary

In Fearnhead and Fryzlewicz (2022), we were concerned with the problem of
testing for, and detecting, a single change-point in the mean of a noisy
univariate data sequence. In this chapter, we deal with the task of detecting
potentially multiple change-points.

The main added challenges of the multiple change-point problem in comparison
with the AMOC problem studied in Fearnhead and Fryzlewicz (2022) are
estimating the number of change-points, and computing the estimate of the
number and locations of change-points in a computationally efficient way. The
increased computational challenge should not surprise, as we now have a
multivariate optimisation problem to solve. In addition, it may be more
difficult to estimate the change-point locations with the same accuracy as in
the AMOC problem, the intuitive reason being that the presence of a change-
point may, depending on the method of detection and the landscape of the
signal, adversely impact location estimation for its neighbour. These aspects
make the multiple change-point problem much harder than the AMOC setting of
Fearnhead and Fryzlewicz (2022).

One way of appreciating the additional difficulty is that in the multiple
change-point problem, the space of all possible piecewise-constant signals is
so large that it is often unclear what sorts of tools are the most appropriate
to perform detection in a given noisy dataset. For example, suppose we were to
localise the problem, in the sense of only considering a subsample of the
data, in the hope of detecting at most a single change-point in it, using a
tool available to us from the AMOC setting (such as the CUSUM statistic). If
we took a subsample that was too large, we might end up with a stretch of the
data containing multiple change-points, which could make our AMOC detection
tool suboptimal. Conversely, if we localised with too small a window, our
detection power might suffer and we could end up not having enough data points
to detect any change-points, even if present. This may suggest that it could
be beneficial to look for an optimal window size to localise with, but this
optimal size may vary from one part of the signal to the next, as the change-
points may be spread unevenly across the domain of the signal.

The multiple change-point detection problem can logically be split into two
sub-problems: estimating the number, and estimating the locations of the
change-points. Some methods available in the existing literature have been
designed to solve these two sub-problems at once, some explicitly separate
them, and some can be represented in either of these ways and return the same
output with either representation. For any given approach, having separate
implementations of solutions to these two sub-problems may be helpful as this
then enables flexible pairing, across different approaches, of modules for
model choice (i.e. estimation of the number of change-points) and location
estimation.

This chapter covers a basic setting in which we assume that we have
observations of a univariate piecewise-constant signal plus independent,
identically distributed, zero-mean Gaussian noise, that is

$X_{i}=f_{i}+\varepsilon_{i},\quad i=1,\ldots,n,$ (1)

with the following properties.

  1. 1.

The signal vector $\mathbf{f}=(f_{1},\ldots,f_{n})^{T}$ is deterministic and
piecewise-constant, with $N_{0}$ ($0\leq N_{0}\leq n-1$) change-points,
denoted $\tau^{0}_{1},\ldots,\tau^{0}_{N_{0}}$, for which
$f_{\tau^{0}_{j}}\neq f_{\tau^{0}_{j}+1}$. For notational completeness, we let
$\tau^{0}_{0}=0$, $\tau^{0}_{{N_{0}}+1}=n$.

  2. 2.

The noise vector
$\bm{\varepsilon}=(\varepsilon_{1},\ldots,\varepsilon_{n})^{T}$ is a
realisation of the $n$-variate normal distribution with
$\mathbb{E}(\bm{\varepsilon})=\mathbf{0}_{n}$ and
$\mbox{Var}(\bm{\varepsilon})=\sigma^{2}\mathbf{I}_{n}$, where
$\mathbf{0}_{n}$ is a vector of zeros of length $n$, the matrix
$\mathbf{I}_{n}$ is the identity matrix of dimension $n\times n$, and
$\sigma^{2}$ is a positive constant.

Occasionally, in this chapter we will be making references to works whose
results apply to larger classes of distributions for $\varepsilon_{i}$, which
include $N(0,\sigma^{2})$ as a special case.

Our main task in this chapter is the estimation of the number $N_{0}$ of
change-points, if unknown, and the estimation of their locations
$\tau^{0}_{1},\ldots,\tau^{0}_{N_{0}}$; these will always be assumed unknown.
There are some pedagogical advantages in beginning our exploration of
estimation in multiple change-point problems with the simple stochastic model
(1). Being arguably the simplest model for multiple change-points, it has seen
the most development in the literature, so it offers the widest possible
background for a discussion of the different methodological approaches to
estimation.

As illustrated in Section 2.5 of Fearnhead and Fryzlewicz (2022), most change-
point detection approaches are sensitive to the degree of serial dependence in
the noise in the input data, and need a good estimate of the dependence
structure to perform well. Under the assumption of the absence of serial
correlation (which we make here), its estimation can obviously be skipped,
which enables us to attribute the performance of each change-point detection
method to the quality of the detection algorithm itself rather than to the
accuracy with which it estimates the dependence structure of the noise. The
variance parameter $\sigma^{2}$ in model (1) has to be assumed unknown in
practical applications. How to estimate $\sigma^{2}$, where this is needed,
will be discussed in the respective sections below; see also Section 7.

Model (1) is a useful testing ground for change-point detection methodologies
in the sense that if a method does not perform well in the setting of (1)
despite its simplicity, it is unlikely to perform well in more complex
stochastic settings (e.g. ones involving correlated or non-Gaussian noise, or
outliers). Even the simple model (1) can easily pose serious challenges for
the state of the art in estimation, especially in situations in which the
number ${N_{0}}$ of change-points is large, the change-points are located
close to each other, or the signal-to-noise ratio is small.

Finally, another way in which model (1) is relevant is that some problems can
be reduced to the setting of (1) by a data transformation, e.g. by performing
a variance-stabilising transformation or local pre-averaging of the input
data, or by fitting a statistical model in a more complex setting and working
with the residuals.

##  1 Impact of Loss Function and Typical Forms of Assumptions and Results

Not every estimation method that produces piecewise-constant estimates of
$\mathbf{f}$ in model (1) is useful for the task of consistently estimating
the change-points in $\mathbf{f}$. In particular, some piecewise-constant
estimation methods are consistent for $\mathbf{f}$ in the mean-square sense,
and can therefore be seen as “signal approximators”, rather than change-point
detectors. Examples of such approaches include the fused lasso (Tibshirani et
al., 2005) and the taut string methods (Davies and Kovac, 2001). Both of these
penalise the least-squares fit to the data with a penalty on the sum of the
absolute first differences of the signal (the $L_{1}$ penalty); we cover
$L_{1}$-penalised methods in more detail in Section 6. Another example is the
wavelet thresholding or shrinkage approach using orthogonal Haar wavelets
(Donoho and Johnstone, 1994); see e.g. Vidakovic (1999) for an accessible
introduction to the use of wavelets in statistics.

In this chapter, we will mainly be focusing on methods that aim to estimate
$N_{0}$ and the locations of all change-points as accurately as possible. Many
of the methods described will come with large-sample consistency results
stating that under certain assumptions on the minimum spacing between change-
points and the minimum magnitudes of the corresponding jumps in $\mathbf{f}$,
for $n$ large enough and on an event set whose probability tends to one with
$n$, the number $N_{0}$ of change-points is estimated exactly (i.e.,
$\hat{N}=N_{0}$) and, if $N_{0}\geq 1$, the change-point location estimators
are guaranteed to lie within certain maximum distances from the corresponding
true locations. As a minimum consistency requirement, the largest of these
distances must be of order $o(n)$.

As argued in Section 1.2 of Fearnhead and Fryzlewicz (2022), in the AMOC
setting with i.i.d. Gaussian noise, if the true change-point location is away
from either edge of the data by a distance of $\asymp n$ (throughout the book,
we use $f\asymp g$ to mean $\exists\,\,C,D>0\,\,\,C|g|<|f|<D|g|$), and if the
jump size is bounded away from zero as $n$ increases, then the error in
estimating the change-point location is $O_{p}(1)$. (Recall that
$V_{n}=O_{p}(a_{n})$ means for all $\delta>0$, there exist $M,N$ such that for
all $n>N$ we have $P(|V_{n}/a_{n}|>M)<\delta$.) This best-possible rate
carries over to the multiple change-point setting if the distances between any
pair of consecutive change-points are $\asymp n$, and if the minimum jump size
is bounded away from zero as $n$ increases. A heuristic argument as to why
this is true is that there are estimation procedures in the multiple change-
point setting that attempt to “isolate” change-points as they detect them,
i.e. to break up the $N_{0}$-dimensional estimation problem into a sequence of
$N_{0}$ one-dimensional problems, and for each of these the error rate for
location is $O_{p}(1)$ as argued in Section 1.2 of Fearnhead and Fryzlewicz
(2022). Some of these estimation methods will be reviewed later in this
chapter.

From the discussion in Fearnhead and Fryzlewicz (2022), it is easy to see that
it is impossible to rule out asymptotic underestimation of $N_{0}$ without
some assumptions on the minimum distance between change-points and the minimum
jump size. This phenomenon happens even in the AMOC setting, not to mention in
the more challenging setting of multiple change-points. A heuristic argument
is as follows. Consider the non-centrality parameter, $\nu$, in formula (3) of
Fearnhead and Fryzlewicz (2022) for the distribution of the test statistic for
a change-point under the alternative of a change, and suppose for simplicity
that the location of a possible change-point is known to be $\tau^{0}_{1}$.
The power of the test, and hence the probability of correctly detecting the
change-point at $\tau^{0}_{1}$ if there is one at that location, is
$P(\chi_{1}^{2}(\nu)>k)$, as described in formula (4) of Fearnhead and
Fryzlewicz (2022). We must have $k\to\infty$ for consistency in case there is
no change-point at $\tau^{0}_{1}$. However, if, for example, $\nu$ is bounded
from above as $n\to\infty$, then $P(\chi_{1}^{2}(\nu)>k)\not\to 1$ as
$n\to\infty$ and hence we do not estimate $N_{0}$ consistently if $N_{0}=1$.
When is $\nu$ bounded from above in $n$? From formula (3) of Fearnhead and
Fryzlewicz (2022), this will happen (assuming a fixed jump size $\Delta$) if
either $\tau^{0}_{1}$ or $n-\tau^{0}_{1}$ does not grow with $n$, or in other
words if the segment size to the left or to the right of the change-point does
not increase with $n$. This illustrates the importance of minimum-distance
assumptions. Similar considerations apply to the minimum jump size, both in
the AMOC and in the multiple change-point setting.

Finally, in addition to convergence in probability, we note that some authors
study almost-sure convergence properties of change-point location estimators,
but we leave a more detailed discussion of this interesting point to one of
the later chapters.

##  2 Relationship to Model Choice in Linear Modelling: Similarities and
Differences

In this section, we interpret the multiple change-point problem as a problem
of model choice in linear regression. We start by representing the signal
vector $\mathbf{f}$ in the Euclidean space $\mathbb{R}^{n}$ as a linear
combination of basis vectors $\mathbf{e}^{j}$, each of which contains only one
change-point and is normalised in such a way that it sums to zero and square-
sums to one. More formally, we define

$\mathbf{e}^{j}_{i}=\left\\{\begin{array}[]{cc}\sqrt{\frac{n-j}{nj}}&i=1,\ldots,j,\\\
-\sqrt{\frac{j}{n(n-j)}}&i=j+1,\ldots,n,\\\
0&\text{otherwise},\end{array}\right.$ (2)

for $j=1,\ldots,n-1$. In addition, define $\mathbf{e}^{0}$ with

$\mathbf{e}^{0}_{i}=\sqrt{\frac{1}{n}},\quad i=1,\ldots,n.$

The reader may wonder why we do not define the $\mathbf{e}^{j}$’s as vectors
of the form $(-1,\ldots,-1,1,\ldots,1)$ or $(0,\ldots,0,1,\ldots,1)$: this is
to ensure that $\sum_{i=1}^{n}\mathbf{e}^{j}_{i}=0$ and
$\sum_{i=1}^{n}(\mathbf{e}^{j}_{i})^{2}=1$, a common requirement in the linear
model selection literature which ensures that the sample correlation between
each pair $\mathbf{e}^{j_{1}}$, $\mathbf{e}^{j_{2}}$ can simply be computed as
$\langle\mathbf{e}^{j_{1}},\mathbf{e}^{j_{2}}\rangle$, their inner product. As
an example, a selection of the basis vectors $\mathbf{e}^{j}$ for $n=250$ is
shown in Figure 1.

Figure 1: Basis vectors $\mathbf{e}^{j}$ for $n=250$ and $j=10,20,\ldots,240$.

With this notation in place, we can represent the signal vector $\mathbf{f}$
as

$\mathbf{f}=\beta_{0}\,\mathbf{e}^{0}+\sum_{j\in
S_{0}}\beta_{j}\,\mathbf{e}^{j},$ (3)

where $S_{0}=\left\\{\tau^{0}_{1},\ldots,\tau^{0}_{N_{0}}\right\\}$ is the set
of change-point locations. This is obviously equivalent to

$\mathbf{f}=\beta_{0}\,\mathbf{e}^{0}+\sum_{j\in
S_{0}}\beta_{j}\,\mathbf{e}^{j}+\sum_{j\in S_{0}^{c}}0\cdot\mathbf{e}^{j},$

with $S_{0}^{c}=\\{1,\ldots,n-1\\}\setminus S_{0}$, or

$\mathbf{f}=\beta_{0}\,\mathbf{e}^{0}+\sum_{j=1}^{n-1}\beta_{j}\,\mathbf{e}^{j},$

where $\beta_{j}=0$ for $j\in S_{0}^{c}$. An illustration of the
representation (3) is in Figure 2.

Figure 2: Example signal of length 250 with change-points at locations 60, 180
(solid black). Basis vectors times the linear model coefficients that make up
the signal: $\beta_{0}\mathbf{e}^{0}$ (dashed red),
$\beta_{60}\mathbf{e}^{60}$ (dashed green), $\beta_{180}\mathbf{e}^{180}$
(dashed blue).

Due to the particular normalisation of the $\mathbf{e}^{j}$ vectors, the
standard setting, in which the jump sizes corresponding to each change-point
are fixed in $n$, translates to the condition that
$c_{1}\leq|\beta_{j}|/\sqrt{n}\leq c_{2}$, $j\in S_{0}$, for two positive
constants $c_{1},c_{2}$. With this representation, our multiple change-point
problem is equivalent to the following: given the linear regression problem

$\mathbf{X}=\beta_{0}\,\mathbf{e}^{0}+\sum_{j=1}^{n-1}\beta_{j}\,\mathbf{e}^{j}+\bm{\varepsilon},$
(4)

identify the set $S_{0}=\\{j\,:\,\beta_{j}\neq 0,j>0\\}$, i.e. the set of
change-point locations
$\left\\{\tau^{0}_{1},\ldots,\tau^{0}_{N_{0}}\right\\}$. Therefore, the
multiple change-point detection problem can be interpreted as a sparse
variable selection problem, in a “high-dimensional” setting in which the
dimensionality of the parameter space equals the number of observations $n$
and the cardinality of the set $S_{0}$ of the relevant variables equals
$N_{0}$, presumed to be much smaller than $n$.

The variable selection problem in sparse high-dimensional linear regression
has been studied extensively in the literature, and various methodologies for
the recovery of the set of non-zero variables have been proposed. We note in
particular such techniques as the lasso (Tibshirani, 1996), the adaptive lasso
(Zou, 2006), the elastic net (Zou and Hastie, 2005) and the Dantzig selector
(Candes and Tao, 2007). We say that a variable selection technique has a sign
consistency property if it recovers, with probability tending to one as the
number of observations increases, the set $S_{0}$, with the correct signs of
the constituent $\beta_{j}$’s.

Any sign consistency result must rely on the design matrix satisfying certain
properties. Indeed, if two columns of the design matrix were identical, no
method would be able to identify which of the two covariates in question (if
either) contributes to the response. This intuition agrees, for example, with
Lounici (2008), whose result shows that the lasso has the sign consistency
property if the maximum correlation between two different columns is less than
$1/7$, while a similar result for the Dantzig selector requires it to be less
than $1/3$.

However, the unusual property of the design matrix in the linear model (4) is
that it is, in a certain sense, asymptotically collinear. As an example, take
two neighbouring columns $\mathbf{e}^{j},\mathbf{e}^{j+1}$ of the design
matrix (where $j\geq 1$). The correlation between them is

$\displaystyle\frac{\langle\mathbf{e}^{j},\mathbf{e}^{j+1}\rangle}{\langle\mathbf{e}^{j},\mathbf{e}^{j}\rangle^{1/2}\langle\mathbf{e}^{j+1},\mathbf{e}^{j+1}\rangle^{1/2}}$
$\displaystyle=$
$\displaystyle\langle\mathbf{e}^{j},\mathbf{e}^{j+1}\rangle=\frac{1}{n}\left(j\sqrt{\frac{(n-j)(n-j-1)}{j(j+1)}}\right.$
$\displaystyle-$
$\displaystyle\left.\sqrt{\frac{j(n-j-1)}{(n-j)(j+1)}}+(n-j)\sqrt{\frac{j(j+1)}{(n-j)(n-j-1)}}\right).$

Taking $j=\alpha n$, where $\alpha$ is a constant in $(0,1)$, the dominant
term in the above is

$\frac{1}{n}\left(\alpha
n\sqrt{\frac{(1-\alpha)^{2}}{\alpha^{2}}}-\sqrt{\frac{\alpha(1-\alpha)}{(1-\alpha)\alpha}}+(1-\alpha)n\sqrt{\frac{\alpha^{2}}{(1-\alpha)^{2}}}\right)=1-\frac{1}{n}.$

Therefore, a large portion of the design matrix exhibits asymptotically
perfect collinearity between neighbouring columns.

Essentially, this asymptotic collinearity property means that our design
matrix $\\{\mathbf{e}^{j}\\}_{j=0}^{n-1}$ violates the assumptions required
for the sign consistency of any variable selection method when applied to
problem (4). This means that unless $\sigma^{2}\to 0$ with $n$, i.e. we are in
the (asymptotically) noiseless setting, or if some of the $\beta_{j}$’s for
$j\in S_{0}$ are larger than $O(\sqrt{n})$, i.e. the corresponding jump sizes
increase in $n$, it is impossible for any method to detect the change-point
locations exactly. This makes sense as we know from Section 1.2 of Fearnhead
and Fryzlewicz (2022) that error-free detection of exact change-point
locations is impossible unless the noise variance $\sigma^{2}$ tends to zero
with the sample size $n$ or the relevant jump size increases.

The “variable selection” point of view on multiple change-point detection is
pursued e.g. in Harchaoui and Lévy-Leduc (2010), where this viewpoint is
enabled by moving from a fused-lasso (Tibshirani et al., 2005) formulation of
a multiple change-point detection problem to a standard lasso formulation. We
describe this in more detail in Section 6.

##  3 Estimators Defined as Optima

One traditional route to solving the multiple change-point detection problem
(1) is to identify it with the problem of fitting a piecewise-constant
function to the data $\mathbf{X}=(X_{1},\ldots,X_{n})^{T}$. If the number
$N_{0}$ of change-points is known, in the i.i.d. Gaussian model (1), such a
fit can be carried out via ordinary least squares, which for a Gaussian model
is equivalent to maximum likelihood estimation. This results in an
$N_{0}$-dimensional optimisation problem. Having solved it and obtained the
fit, the change-points in the piecewise-constant estimate serve as estimates
of the change-point locations $\tau^{0}_{1},\ldots,\tau^{0}_{N_{0}}$.

This optimisation-based approach to the multiple change-point problem extends
to the case of unknown number $N_{0}$ of change-points via penalised cost
approaches. In these, the estimator is defined as the minimum of a criterion
which can be decomposed as a measure of fit to the data plus a penalty. The
role of the penalty is to discourage solutions with particularly many change-
points (the exact fit to the data is given by the solution
$\hat{\mathbf{f}}=\mathbf{X}$, which is of no use as a statistical estimator).

The theoretical consistency analysis is then typically done for this
theoretical minimum. Computing the minimum is a different matter entirely.
This section covers the conceptual, theoretical and computational aspects of
unpenalised (for $N_{0}$ known) and penalised (for $N_{0}$ unknown) cost
approaches to multiple change-point detection.

###  3.1 Estimating Change-point Locations When $N_{0}$ Is Known

In this section, we assume that $N_{0}>0$ is known: if not, then we can
estimate it consistently as outlined in the following sections. For any
$\\{k_{1},\ldots,k_{m}\\}\subseteq\\{1,\ldots,n-1\\}$, define

$S_{n}(k_{1},\ldots,k_{m})=\sum_{r=1}^{m+1}\sum_{i=k_{(r-1)}+1}^{k_{(r)}}(X_{i}-\bar{X}_{(k_{(r-1)}+1):k_{(r)}})^{2},$
(5)

where $k_{(0)}=0$, $k_{(m+1)}=n$, $0<k_{(1)}<\ldots k_{(m)}<n$ are the ordered
versions of $k_{1},\ldots,k_{m}$, and
$\bar{X}_{a:b}=\frac{1}{b-a+1}\sum_{i=a}^{b}X_{i}$, the mean of the data
$X_{a},\ldots,X_{b}$. Consider

$(\hat{\tau}_{1},\ldots,\hat{\tau}_{N_{0}})=\arg\min_{l_{1}<\ldots<l_{N_{0}}}S_{n}(l_{1},\ldots,l_{N_{0}}),$
(6)

with the convention $\hat{\tau}_{0}=0$, $\hat{\tau}_{N_{0}+1}=n$. Let
$\hat{q}_{j}=\hat{\tau}_{j}/n$.

The estimators $(\hat{\tau}_{1},\ldots,\hat{\tau}_{N_{0}})$ are the least-
squares estimators of the change-point locations. As mentioned above, they are
also the maximum likelihood estimators if the noise
$\\{\varepsilon_{i}\\}_{i=1}^{n}$ is i.i.d. standard normal. In this section,
we do not discuss the computation of
$(\hat{\tau}_{1},\ldots,\hat{\tau}_{N_{0}})$; this will be covered in Section
3.6. On the face of it, the computation of
$(\hat{\tau}_{1},\ldots,\hat{\tau}_{N_{0}})$ may appear to be a difficult
problem: the parameter space is discrete and constrained, and the minimisation
is $N_{0}$-dimensional. Fortunately, as we will discover in Section 3.6, there
are fast algorithms for computing this multivariate minimum. In this section,
we are concerned with the stochastic behaviour of
$(\hat{\tau}_{1},\ldots,\hat{\tau}_{N_{0}})$ as a theoretically formulated
argument-minimum, regardless of how easy or fast it is to compute it.

Yao and Au (1989) describe the behaviour of
$(\hat{q}_{1},\ldots,\hat{q}_{N_{0}})$ for continuously distributed, i.i.d.
noise $\\{\varepsilon_{i}\\}_{i=1}^{n}$ such that
$\mathbb{E}(\varepsilon_{i})=0$ and $\mathbb{E}(\varepsilon_{i}^{6})<\infty$,
which includes the i.i.d. standard normal case, the main focus of this
chapter. This section summarises this discussion. We first quote and discuss a
result that states that the estimators $(\hat{q}_{1},\ldots,\hat{q}_{N_{0}})$
are consistent for $(q_{1}^{0},\ldots,q_{N_{0}}^{0})$ (where $q_{j}^{0}$ is a
rescaled-time version of $\tau_{j}^{0}$ and is formally defined in Theorem 3.1
below), but with no mention of the rates of convergence just yet. We have the
following result, adapted from Yao and Au (1989).

###### Theorem 3.1

Let the i.i.d. noise $\\{\varepsilon_{i}\\}_{i=1}^{n}$ be such that
$\mathbb{E}(\varepsilon_{i})=0$ and $\mathbb{E}(\varepsilon_{i}^{6})<\infty$.
Assume $\tau_{j}^{0}=[nq_{j}^{0}]$, where $[\cdot]$ denotes the nearest
integer to its argument, and
$0=q_{0}^{0}<q_{1}^{0}<\ldots<q_{N_{0}}^{0}<q_{N_{0}+1}^{0}=1$ are fixed in
the rescaled-time interval $[0,1]$ as $n\to\infty$. Assume also that
$|f_{\tau^{0}_{j}}-f_{\tau^{0}_{j}+1}|>0$ remains constant with $n$ for each
$j=1,\ldots,{N_{0}}$. Then

$|\hat{q}_{j}-q_{j}^{0}|=o_{p}(1),$ (7)

for $1\leq j\leq N_{0}$, as $n\to\infty$.

We recall that $V_{n}=o_{p}(a_{n})$ means that for all $\delta>0$, we have
$\lim_{n\to\infty}P(|V_{n}/a_{n}|>\delta)=0$.

In the statement of Theorem 3.1, the concept of “rescaled time” refers to the
fact that we are scaling the true change-point locations $\tau_{j}^{0}$ by $n$
to create $q_{j}^{0}$, thereby embedding the true domain $\\{1,2,\ldots,n\\}$
in $(0,1]$, and examining the convergence of the change-point location
estimators $\hat{q}_{j}$ in the new domain $(0,1]$. This is because apart from
trivial situations (such as the variance $\sigma^{2}$ of the innovations
converging to zero with $n$), it is impossible for $\hat{\tau}_{j}$ to
converge to $\tau_{j}^{0}$: we saw in Section 1.2 of Fearnhead and Fryzlewicz
(2022) that the error of $\hat{\tau}_{j}$ as a location estimator is at least
$O_{p}(1)$.

Theorem 3.1 is a bare-bones consistency result for the change-point location
estimators $\hat{q}_{j}$. Indeed, (7) can be viewed as a basic sanity check
and if it were not to hold, the estimators $\hat{q}_{j}$ could not be regarded
as useful. The key ideas of the proofs of Theorem 3.1 and the following
Theorem 3.2 are in Section 3.2, intended for the more theoretically-minded
reader.

We now improve on the result of Theorem 3.1 with the following result, also
from Yao and Au (1989).

###### Theorem 3.2

Let the i.i.d. noise $\\{\varepsilon_{i}\\}_{i=1}^{n}$ be such that
$\mathbb{E}(\varepsilon_{i})=0$ and $\mathbb{E}(\varepsilon_{i}^{6})<\infty$.
Assume $\tau_{j}^{0}=[nq_{j}^{0}]$, where $[\cdot]$ denote the nearest integer
to its argument, and
$0=q_{0}^{0}<q_{1}^{0}<\ldots<q_{N_{0}}^{0}<q_{N_{0}+1}^{0}=1$ are fixed in
the rescaled-time interval $[0,1]$ as $n\to\infty$. Assume also that
$|f_{\tau^{0}_{j}}-f_{\tau^{0}_{j}+1}|>0$ remains constant with $n$ for each
$j=1,\ldots,{N_{0}}$. Then

$n|\hat{q}_{j}-q_{j}^{0}|=O_{p}(1),$ (8)

for $1\leq j\leq N_{0}$, as $n\to\infty$.

We remark that the result of Theorem 3.2 is optimal in the sense that the rate
$n|\hat{q}_{j}-q_{j}^{0}|=O_{p}(1)$ could not be improved even if the noise
$\\{\varepsilon_{i}\\}_{i=1}^{n}$ were i.i.d. Gaussian, and even if we were in
the AMOC setting (as we saw in Section 1.2 of Fearnhead and Fryzlewicz
(2022)). Therefore, we can see that Gaussianity is not needed for optimality
in this sense and that the finiteness of the sixth moment suffices.

To better appreciate the result of Theorem 3.2, we consider an example
illustrating how the error in estimating the change-point locations changes as
we vary $n$.

Figure 3 shows the `blocks` signal from Donoho and Johnstone (1994), sampled
at $n=500$ points. The signal has 11 change-points. It is contaminated with
i.i.d. Gaussian noise with mean zero and variance $\sigma^{2}=100$. We up-
sample the signal to lengths $n=500,1000,2000,4000,8000,16000$ and for each
length, we estimate the $N_{0}=11$ change-point locations via the least-
squares fit in formula (6). We repeat over 250 realisations of the Gaussian
noise. As is evident from the right-hand box plots in Figure 3, the empirical
distribution of the quantity
$\max_{j=1,\ldots,11}|\tau_{j}^{0}-\hat{\tau}_{j}|$ appears to settle on a
limiting distribution. This is consistent with the message of Theorem 3.2,
which states that this distance should behave like $n^{-1}$ in rescaled time,
and therefore be of a constant order (with respect to $n$) on the original
time scale, as $n$ gets larger. The exact limiting distribution of the maximum
absolute distance between the estimated change-points and the truth is a
complicated random function of the jump sizes and jump signs, and can be
obtained from Theorem 1 in Yao and Au (1989). We omit this result here. An
important observation, though, is that increasing the number of data points
between successive change-points increases the power of detecting that there
is a change, but does not impact the accuracy of the estimate of the location.
The latter is primarily determined by the pattern of the data close to the
change-point.

Figure 3: Left: The blocks signal for $n=500$ (black), its noisy realisation
(grey) and the 11 best change-point estimates via the least-squares criterion
(red). Right: box plots of the distribution of
$\max_{j=1,\ldots,11}|\tau_{j}^{0}-\hat{\tau}_{j}|$ over 250 simulated sample
paths, for lengths $n=500,1000,2000,4000,8000,16000$.

In the above example, we computed the argument minimum in (6) via the Segment
Neighbourhood approach (Auger and Lawrence, 1989). In R, for a data vector
stored in `x`, this can be carried out via the calls

    library(changepoint)
    cpts(cpt.mean(x, "None", method = "SegNeigh", Q = 12))

In the above, `Q` should be set to the number of segments, $N_{0}+1$.
Computation for change-point estimators defined as optima is covered in more
detail in Section 3.6.

###  3.2 (*) Estimating Change-point Locations When $N_{0}$ Is Known –
Theoretical Details

This section collects some key fact and arguments used in the proofs of
Theorems 3.1 and 3.2.

Key ideas of the proof of Theorem 3.1. The main idea of the proof is simple,
and the proof is, in essence, by contradiction. It is possible to get lower
bounds on $S_{n}(l_{1},\ldots,l_{N_{0}})$ in terms of other values of
$S_{n}(\cdots)$ by allowing for more than $N_{0}$ change-points: as adding
change-points to the set $(l_{1},\ldots,l_{N_{0}})$ can only reduce
$S_{n}(\cdots)$. If (7) did not hold, then the value of
$S_{n}(l_{1},\ldots,l_{N_{0}})$ at the estimated change-points could be lower-
bounded by such an other value of $S_{n}(\cdots)$, which included as one of
its summands the sum of squares

$\sum_{i=r_{s}}^{r_{e}}(X_{i}-\bar{X}_{r_{s}:r_{e}})^{2},$

where $r_{s}=[n(q_{j}^{0}-\nu)]+1$ and $r_{e}=[n(q_{j}^{0}+\nu)]$, for some
$j=1,\ldots,N_{0}$ and a (small) fixed constant $\nu\in(0,1)$. The fact that
the interval $[r_{s},r_{e}]$ included the change-point $q_{j}^{0}$ would cause
this sum of squares to be “large”: of order $\asymp n$. This in turn would
cause this lower bound for $S_{n}(l_{1},\ldots,l_{N_{0}})$ to be substantially
larger than $S_{n}(\tau^{0}_{1},\ldots,\tau^{0}_{N_{0}})$. This demonstrates
the impossibility of $S_{n}(l_{1},\ldots,l_{N_{0}})$ being minimised at points
not satisfying (7), and therefore leads us to deduce that (7) must hold.

There is one more aspect of this result that is worth discussing here. Namely,
in the reasoning outlined above, how can we be sure that the behaviour of the
noise $\\{\varepsilon_{i}\\}_{i=1}^{n}$ will not be such that the lower bound
for $S_{n}(l_{1},\ldots,l_{N_{0}})$ considered above is spuriously smaller
than $S_{n}(\tau^{0}_{1},\ldots,\tau^{0}_{N_{0}})$, or at least approaches it,
which would cause this reasoning to break down? The answer lies in the
following lemma shown in Yao and Au (1989), which is a weaker version of the
sharp Gaussian Lemma 3.2 but applying to a wider class of distributions:
centred with a finite $2m$th moment.

###### Lemma 3.1

Suppose $\epsilon_{1},\ldots,\epsilon_{n}$ are i.i.d. with common mean $0$ and
such that $\mathbb{E}(\epsilon_{i}^{2m})<\infty$ for some positive integer
$m$. Then as $n\to\infty$,

$\max_{0\leq i<j\leq
n}(\epsilon_{i+1}+\ldots+\epsilon_{j})^{2}/(j-i)=O_{p}(n^{2/m}).$

Lemma 3.1 ensures that with a high probability, in the a contrario sketch
proof outlined above, the lower bound on $S_{n}(l_{1},\ldots,l_{N_{0}})$ does
not fluctuate “too much” and therefore still remains significantly above
$S_{n}(\tau^{0}_{1},\ldots,\tau^{0}_{N_{0}})$. Lemma 3.1 is fairly crude and
the main tool in its proof is Markov’s inequality. What may come as a surprise
is that even the crude Lemma 3.1 is sufficient to ensure consistency of the
change-point location estimators (not only for noise with more general
distributions than the Gaussian, but also, as we will see next, with optimal
rates).

Looking ahead, the sharp result in Lemma 3.2, whose proof is fairly
complicated, will later improve the bound from Lemma 3.1 for the Gaussian
distribution. We will see that we will need the sharpness of Lemma 3.2 to
correctly estimate the number $N_{0}$ of change-points via Schwarz Information
Criterion (SIC). The sharpness of Lemma 3.2 needed to estimate $N_{0}$ well
for Gaussian noise, contrasted with the relative looseness of Lemma 3.1
sufficient to estimate the change-point locations optimally, and for more
general noise distributions, may signal that estimating the number of change-
points is a more difficult problem than estimating their locations.

Key ideas of the proof of Theorem 3.2. For the result to hold, it is enough to
show, for any fixed $\delta>0$, that there exists $M>0$ such that
$P(n|\hat{q}_{j}-q_{j}^{0}|>M)<\delta$ for $n$ large enough. There are in
essence two sets of results that are used in the proof of Theorem 3.2:

  1. 1.

The result of Theorem 3.1.

  2. 2.

What Yao and Au (1989) refer to as the strong law of large numbers, which
ensures that there exist $C>0$ and integer $M>0$ such that the event
$E_{j}(n,C,M)$:

$\displaystyle|\bar{X}_{[nq^{0}_{s}]:([nq^{0}_{s}]+i)}|$ $\displaystyle<$
$\displaystyle C\quad\mbox{for}\quad
i=1,\ldots,[nq^{0}_{s+1}]-[nq^{0}_{s}],\quad s=0,\ldots,N_{0},$
$\displaystyle|\bar{X}_{([nq^{0}_{s+1}]-i):[nq^{0}_{s+1}]}|$ $\displaystyle<$
$\displaystyle C\quad\mbox{for}\quad
i=1,\ldots,[nq^{0}_{s+1}]-[nq^{0}_{s}],\quad s=0,\ldots,N_{0},$
$\displaystyle|\bar{X}_{([nq^{0}_{j}]-i):[nq^{0}_{j}]}-f_{\tau^{0}_{j}}|$
$\displaystyle<$ $\displaystyle\delta_{1}\quad\mbox{for}\quad
i=M,\ldots,[nq^{0}_{j}]-[nq^{0}_{j-1}],$
$\displaystyle|\bar{X}_{[nq^{0}_{j}]:([nq^{0}_{j}]+i)}-f_{\tau^{0}_{j+1}}|$
$\displaystyle<$ $\displaystyle\delta_{1}\quad\mbox{for}\quad
i=M,\ldots,[nq^{0}_{j+1}]-[nq^{0}_{j}],$ (9)

where $\delta_{1}$ is a particular small (fixed, signal-dependent) constant,
has a probability larger than $1-\delta/4$ for $n$ large enough. The
conditions (9) collectively ensure that the input data, for $n$ sufficiently
large, is “well-behaved” with an arbitrarily high probability, in a certain
uniform sense. More specifically, the first two conditions in (9) ensure that
the averages of the data to the left and to the right of each change-point are
uniformly bounded by a constant, while the last two conditions ask that the
averages of the data to the left and to the right of the $j$th change-point
(the change-point singled out in the statement of the theorem) are within a
certain specific signal-dependent distance of the corresponding population
means, provided that the spans of these averages are at least $M$; the
constant $M$ can be arbitrary here, as long as it exists and is fixed (i.e.
independent of $n$). Similar conditions ensuring that the input data is well-
behaved are commonplace in consistency proofs in the change-point literature.
This should not surprise: for example, if the data were not to be well-behaved
in the sense of having a cluster of, say, particularly high values, any
change-point procedure might be inclined to spuriously detect extra change-
points delineating the cluster. Of course, clusters of particularly high
values can occur even under Gaussianity, but the probability of such an
extreme event is low – hence the requirement that (9) should hold with a
probability that is “high” but less than one. It would be unrealistic to
request that (9), or a similar set of conditions, should hold with probability
one.

It is insightful to examine, at this point, what exactly is meant by an
increasing sample size $n$ in the results of this section. From the “rescaled
time” discussion in Section 3.1, we can unambiguously deduce that we operate
in a stochastic regime in which $\asymp n$ data points are collected between
each pair of consecutive change-points. However, one theoretically important
question to ask is whether, as $n$ increases, we “throw away” all existing
data points between each pair of consecutive change-points and replace them by
new, independent data with length increased by one, or whether we simply add
one more new data point to the existing sequence of observations between each
pair of consecutive change-points. Although this is not an issue that is
discussed in full in Yao and Au (1989), the steps taken in the proof suggests
a setting in which as $n$ increases, a new data point gets added to the
existing data between each pair of consecutive change-points.

The main logic of the proof is that from the constraints that describe the set
$E_{j}(n,C,M)$ (formula (9)), it is possible to show the following. Let

$B(n,\delta_{2})=\\{(l_{1}/n,\cdots,l_{N_{0}}/n)\quad:\quad
0<l_{1}<\ldots<l_{N_{0}}<n,\,\,|l_{j}/n-q_{j}^{0}|<\delta_{2}\\},$

for a suitably chosen, sufficiently small, $\delta_{2}$, and

$B_{j}(n,\delta_{2},M)=\\{(l_{1}/n,\cdots,l_{N_{0}}/n)\in
B(n,\delta_{2})\quad:\quad l_{j}-nq_{j}^{0}<-M\\}.$

By Theorem 3.1, $(\hat{q}_{1},\ldots,\hat{q}_{N_{0}})\in B(n,\delta_{2})$ with
probability larger than $1-\delta/4$ for large $n$. It is claimed that
$(\hat{q}_{1},\ldots,\hat{q}_{N_{0}})\in B_{j}(n,\delta_{2},M)$ with
probability less than $\delta/4$ for large $n$. Assuming that this claim
holds, we have, for large $n$,

$\displaystyle P(n(\hat{q}_{j}-q^{0}_{j})<-M)$ $\displaystyle\leq$
$\displaystyle P((\hat{q}_{1},\ldots,\hat{q}_{N_{0}})\not\in B(n,\delta_{2}))$
$\displaystyle+$ $\displaystyle P((\hat{q}_{1},\ldots,\hat{q}_{N_{0}})\in
B_{j}(n,\delta_{2},M))$ $\displaystyle<$
$\displaystyle\delta/4+\delta/4=\delta/2.$

Similarly, it can be shown that $P(n(\hat{q}_{j}-q^{0}_{j})>M)<\delta/2$ for
large $n$, and so $P(n|\hat{q}_{j}-q^{0}_{j}|>M)<\delta$ for large $n$. To
prove the claim, for every $(l_{1}/n,\cdots,l_{N_{0}}/n)\in
B_{j}(n,\delta_{2},M)$, let $(l_{1}^{\prime}/n,\cdots,l_{N_{0}}^{\prime}/n)\in
B(n,\delta_{2})$ be such that $l_{s}^{\prime}=l_{s}$ for $s\neq j$ and
$l_{j}^{\prime}=[nq_{j}^{0}]$.

After some fairly elementary algebra (we omit the details here), the
constraints in the set $E_{j}(n,C,M)$ permit us to establish that

$\displaystyle\min_{(l_{1}/n,\cdots,l_{N_{0}}/n)\in
B_{j}(n,\delta_{2},M)}\\{S_{n}(l_{1},\ldots,l_{N_{0}})-S_{n}(l_{1}^{\prime},\ldots,l_{N_{0}}^{\prime})\\}$
$\displaystyle\geq$
$\displaystyle\frac{1}{16}|f_{\tau^{0}_{j}}-f_{\tau^{0}_{j}+1}|^{2}\min_{(l_{1}/n,\cdots,l_{N_{0}}/n)\in
B_{j}(n,\delta_{2},M)}([nq^{0}_{j}]-l_{j})$ $\displaystyle\geq$
$\displaystyle\frac{1}{16}|f_{\tau^{0}_{j}}-f_{\tau^{0}_{j}+1}|^{2}M>0,$

which implies that $(\hat{q}_{1},\ldots,\hat{q}_{N_{0}})\not\in
B_{j}(n,\delta_{2},M)$. Therefore, for large $n$,

$P((\hat{q}_{1},\ldots,\hat{q}_{N_{0}})\in B_{j}(n,\delta_{2},M))\leq
1-P(E_{j}(n,C,M))<\delta/4.$

This proves the claim and thus the Theorem.

###  3.3 When $N_{0}$ Is Unknown: General Principles of Penalised Estimation

In many cases of practical interest, the number $N_{0}$ of change-points will
not be known to the analyst and will need to be estimated from the data along
with the change-point locations. In this section, we discuss general
principles behind constructing criteria, penalised for model complexity, whose
minimisation will yield estimators of both $N_{0}$ and
$(\tau^{0}_{1},\ldots,\tau^{0}_{N_{0}})$.

We saw in Section 3.1 that for a known $N_{0}$, a valid way to obtain
estimators of the change-point locations
$(\tau^{0}_{1},\ldots,\tau^{0}_{N_{0}})$ was to minimise the least-squares
criterion $S_{n}(l_{1},\ldots,l_{N_{0}})$. Some of the criteria for unknown
$N_{0}$ outlined in the next sections also use $S_{n}(\cdot)$ as a building
block, while some others use the likelihood of the unknown true model
parameters
$\theta_{0}:=(N_{0},(\tau^{0}_{1},\ldots,\tau^{0}_{N_{0}}),\sigma^{2},(f_{\tau^{0}_{1}},\ldots,f_{\tau^{0}_{N_{0}}},f_{n}))$.
Let a candidate model parameter be denoted by
$\theta=(N^{\prime},(l_{1},\ldots,l_{N^{\prime}}),\tilde{\sigma}^{2},(\tilde{f}_{l_{1}},\ldots,\tilde{f}_{l_{N^{\prime}}},\tilde{f}_{n}))$.
The likelihood of $\theta$ given the data in model (1) is

$L(\theta|X_{1},\ldots,X_{n})=(2\pi\tilde{\sigma}^{2})^{-n/2}\exp\left\\{-\frac{\sum_{j=1}^{N^{\prime}+1}\sum_{i=l_{j-1}+1}^{l_{j}}(X_{i}-\tilde{f}_{l_{j}})^{2}}{2\tilde{\sigma}^{2}}\right\\}.$

Taking the logarithm, multiplying by $-1$ and removing any additive terms that
do not depend on $\theta$, we obtain the negative log-likelihood $-l(\cdot)$,
defined by

$-l(\theta|X_{1},\ldots,X_{n})=\frac{n}{2}\log\,\tilde{\sigma}^{2}+\frac{\sum_{j=1}^{N^{\prime}+1}\sum_{i=l_{j-1}+1}^{l_{j}}(X_{i}-\tilde{f}_{l_{j}})^{2}}{2\tilde{\sigma}^{2}}.$
(10)

Both $S_{n}(\cdot)$ and $-l(\cdot)$ are measures of fit to the data of a
particular candidate model. But while their minimisation can be used to
estimate the change-point locations for a fixed $N^{\prime}$, an estimate of
$N_{0}$ cannot be obtained by minimising these criteria alone. This is
because, heuristically speaking, the fit of a model to the data improves as
its putative number of change-point increases. For example, for the least-
squares measure $S_{n}(\cdot)$, we generally have (except some trivial cases
in which we have equality)

$\min_{l_{1}<\ldots<l_{N^{\prime}}}S_{n}(l_{1},\cdots,l_{N^{\prime}})>\min_{l_{1}<\ldots<l_{N^{\prime\prime}}}S_{n}(l_{1},\cdots,l_{N^{\prime\prime}})$

for $N^{\prime}<N^{\prime\prime}$, as on the right-hand side we are taking the
minimum over a larger model class. Therefore, the absolute minimum of
$S_{n}(\cdot)$ is achieved for $N^{\prime}=n-1$ and $l_{j}=j$, or in other
words when “every point is a change-point”, in which case
$S_{n}(l_{1},\ldots,l_{N^{\prime}})=0$. Clearly, this is not a sensible
statistical model, as it in no way generates additional insight. A very
similar argument can be made for the negative log-likelihood, $-l(\cdot)$.

For this reason, for the analyst wishing to estimate both $N_{0}$ and
$(\tau^{0}_{1},\ldots,\tau^{0}_{N_{0}})$ via minimising a model cost function,
a natural and valid approach is to consider model cost functions of the form

$\mbox{fit to data}+\mbox{penalty},$ (11)

where the “fit to data” part measures the quality of fit of a model candidate
(the reader is invited to think of $S_{n}(\cdot)$ or $-l(\cdot)$ here) and the
“penalty” term penalises more complex models, e.g. ones containing more
change-points, or containing change-points in particular configurations which
are deemed unlikely (e.g. being located close to one another). In this
approach, the preferred model is one that minimises a given criterion of the
form (11).

Many readers will be familiar with the penalised approach to model selection
in other statistical contexts (e.g. model selection in linear regression). In
change-point problems, the most frequently used penalties are proportional to
the number of estimated change-points, that is they are of the form
$N^{\prime}\rho(n)$, where $\rho(\cdot)$ is a particular function of the
sample size (and possibly $\sigma^{2}$, if known). However, more complex
penalties exist and we review a selection of the most popular ones later in
this chapter.

All methods of this form reviewed in this book will use either $S_{n}(\cdot)$
or $-l(\cdot)$ as the “fit to data” term. We now briefly discuss their
particular features as well as the similarities and differences between them.
If a particular component of the parameter vector $\theta$ does not appear in
the penalty, it can be replaced in $-l(\cdot)$ by its unconstrained maximum
likelihood estimator conditioned on those elements of the parameter vector
$\theta$ that do appear in the penalty. For example, in Schwarz Information
Criterion reviewed in detail in Section 3.4 below, the penalty is of the form
$N^{\prime}\log(n)$, and we therefore make the following replacements in the
criterion $-l(\cdot)+N^{\prime}\log(n)$:

$\displaystyle(l_{1},\ldots,l_{N^{\prime}})$ $\displaystyle\rightsquigarrow$
$\displaystyle\arg\min_{l_{1}<\ldots<l_{N^{\prime}}}S_{n}(l_{1},\ldots,l_{N^{\prime}}),$
$\displaystyle\tilde{\sigma}^{2}$ $\displaystyle\rightsquigarrow$
$\displaystyle\hat{\sigma}_{N^{\prime}}^{2},$ $\displaystyle\tilde{f}_{l_{j}}$
$\displaystyle\rightsquigarrow$ $\displaystyle\bar{X}_{(l_{j-1}+1):l_{j}},$

where the quantities on the right-hand sides of these expressions are the MLEs
of the quantities on the left under the assumption that the number of change-
points equals $N^{\prime}$; in particular,

$\hat{\sigma}_{N^{\prime}}^{2}=n^{-1}\min_{l_{1}<\ldots<l_{N^{\prime}}}S_{n}(l_{1},\ldots,l_{N^{\prime}}).$

This replacement, and the removal of an additive constant which makes no
difference to the minimisation, reduce the $-l(\cdot)$ function in this
particular case to $\frac{n}{2}\log(\hat{\sigma}_{N^{\prime}}^{2})$. This
causes Schwarz Information Criterion to have the form
$\frac{n}{2}\log(\hat{\sigma}_{N^{\prime}}^{2})+N^{\prime}\log(n)$. There is
more material on under what conditions Schwarz Information Criterion leads to
consistent model selection in Section 3.4.

On the face of it, Schwarz Information Criterion is a criterion for the
estimation of $N_{0}$ only, rather than the change-point locations, as these
do not explicitly feature in the expression
$\frac{n}{2}\log(\hat{\sigma}_{N^{\prime}}^{2})+N^{\prime}\log(n)$. However,
to obtain this expression, we implicitly had to estimate the change-point
locations for each $N^{\prime}$ separately by
$\arg\min_{l_{1}<\ldots<l_{N^{\prime}}}S_{n}(l_{1},\ldots,l_{N^{\prime}})$.
This specific observation leads us to make the following more general
statement.

###### Note 3.1

In penalised problems of the form $-l(\cdot)+\mbox{penalty}$ (or
$S_{n}(\cdot)+\mbox{penalty}$, respectively), if the penalty depends on the
number of change-points, but not on their locations, joint estimation of the
number and locations of change-points can be achieved by following the
workflow below:

  1. 1.

estimate the change-point locations for each postulated number $N^{\prime}$ of
change-points separately via conditional MLE, by minimising $-l(\cdot)$,
without a penalty, under the assumption that the number of change-points
equals $N^{\prime}$ (which is equivalent to the minimisation of $S_{n}(\cdot)$
in the signal + i.i.d. Gaussian noise setting, under the assumption that the
number of change-points equals $N^{\prime}$);

  2. 2.

estimate the number $N_{0}$ of change-points by minimising
$-l(\cdot)+\mbox{penalty}$ (or $S_{n}(\cdot)+penalty$, respectively); in these
expressions, substitute the unknown change-point locations by their MLEs
obtained in the previous step.

This provides a convenient separation of the joint estimation of the number
and locations of change-points, for penalties which only depend on the former
but not the latter, into unpenalised estimation of the locations for each
number first, and then penalised estimation of the number.

The above discussion concerns the case of
$\sigma^{2}=\mbox{Var}(\varepsilon_{i})$ being unknown. It is also useful to
study the change-point problem (1) in the case of $\sigma^{2}$ being known,
the main reason being that even if it is unknown, it can typically be
estimated fairly accurately, at least when the number $N_{0}$ of change-points
is not very large in relation to the sample size $n$. This can be achieved,
for example, by using the Median Absolute Deviation (Hampel, 1974) or Inter-
Quartile Range (Rousseeuw and Croux, 1993) estimators suitable for Gaussian
data on the input sequence $\\{2^{-1/2}(X_{i+1}-X_{i})\\}_{i=1}^{n-1}$. Here,
the rationale is that most elements of this sequence (more precisely, those
for which both $X_{i}$ and $X_{i+1}$ are taken from the same segment of
constancy) will have mean zero and variance $\sigma^{2}$, and therefore robust
estimators of $\sigma^{2}$ should be accurate. There is more on the estimation
of $\sigma^{2}$ in Section 7.

Let us now investigate how the minimisation of $-l(\cdot)+\mbox{penalty}$ will
simplify in the case of $\sigma^{2}$ being known. From formula (10),
disregarding the (now assumed known) additive term
$\frac{n}{2}\log\,\tilde{\sigma}^{2}=\frac{n}{2}\log\,\sigma^{2}$, we can see
that the minimisation of $-l(\cdot)+\mbox{penalty}$ is equivalent to the
minimisation of $S_{n}(\cdot)+2\sigma^{2}\mbox{penalty}$. This establishes the
equivalence, up to the multiplicative factor of $2\sigma^{2}$, of the
likelihood and least-squares approaches in the case of $\sigma^{2}$ known.

We conclude this section with definitions of a solution path in change-point
problems. Note 3.1 implies that it is of interest to store sets of change-
point estimates
$(\hat{\tau}_{1}^{(N^{\prime})},\ldots,\hat{\tau}_{N^{\prime}}^{(N^{\prime})})$
for $N^{\prime}=1,\ldots,N_{U}$, where $N_{U}$ is a certain upper bound. We
refer to the list
$(\hat{\tau}_{1}^{(N^{\prime})},\ldots,\hat{\tau}_{N^{\prime}}^{(N^{\prime})})_{N^{\prime}=1}^{N_{U}}$
as the solution path to a given change-point problem associated with a
particular estimation method. We refer to a solution path as complete if
$N_{U}$ is the largest possible, i.e. $N_{U}=n-1$.

We refer to an estimation method as nested if for each
$N^{\prime}<N^{\prime\prime}$ we have
$(\hat{\tau}_{1}^{(N^{\prime})},\ldots,\hat{\tau}_{N^{\prime}}^{(N^{\prime})})\subset(\hat{\tau}_{1}^{(N^{\prime\prime})},\ldots,\hat{\tau}_{N^{\prime\prime}}^{(N^{\prime\prime})})$.
Estimation methods based on penalised optimisation are typically not nested;
heuristically, the reason for this is that estimation (via optimisation) for
$N^{\prime}$ postulated change-points is conducted separately from estimation
for $N^{\prime\prime}$ postulated change-points, so there is no algorithmic
link between the two solution sets and hence no reason for one to be a subset
of the other. However, some other methods introduced later in the chapter will
have the nested property. For a nested estimation method, denote for
$N^{\prime}=1,\ldots,N_{U}$,

$\hat{\tau}^{(N^{\prime})}_{j(N^{\prime})}=(\hat{\tau}_{1}^{(N^{\prime})},\ldots,\hat{\tau}_{N^{\prime}}^{(N^{\prime})})\setminus(\hat{\tau}_{1}^{(N^{\prime}-1)},\ldots,\hat{\tau}_{N^{\prime}-1}^{(N^{\prime}-1)}),$

the additional estimated change-point when we go from $N^{\prime}-1$ to
$N^{\prime}$ estimated change-points. The ordered set
$\\{\hat{\tau}^{(N^{\prime})}_{j(N^{\prime})}\\}_{N^{\prime}=0}^{N_{U}}$ is
then referred to as the nested solution path of the nested estimation method
under consideration, and it is said to be complete if $N_{U}=n-1$.

###  3.4 When $N_{0}$ Is Unknown: Schwarz Information Criterion

Yao (1988) proposes to estimate the number ${N_{0}}$ of change-points via
Schwarz (a.k.a. the Bayesian) Information Criterion (Schwarz, 1978). As
motivated in Section 3.3, the criterion is

$\mbox{SIC}(N^{\prime})=\frac{n}{2}\log(\hat{\sigma}_{N^{\prime}}^{2})+N^{\prime}\log(n),$
(12)

where $\hat{\sigma}_{N^{\prime}}^{2}$ is the Maximum Likelihood Estimator of
$\sigma^{2}$ (see formula (13) below), and the estimator is

$\hat{N}_{\mathrm{SIC}}=\arg\min_{N^{\prime}\in\\{0,\ldots,N_{U}\\}}\mbox{SIC}(N^{\prime}),$

where $N_{U}$ is a given fixed upper bound for ${N_{0}}$. It is instructive to
understand the properties of this estimator and how to demonstrate them. For
each $N^{\prime}\leq N_{U}$, let
$\hat{\tau}_{1,N^{\prime}},\ldots,\hat{\tau}_{N^{\prime},N^{\prime}}$ be the
values of $\tau_{1},\ldots,\tau_{N^{\prime}}$ that minimise
$S_{n}(\tau_{1},\ldots,\tau_{N^{\prime}})$. We have

$\hat{\sigma}_{N^{\prime}}^{2}=\frac{1}{n}S_{n}(\hat{\tau}_{1,N^{\prime}},\ldots,\hat{\tau}_{N^{\prime},N^{\prime}}).$
(13)

Before proceeding, we remark that if $\hat{N}_{\mathrm{SIC}}=N_{0}$ (which,
under certain assumptions, is the case with probability tending to 1 with $n$,
as the rest of this section demonstrates), the change-point location
estimators
$\hat{\tau}_{1,\hat{N}_{\mathrm{SIC}}},\ldots,\hat{\tau}_{\hat{N}_{\mathrm{SIC}},\hat{N}_{\mathrm{SIC}}}$
are exactly the same as the analogous estimators obtained for $N_{0}$ known,
defined in formula (6). This illustrates more concretely the statement made in
Section 3.3 that Schwarz Information Criterion is a method for estimating not
only $N_{0}$ but also, implicitly, the change-point locations. Yao (1988)
proves the following result.

###### Theorem 3.3

Assume ${N_{0}}\leq N_{U}$ and $\tau^{0}_{j}/n\to q_{j}\in(0,1)$ for
$j=1,\ldots,{N_{0}}$ as $n\to\infty$. Assume also that
$|f_{\tau^{0}_{j}}-f_{\tau^{0}_{j}+1}|>0$ remains constant with $n$ for each
$j=1,\ldots,{N_{0}}$. Then

$P(\hat{N}_{SIC}={N_{0}})\to 1$

as $n\to\infty$.

Before outlining the proof of Theorem 3.3, we discuss a few key aspects of its
formulation. The asymptotic regime assumed in Theorem 3.3 is that the change-
point locations $\tau^{0}_{j}$ are asymptotically fixed in rescaled time,
which means that each sequence $\tau^{0}_{j}/n$ has a fixed limit in the
interval $(0,1)$. The significance of this assumption is twofold: firstly,
that the spacing between each pair of consecutive true change-points is of
asymptotic order $\asymp n$, which is the largest possible, since the length
of the data sample is $n$. Secondly, this trivially implies that the true
number ${N_{0}}$ of change-points is bounded as a function of $n$. The
assumption of a “known” upper bound (here: $N_{U}$) on the size of the model
is common in model selectors that use SIC, as, classically, it is a criterion
for model choice from among a finite set of models (Schwarz, 1978). Finally,
Theorem 3.3 requires that the size of each jump remains constant with $n$,
rather than being allowed to decrease with $n$ at a certain speed.

We now digress slightly to discuss the nature and usefulness of asymptotic
results in change-point problems, such as that in Theorem 3.3. Of course,
asymptotic considerations of this type have to be taken with a pinch of salt
by the applied user. To start with, in most scenarios involving a posteriori
change-point detection, there is only one dataset (of length $n$) available
for analysis, so there is no increasing $n$ regime to speak of in the first
place. And even if there were, it would be unrealistic to expect that any
change-points remain fixed in rescaled time – this would require data sampling
between each pair of consecutive change-points each time $n$ increased, and
would therefore preclude typical time series scenarios in which new data were
added only at the end of the sample as $n$ increased.

Still, we would argue that asymptotic results of this type can be of some use.
Firstly, they offer a basic minimum standard (some readers may like to think
of a basic “sanity check” here): we would likely be wary of an estimator that
was not able to achieve consistency even in this idealised setting, in which
results of this type are possible and achievable. Secondly, consistency
results of this type may be useful for comparing estimators; for example, if
estimator A of change-point locations were shown theoretically to be more
accurate, or work under less strict assumptions, than estimator B, we might
find it easier to justify our preference for A in practice.

Finally, asymptotic results in change-point problems can provide a practically
useful mathematical abstraction of reality. For example, one possible lead
from Theorem 3.3 is that SIC may be an accurate model selector in change-point
problems in which the data at hand exhibit relatively infrequent change-
points. Section 9 shows how this observation can manifest itself in practice.

The statement of Theorem 3.3, that $P(\hat{N}_{SIC}={N_{0}})\to 1$ as
$n\to\infty$, is typical of results on the accuracy of model selectors in
change-point problems. In fact, as will become apparent from the proof, even
more can be said in this and typically in other cases: the proof will identify
a specific subset of $\Omega$, the sample space, on which
$\hat{N}_{SIC}={N_{0}}$. The subset will be one in which the noise
$\bm{\varepsilon}$ behaves “reasonably” in the sense that a certain suitably
rich collection of linear transformations of $\bm{\varepsilon}$ uniformly does
not exhibit extreme behaviour in a certain precisely defined sense.

Key ideas of the proof of Theorem 3.3. As with many other proofs of the
validity of Schwarz Information Criterion, the argument here is split into two
parts, consisting in showing why the events (14) and (15) below are
impossible, other than on an event set whose probability tends to zero with
$n$:

$\displaystyle\hat{N}_{SIC}$ $\displaystyle<$ $\displaystyle N_{0},$ (14)
$\displaystyle\hat{N}_{SIC}$ $\displaystyle>$ $\displaystyle N_{0}.$ (15)

Before we discuss the above two parts of the proof, we firstly explain how to
understand the low-probability event on which $\hat{N}_{SIC}\neq{N_{0}}$.
Clearly, if a particular realisation of the noise vector $\bm{\varepsilon}$ is
in some sense extreme, then we cannot expect most model selection criterions
including SIC to identify the correct number of change-points. For example,
imagine that the values of $\varepsilon_{1},\ldots,\varepsilon_{m}$ for a
certain $m<n$ randomly happen to be extremely high, much higher than the
values of $\varepsilon_{m+1},\ldots,\varepsilon_{n}$. Take
$\mathbf{f}=\mathbf{0}_{n}$. Chances are that many model selection criterions,
SIC included, will mistake this random fluctuation for signal, and will
spuriously identify a change-point at location $m$, even though the true
signal $\mathbf{f}$ contains no change-points.

Fortunately, as part of the proof of Theorem 3.3, Yao (1988) shows a key
lemma, which is important in two ways: (a) it shows that the probability of
“extreme” behaviour of the noise $\bm{\varepsilon}$ (such as that described in
the toy example of the previous paragraph) tends to zero with $n$ uniformly in
a certain strong sense (more on this below), and (b) it is (practically)
sufficient for us to be outside this set of extreme noise behaviour to
guarantee $\hat{N}_{SIC}=N_{0}$ as $n\to\infty$. This key lemma appears as
Lemma 3.2 below.

###### Lemma 3.2

Suppose $\epsilon_{1},\ldots,\epsilon_{n}$ are i.i.d. normal with common mean
$0$ and variance $\sigma^{2}$. Then for any $\delta>0$, as $n\to\infty$,

$P\left(\max_{0\leq i<j\leq
n}(\epsilon_{i+1}+\ldots+\epsilon_{j})^{2}/(j-i)>2(1+\delta)\sigma^{2}\log\,n\right)\to
0.$

This is a remarkable result, and we now explain why. Lemma 3.2 shows that
outside a set whose probability tends to zero with $n$, all normalised partial
sums of an i.i.d. $N(0,\sigma^{2})$ sequence
$\epsilon_{1},\ldots,\epsilon_{n}$ of the form

$z_{i,j}=\frac{\epsilon_{i+1}+\ldots+\epsilon_{j}}{(j-i)^{1/2}}$

are bounded in absolute value by $\sigma\sqrt{2(1+\delta)\log\,n}$, for any
positive $\delta$ (obviously, the smaller the $\delta$, the higher the
probability that the bound does not hold). Each variable $z_{i,j}$ is again
distributed as $N(0,\sigma^{2})$, and there are $n^{2}$ of them. To understand
why the bound of $\sigma\sqrt{2(1+\delta)\log\,n}$ may come as a surprise, let
us step back and understand the lowest possible constant $a$ for which

$\displaystyle P(\max_{1\leq i\leq
n}|\epsilon_{i}|>\sigma\sqrt{a\,\log\,n})\to 0.$ (16)

Setting $a=2$ clearly works: using the simple Bonferroni bound and a well-
known Mills ratio inequality for the Gaussian distribution, we have

$P(\max_{1\leq i\leq
n}|\epsilon_{i}|>\sigma\sqrt{2\,\log\,n})\leq\sum_{i=1}^{n}P(|\epsilon_{i}|>\sigma\sqrt{2\,\log\,n})\leq\frac{2n}{\sigma\sqrt{2\,\log\,n}}\phi(\sqrt{2\,\log\,n}),$

where $\phi()$ is the pdf of the standard normal. This simplifies to
$\frac{2}{\sigma\sqrt{2\,\log\,n}}$, which tends to 0 with $n$. Using more
advanced techniques (see e.g. Example 2.29 and Exercise 2.11 in Wainwright
(2019)), it can be shown that no constant $a<2$ leads to (16) being satisfied,
and therefore $a=2$ is the lowest constant for which (16) holds. Now, Lemma
3.2 tells us that the bound $\sigma\sqrt{2(1+\delta)\log\,n}$ for all $n^{2}$
variables $|z_{i,j}|$ is only very slightly higher than the lowest possible
bound $\sigma\sqrt{2\log\,n}$ for the $n$ variables $|\epsilon_{i}|$. Clearly,
if all $z_{i,j}$ were independent, then by the above discussion, the lowest
possible bound for all the $|z_{i,j}|$’s would have to be
$\sigma\sqrt{2\log\,(n^{2})}=\sigma\sqrt{4\log\,n}\gg\sigma\sqrt{2(1+\delta)\log\,n}$
(bear in mind that any positive $\delta$ is permitted, and we therefore think
of $\delta$ being small). However, the particular dependence structure in the
ensemble of variables $\\{z_{i,j}\\}_{i,j=1}^{n}$ causes the bound for this
ensemble to be “only” $\sigma\sqrt{2(1+\delta)\log\,n}$. For completeness, we
note that similar bounds for more general sub-Gaussian distributions were
independently obtained later by Shao (1995).

We now outline the two parts of the proof of Theorem 3.3: showing the
impossibility of (14) and of (15).

Showing the impossibility of $\hat{N}_{SIC}<N_{0}$. If we were to have
$\hat{N}_{SIC}<N_{0}$, there would be at least one true change-point, say
$\tau^{0}_{r}$ for $1\leq r\leq N_{0}$, with no estimated change-points in its
neighbourhood $[r_{s},r_{e}]$. We could choose $r_{s},r_{e}$ such that
$r_{e}-r_{s}\asymp n$ due to the assumption that $\tau^{0}_{j}/n\to q_{j}$.
Consider the quantity $\hat{\sigma}^{2}_{\hat{N}_{SIC}}$, defined in (13) and
needed in the evaluation of $\mbox{SIC}(\hat{N}_{SIC})$, defined in (12). The
presence of the change-point $\tau^{0}_{r}\in[r_{s},r_{e}]$ would cause the
contribution of the sum of squares

$\sum_{i=r_{s}}^{r_{e}}(X_{i}-\bar{X}_{r_{s}:r_{e}})^{2}$

to $\hat{\sigma}^{2}_{\hat{N}_{SIC}}$, see formulae (5) and (13), to be,
stochastically, much larger than $\sigma^{2}(r_{e}-r_{s}+1)$. As
$r_{e}-r_{s}\asymp n$, this increased contribution would take place over a
relatively long interval, and therefore the remaining components of the sum of
squares (5) would not be able to compensate for it and as a result,
$\hat{\sigma}^{2}_{\hat{N}_{SIC}}$ would end up much larger than $\sigma^{2}$.
This would in turn mean that $\mbox{SIC}(N^{\prime})$ could not be minimised
at $\hat{N}_{SIC}$, which leads to a contradiction.

To formalise this reasoning, Yao (1988) shows two lemmas, which we give below
without proofs.

###### Lemma 3.3

As $n\to\infty$, we have
$0\leq\sum_{i=1}^{n}\varepsilon_{i}^{2}-n\hat{\sigma}^{2}_{N_{0}}=O_{p}(\log\,n)$.

This lemma says that at the true number of change-points $N_{0}$, the
estimator $\hat{\sigma}^{2}_{N_{0}}$ is close to and bounded from above by
$\frac{1}{n}\sum_{i=1}^{n}\varepsilon_{i}^{2}$ in probability.

###### Lemma 3.4

For every $N<N_{0}$, there exists $\nu>0$ such that
$P(\hat{\sigma}^{2}_{N}>\sigma^{2}+\nu)\to 1$ as $n\to\infty$.

This lemma says that if $N<N_{0}$, the estimator $\hat{\sigma}^{2}_{N}$ is
significantly larger than $\sigma^{2}$ in probability. Yao (1988) then argues
that by Lemma 3.3, $\hat{\sigma}^{2}_{N_{0}}\to\sigma^{2}$ in probability, so
by Lemma 3.4 and the definition of $\mbox{SIC}(N)$, $P(\hat{N}_{SIC}\geq
N_{0})\to 1$.

Even though this is not discussed in Yao (1988), it is worth taking a look at
other penalty values that would lead to similar conclusions. Suppose that we
extend SIC to the more general form

$\mbox{IC}_{p(n)}(N^{\prime})=\frac{n}{2}\log(\hat{\sigma}_{N^{\prime}}^{2})+N^{\prime}p(n),$
(17)

where $p(n)$ is a general function of $n$ (not necessarily $\log\,n$). With
$\mbox{IC}_{p(n)}$ in place of SIC, by the above discussion, for $N<N_{0}$
(for the $\nu$ from Lemma 3.4 and with probability approaching 1), we would
have

$\mbox{IC}_{p(n)}(N)-\mbox{IC}_{p(n)}(N_{0})\geq\frac{n}{2}\log\left(1+\frac{\nu}{2\sigma^{2}}\right)+(N-N_{0})p(n),$
(18)

which is guaranteed to be positive (uniformly over $N$, since $N_{0}$ is
finite) as long as $p(n)=o(n)$.

Showing the impossibility of $\hat{N}_{SIC}>N_{0}$. In typical proofs of the
consistency of SIC, showing the impossibility of this side inequality is
usually more involved than the other side. This is fundamentally because it
can be difficult to exclude the possibility that $\hat{N}_{SIC}$ is “close to
but larger than $N_{0}$”. This could potentially be caused by
$\hat{\sigma}^{2}_{N}$ being so much below $\sigma^{2}$ for $N>N_{0}$ that the
penalty (which increases in $N$) would not be able to compensate for this
drop. Fortunately, with a high probability, this does not happen: if $N$ is
larger than but close to $N_{0}$, $\hat{\sigma}^{2}_{N}$ still does not
underestimate $\sigma^{2}$ by too much. Heuristically, this is for the
following reason. Suppose we have already estimated the $N_{0}$ true change-
points accurately, but we also add a small number of spuriously estimated
change-points, and denote the overall number by $N$. In these circumstances,
$\hat{\sigma}^{2}_{N}$ would be slightly smaller than
$\hat{\sigma}^{2}_{N_{0}}$ (the latter being a good estimator of
$\sigma^{2}$), but not by much, as the extra spurious estimated change-points
(being located, by definition, over a section of the data where the signal is
constant) would not induce much difference in the residual sum of squares, a
key ingredient of $\hat{\sigma}^{2}_{N}$.

More formally, to quantify the magnitude of $\hat{\sigma}^{2}_{N}$ for
$N>N_{0}$ with the precision required to show the impossibility of
$\hat{N}_{SIC}>N_{0}$, Yao (1988) introduces the Lemma below, whose proof is
considerably more involved than the proof of Lemmas 3.3 and 3.4.

###### Lemma 3.5

For every $N$ ($N_{0}<N\leq N_{U}$) and for any $\nu>0$, with probability
approaching 1,

$0\leq\sum_{i=1}^{n}\varepsilon_{i}^{2}-n\hat{\sigma}^{2}_{N}\leq\\{\nu+(N-N_{0}-1)2(1+\nu)\\}\sigma^{2}\log\,n.$

This lemma says that for $N>N_{0}$, as long as $N\leq N_{U}$, a pre-set upper
bound, the estimator $\hat{\sigma}^{2}_{N}$ is lower bounded by the quantity

$\frac{1}{n}\sum_{i=1}^{n}\varepsilon_{i}^{2}-\\{\nu+(N-N_{0}-1)2(1+\nu)\\}\frac{\sigma^{2}}{n}\log\,n.$

In other words, Lemma 3.5 ensures that $\hat{\sigma}^{2}_{N}$, for
$N_{0}<N\leq N_{U}$, is not too small, and therefore that there is no danger
of $\mbox{SIC}(N)$ being minimised at $N>N_{0}$. More specifically, with Lemma
3.5 in place, Yao (1988) shows the impossibility of $\hat{N}_{SIC}>N_{0}$ as
follows. Note that $\sum_{i=1}^{n}\varepsilon_{i}^{2}>n(\sigma^{2}-\nu)$ with
probability approaching 1 for any $\nu>0$. By Lemmas 3.3 and 3.5, for
$N_{0}<N\leq N_{U}$, for any $\nu>0$, with probability approaching 1,

$\displaystyle 2\\{\mbox{SIC}(N)-\mbox{SIC}(N_{0})\\}$ $\displaystyle\geq$
$\displaystyle
n\,\log\,\hat{\sigma}_{N}^{2}-n\log\left\\{\sum_{i=1}^{n}\varepsilon_{i}^{2}/n\right\\}+2(N-N_{0})\log\,n$
$\displaystyle=$ $\displaystyle
n\log\left\\{1-\frac{\sum_{i=1}^{n}\varepsilon_{i}^{2}-n\hat{\sigma}^{2}_{N}}{\sum_{i=1}^{n}\varepsilon_{i}^{2}}\right\\}+2(N-N_{0})\log\,n$
$\displaystyle\geq$ $\displaystyle
n\log\left\\{1-\frac{\\{\nu+(N-N_{0}-1)2(1+\nu)\\}\sigma^{2}\log\,n}{n(\sigma^{2}-\nu)}\right\\}$
$\displaystyle+$ $\displaystyle 2(N-N_{0})\log\,n,$

which, using $\log(1-x)>(1+\nu)(-x)$ for small $x>0$, is greater than

$-(1+\nu)\frac{\\{\nu+(N-N_{0}-1)2(1+\nu)\\}\sigma^{2}\log\,n}{\sigma^{2}-\nu}+2(N-N_{0})\log\,n,$
(19)

for large $n$. Since (19) is positive for sufficiently small $\nu>0$ (as it is
positive for $\nu=0$ and continuous at that point), we have
$P(\mbox{SIC}(N)-\mbox{SIC}(N_{0})>0)\to 1$, completing the proof.

Let us again pick up the discussion of the more general penalty
$\mbox{IC}_{p(n)}$ defined in (18). For what other values of $p(n)$ could we
obtain the conclusion that the equivalent of (19) is positive, and therefore
that $\hat{N}_{SIC_{p(n)}}>N_{0}$ is impossible?

Clearly, from (19), $p(n)=\log\,n$ is practically the smallest penalty that
works, but it is interesting to observe that anything larger than that would
also do the job. Therefore, combining this discussion with the restriction on
$p(n)$ obtained through formula (18), we can see that $p(n)$ in the region

$p(n)\in[\log\,n,\bar{p}(n)],$ (20)

where $\bar{p}(n)$ is any sequence satisfying $\bar{p}(n)=o(n)$, permits the
same proof of consistency for $N_{0}$ of the $\mbox{IC}_{p(n)}$ criterion as
the one we saw in this section for the classical SIC criterion. From this
discussion, we can also see that the classical SIC penalty $p(n)=\log n$ is
fundamentally the smallest penalty for which this argument for consistency
works, and, from this point of view, the classical SIC leads to the smallest
permitted amount of penalisation.

It is interesting to observe that the smallest permitted penalty (here:
$p(n)=\log n$) is also the most desirable, because it gives us the highest
power in detecting change-points. This is despite the fact that often,
consistency proofs for the number of change-points can be substantially
simplified for penalties higher than the lowest permitted one. For an example
of a simpler consistency proof than that of Theorem 3.3 – because it uses a
larger penalty – see Fryzlewicz (2014), Theorem 3.3.

Yao and Au (1989) show how to adapt the penalty to estimate $N_{0}$
consistently for noise for which only a finite number of moments are finite.
Lee (1995) generalises SIC to situations in which the minimum spacing between
consecutive change-points is smaller than in Yao (1988).

###  3.5 When $N_{0}$ Is Unknown: Penalised Approaches Derived from Other
Information Criteria

Unless specified otherwise, estimation of the relevant quantities throughout
this section is performed via the minimisation of the criteria with respect to
their parameters.

Akaike Information Criterion. Besides the Bayesian (Schwarz) Information
Criterion, the other arguably the most commonly used model selection criterion
in statistics generally is Akaike Information Criterion (AIC). To estimate the
number $N_{0}$ of change-points via AIC, define

$\mbox{AIC}(N^{\prime})=\frac{n}{2}\log(\hat{\sigma}^{2}_{N^{\prime}})+2N^{\prime},$
(21)

which leads to the AIC estimator

$\hat{N}_{AIC}=\arg\min_{N^{\prime}\in\\{0,\ldots,N_{U}\\}}\mbox{AIC}(N^{\prime}).$

For the purpose of model selection in change-point problems, AIC appears to
have been used less extensively than SIC in the literature. A contributing
factor to this relative lack of popularity may be the well-documented lack of
statistical consistency of the AIC criterion in selecting the true model in
other problems, such as linear model choice in regression. Indeed, condition
(20) excludes the AIC penalty from the range of penalties that lead to
consistency for $N_{0}$. In practice, since AIC penalises models with more
change-points less heavily than SIC, we can expect it to overestimate the
number of changes – and substantial over-estimation is often seen in practice.
For Gaussian models, AIC is equivalent to Mallows’ $C_{p}$ (Birgé and Massart,
2001; Boisbunon et al., 2014).

Modified Bayesian Information Criterion. Zhang and Siegmund (2007) propose the
so-called modified BIC (mBIC), constructed as follows. Assume a Bayesian set-
up in which we place the uniform prior on the number of change-points, on
their locations, on their sizes, and on the overall mean of the signal
$\mathbf{f}$. Assume further that $\mbox{Var}(\varepsilon_{i})$ is known to be
1. Then Theorem 1 in Zhang and Siegmund (2007) shows that the log Bayes factor
of a model with $N^{\prime}$ change-point versus the no-change-point model is
asymptotically equal to the maximised log-likelihood under a model with
$N^{\prime}$ change-points (that is, the log-likelihood of the most likely
model with $N^{\prime}$ change-points), minus the following term, which can be
interpreted as the mBIC penalty.

$\mbox{mBIC}_{\mbox{pen}}(N^{\prime},q_{1},\ldots,q_{N^{\prime}})=\frac{1}{2}\left\\{3N^{\prime}\log\,n+\sum_{j=1}^{N^{\prime}+1}\log(q_{j}-q_{j-1})\right\\},$
(22)

where $0=q_{0}<q_{1}<\ldots<q_{N^{\prime}}<q_{N^{\prime}+1}=1$ are the
postulated change-point locations rescaled to lie in the interval $[0,1]$,
i.e. $q_{j}=\tau_{j}/n$. Therefore, unlike the standard BIC (SIC) penalty
discussed in Section 3.4, the mBIC penalty penalises not just the number of
change-points, but also their relative locations.

As done in Zhang and Siegmund (2007), it is interesting to investigate how the
second, location-dependent term behaves in two extreme cases: for rescaled
change-point locations placed at equal intervals over $[0,1]$, and clustered
together as closely as possible (e.g.
$\tau_{1}=nq_{1}=1,\ldots,\tau_{N^{\prime}}=nq_{N^{\prime}}=N^{\prime}$).
Simple algebra shows that in the former case, we have

$\sum_{j=1}^{N^{\prime}+1}\log(q_{j}-q_{j-1})=-(N^{\prime}+1)\log(N^{\prime}+1),$
(23)

whereas in the latter case, we have

$\sum_{j=1}^{N^{\prime}+1}\log(q_{j}-q_{j-1})\approx-N^{\prime}\log\,n,$ (24)

which leads to, in this special case,
$\mbox{mBIC}_{\mbox{pen}}(N^{\prime},q_{1},\ldots,q_{N^{\prime}})=\mbox{SIC}(N^{\prime})=N^{\prime}\log(n)$.
As, typically, $N^{\prime}\log\,n>(N^{\prime}+1)\log(N^{\prime}+1)$, the mBIC
penalty penalises change-points located close to each other less severely than
it penalises uniformly spread ones. This is somehow counter-intuitive. This
particular aspect of the penalty is not discussed in Zhang and Siegmund
(2007). However, the authors mention that change-points located too close to
each other are precluded by their assumption regarding the minimum distance
between each pair of consecutive change-points, which is assumed to be $\asymp
n$. This leads them to conclude that the case described in formula (23) is the
prevalent one, and therefore it is the first term in (22) that typically
dominates, which makes the mBIC penalty close to
$\frac{3}{2}N^{\prime}\log\,n$. However, they also remark that $\log\,n$ is
such a slowly-growing function that the second term in (22) “can still make a
significant contribution to the penalty”. An extension to the case of an
unknown $\mbox{Var}(\varepsilon_{i})$ is also described. Du et al. (2016)
propose a related, empirical Bayes approach to change-point estimation in a
marginal likelihood framework, which, disregarding any information supplied by
the prior, also leads to a penalty which prefers change-point estimates
located close to each other.

In practice, Zhang and Siegmund (2007) propose to estimate a change-point
solution path via an algorithm related to Circular Binary Segmentation (Olshen
et al., 2004), and only use mBIC to choose from among the candidate models
along the thus-constructed solution path. This way of proceeding goes some way
towards mitigating the preference of mBIC for models with change-points
located close to each other (among models with the same numbers of change-
points): this potential issue will be of no consequence if the solution path
includes no such models. Our experience is that, if we directly use the mBIC,
that is we use it also to find the best segmentation with a given number of
change-points, then it performs less well due to its preference to fit shorter
segments to the data.

Minimum Description Length. Davis et al. (2006) propose to use the Minimum
Description Length (MDL) principle to find segments of homogeneity within a
piecewise-stationary autoregressive process. Based on this work, we examine
how this model selection principle simplifies in our signal + Gaussian i.i.d.
noise setting. Heuristically speaking, the code length of a model is the
number of bits used to encode the parameters, plus the number of bits used to
encode the residuals. We refer to Rissanen (1989, 2007) for the basics of
information-theoretic encoding of random variables and parameters, which we do
not reproduce here. To encode the parameters of the model, we need to include
the following components.

  * •

The postulated number $N^{\prime}$ of change-points. To encode an unbounded
integer $N^{\prime}$, approximately $\log_{2}(N^{\prime}+1)$ bits are
required. We note that Davis et al. (2006) propose to use
$\log_{2}N^{\prime}$, but this technically cannot be correct if the search for
the number of models is to permit a model with $N^{\prime}=0$ (which it does
in Davis et al. (2006)), as $\log_{2}0=-\infty$, so the procedure would always
end up selecting the null model. Davis (2019, personal communication)
clarifies that this term should equal, in our notation,
$\log_{2}(N^{\prime}+1)$. We also note that despite being naturally bounded by
$n-1$, the postulated number $N^{\prime}$ of change-points is treated as
unbounded here as in practice we usually have $N^{\prime}\ll n$, so typical
values for $N^{\prime}$ are far from this upper bound.

  * •

The postulated segment lengths. Each of the $N^{\prime}+1$ segment lengths is
naturally bounded by $n$ (and this upper bound can easily be approached in
practice). Therefore, each of the first $N^{\prime}$ segments requires
$\log_{2}n$ bits to be encoded. We do not need to encode the length of the
final segment as the segment lengths add up to $n$. Therefore encoding the
segment lengths requires $N^{\prime}\log_{2}n$ bits.

  * •

The maximum likelihood estimates of $\sigma^{2}$ and $f_{\tau^{0}_{j}+1}$. To
encode the maximum likelihood estimates of the signal variance $\sigma^{2}$
and the signal levels $f_{\tau^{0}_{j}+1}$ for $j=0,\ldots,N^{\prime}$, we use
a result by Rissanen (1989) which states that maximum likelihood estimates of
a real parameter computed from $n^{\prime}$ observations can be effectively
encoded with $\frac{1}{2}\log_{2}n^{\prime}$ bits. The signal level
$f_{\tau^{0}_{j}+1}$ is estimated on a data stretch of length $l_{j+1}-l_{j}$,
where $l_{j}$ are the postulated change-point locations. The signal variance
$\sigma^{2}$ is estimated on the entire dataset. Therefore to encode these
parameters, we require
$\frac{1}{2}\sum_{j=0}^{N^{\prime}}\log_{2}(l_{j+1}-l_{j})+\frac{1}{2}\log_{2}n$
bits.

  * •

The residuals from each fitted model. From Shannon’s results in information
theory, Rissanen (1989) shows that the code length of the residuals is given
by the negative of the log-likelihood (where the base 2 logarithm is used) of
the corresponding fitted model.

Putting together the above pieces, ignoring the additive constants in the same
way as it is done in Schwarz Information Criterion (Yao, 1988), and switching
to base $e$ for the logarithm, which is equivalent to just multiplying the
criterion by a constant, we obtain the MDL criterion as

$\displaystyle\mbox{MDL}(N^{\prime},l_{1}<\ldots<l_{N^{\prime}})$
$\displaystyle=$
$\displaystyle\frac{n}{2}\log(\hat{\sigma}^{2}_{N^{\prime}})+\log(N^{\prime}+1)+N^{\prime}\log\,n$
$\displaystyle+$
$\displaystyle\frac{1}{2}\sum_{j=0}^{N^{\prime}}\log(l_{j+1}-l_{j})$
$\displaystyle=:$
$\displaystyle\frac{n}{2}\log(\hat{\sigma}^{2}_{N^{\prime}})+\mbox{MDL}_{\mbox{pen}}(N^{\prime},l_{1}<\ldots<l_{N^{\prime}}).$

As with the other penalties described in this section, we note that the MDL
criterion can also be used for the selection of $N^{\prime}$ only, for any
(other) solution path algorithm that produces change-point location estimators
for any given $N^{\prime}$.

If the MDL criterion is used to select $N^{\prime}$ as well as the change-
point locations, it is interesting to observe that from the concavity of the
logarithm, as with the mBIC criterion described earlier, the MDL also
encourages estimated change-points that cluster together, rather than being
spread out.

Finally, we note that the thus-defined MDL penalty is uniformly larger than
the SIC penalty defined earlier. If we compare it to mBIC, using
$q_{j}=l_{j}/n$, it is easy to see that

$\mbox{mBIC}_{\mbox{pen}}(N^{\prime},q_{1},\ldots,q_{N^{\prime}})-\mbox{MDL}_{\mbox{pen}}(N^{\prime},l_{1}<\ldots<l_{N^{\prime}})=-\log(N^{\prime}+1)+f(n),$

where the $f(n)$ term does not depend on the postulated number of change-
points $N^{\prime}$ – so from the point of view of minimisation, the
$\log(N^{\prime}+1)$ term is the only term that separates the two penalties.

Potts functional. Boysen et al. (2009) consider minimisation of what they
refer to as the Potts functional, defined by

$H_{\gamma}(\mathbf{f}^{\prime},\mathbf{X})=\frac{1}{n}\sum_{i=1}^{n}(X_{i}-f^{\prime}_{i})^{2}+\gamma|J(\mathbf{f}^{\prime})|,$
(25)

where $J(\mathbf{f}^{\prime})=\\{i\,:\,1\leq i\leq n-1,\,f^{\prime}_{i}\neq
f^{\prime}_{i+1}\\}$ is the set of change-points in the candidate fit
$\mathbf{f}^{\prime}$, and $|\cdot|$ is the cardinality of its argument
(denoted by $N^{\prime}$ earlier in this chapter). Theoretically, one of the
requirements on $\gamma$ is that $\gamma n/\log\,n\to\infty$; however, in
practical terms, the authors recommend that $\gamma$ be set to
$2.5\,\hat{\sigma}^{2}\log\,n/n$, where $\hat{\sigma}^{2}$ is a consistent
estimator of $\sigma$. This, however, means that this approach is very closely
related to SIC: by the arguments outlined in Section 3.3, if $\sigma^{2}$ were
assumed known, in the notation of Section 3.1, SIC for the estimation of the
number and the locations of change-points would be

$\mathrm{SIC}(N^{\prime},l_{1}<\ldots<l_{N^{\prime}})=\frac{1}{n}S_{n}(l_{1},\ldots,l_{N^{\prime}})+2\sigma^{2}N^{\prime}\log\,n/n,$
(26)

which agrees with (25) except the penalty uses the multiplicative constant of
2, rather than $2.5$ (if $\sigma^{2}$ were unknown, it would be advisable to
use the likelihood form of SIC as in formula (12), which does not require the
knowledge of $\sigma^{2}$, or alternatively use a consistent estimator of
$\sigma^{2}$ in (26), which would bring (25) and (26) in even closer
agreement).

We conclude this section with pointers to other material on penalised model
selection in this monograph. A future chapter of this book will describe
penalised model selection in parametric multiple change-point problems
involving settings beyond the Gaussian i.i.d. model. Section 6 later in this
chapter describes $L_{1}$-penalised methods for model selection in multiple
change-point problems; the reason for the separate treatment of this material
is the distinct computational nature of these methods, which tend to be
closely linked to the lasso problem. Finally, Section 8 covers methods in
which the penalties themselves are chosen from the data.

###  3.6 Computation for Estimators Defined as Optima

Segment Neighbourhood

We now turn to the problem of how we can calculate the estimators we have
described. First we return to estimators that assume $N_{0}$ is known, such as
the one based on minimising the residual sum of squares

$(\hat{\tau}_{1},\ldots,\hat{\tau}_{N_{0}})=\arg\min_{l_{1}<\ldots<l_{N_{0}}}S_{n}(l_{1},\ldots,l_{N_{0}}),$
(27)

where the residual sum of squares function $S_{n}(l_{1},\ldots,l_{m})$ is
defined in formula (5).

A naive approach to calculating $(\hat{\tau}_{1},\ldots,\hat{\tau}_{N_{0}})$
would be to evaluate $S_{n}(l_{1},\ldots,l_{N_{0}})$ for all possible choices
of $l_{1}<l_{2}<\cdots<l_{N_{0}}$. However, the computation involved is
$O(n^{N_{0}})$, which is impracticable for moderate to large values of
$N_{0}$. Fortunately there are dynamic programming algorithms that can solve
the minimisation problem (27) with a cost that is $O(N_{0}n^{2})$.

These algorithms can be applied to a broader class of optimisation problems
than just minimising the residual sum of squares (27), including many change-
point estimators based on maximising some form of log-likelihood. To make this
generality clear, we will introduce segment cost functions,
$C({X}_{1:n},k,l)$, for $l>k$, which can be viewed as some measure of fit,
such as based on a negative log-likelihood, for fitting data from time $k+1$
to time $l$ as a single segment. For the residual sum of squares,

$C({X}_{1:n},k,l)=\sum_{i=k+1}^{l}(X_{i}-\bar{X}_{(k+1):l})^{2}.$

In most cases these cost functions will depend on the data $X_{1:n}$ just
through the data within the segment, $X_{(k+1):l}$; though this is not needed
for the algorithms we will describe to work.

We require that the minimisation problem we wish to solve involves minimising
a sum of the segment costs. That is, we wish to minimise, for some specified
$m$,

$Q_{n}(l_{1},\ldots,l_{m})=\sum_{j=1}^{m+1}C({X}_{1:n},l_{j-1},l_{j}),$

where as above we require $0<l_{1}<l_{2},\ldots<l_{m}<n$ and we set $l_{0}=0$
and $l_{m+1}=n$.

It is helpful to consider a wider class of minimisation problems,

$Q_{t}(l_{1},\ldots,l_{m})=\sum_{j=1}^{m}C({X}_{1:n},l_{j-1},l_{j})+C({X}_{1:n},l_{m},t),$

where we now consider only segmenting the data up to time $t$ and require
$l_{m}<t$. Denote the minimum of this function as

$Q_{t,m}=\min_{l_{1}<\ldots<l_{m}<t}Q_{t}(l_{1},\ldots,l_{m}).$

If $m\geq t$ then, as there is no possible segmentation of $t$ data points
into $m+1$ segments, we set $Q_{t,m}=\infty$.

Define these constants, $Q_{t,m}$, for all $t=1,\ldots,n$ and
$m=1,\ldots,N_{0}$. The Segment Neighbourhood algorithm (Auger and Lawrence,
1989) is based on the idea that if we are forced to place the last change-
point before $t$ at some location ($k$, say), then there is a simple recursion
for the minimum cost of segmenting the data with this additional constraint.
To see this let

$Q_{t,m,k}=\min_{l_{1}<\ldots<l_{m-1}<k}Q_{t}(l_{1},\ldots,l_{m-1},k).$

Then we have for $m\geq 1$

$Q_{t,m,k}=Q_{k,m-1}+C(X_{1:n},k,t),$

the minimum cost of fitting the data ${X}_{1:k}$ with $m-1$ change-points,
plus the segment cost for the data ${X}_{(k+1):t}$. Thus we can calculate
$Q_{t,m}$ by considering each possible value for $k$ and taking the minimum of
the costs

$Q_{t,m}=\min_{k}Q_{t,m,k}=\min_{k}\\{Q_{k,m-1}+C(X_{1:n},k,t)\\}.$ (28)

Finally, it is trivial to calculate $Q_{t,m}$ when $m=0$ as
$Q_{t,0}=C(X_{1:n},0,t)$. Thus we can calculate $Q_{t,0}$ for $t=1,\ldots,n$,
and then solve (28) in turn for increasing $m=1,\ldots,N_{0}$ and, for each
$m$, for $t=m+1,\ldots,n$.

Once we have solved the recursion it is simple to extract the optimal set of
change-points. The optimal location of the last change-point for a
segmentation with $m$ change-points is just

$\hat{\tau}_{m}=\arg\min_{k}\\{Q_{k,m-1}+C(X_{1:n},k,n)\\}.$

Then we can recurse backwards in time, so that if the optimal location of the
$j$th change-point is $\hat{\tau}_{j}$ then, for $j>1$, the optimal location
of the $(j-1)$th change-point is

$\hat{\tau}_{j-1}=\arg\min_{k}\\{Q_{k,j-2}+C(X_{1:n},k,\hat{\tau}_{j})\\}.$

It is straightforward to calculate the computational complexity of this
algorithm. One implementation is to calculate and store $C(X_{1:n},k,l)$ for
all $0<k<l\leq n$ – which for most models will have an $O(n^{2})$ cost for
both space and time. For each value of $t$ and $m$, solving (28) involves
$O(t)$ calculations, thus solving the recursion for any $m$ will have an
$O(n^{2})$ computational cost. As the recursions need to be solved for
$m=1,\ldots,N_{0}$ this leads to an overall computational cost of
$O(N_{0}n^{2})$. Furthermore it is often possible to implement the algorithm
with the same time complexity, but without storing all segment costs if needed
(see Maidstone et al., 2017).

One nice feature of the Segment Neighbourhood algorithm is that it not only
gives the optimal location of the change-points for a segmentation with
$N_{0}$ change-points, but also for those with $1,2,\ldots,N_{0}-1$. Thus, if
the number of change-points is unknown, we can run Segment Neighbourhood to
find the best segmentations with $1,2,\ldots,N_{U}$ change-points for some
suitable chosen upper bound, $N_{U}$, on the number of change-points we think
are present. We can then estimate $N_{0}$ from this solution path of best
segmentations that minimise our cost for differing numbers of change-points,
using an appropriate criterion such as one of the information criteria
discussed in Section 3.3, or one of the adaptive procedures that we will
introduce in Section 8. However, if we intend to estimate $N_{0}$ using an
information criterion that is linear in the number of change-points, there is
a more efficient dynamic programming approach that can jointly estimate the
number and locations of the change-points under such a criterion: Optimal
Partitioning.

Optimal Partitioning

As described in Section 3.3, when we do not know the number of change-points,
$N_{0}$, we can estimate it by minimising an information criterion. For
example, a common information criterion for the change-in-mean problem is
Schwarz Information Criterion under a model where the variance is assumed
known (which is equivalent to using the criterion defined in Equation 26).
This is

$\mathrm{SIC}(m;l_{1},\ldots,l_{m})=\frac{1}{\sigma^{2}}S_{n}(l_{1},\ldots,l_{m})+m(2\log
n),$ (29)

where $S_{n}(l_{1},\ldots,l_{m})$ is the residual sum of squares for the
segmentation with change-points at $l_{1}<\cdots<l_{m}$.

Using the notation introduced earlier in this section, $Q_{n,m}=\min
S_{n}(l_{1},\ldots,l_{m})$, and we can re-write the criterion as estimating
the number of change-points by the value that minimises
$(1/\sigma^{2})Q_{n,m}+m(2\log n)$. As discussed above, we can calculate this
by using the Segment Neighbourhood algorithm to find $Q_{n,m}$ for
$m=0,1,\ldots,N_{\max}$ for some suitable choice of $N_{\max}$. However, it
turns out that there is a faster algorithm that we can use to directly
minimise (29) with respect to both the number and location of the change-
points.

This algorithm, called Optimal Partitioning (Jackson et al., 2005), is a
dynamic programming algorithm that can be used to minimise functions of the
form

$Q_{n,\lambda}(m;l_{1},\ldots,l_{m})=\sum_{j=1}^{m+1}C({X}_{1:n},l_{j-1},l_{j})+m\lambda$
(30)

with respect to both the number of change-points, $m$, and their locations,
$l_{1}<\cdots<l_{m}$. This is often called a penalised cost, and the key is
that it involves a term that is the sum of segment costs plus a penalty that
is linear in the number of change-points.

For any $t=1,\ldots,n$ and $\lambda>0$ it is helpful to define

$Q_{t,\lambda}=\min_{m;l_{1}<\cdots<l_{m}<t}\left\\{\sum_{j=1}^{m}C({X}_{1:n},l_{j-1},l_{j})+C({X}_{1:n},l_{m},t)+m\lambda\right\\},$

the minimum value of the penalised cost for segmenting the data up to time
$t$. This definition is a slight abuse of notation, given our earlier
definition of $Q_{n,m}$, but the context should mean that this causes no
confusion.

The idea of Optimal Partitioning is that we can derive, and solve, a recursion
for $Q_{t,\lambda}$. As with Segment Neighbourhood, we first consider fixing
the location of the most recent change-point prior to $t$. If this is fixed at
$k$, it is immediate that the minimum penalised cost for segmenting the data
to time $t$ given the last change is at $k$ will be
$Q_{k,\lambda}+C({X}_{1:n},k,t)+\lambda$. Thus we have

$Q_{t,\lambda}=\min\left\\{C({X}_{1:n},0,t),\min_{k=1,\ldots,t-1}\left\\{Q_{k,\lambda}+C({X}_{1:n},k,t)+\lambda\right\\}\right\\}.$

Here the first term on the right-hand side corresponds to no change-points,
and the remaining terms correspond to one or more change-points with the most
recent change being at time $k$. If we define $Q_{0,\lambda}=-\lambda$ this
can be simplified to

$Q_{t,\lambda}=\min_{k=0,\ldots,t-1}\left\\{Q_{k,\lambda}+C({X}_{1:n},k,t)+\lambda\right\\}.$
(31)

This recursion can be solved for $t=1,\ldots,n$. Once it has been solved it is
straightforward to obtain the estimated change-point locations. To describe
this it is simplest to denote the change-points in reverse order, so
$\hat{\tau}_{1}>\hat{\tau}_{2}>\cdots>\hat{\tau}_{m}$. The location of the
last change-point will be

$\hat{\tau}_{1}=\arg\min\left\\{Q_{k,\lambda}+C({X}_{1:n},k,t)+\lambda\right\\}.$

Then if we have estimated the $j$th to last change-point, $\hat{\tau}_{j}>0$
we set

$\hat{\tau}_{j+1}=\arg\min_{k=0,\ldots,\hat{\tau}_{j}-1}\left\\{Q_{k,\lambda}+C({X}_{1:n},k,\hat{\tau}_{j})+\lambda\right\\}.$

We repeat this until we obtain $\hat{\tau}_{m+1}=0$, which corresponds to no
earlier change-points.

The advantage of Optimal Partitioning over Segment Neighbourhood is that it
has a lower computational complexity. For most models calculating
$C({X}_{1:n},k,t)$ can be done with cost $O(1)$ regardless of $k$ and $t$ (for
example after storing appropriate summaries of the data as discussed in
Section 1.2 of Fearnhead and Fryzlewicz (2022)). Assuming this is the case,
then solving the Optimal Partitioning recursion (31) for one value of $t$ has
complexity that is linear in $t$. As we need to solve the recursion for
$t=1,\ldots,n$ this leads to an $O(n^{2})$ computational complexity.

When minimising the penalised cost, the choice of $\lambda$ can greatly affect
the number of estimated change-points and hence how accurate the estimated
segmentation is. Whilst there are default choice for $\lambda$, for example
based on various information criteria described in Sections 3.3–3.5, these can
perform poorly in real applications when the strong assumptions that are
required by the theory supporting the information criteria does not hold. In
such situations it is recommended to compare how the segmentations vary as we
vary $\lambda$. This can be done efficiently using the CROPS algorithm of
Haynes et al. (2017b).

Pruning

Figure 4:  Plot of simulated data and underlying mean (left) and plot of
conditional cost (see the text for the definition of this term) in the Optimal
Partitioning recursion (right). The conditional cost has been shifted so its
minimum is at 0, and the horizontal red line is at $\lambda$ and determines
the condition for the PELT pruning.

If we consider the Optimal Partitioning recursion, then at time $t$ we are
minimising over terms, $Q_{k,\lambda}+C({X}_{1:n},k,t)+\lambda$, which we will
call the conditional cost for the most recent change at $k$. Figure 4 shows
example data and a plot of this conditional cost minus its minimum value. For
this example, $t=200$ and the data has had a change at time 160. We notice
that for locations of the most recent change-point, $k$, prior to 160 the
conditional cost tends to increase as $k$ decreases. This is intuitive as the
conditional cost includes a segment from $k$ to time 200, and these will
ignore the change in the data. The increasing conditional cost is indicative
of the increasing evidence that there has been a change-point since time $k$.

A natural question is whether we could use such information to speed up the
Optimal Partitioning algorithm. That is, if the conditional cost at some
location $k$ is sufficiently high, does that mean that $k$ can never be the
best choice of the most recent change-point for future time-steps, regardless
of what the future data is? If this was the case we could “prune” $k$ and
reduce the set of possible most recent change-points we search over at future
steps of the algorithm.

Killick et al. (2012) show that in many instances this is the case. In
particular if the segment cost is such that adding change-points can only ever
reduce the unpenalised cost, then we can prune any $k$ for which the
conditional cost is greater than $\lambda$ above the minimum conditional cost.
This corresponds to all time-points, $k$, for which the conditional cost is
above the red line in the right-hand plot of Figure 4.

The intuition behind this pruning rule is that at any future time-point,
regardless of the future data, there would be a smaller cost of adding the
most recent change at the current time-point than having the most recent
change at $k$. To prove this formally, let the current step of the algorithm
be $t$, the value that we are considering pruning being $k$, and $T$ being any
future time-point. Then the conditional cost at time $T$ with a most recent
change-point at $k$ is

$Q_{k,\lambda}+C({X}_{1:n},k,T)+\lambda.$

Now our assumption on the cost means that $C({X}_{1:n},k,T)\geq
C({X}_{1:n},k,t)+C({X}_{1:n},t,T)$ as adding any change-point, in this case at
time $t$, can only reduce the unpenalised cost. Thus we have

$Q_{k,\lambda}+C({X}_{1:n},k,T)+\lambda\geq
Q_{k,\lambda}+C({X}_{1:n},k,t)+C({X}_{1:n},t,T)+\lambda.$

The conditional cost at time $T$ of a most recent change at $t$ is just

$Q_{t,\lambda}+C({X}_{1:n},t,T)+\lambda=\min_{l}\left\\{Q_{l,\lambda}+C({X}_{1:n},l,t)\right\\}+\lambda+C({X}_{1:n},t,T)+\lambda,$

where we have just plugged in the expresssion for $Q_{t,\lambda}$. This will
be smaller than the conditional cost at time $T$ of a change at $k$ if and
only if

$\min_{l}\left\\{Q_{l,\lambda}+C({X}_{1:n},l,t)\right\\}+\lambda\leq
Q_{k,\lambda}+C({X}_{1:n},k,t).$

The left-hand side is the minimum conditional cost at time $t$ plus $\lambda$.
Thus if the conditional cost at time $t$ of a most recent change at $k$ is
above this value, $k$ can never be the best choice of most recent change-point
at future time points, hence it can be pruned.

Killick et al. (2012) introduce an algorithm, PELT, that implements this
pruning rule. They show that it can give substantial improvements in
computational speed over Optimal Partitioning, particular if the data has
segments that are short relative to the length of data. In fact they show that
if we consider a data generating mechanism where, on average, the segment
lengths are constant as we increase $n$ (and hence the number of change-points
increases linearly with $n$) then PELT can have an expected computational
complexity that is linear in $n$.

The intuition behind this result comes from Figure 4, as we consider most
recent changes that are older than the actual most recent change-point, the
conditional cost tends to increase linearly. Thus the PELT pruning step will
often only keep candidate values for the most recent change-point that are
close to or more recent than the actual most recent change-point. Thus, on
average the computational cost of one step of the algorithm will be of the
order of the length of a segment rather than proportional to $t$.

##  4 Binary Segmentation and Other Hierarchical Methods

In the preceding sections, we described an approach to change-point detection
in which the estimator was formulated as the theoretical minimum of a
criterion. An algorithm was then needed to compute the minimum. By contrast,
the estimators described in this section will already be defined as algorithms
to begin with.

This section focuses on hierarchical methods for change-point detection, i.e.
methods for which the corresponding algorithms either:

  * •

first look at the data as a whole to identify change-point candidate(s)
considered to be the most prominent, and then gradually narrow their focus to
smaller and smaller sections of the data to propose less prominent change-
point candidates – we will refer to these types of techniques as top-down; or

  * •

start at the finest level of resolution of the data, where they remove from
consideration locations considered to be the least likely change-points, and
then gradually broaden their focus to bigger and bigger sections of the data
to propose ever more prominent change-point candidates – these types of
techniques will be referred to as bottom-up.

In either case, hierarchical methods gradually add ever less likely change-
point candidates (top-down methods), or remove ever more likely ones (bottom-
up methods), and therefore naturally arrange the change-point candidates in
their implied order of importance. This means that they tend to yield nested
solution paths to the multiple change-point problem, in the sense defined in
Section 3.3.

Model selection, i.e. deciding which set of change-point candidates along the
solution path should be preferred, will also be discussed for each method
separately.

###  4.1 Binary Segmentation

Arguably the simplest hierarchical, top-down method for multiple change-point
detection is binary segmentation. Before we describe it in simple algorithmic
terms, one necessary building block that needs introducing is a more general
definition of the CUSUM statistics compared to that introduced in Fearnhead
and Fryzlewicz (2022). The reason for the generalisation is that the CUSUM
statistic $C_{t}$ (and the signed CUSUM $C_{t}^{*}$) operated on the entire
data sequence, i.e. on the time domain $[1,n]$. In this section, we will need
CUSUMs that can be applied on arbitrary intervals $[s,e]\subseteq[1,n]$ (we
will use the term interval to denote an interval subset of the integers
$[1,n]$).

For any input sequence $\mathbf{X}=(X_{1},\ldots,X_{n})^{\prime}$ and $e>b\geq
s$ with $m=|[s,e]|=e-s+1$, we define the signed CUSUM statistic
$C^{*}_{s,b,e}(\mathbf{X})$, and the CUSUM statistic $C_{s,b,e}(\mathbf{X})$
for $\mathbf{X}$ on the interval $[s,e]$ by

$\displaystyle C^{*}_{s,b,e}(\mathbf{X})$ $\displaystyle=$
$\displaystyle\sqrt{\frac{e-b}{m(b-s+1)}}\sum_{t=s}^{b}X_{t}-\sqrt{\frac{b-s+1}{m(e-b)}}\sum_{t=b+1}^{e}X_{t}$
(32) $\displaystyle=$
$\displaystyle\sqrt{\frac{(b-s+1)(e-b)}{m}}(\bar{X}_{s:b}-\bar{X}_{(b+1):e}),$
$\displaystyle C_{s,b,e}(\mathbf{X})$ $\displaystyle=$
$\displaystyle|C^{*}_{s,b,e}(\mathbf{X})|.$ (33)

Remember that the CUSUM statistic is giving the evidence for a change at $b$,
in this case based just on the data $X_{s:e}$. This evidence is in terms of a
suitably re-scaled difference of the empirical mean of the data before and
after $b$. Larger values correspond to greater evidence for a change.

The binary segmentation algorithm is best defined recursively and hence
described by pseudocode. The main function is described as follows.

1:function Binary Segmentation($s$, $e$, $\zeta_{n}$)

2: if $e-s<1$ then

3: STOP

4: else

5: $b_{0}:=\arg\max_{b\in\\{s,\ldots,e-1\\}}C_{s,b,e}(\mathbf{X})$

6: if $C_{s,b_{0},e}(\mathbf{X})\geq\zeta_{n}$ then

7: add $b_{0}$ to the set of change-point candidates

8: Binary Segmentation($s$, $b_{0}$, $\zeta_{n}$)

9: Binary Segmentation($b_{0}+1$, $e$, $\zeta_{n}$)

10: else

11: STOP

12: end if

13: end if

14:end function

With the above definition in place, the standard binary segmentation procedure
is launched by the call Binary Segmentation($1$, $n$, $\zeta_{n}$), where
$\zeta_{n}$ is a threshold parameter. If $\zeta_{n}=0$, then the above
function acts in the “solution path” mode: no judgment is made as to whether
the successively proposed change-points are significant or not. For this mode
to be useful in change-point detection, a model selection procedure needs to
be applied at a later stage. The same principle applies in the case in which
$\zeta_{n}$ is positive but small – smaller than the level required for
consistent model selection. If $\zeta_{n}>0$ and has been set at an
appropriate level, the binary segmentation routine can combine solution path
generation and model selection. We discuss model selection for binary
segmentation and suitable choices of $\zeta_{n}$ later on.

The main idea of binary segmentation is very simple: at the start, we search
for the single most prominent change-point, by proceeding in the same way as
in the AMOC setting described in Fearnhead and Fryzlewicz (2022). If
$\zeta_{n}=0$, or if $C_{1,b_{0},n}(\mathbf{X})\geq\zeta_{n}$, we save the
location of this change-point candidate and proceed in the same way to the
left and to the right of it, thereby “zooming in” on smaller and smaller
sections of the data. On any given current interval $[s,e]$, we stop if the
interval is too short for us to proceed ($e-s<1$) or if we do not discover
significant change-points ($C_{s,b_{0},e}(\mathbf{X})<\zeta_{n}$). On
completion, we execute one of the following steps.

  * •

In the solution path mode ($\zeta_{n}=0$ or small), we sort the thus-obtained
set of change-point candidates in decreasing order of their corresponding
detection CUSUMs to create the solution path

$\mathcal{P}_{BS}=\left((b_{i}),C_{(s_{i}),(b_{i}),(e_{i})}(\mathbf{X})\right)_{i=1}^{n-1},$

where the ordering $(b_{i})$ is such that

$C_{(s_{1}),(b_{1}),(e_{1})}(\mathbf{X})\geq\ldots\geq
C_{(s_{n-1}),(b_{n-1}),(e_{n-1})}(\mathbf{X}).$

The solution path $\mathcal{P}_{BS}$ is then passed onto a model selection
routine (more details below).

  * •

In the model selection mode ($\zeta_{n}>0$ and set appropriately), we return
the obtained set of change-point candidates as the final set of change-point
location estimates; it is good practice to sort it in increasing order
beforehand.

For the solution path mode, any existing model selection tool can, in
principle, be used to choose the model based on the solution path
$\mathcal{P}_{BS}$. This is done by choosing from among the nested sequence of
models

$\displaystyle\mathcal{M}_{0}=\emptyset$
$\displaystyle\mathcal{M}_{1}=\\{(b_{1})\\}$ $\displaystyle\cdots$
$\displaystyle\mathcal{M}_{n-1}=\\{(b_{1}),\ldots,(b_{n-1})\\}$

One natural option would be to use SIC, i.e. choose the model that minimises

$\mbox{SIC}(N^{\prime},(b_{1}),\ldots,(b_{N^{\prime}}))=\frac{n}{2}\log\left\\{\frac{1}{n}S_{n}((b_{1}),\ldots,(b_{N^{\prime}}))\right\\}+N^{\prime}\log(n),$
(34)

over $N^{\prime}\in\\{0,\ldots,N_{U}\\}$, where $N_{U}$ is a pre-determined
upper bound for the number of change-points.

It is easy to think of binary segmentation as an appealing algorithm: to start
with, it is conceptually simple and easy to explain and to code – its
naturally recursive structure suggests recursion as an appropriate programming
device. Binary segmentation breaks up the difficult multivariate problem of
finding possibly multiple change-points into a series of one-dimensional (and
hence easier) problems: observe that the only task of the Binary Segmentation
routine is to find the single most prominent change-point on its current
domain of execution, before calling itself on the two implied children
intervals. (The property of binary segmentation whereby it fits the “best”
single-change-point model on each current interval $[s,e]$, i.e. the model
that maximises the CUSUM $C_{s,b,e}(\mathbf{X})$ over $b$, can be describes as
its greediness, a term used in algorithmic science to describe algorithms that
make locally optimal choices at each stage.) Another advantage of binary
segmentation is that it is relatively fast: a typical execution, in which
successive change-point candidates are found in the middle part of each
current interval (i.e. away from its edges) passes the entire dataset
$O(\log(n))$ times and hence takes $O(n\log(n))$ operations. Finally, an
appealing feature of binary segmentation is the generality of its principle –
it is possible to replace the CUSUM test with any other test appropriate for a
given type of change. However, care is needed as early iterations of binary
segmentation involved detecting a single change for intervals of data which
may contain multiple changes. As we saw, for example, in Section 2.4 of
Fearnhead and Fryzlewicz (2022) for detecting a change in slope, some methods
for detecting a single change can perform poorly when there are actually
multiple changes – for example, detecting the change in the wrong position
even for data where the signal for the true change locations is strong. A
future chapter of this book will return to this issue, and describe an
approach that overcomes this drawback of binary segmentation.

This latter issue shows that binary segmentation should only be used for
change-point detection methods that if applied to data without noise but with
multiple change-points, would detect one of these change-points. That is that
the maximum of the CUSUM statistic will be at a value $b$ that corresponds to
a true change-point. This can be shown to be the case for detecting changes in
mean with the CUSUM statistic. However even in this case, the fact the binary
segmentation can involve detecting a single change-point on intervals of data
which may contain multiple changes can affect its accuracy. We now illustrate
this.

Figure 5: Both panels: signal (red), its noisy realisation $\mathbf{X}$
(grey), corresponding CUSUM statistic $C_{1,b,200}(\mathbf{X})$ as a function
of $b$ (black).

To begin with, consider the left panel of Figure 5. Starting with the whole
data domain $[1,200]$, the CUSUM fitted by binary segmentation on $[1,200]$
clearly indicates $b=100$ as the first change-point candidate. This is
unsurprising as the two change-points at times $20$ and $180$ appear to be
secondary in importance to the change-point at time 100, so overall the signal
is, in a way, not far from a signal with just one change-point at 100; hence
fitting the single-change-point model to the entire data would not be “very
wrong” in this particular instance.

Let us contrast this with the situation in the right panel of Figure 5. Here,
visually, there does not appear to be one dominant change-point, and it is
unclear how to approximate the signal well by a signal with only a single
change-point. This, coupled with the fact that the noise level is relatively
large, causes the CUSUM statistic to miss all of the change-points: the
maximum of the CUSUM is attained around time $t=30$, far from any of the
actual change-points. In this instance, fitting the best single-change-point
model to the entire data is not a good idea, and provides an example of how
binary segmentation can be led onto the wrong track. As we will see in Section
4.2, the detection failure here is not because the problem may be inherently
difficult for any method to detect the change-points correctly (it is not the
case); the root of the failure of binary segmentation is that it considers the
whole data sample at once, which “averages away” the interesting features of
the data. As we will see in Section 4.2, in this example it is a better idea
to focus on smaller sections of the data with more prominent presence of the
features of interest, for example $X_{1:100}$ or $X_{81:120}$, which binary
segmentation is unable to do.

Authors such as Vostrikova (1981) and Venkatraman (1993) provide theoretically
interesting consistency results regarding binary segmentation, but we do not
reproduce them here: the most important messages that we want to convey are
that it is a particularly simple and fast algorithm, but that it can fail
unless the change-points have a natural ordering in terms of
prominence/importance, i.e. unless there is one dominant change-point on every
current interval under consideration. However, our practical experience with
problems of this type suggests that whether or not this holds for a given
dataset may be difficult to determine a priori. Methods described in the
following sections have been designed to mitigate, in one way or another, this
weakness of binary segmentation.

The fact that binary segmentation discovers change-points one by one has one
other important consequence: the best binary segmentation fit with
$N^{\prime}$ change-points (that is, the fit corresponding to the first
$N^{\prime}$ elements of the binary segmentation solution path) is not
guaranteed to be the same as the least-squares fit to the data of a piecewise-
constant function with $N^{\prime}$ change-points. However, the least-squares
fits to the data are not nested in the sense that the best $N^{\prime}+1$
change-point locations are not guaranteed to be a superset of the best
$N^{\prime}$ change-point locations, which may have computational speed
consequences when one wishes to compare least-squares model fits for a range
of proposed change-point numbers. By contrast, most algorithms that add
change-points one by one (such as binary segmentation) produce nested solution
paths which can by construction be traversed quickly when one wishes to
compare different model fits.

###  4.2 Wild Binary Segmentation

In the situation from the right-hand plot of Figure 5, rather than starting
with a CUSUM taken over the entire data sample, it would have been more
advisable to “zoom in” on the three-change-point feature of interest in the
middle part of the signal, by taking the CUSUM statistic over a shorter
interval containing this feature. Figure 6 shows the CUSUM statistic taken
over the interval $[66,135]$, its endpoints being much closer to the change-
points at 80, 100 and 120 than those of the original CUSUM over $[1,200]$. As
the section of the signal on the interval $[66,135]$ can be approximated
relatively well by a signal with a single change-point, it is not surprising
to see that the CUSUM computed over $[66,135]$ brings out the location of one
of the change-points, the middle one, very well.

Figure 6: Signal (red), its noisy realisation $\mathbf{X}$ (grey),
corresponding CUSUM statistic $C_{1,b,200}(\mathbf{X})$ as a function of $b$
(black) and the CUSUM statistic $C_{66,b,135}(\mathbf{X})$ as a function of
$b$ (blue).

If this observation were to be built into a change-point detection algorithm,
the obvious first issue to face is that the analyst (more often than not) has
little knowledge of where any features of interest may be located, so it is
unclear which parts of the signal to zoom in on. The main idea underlying Wild
Binary Segmentation (WBS; Fryzlewicz (2014)) is to take the CUSUM statistic
over each of a “representative sample” of intervals of varying lengths and
locations within $[1,n]$, in the hope that some of these intervals will align
well with the features of interest in the signal (as illustrated in Figure 6)
and therefore the CUSUMs computed on those will be informative for change-
point locations. At each stage of the procedure, WBS works with the largest
CUSUM (over the collection of all intervals drawn) and chooses the next
change-point candidate as the argument-maximum of this largest CUSUM.

To define WBS, we denote by $F_{n}^{M_{1}}$ a set of $M_{1}\geq M$ intervals
$[s_{m},e_{m}]$, $m=1,\ldots,M_{1}$, produced in either of the two following
ways.

  1. 1.

(Random grid.) Their start- and endpoints are drawn uniformly at random,
independently with replacement, from the set $\\{1,\ldots,n\\}$. In this
random setting we draw $M$ intervals and therefore $M_{1}=M$.

  2. 2.

(Deterministic grid.) Construct an equispaced grid
$\mathcal{K}=\\{i/(K-1)\,\,:\,\,i=0,\ldots,K-1\\}$ on $[0,1]$ such that the
number of gridpoints is the smallest integer $K$ for which $M_{1}=K(K-1)/2$,
the number of all intervals with start- and endpoints in $\mathcal{K}$, is
such that $M_{1}\geq M$. Transform the grid $K$ to an integer grid on $[1,n]$
through the transformation $i/(K-1)\to[(n-1)i/(K-1)+1]$, where $[\cdot]$ is
rounding to the nearest integer. Take $[s_{m},e_{m}]$ to be all intervals with
start- and endpoints on the thus-transformed grid.

Other ways of choosing the intervals – mixed random plus deterministic, random
but drawn from a distribution different than uniform, a different
deterministic scheme such as that of Seeded Binary Segmentation (Kovacs et
al., 2023) – are also acceptable as long as they ensure the property that,
with a high probability, each true change-point is well covered by at least
one interval which does not contain any other true change-points. This can
typically be achieved by (a) ensuring that the number of intervals is
sufficiently large and (b) that the start- and end-points of the intervals
cover the domain of the data sufficiently densely.

In either of the above cases (random or deterministic grid), we are guaranteed
at least $M$ intervals, where $M$ is decided by the user. Version 2.2 of the R
package `breakfast` (Anastasiou et al., 2022) uses $M=10000$ by default. An
appropriate choice of $M$ is key: on the whole, the larger the value of $M$,
the more accurate but slower the procedure. Theoretical lower bounds on $M$
are available, but, unsurprisingly, they depend on the minimum spacing between
neighbouring change-points, which the analyst typically does not have access
to. Naturally, if $n$ is so small that the $M$ requested by the user is larger
than or equal to the number of all intervals on $[1,n]$, then WBS simply uses
all those.

Using pseudocode, the main function of the WBS algorithm, with model selection
performed via thresholding with threshold $\zeta_{n}$ (as with binary
segmentation, if $\zeta_{n}=0$ or is positive but small, then the procedure
acts in the solution path mode), is defined in the following way.

1:function WBS($s$, $e$, $\zeta_{n}$)

2: if $e-s<1$ then

3: STOP

4: else

5: ${\mathcal{M}}_{s,e}:=$ set of those indices $m$ for which
$[s_{m},e_{m}]\in F_{n}^{M_{1}}$ is such that $[s_{m},e_{m}]\subseteq[s,e]$

6:
$(m_{0},b_{m_{0}}):=\arg\max_{m\in{\mathcal{M}}_{s,e},b\in\\{s_{m},\ldots,e_{m}-1\\}}C_{s_{m},b,e_{m}}(\mathbf{X})$

7: if $C_{s_{m_{0}},b_{m_{0}},e_{m_{0}}}(\mathbf{X})\geq\zeta_{n}$ then

8: add $b_{m_{0}}$ to the set of change-point candidates [with the additional
side information
$(m_{0},s_{m_{0}},e_{m_{0}},C_{s_{m_{0}},b_{m_{0}},e_{m_{0}}}(\mathbf{X}))$]

9: WBS($s$, $b_{m_{0}}$, $\zeta_{n}$)

10: WBS($b_{m_{0}}+1$, $e$, $\zeta_{n}$)

11: else

12: STOP

13: end if

14: end if

15:end function

The WBS procedure is launched by the call WBS($1$, $n$, $\zeta_{n}$).

The recursive calls in lines 9 and 10 of the WBS algorithm may create an
impression that the intervals on which WBS performs computation have to be re-
drawn, but this is not the case: all recursive calls use the same set of
intervals drawn prior to the launch of the WBS routine, and the same set of
CUSUM statistics computed (once, prior to the start of the recursive
mechanism) on these interevals. Therefore, WBS can equally be implemented non-
recursively as follows.

1:function WBS.NONREC($\zeta_{n}$)

2: while $F_{n}^{M_{1}}$ non-empty do

3: ${\mathcal{M}}:=$ set of those indices $m$ for which $[s_{m},e_{m}]\in
F_{n}^{M_{1}}$

4:
$(m_{0},b_{m_{0}}):=\arg\max_{m\in{\mathcal{M}},b\in\\{s_{m},\ldots,e_{m}-1\\}}C_{s_{m},b,e_{m}}(\mathbf{X})$

5: if $C_{s_{m_{0}},b_{m_{0}},e_{m_{0}}}(\mathbf{X})\geq\zeta_{n}$ then

6: add $b_{m_{0}}$ to the set of change-point candidates [with the additional
side information
$(m_{0},s_{m_{0}},e_{m_{0}},C_{s_{m_{0}},b_{m_{0}},e_{m_{0}}}(\mathbf{X}))$]

7: remove from $F_{n}^{M_{1}}$ those $m$ for which
$b_{m_{0}}\in\\{s_{m},\ldots,e_{m}-1\\}$

8: else

9: STOP

10: end if

11: end while

12:end function

Suppose we wish to use WBS to estimate the locations of significant change-
points, as opposed to merely producing a solution path of all possible models.
It is clear that we should use $\zeta_{n}>0$, but what values of $\zeta_{n}$
are appropriate? One natural requirement is that $\zeta_{n}$ is large enough
to control the probability of false alarms, in the sense that if the signal
$\mathbf{f}$ contains no change-points, the WBS estimator also returns no
change-points with a suitably high probability. Let us examine how high
$C_{s_{m_{0}},b_{m_{0}},e_{m_{0}}}(\mathbf{X})$ can possibly get (with a high
probability) given that $\mathbf{f}\equiv 0$ (the same bound will apply in the
more general case $\mathbf{f}\equiv$ const, as the CUSUM statistic ignores
additive constants). Trivially, it is clear that

$C_{s_{m_{0}},b_{m_{0}},e_{m_{0}}}(\mathbf{X})\leq\max_{1\leq s\leq b<e\leq
n}C_{s,b,e}(\mathbf{X}).$

We wish to choose $\zeta_{n}$ such that

$P\left(\max_{1\leq s\leq b<e\leq n}C_{s,b,e}(\mathbf{X})>\zeta_{n}\right)$

is small. Perhaps the simplest, albeit suboptimal, bound is provided by the
Bonferroni correction. We first observe that for a fixed triple $(s,b,e)$, the
quantity $C_{s,b,e}(\mathbf{X})$ is distributed as $|N(0,\sigma^{2})|$. This
is easy to establish directly: $C^{*}_{s,b,e}(\mathbf{X})$ is clearly normally
distributed as a linear combination of the data points, and considering its
expectation and variance yields the above fact. Using the Bonferroni
correction and Mill’s ratio inequality for the Gaussian distribution, we can
therefore bound

$\displaystyle P\left(\max_{1\leq s\leq b<e\leq
n}C_{s,b,e}(\mathbf{X})>\zeta_{n}\right)$ $\displaystyle\leq$
$\displaystyle\sum_{1\leq s\leq b<e\leq n}P(C_{s,b,e}(\mathbf{X})>\zeta_{n})$
$\displaystyle\leq$ $\displaystyle n^{3}2(1-\Phi(\zeta_{n}/\sigma))$
$\displaystyle\leq$ $\displaystyle
Cn^{3}\exp\\{-\zeta_{n}^{2}/(2\sigma^{2})\\},$

where $\Phi()$ is the standard normal cdf, and $C$ is a certain positive
constant. This upper bound converges to zero as $n\to\infty$ as long as
$\zeta_{n}=\sigma\sqrt{C_{0}\log\,n}$ with $C_{0}>6$. As the Bonferroni
correction is, in this case, very loose (because of the strong dependence
between $C_{s,b,e}(\mathbf{X})$ for neighbouring values of $(s,b,e)$), the
constant $C_{0}=6$ can be expected to be too high in practice. On the other
hand, the discussion in Section 3.4 implies that the $C_{0}$ in
$\zeta_{n}=\sigma\sqrt{C_{0}\log\,n}$ must exceed 2 for $P\left(\max_{1\leq
s\leq b<e\leq n}C_{s,b,e}(\mathbf{X})\right)$ to have any hope of tending to
0. To summarise, it is a good strategy to look for a suitable value of $C_{0}$
in the region $(2,6]$. In Fryzlewicz (2014), in which $\sigma$ is estimated
via the Median Absolute Deviation (MAD) estimator with
$\\{2^{-1/2}(X_{t+1}-X_{t})\\}_{t=1}^{T-1}$ on input (see Section 7 for more
on the MAD estimator and on estimating the noise variance in general), two
particular values of $C_{0}$ are considered and tested numerically on a
variety of signals: $C_{0}=2$ (it is concluded that this constant value is too
low in practice and tends to lead to the overestimation of the number of
change-points) and $C_{0}=3.38$ (which appears too high in practice and tends
to underestimate the number of change-points). A good universal choice of
$C_{0}$ is, at the time of writing, still an open question. For more details
on this problem, and partial solutions, we refer the reader to Fang et al.
(2020), and to Corollary 1 in Anastasiou and Fryzlewicz (2022).

Fryzlewicz (2014) also explores the use of WBS with an information criterion
referred to as the strengthened SIC. This is close to the standard SIC, and we
do not elaborate on its form here. Conceptually, WBS with an information
criterion as the model selector proceeds similarly to binary segmentation (see
Section 4.1): first, a solution path is executed via the call WBS($1$, $n$,
0). Next, the change-point candidates are sorted in decreasing order of their
corresponding detection CUSUMs to create the solution path

$\mathcal{P}_{BS}=\left((b_{i}),C_{(s_{i}),(b_{i}),(e_{i})}(\mathbf{X})\right)_{i=1}^{M_{2}},$

where $M_{2}\leq M$. Finally, SIC is applied on the resulting nested sequence
of model candidates in the same way as in formula (34). A different model
selection procedure for WBS, termed Steepest Drop to Low Levels (Fryzlewicz,
2020), directly processes the WBS solution path and is covered in Section 8.2.

In WBS, if $M$ is large enough, each true change-point is guaranteed (with
probability approaching one in $n$) to be contained within at least one of the
intervals $[s_{m},e_{m}]$ which contains only that single change-point, and to
be sufficiently distanced from either $s_{m}$ or $e_{m}$, i.e. to sit
sufficiently deeply within the interior of $[s_{m},e_{m}]$. This alone
provides suitable guarantees for the size of the largest maximum of the CUSUM
statistics, and is enough to yield theoretical and practical advantages of WBS
over binary segmentation, as similar guarantees are unavailable in the latter
method. However, WBS, as defined here, is still not rate-optimal in terms of
the permitted minimum distance between change-points. The fundamental reason
for this is the fact that the interval of the largest CUSUM, selected by WBS
at each stage, is not necessarily guaranteed to contain a single change-point
only, and therefore the CUSUM spanned on that interval is not necessarily the
best possible estimator for the change-point location. The Narrowest-Over-
Threshold (NOT) localisation approach described in Baranowski et al. (2019)
overcomes this issue by preferring the shortest intervals on which the
contrast statistics (such as CUSUMs) are deemed to be significant. However,
NOT applies to more general models than the multiple mean shift model
considered in this chapter and will therefore be covered in a future chapter.

In WBS, the collection $F_{n}^{M_{1}}$ is drawn only once at the start of the
procedure. In a signal with many change-points, this may potentially mean that
WBS can “run out of intervals”, if there are still undetected change-points
but no intervals left to operate on. This can be referred to as data-
nonadaptivity of WBS with respect to interval choice, and there is no
mechanism in WBS itself to remedy this.

A solution to this problem is proposed in Fryzlewicz (2020) by the method
termed “Wild Binary Segmentation 2” (WBS2), in which intervals are drawn
afresh on each interval on which the procedure currently operates (which can
be seen as a data-adaptive interval drawing scheme). This ensures that there
is always a supply of intervals to work with. While this can in principle lead
to an increase in the computational cost over WBS, Fryzlewicz (2020) notes
that each interval draw in WBS2 requires of the order of 100 or 1000
intervals, which is 100-10 times fewer than the corresponding default value in
WBS used in the R package `breakfast`.

Conceptually, WBS is a modular procedure in the sense that it is relatively
straightforward to see how it can be extended to the multiple mean-shift
problem in noise settings other than i.i.d. Gaussian: the CUSUM statistic
$C_{s,b,e}$ can be replaced, for example, by a suitable generalised likelihood
ratio test. This preserves the construction principle of WBS whereby it fits
simple single-change-point models on subsamples of the data, and combines the
results to create a multi-change-point fit on the entire sample. Plain binary
segmentation is equally easy to extend conceptually to more complex settings;
therefore, part of the appeal of the methods belonging to the binary
segmentation family is that they can potentially be applied to multiple change
detection e.g. in high-dimensional regression models or network data – topics
we will return to in a future chapter of this book.

###  4.3 From the Bottom Up: Tail-Greedy Unbalanced Haar

Both binary segmentation and WBS are top-down hierarchical methodologies,
successively subdividing the data by the next most likely change-point
candidate. An alternative direction may be to proceed from the bottom up, by
successively merging regions of the data considered the least likely to
contain change-points. This principle is formalised in the Tail-Greedy
Unbalanced Haar (TGUH) approach of Fryzlewicz (2018), which can be described
as a hierarchical agglomeration technique. As the generic notation for TGUH
can get a little involved, we introduce this algorithm on an example, rather
than in its generality.

Suppose that $n=6$, and our input data vector is $(X_{1},X_{2},\ldots,X_{6})$.
We initially look at all pairwise differences between consecutive data points,