Optimal Covariance Change Point Detection in High Dimension
Daren Wang1 , Yi Yu2 , and Alessandro Rinaldo1

arXiv:1712.09912v1 [math.ST] 28 Dec 2017

1

Department of Statistics and Data Science, Carnegie Mellon University
2
School of Mathematics, University of Bristol

Abstract
We study the problem of change point detection for covariance matrices in high dimensions.
We assume that {Xi }i=1,...,n is a sequence of independent, centered p-dimensional sub-Gaussian
random vectors is observed whose covariance matrices are piece-wise constant. Our task is to
recover with high accuracy the number and locations the of change points, which are unknown.
Our generic model setting allows for all the model parameters to change with n, including the
dimension p, the minimal spacing between consecutive change points, the magnitude of smallest
change and the maximal operator norm of the covariance matrices of the sample points. We
introduce two procedures, one based on the binary segmentation algorithm (e.g. Vostrikova,
1981) and the other on its extension known as wild binary segmentation of Fryzlewicz (2014),
and demonstrate that, under suitable conditions, both are able to consistently estimate the
number and locations of change points. Our second algorithm, called Wild Binary Segmentation
through Independent Projection (WBSIP) is shown to be minimax optimal in the in terms of
all the relevant parameters. Our minimax analysis also reveals a phase transition effect based
on our generic model setting. To the best of our knowledge, this type of results has not been
established elsewhere in the change point detection literature.
Keywords: Change point detection; High dimensional covariance estimation; High dimensional covariance testing; Binary segmentation; Wild binary segmentation; Random projection;
Minimax optimal.

1

Introduction

Change point detection in time series data has a long history and can be traced back to World War
II, during which the change point detection was specifically demonstrated in a sequential fashion
and heavily used in quality control. Wald (1945) was published shortly after the war and formally
introduced sequential probability ratio test (SPRT), which consists of simple null and alternative
hypotheses and is based on sequentially calculating the cumulative sum of log-likelihood ratios.
Page (1954) relaxes the statistic to the cumulative sum of weights, which are not restricted to loglikelihood ratios, namely cumulative sum (CUSUM). Since then, a tremendous amount of efforts
have been dedicated to this problem (see, e.g. James et al., 1987; Siegmund and Venkatraman,
1995) with applications ranging from clinical trials to education testing, signal processing and
cyber security.
In the last few decades, with the advancement of technology, an increasing demand from emerging application areas, e.g. finance, genetics, neuroscience, climatology, has fostered the development
of statistical theorties and methods for change point detection also in an offline fashion. Given a
1

time series data set, instead of assuming it is stationary over the whole time course, it is more robust
and realistic to assume that it is stationary only within time segments, and the underlying models
may change at certain time points. Formally, we assume that we observed n variables X1 . . . , Xn
such that,


F1 , i = 1, . . . , t1 − 1,



F , i = t , . . . , t − 1,
2
1
2
Xi ∼
(1)

.
.
.
,



F , i = t
K
K−1 , . . . , n,

where F1 , . . . , FK , K ≥ 2, are distribution functions, Fk−1 6= Fk , k = 2, . . . , K, and {t1 , . . . tK−1 } ⊂
{2, . . . , n − 1} are unknown change point locations. Testing the existence of change points and
estimating the locations of the change points are of primary interest.
The simplest and best studied scenario of model (1) is that {Xi = fi + εi }ni=1 is a univariate
time series, and {fi }ni=1 is piecewise constant signal. Change point detection here means detecting
changes in the mean of a univariate time series. An abundance of literature exists on this model
and variants thereof. See Section 1.1 below for a literature review.
The properties of the covariance are also of key importance in statistics, and detecting the
covariance changes in a time series data set also has a long history. See Section 1.1 below for
some literature review. In this paper we investigate the problem of change point detection and
localization for covariance matrices. We consider what is arguably the most basic setting for the
problem in (1): we observed a sequence (X1 . . . , Xn ) of independent and centered random vectors
in Rd wih respective covariance matrices (Σ1 , . . . , Σn ) with that maxi kΣi kop ≤ B. We assume that
Σi = Σi+1 for all coordinates i except for a few, which we refer to as the change points. We denote
with ∆ the minimal distance between two consecutive change points and with κ the magnitude
of the smallest change, measured by the operator norm of the difference between the covariance
matrices at two consecutive change points. The parameters p, ∆, B and κ completely characterize
the difficulty of change point localization problem, which, for a given n, is increasing in p and B
and decreasing in ∆ and κ. In fact, for reason that will become clear later on, it will be convenient
to aggregate the parameter B and κ into one parameter
κ
,
B2

which effectively plays the role of the signal to noise ratio. Our goal to identify all the change points
and to localize them, i.e. estimate their values accurately. In our analysis, we allow all the relevant
parameters p, ∆ and Bκ2 to change with n and derive upper and lower bounds for the localization
rate in terms of such parameters.
We propse two algorithms for covariance change point localization. The first one, called BSOP
(Binary Segmentation in Operator Norm), is based on an adaptation to the covariance setting of
the popular binary segmentaiton (BS) algorithm (e.g. Vostrikova, 1981) for change point detection
of a univariate signal. The BSOP algorithm, described below in Algorithm 1, is rather simple to
implement. Under appropriate assumptions we show in Theorem 1 that BSOP can consistently
estimate all the change point, with a localization rate that is sub-optimal, especially, when the
dimension p is allowed to be increasing with n. Our second algorithm, called WBSIP (Wilde Binary
Segmentation through Independent Projections) is significantly more refined, as it combines data
splitting, random projections and the wild binary segmentation (WBS) procedure of Fryzlewicz
(2014); see Algorithm 3. For all its complexity, WBSIP yields much sharper localization rates
2

than WBOP under a different, milder set of assumptions; see Theorem 2. We further obtain lower
bounds on the localization rate and demonstrate that WBSIP is in fact minimax rate optimal.
Interestingly, our analysis reveals
 a phase transition for the problem of covariance change point
localization: if ∆ = Θ B 4 pκ−2 then no estimator will be able to estimate the location of the
change point at a vanishing rate; see Lemma 3. On the other hand, if ∆ = Ω B 4 pκ−2 log n then
WBSIP will guarantee a localization rate of the order B 4 κ−2 log n, which, up to a log n term, is a
minimax optimal rate; see Lemma 4. It appears that this type of phase transition is unique to our
generic model setting. In fact, many authors (e.g. Aue et al., 2009; Avanesov and Buzun, 2016)
have already studied the consistency of change point estimation in covariate structure. To the best
of our knowledge, this phase transition based on the spacing parameter ∆ and the dimension p
has not been established elsewhere. Overall, our results provide a complete characterization of the
problem of change point localization in the covariance setting described above.

1.1

Relevant and Related Literature

For the classical problem of change point detection and localization of the mean of a univariate time
series, least squares estimation is a natural choice. For example, Yao and Au (1989) uses the least
squares estimators to show that the distance between the estimated change points and the truth
are within Op (1). Lavielle (1999) and Lavielle and Moulines (2000) consider penalized least-squares
estimation and prove consistency in the presence of strong mixing and long-range dependence of
the error terms εi ’s. Harchaoui and Lévy-Leduc (2010) consider the least squares criterion with a
total variation.
Besides least squares estimation based methods, other attempts have also been made. For
instance, Davis et al. (2006) propose a model based criterion on an autoregressive model. Wang
(1995) detects the change points using fast discrete wavelet transform. Harchaoui and Lvy-Leduc
(2010), Qian and Jia (2012), Rojas and Wahlberg (2014) and Lin et al. (2017) consider the fused
lasso (Tibshirani et al., 2005) to penalise the differences between the signals at two consecutive
points to detect the change points. Frick et al. (2014) introduce a multiscale statistic based procedure, namely simultaneous multiscale change point estimator (SMUCE), to detect change points
in piecewise constant signals with errors from exponential families. More recently, Li et al. (2017)
extends the multiscale method to a class of misspecified signal functions which are not necessarily
piecewise constants. The method of Davies and Kovac (2001) can also be used to determine the
location of change points in a piece-wise constant signal. Chan and Walther (2013) study the power
of the likelihood ratio test statistics to detect the presence of change points.
Among all the methods, binary segmentation (BS) (e.g. Vostrikova, 1981) is ‘arguably the most
widely used change-point search method’ (Killick et al., 2012). It goes through the whole time
course and searches for a change point. If a change point is detected, then the whole time course
is split into two, and the same procedure is conducted separately on the data sets before and
after the detected change point. The procedure is carried on until no change point is detected,
or the remaining time course consists of too few time points to continue the test. Venkatraman
(1992) proves the consistency results of the BS method in univariate time series mean change point
detection, with the number of change points allowed to increase with the number of time points.
It is worth to mention a variant of BS, namely wild binary segmentation (WBS), which is
proposed in Fryzlewicz (2014) and which can be viewed as a flexible moving window techniques, or
a hybrid of moving window and BS. Instead of starting with the whole data set and doing binary
segmentation, WBS randomly draws a collection of intervals under certain conditions, conducts
3

BS on each interval, and return the one which has the most extreme criterion value among all
the intervals. Compared to BS, under certain conditions, WBS is more preferable when multiple
change points are present. In the univariate time series mean change point detection problem,
Venkatraman (1992) shows that in order to achieve the estimating consistency using BS algorithm,
the minimum gap between two consecutive change points should be at least of order n1−β , where
n is the number of time points, and 0 ≤ β < 1/8; as claimed in Fryzlewicz (2014), by using WBS
algorithm, this rate can be reduced to log(n).
All the literature mentioned above tackles univariate time series models, however, in the big data
era, data sets are now routinely more complex and often appear to be multi- or high-dimensional,
i.e. Xi ∈ Rp , where p is allowed to grow with the number of data points n. Horváth and Hušková
(2012) propose a variant of the CUSUM statistic by summing up the square of the CUSUM statistic
in each coordinate. Cho and Fryzlewicz (2012) transform a univariate non-stationary time series
into multi-scale wavelet regime, and conduct BS at each scale in the wavelet context. Jirak (2015)
allows p to tend to infinity together with n, by taking maxima statistics across panels coordinatewise. Cho and Fryzlewicz (2015) propose sparsified binary segmentation method which aggregates
the CUSUM statistics across the panel by adding those which exceed a certain threshold. Cho
(2015) proposes the double CUSUM statistics which, at each time point, picks the coordinate which
maximizes the CUSUM statistic, and de facto transfers the high-dimensional data to a univariate
CUSUM statistics sequence. Aston and Kirch (2014) introduces the asymptotic concept of highdimensional efficiency which quantifies the detection power of different statistics in this setting.
Wang and Samworth (2016) study the problem of estimating the location of the change points of a
multivariate piecewise-constant vector-valued function under appropriate sparsity assumptions on
the number of changes.
As for change point detection in more general scenarios, the SPRT procedure (Wald, 1945) can
be easily used for the variance change point detection. Based on a generalized likelihood ratio
statistic, Baranowski et al. (2016) tackle a range of univariate time series change point scenarios,
including the variance change situations, although theoretical results are missing. Picard (1985)
proposes tests on the existence of change points in terms of spectrum and variance. Inclan and Tiao
(1994) develop an iterative cumulative sums of squares algorithm to detect the variance changes.
Gombay et al. (1996) propose some tests on detection of possible changes in the variance of independent observations and obtained the asymptotic properties under the non-change null hypothesis.
Berkes et al. (2009), among others, extend the tests and corresponding results to linear processes, as
well as ARCH and GARCH processes. Aue et al. (2009) considers the problem of variance change
point detection in a multivariate time series model, allowing the observations to have m-dependent
structures. Note that the consistency results in Aue et al. (2009) are in the asymptotic sense that
the number of time points diverges and the dimension of the time series remains fixed. Aue et al.
(2009) also require the existence of good estimators of the covariance and precision matrices, and
the conditions thereof are left implicit. Barigozzi et al. (2016) deal with a factor model, which
is potentially of high dimension p/n = O(log2 (n)), and use the wavelet transforms to make the
data possibly dependent across the timeline. Note that the model in Barigozzi et al. (2016) can be
viewed as a specific covariance change point problem, where the additional structural assumption
allows the dimensionality to go beyond the sample size.
As for the problem of hypothesis testing for high dimensional covariance matrices, which corresponds to the problem of change point detection, the literature is also abundant. In the cases where
p fixed and n → ∞, the likelihood ratio test (Anderson, 2003) has a χ2p(p+1)/2 limiting distribution
4

under H0 : Σ = I. When both n, p → ∞ and p/n → c ∈ (0, ∞), Johnstone (2001) extended Roy’s
largest root test (Roy, 1957) and derived the Tracy–Widom limit of its null distribution, to name
but a few. In the case where n, p → ∞ and p/n → ∞, Birke and Dette (2005) derived the asymptotic null distribution of the Ledoit–Wolf test (Ledoit and Wolf, 2002). More recently, Cai and Ma
(2013) studied the testing problem in the setting p → ∞ from
p a minimax point of view and derived
the testable region in terms of Frobenius norm is of order p/n.

Organization of the paper

The rest of the paper is organized as follows. In Section 2, we propose two different methods
for covariance change point detection problem. All three methods are shown to be consistent,
but under different conditions, and the wild binary segmentation with random projection has the
location error rate being log(n). In Section 3, we show the lower bound of the location error rate
in the covariance change point detection problem is log(n), which implies that the wild binary
segmentation with random projection is minimax optimal. Further discussion and future work
directions can be found in Section 4.
Notation
For a vector v ∈ Rp and matrix Σ ∈ Rp×p , kvk and kΣkop = maxkvk=1 |v ⊤ Σv| indicate the Euclidean
Rand the operator norm, respectively; for any integrable function f (·) : R 7→ R, denote kf k1 =
x∈R |f (x)|dx as the ℓ1 -norm of f (·). For functions f (·) and g(·) satisfying f (n)/g(n) → a, if
a ∈ (0, ∞), we denote by f (n) = O(g(n)).

2

Main Results

In this paper, we study the covariance change point detection in high dimension. To be specific, we
consider a centered and independent time series {Xi }ni=1 ⊂ Rp⊗n . Let {ηk }K
k=1 ⊂ {1, . . . , n} be the
collection of time points at which the covariance matrices of Xi ’s change. The model is formally
summarized as follows with Orlicz norm k · kψ2 , which is defined in Definition 2 in Appendix A.
Assumption 1. Let X1 , . . . , Xn ∈ Rp be independent sub-Gaussian random vectors such that
E(Xi ) = 0 , E(Xi Xi′ ) = Σi and kXi kψ2 ≤ B < ∞ for all i. Let {ηk }K+1
k=0 ⊂ {0, . . . , n} be a
collection of change points, such that η0 = 0 and ηK+1 = n and that
Σηk +1 = Σηk +2 = . . . = Σηk+1 , for any k = 0, . . . , K.
Assume the jump size κ = κ(n) and the spacing ∆ = ∆(n) satisfy that
inf

{ηk − ηk−1 } ≥ ∆ > 0,

k=1,...,K+1

and
kΣηk − Σηk−1 kop = κk ≥ κ > 0, for any k = 1, . . . , K + 1.
Remark 1. The parameter B. The assumption that maxi kXi kψ2 ≤ B is imposed in order to
control the order of magnitude of the random fluctuations of the CUSUM covariance statistics, given
below in Definition 1, around its means. See Lemma 5 in the Appendix for details. At the same
5

time, the chain of inequalities (2) shows also that maxi kΣi kop ≤ 2B 2 , so that the same assumptions
amounts also to a uniform upper bound on the operator norm of the covariance matrices of the date
points. In this regard, B may be reminiscent of the assumption that the signal be of bounded
magnitude sometimes used in the litarature on mean change point detection: see, e.g. Fryzlewicz
(2014, Assumption 3.1(ii)) and Venkatraman (1992, Condition (iii) on page 11). However, while
in the mean change point detection problem such boundedness assumption can be in fact removed
because the optimal solution is translation invariant (see, e.g., Wang and Samworth, 2016), this
is not the case in the present setting. Indeed, a constant shift in the largest eigenvalue of the
population covariance matrices also implies a change in the precision with which such matrices can
be estimated. Thus, it is helpful to think of B 2 as some form of variance term.
Remark 2. Relationship between κ and B. The parameters κ and B are not variation independent, as they satisfy the inequality κ ≤ B 2 /4. In fact,


K
n
n
n
κ ≤ max kΣηk − Σηk−1 kop ≤ 2 max kΣi kop = 2 max sup E (v ⊤ Xi )2 ≤ 4 max kXi k2ψ2 ≤ 4B 2 , (2)
i=1

k=1

i=1 v∈S p−1

i=1

where the second-to-last inequality follows from Equation (16) in Appendix A. For ease of readability, we will instead use the weaker bound κ ≤ B 2 throughout. In fact, in our analysis we will
quantify the combined effect of both κ and B with their ratio Bκ2 , which we will refer to as the
signal-to-noise ratio. Larger values of such quantities lead to better performance of our algorithm.
It is important to notice that the signal-to-noise ratio, and the task of change-point detection itself,
remains invariant with respect to any multiplictive rescaling of the data by an arbitsraty non-zero
constant.
In our model Assumption 1, all the relevant parameters p, ∆, K, B and κ are allowed to be
functions of the sample size n, although we do make this dependence explicit in our notation for
ease of readability. This generic setting allows us to study the covariance change point problem
with potential high dimensional data and with growing number of change points and decreasing
jump sizes.
Motivated by the univariate CUSUM statistic for mean change point detection, we define the
CUSUM statistic in the covariance context.
Definition 1 (Covariance CUSUM). For X1 , . . . , Xn ∈ Rp , a pair of integers (s, e) such that
0 ≤ s < e − 1 < n, and any t ∈ {s + 1, . . . , e − 1}, the covariance CUSUM statistic is defined as
s
s
t
e
X
X
e
−
t
t−s
s,e
⊤
e
Xi Xi −
Xi Xi⊤ .
St =
(e − s)(t − s)
(e − s)(e − t)
i=s+1

i=t+1

Its expected value is

e s,e
Σ
t

=

s

t
X
e−t
Σi −
(e − s)(t − s)
i=s+1

s

e
X
t−s
Σi .
(e − s)(e − t)
i=t+1

In the rest of this section, we propose two algorithms to detect the covariance change points.
Detailed algorithms and the main consistency results are presented in this section, with the proofs
provided in the Appendices. The advantages of each algorithm will be discussed later in the section.
6

2.1

Consistency of the BSOP algorithm

We begin our study by analyzing the performance of a direct adaptation of the binary segmentation
(BS) algorithm to the matrix setting based on the distance induced by the operator norm. The
resulting algorithm, which we call BSOP, is given in Algorithm 1. The BSOP procedure works
as follows: given any time interval (s, e), BSOP first computes the maximal operator norm of the
covariance CUSUM statistics over the time points in (s + ⌈p log n⌉, e − ⌊p log n⌋); if such maximal
value exceeds a predetermined threshold τ , then BSOP will identify the location b of the maximum
as a change point. The interval (s, e) is then split into two subintervals at b and the procedure is
then iterated separately on each of the resulting sub-intervals (s, b) and (b, e) until an appropriate
stopping condition is met.
The BSOP algorithm differs from the standard BS implementation in one aspect: the maximization of the of the operator norm of the CUSUM covariance operator is carried out only over
the time points in (s, e) that are away by at least p log n from the endpoints of the interval. Such restriction is needed to obtain adequate tail bounds for the operator norm of the covariance CUSUM
statistics S̃ts,e given in Definition 1 and of the centered and weighted empirical covariance matrices.
See Lemma 5 in Appendix A.
Algorithm 1 Binary Segmentation through Operator Norm. BSOP((s, e), τ )
INPUT: {Xi }ei=s+1 ⊂ Rp⊗(e−s), τ > 0.
Initial FLAG ← 0,
while e − s > 2p log(n) + 1 and FLAG = 0 do
a ← max⌈s+p log(n)⌉≤t≤⌊e−p log(n)⌋ Sets,e op
if a ≤ τ then
FLAG ← 1
else
b ← arg max⌈s+p log(n)⌉≤t≤⌊e−p log(n)⌋ Sets,e op
add b to the collection of estimated change points
BSOP((s, b − 1), τ )
BSOP((b, e), τ )
end if
end while
OUTPUT: The collection of estimated change points.

Remark 3 (Uniqueness of the maxima of the covariance CUSUM statistic.). The probability that
the maximiser of kSets,e kop over multiple time points is non unique is zero.

To analyze the performance of the BSOP algorithm we will impose the following assumption,
which is, for the most part, modeled after Assumption 3.2 in Fryzlewicz (2014).
Assumption 2. For a sufficiently large constant Cα > 0 and sufficient small constant cα , cκ > 0,
assume that ∆κB −2 ≥ Cα nΘ , p ≤ cα n8Θ−7 / log(n), where Θ ∈ (7/8, 1].
When the parameters κ and B are fixed, the above assumption requires ∆, the minimal spacing
between consecutive change point, to be of slightly smaller order than n, the size of the time series.
This is precisely Assumption 3.3 in Fryzlewicz (2014) (see also Cho and Fryzlewicz, 2015). The
7

fact that ∆ cannot be too small compared to n in order for the BS algorithm to exhibit good
performance is well known: see, e.g., Olshen et al. (2004). In Assumption 2, we require also the
8Θ−7
dimension p to be upper bounded by nlog n , which means that p is allowed to diverge as n → ∞.
Remark 4 (Generalizing Assumption 2). In Assumption 2 we impose certain constraints on the
scaling of the quantities B, κ, ∆ and p in relation to n that are captured by a single parameter Θ,
whose admissible values lye in (7/8, 1]. The strict lower bound of 7/8 on Θ is determined by the
calculations outlined below on (5) and (6), which are needed to ensure the existence of a non-empty
range of value for the input parameter τ to the BSOP algorithm. In fact, Assumption 2 may be
generalized by by allowing for different types of scaling in n of the signal-to-noise ratio κB −2 , the
minimal distance ∆ between consecutive change points and the dimension p. In detail, we may
require that κB −2  nΘ1 , ∆  nΘ2 and p log n  nΘ3 for a given triplet of parameters (Θ1 , Θ2 , Θ3 )
in an appropriate subset of [0, 1]3 . Such a generalization would then lead to consistency rates in
n that depend on all these parameters simultaneously. However, the range of allowable values of
(Θ1 , Θ2 , Θ3 ) is not a product set due to non-trivial constraints among them. We will refrain from
providing details and instead rely on the simpler formulation given in Assumption 2.
b

Theorem 1 (Consistency of BSOP). Under Assumptions 1 and 2, let B = {η̂k }K
k=1 be the collection
of the estimated change points from the BSOP((0, n), τ ) algorithm, where the parameter τ satisfies
p
√
(3)
B 2 p log(n) + 2 ǫn B 2 < τ < C1 κ∆(e − s)−1/2 ,
for some constant C1 ∈ (0, 1) and where

ǫn = C2 B 2 κ−1 n5/2 ∆−2
for some C2 > 0. Then,

b =K
P K

and

p

p log(n)


max |ηk − η̂k | ≤ ǫn ≥ 1 − n3 9p 2n−cp ,

k=1,...,K

(4)

for some absolute constant c > 0.
Remark 5. The condition (3) on the admissible values of the input parameter τ to the BSOP
algorithm is well defined. Indeed, by Assumption 2, for all pairs (s, e) such that e − s > 2p log n,
we have that
p
cα
κ∆n−1/2 ≤ (1/8)κ∆(e − s)−1/2
(5)
B 2 p log(n) ≤ B 2 cα n4Θ−7/2 ≤ B 2 cα nΘ n−1/2 ≤
Cα
and

√
1/2
1/2
1/2 5/4+Θ−7/4
2 ǫn B 2 = 2C2 B 3 κ−1/2 n5/4 ∆−1 (p log(n))1/4 ≤ (2C2 Cα−1 c1/4
n
α )Bκ
≤ (C1 /8)κ∆n−1/2 B −1 κ1/2 ≤ (C1 /8)κ∆(e − s)−1/2 ,

(6)

where in the chain of inequalities we have used Assumption 2 repeatedly. It is also worth noting
that the difference between the right hand side and the left hand side of (3) increase as Θ increases
to 1. Finally we remark that in the proof of Theorem 1, we actually let C1 = 1/8, but this is an
arbitrary choice and it essentially depends on the constants Cα and cα from Assumption 2.
8

Proof of Theorem 1. By induction, it suffices to consider any pair of integers s and e such that
(s, e) ⊂ (0, T ) and satisfying
ηr−1 ≤ s ≤ ηr ≤ . . . ≤ ηr+q ≤ e ≤ ηr+q+1 ,

q ≥ −1,

max{min{ηr − s, s − ηr−1 }, min{ηr+q+1 − e, e − ηr+q }} ≤ ǫn ,

where q = −1 indicates that there is no change point contained in (s, e). It follows that, for
sufficiently small cα and sufficiently large Cα ,
p
C2 B 2 κ−1 n5/2 p log(n)∆−2
ǫn
≤
∆/4
∆/4
1/2

≤ 4C2 B 2 κ−1 n5/2

cα n4Θ−7/2
Cα3 κ−3 B 6 n3Θ

−3
2 −4 Θ−1
≤ (4C2 c1/2
α Cα )(κ B )n

≤ (κ2 B −4 )nΘ−1
≤1

where the second inequality stems from Assumption 2, the third inequality holds by choosing
sufficiently small cα and sufficiently large Cα and the last inequality follows from the fact that
κ ≤ B 2 . Then, for any change point should be (s, e), it is either the case that
|ηp − s| ≤ ǫn ,
or that
|ηp − s| ≥ ∆ − ǫn ≥ ∆ − ∆/4 = 3∆/4.
Similar considerations apply to the other endpoint e. As a consequence, the fact that min{|ηp −
e|, |ηp − s|} ≤ ǫn implies that ηp is a detected change point found in the previous induction step,
while if min{ηp − s, ηp − e} ≥ 3∆/4 we can conclude that ηp ∈ (s, e) is an undetected change point.
In order to complete the induction step, it suffices to show that BSOP([s, e], τ ) (i) will not find
any new change point in the interval (s, e) if there is none, or if all the change points in (s, e) have
been already detected and (ii) will identify a location b such that |ηp − b| ≤ ǫn if there exists at
least one undetected
p change point in (s, e).
Set λ = B 2 p log(n). Then, the event A1 ({Xi }ni=1 , λ) defined in Equation (17) holds with
probbability at least 1−2×9p n3 ncp , for some universal constant c > 0. The proof will be completed
in two steps.
Step 1. First we will show that on the event A1 ({Xi }ni=1 , λ), BSOP((s, e), τ ) can consistently
detect or reject the existence of undetected change points within (s, e).
Suppose there exists ηp ∈ (s, e) such that min{ηp − s, ηp − e} ≥ 3∆/4. Set δ = p log n. Then
3
∆, since
δ ≤ 32
p log(n) ≤ cα n8Θ−7 ≤ cα nΘ ≤ cα Cα−1 ∆B −2 κ1 ≤ 3∆/32,
where the last inequality follows from Assumption 2. With this choice of δ, we apply Lemma 20 in
Appendix D (where we set c1 = 3/4) and obtain that
max

t=⌈s+δ⌉,...,⌊e−δ⌋

e s,e kop ≥ (3/8)κ∆(e − s)−1/2 .
kΣ
t
9

On the event A1 ({Xi }ni=1 , λ),
max

t=⌈s+δ⌉,...,⌊e−δ⌋

kSets,e kop ≥

max

t=⌈s+δ⌉,...,⌊e−δ⌋

e s,e kop − λ ≥ (3/8)κ∆(e − s)−1/2 − λ ≥ (1/8)κ∆(e − s)−1/2 ,
kΣ
t
(7)

where the last inequality follows from (5) (in the last step we have set C1 = 1/8). If (3) holds, then,
on the event A1 ({Xi }ni=1 , λ), BSOP((s, e), τ ) detects the existence of undetected change points if
there are any.
Next, suppose there does not exist any undetected change point within (s, e). Then, one of the
following cases must occur.
(a) There is no change point within (s, e);
(b) there exists only one change point ηr within (s, e) and min{ηr − s, e − ηr } ≤ ǫn ;
(c) there exist two change points ηr , ηr+1 within (s, e) and that max{ηr − s, e − ηr+1 } ≤ ǫn .
Observe that if case (a) holds then, on the event A1 ({Xi }ni=1 , λ), we have that
max

t=⌈s+δ⌉,...,⌊e−δ⌋

kSets,e kop ≤

max

t=⌈s+δ⌉,...,⌊e−δ⌋

e s,e kop + λ = λ < τ,
kΣ
t

where the last inequality follows from (3). If situation (c) holds, then on, the event A1 ({Xi }ni=1 , λ),
we have
max

t=⌈s+δ⌉,...,⌊e−δ⌋

kSets,e kop ≤

max

t=⌈s+δ⌉,...,⌊e−δ⌋

e s,e kop +λ ≤ max{kΣ
e s,e kop , kΣ
e s,e kop }+λ ≤ 2√ǫn B 2 +λ,
kΣ
ηr
ηr+1
t

where the first inequality follows from A1 ({Xi }ni=1 , λ), the second inequality from Lemma 19 and
the third inequality from Lemma 23. (Both Lemmas are in Appendix D.2.) Case (b) can be handled
in a similar manner. Thus, if (3) holds, then on the event A1 ({Xi }ni=1 , λ), BSOP((s, e), τ ) has no
false positives when there are no undetected change points in (s, e).
Step 2. Assume now that there exists a change point ηp ∈ (s, e) such that min{ηp −s, ηp −e} ≥ 3∆/4
and let
b ∈ arg
max
Sets,e op .
t=⌈s+δ⌉,...,⌊e−δ⌋

To complete the proof it suffices to show that |b − ηk | ≤ ǫn .
Let v be such that
v ∈ arg max |u⊤ Sebs,e u|.
kuk=1

Consider the univariate time series {Yi (v)}ni=1 and {fi (v)}ni=1 defined in (53) and (54) in Appendix D.2. By Lemma 21 there, b ∈ arg maxs≤t≤e |Yet (v)|. Next, we wish to apply Corollary 10
to the time series {Yi (v)}ei=s and {fi (v)}ei=s . Towards that end, we first need to ensure that the
conditions required for that result to hold are verified. (Notice that in the statement of Corollary 10, the fi ’s are assumed to be uniformly bounded by B1 , while in this proof the fi (v)’s defined
in (54) are assumed to be bounded by 2B 2 .) First, the collection of the change points of the
time series {fi (v)}ei=s+1 is a subset of {ηk }K+1
k=0 ∩ (s, e). The condition (26) and the inequality
10

√
2 δB 2 ≤ (3c1 /4)κ∆(e − s)−1/2 are straightforward consequences of Assumption 2, while (34)
follows from the fact that
e s,e
|fets,e (v) − Yets,e (v)| ≤ Sets,e − Σ
t

op

≤ λ.

Similarly, (33) stems from the relationships
max

t=⌈s+δ⌉,...,⌊e−δ⌋

|Yets,e (v)| =

max

t=⌈s+δ⌉,...,⌊e−δ⌋

kSets,e kop ≥

max

t=⌈s+δ⌉,...,⌊e−δ⌋

e s,e kop −λ ≥ (1/8)κ∆(e−s)−1/2
kΣ
t

where the first inequality holds under the event A1 ({Xi }ni=1 , λ) and the second inequality is due to
(7) and Assumption 2. Thus, all the assumptions of Corollary 10 are met. An application of that
result yields that there exists ηk , a change point of {fi (v)}ei=s satisfying (29), such that
|b − ηk | ≤ C2 λ(e − s)5/2 ∆−2 κ−1 ≤ ǫn .
The proof is complete by observing that (29) implies min{ηk − s, ηk − e} ≥ 3∆/4, as discussed
in the argument before Step 1.
The proof of the theorem relies on a non-trivial extension of the arguments for proving consistency of the BS algorithm in one-dimensional mean change-point detection problems, as done
in Venkatraman (1992) (see also Fryzlewicz, 2014). The main difficulty that prevents a direct
e s,e
application of those results is the fact the regions of monotonicity of the function t 7→ Σ
t
op
are hard to derive. Instead, for each par of integers 1 ≤ s < e ≤ n with e − s > 2p log n, we
study the one-dimensinal time series {(v ⊤ Xi )2 }i=1,...,n of the squared coefficients of the projeciton
of the data along a one-dimensional linear subspace spanned by a distinguished unit vector v,
e s,e , where
which we term the shadow vector. This is simply the leading singular vector of Σ
b
op
e s,e . As it turns out, with such a choice of the shadow vector,
b = arg maxt∈(s+p log n,e−p log n) Σ
t
op
the local maxima of CUSUM statistic applied to the corresponding one-dimensional time series
coincide with the local maxima of the time series of the values of the operator norm of the CUSUM
covarance statistics. As a result, for the purpose of detecting local maxima of the CUSUM covariance statistic, it is enough and in fact much simpler to study the univariate times series of the
squared projections onto the appropriate ghost vector. Of course, at each iteration of the BSOP
algorithm a new ghost vector and a new univarite time series is obtained and a new local maximum
is found. See Appendix D.2 for further comments on the uses and intepretation of the shadow
vector.
Theorem 1 implies, that with high probability, the BSOP algorithm will identify all the change
points and estimate their locations with an error that is bounded by
ǫn 

B 2 −2 5/2 p
∆ n
p log n.
κ

Noitce that, as expected, the performance of BSOP is decreasing in the inverse of the signal-tonoise ratio parmeter Bκ2 , the inverse of the minima distance ∆ between change points and the
dimension p. The above bound yields a family of rates of consistency for BSOP, depending on the
scaling of each of the quantities involved in it. For example, in the simplest and most favorable

11

scenario whereby B, κ and the dimension p are constants, the bound implies a rate for change point
localization of the order
p
ǫn
 n−2Θ+3/2 log n,
n
which is decreasing in the Θ ∈ (7/8, 1]. In particular, when the number
q of change points is also

kept constant, we have that Θ = 1, yielding a localization rate of order logn n .
As we will see in the next section, the dependence on the parameters B, κ and ∆ is suboptimal. The advantage of BSOP over the rate-optimal algorithm we introduce next is that BSOP
only requires one input parameter, the threshold value τ . Furthermore, when the spacing parameter
∆ is comparable with n and the dimension p of the data grows slowly with respect with n, then
BSOP can still deliver good consistency rates. Therefore, despite its suboptimality in general,
BSOP is a simple and convenient algorithm which may serve as a competitive benchmark for other
proceudres.

2.2

Consistency of the WBSIP algorithm

In this section we describe and analyze the performance of a new algorithm for covariance changepoint detection, which we term WBSIP for Wild Binary Segmentation through Independent Projections. The WBSIP algorithm is a generalization of the wild binary segmentation or WBS procedure
of Fryzlewicz (2014) for mean change-point detection and further exploits the properties shadow
vectors. The WBSIP procedure begins by splitting the data into halves and by selecting at random
a collection of M pairs of integers (s, e) such that 1 ≤ s < e ≤ n and e−s > p log n+1. In its second
step, WBSIP computes, for each of the M random integer intervals previously generated, a shadow
vector using one half od the data and its corresponding one-dimensional time series using the other
half. The final step of the procedure is to apply the WBS algorithm over the resulting univariate
time series. The details of the algorithm are given in Algorithm 2, which describes the computation
of the shadow vectors by principal component methods, and Algorithm 3, which implements WBS
to the resulting one dimensional time series.
Remark 6. The idea of combining the WBS algorithm with sample splitting is due to Wang and Samworth
(2016), who applied it to the problem of mean change-point detection in multivariate settings. Like
us, Wang and Samworth (2016) employ one half of the sample to compute an optimal direction
and then project the second half of the data along such direction to obtain a univariate time series
which can then be handled using WBS. One of the main differences between our approach and that
of Wang and Samworth (2016) is the fact that we do no require the shadow vectors to be consistent estimators of subspaces related the true covariance matrices. In particular, our analysis holds
without any eigengap assumptions.
In order to analyze the performance of the WBSIP procedure, we will impose the following
assumption, which is significantly weaker then Assumption 2.
Assumption 3. There exists a sufficiently large absolute constant C > 0 such that ∆κ2 ≥
Cp log(n)B 4 .
Remark 7. We recall that all the parameters ∆, κ, p, B are allowed to depend on n. Since κ ≤ B 2 ,
and assuming without loss of generality that the constant C in the previous assumption is larger
than 8, we further have that
p log(n) ≤ ∆κ2 B −4 C −1 ≤ ∆/8,
12

√
which is used repeatedly below. In fact, in the proof we will set C = 32 2. This choice is of course
arbitrary and is only made for convenience in carrying out the calculations below. We also recall
that the ratio B 4 κ−2 is invariant of multiplicative scaling, i.e. if Xi′ = αXi for all i, then the
corresponding ratio B 4 κ−2 stays the same.
Algorithm 2 Principal Component Estimation PC({Xi }ni=1 , {(αm , βm )}M
m=1 )

INPUT: {Xi }ni=1 , {(αm , βm )}M
m=1
for m = 1, . . . , M do
if βm − αm > 2p log(n) + 1 then
dm ← arg max⌈αm +p log(n)⌉≤t≤⌊βm −p log(n)⌋ kSetαm ,βm kop
um ← arg maxkvk2 =1 v ⊤ Sedαmm ,βm v|
else
um ← 0
end if
end for
OUTPUT: {um }M
m=1 .

Algorithm 3 Wild Binary Segmentation through Independent Projection.
{(αm , βm )}M
m=1 , τ, δ)
INPUT: Two independent samples {Wi }ni=1 , {Xi }ni=1 , τ , δ.
n
M
{um }M
m=1 = P C({Wi }i=1 , {(αm , βm )}m=1 )
for i ∈ {s, . . . , e} do
for m = 1, . . . , M do
2
Yi (um ) ← u⊤
m Xi
end for
end for
for m = 1, . . . , M do
(s′m , e′m ) ← [s, e] ∩ [αm , βm ] and (sm , em ) ← (⌈s′m + δ⌉, ⌊e′m − δ⌋)
if em − sm ≥ 2 log(n) + 1 then
bm ← arg maxsm +log(n)≤t≤em −log(n) |Yetsm ,em (um )|
am ← Yebsmm ,em (um )
else
am ← −1
end if
end for
m∗ ← arg maxm=1,...,M am
if am∗ > τ then
add bm∗ to the set of estimated change points
WBSRP((s, bm∗ ), {(αm , βm )}M
m=1 , τ )
WBSRP ((bm∗ + 1, e), {(αm , βm )}M
m=1 , τ )
end if
OUTPUT: The set of estimated change points.

13

WBSIP((s, e),

Similarly to the BSOP algoithm, WBSIP also applies a slight modification to the WBS algorithm
as originally proposed in Fryzlewicz (2014). When computing the shadow vectors in Algorithm 2,
the search for the optimal direction onto which to project the data is restricted, for any given
candidate interval, only to the time points that are at least p log n away from the endpoints of the
interval. As remarked in the previous section, this ensures good tail bounds for the operator norms
of the matrices involved.
A second, more substantial, adaptation of WBS is used in Algorithm 3: when searching for
candidate change points inside a given interval, the algorithm only consider time points that are
δ-away from the endpoints of the interval, where δ is an upper bound on the localization error –
the term ǫn in Theorem 2 below. The reasons for such restriction are somewhat subtle: once an
estimated change point is found near a true change point, in its next iteration the algorithm can
no longer look for change points in that proximity, since this will result, with high probability,
in spurious detections. This phenomenon is due to the fact that the behavior of the CUSUM
statistic of the projected data is not uniform around its local maxima. Therefore, after a true
change point has been detected, the algorithm must scan only nearby regions of low signal-to-noise
ratio – so that the probability of false positives can be adequately controlled with a proper choice
of the thresholding parameter τ . Thus the need to stay away by at leat the localization error
form the detected chaing points. A very similar condition is imposed in the main algorithm of
Wang and Samworth (2016) and implicitly in Korkas and Fryzlewicz (2017). The value of δ is left
as an input paramerer (as in Wang and Samworth, 2016), but any value between the localization
error ǫn given in the statement of Theorem 2 and the minimal distance ∆ between change points
will do. In the proof of we Theorem 2 we set δ ≤ 3∆/32.
M
Theorem 2 (Consistency of WBSIP). Let Assumptions 1 and 3 hold let {(αm , βm )}M
m=1 ⊂ (0, T )
be a collection of intervals whose end points are drawn independently and uniformly from {1, . . . , T }
and such that max1≤m≤M (βm − αm ) ≤ C∆ for an absolute constant C > 0. Set

ǫn = C1 B 4 log(n)κ−2 ,
for a C1 > 0. Suppose there exist c2 , c3 > 0, sufficiently small, such that the input parameters τ
and δ satisfy
p
√
B 2 log(n) < τ < c2 κ ∆,
(8)
ǫn < δ ≤ c3 ∆.

b

Then the collection of the estimated change points B = {η̂k }K
k=1 returned by WBSIP with input
,
τ
and
δ
satisfies
parameters of (0, n), {(αm , βm )}M
m=1
n
o
b = K;
P K
max |ηk − η̂k | ≤ ǫn
k=1,...,K

2
−c
≥1 − 2n M n − n3 9p 2n−cp − exp log(n/∆) − M ∆2 /(16n2 )

for some absolute constants c > 0.

Remark 8 (The relationships among the constants in Theorem 2). The choice of the
constant C is essentially arbitrary but will affect the choice of the constants C1 , c2 and c3 , where
c2 and c3 in particular have to be picked small enough. This dependence can be tracked in the proof
but we refrain from giving further details.
14

The above theorem implies that the WBSIP algorithm can estimate the change points perfectly
well, with high probability, with a localization rate upper bounded by
B 4 log n
ǫn
 2
.
n
κ n
The fact that the dimension p does not appear explicitly in the localization rate is an interesting, if
not perhaps surprising, finding. Of course, the dimension does affect (negatively) the performance
of the algorithm through Assumption 3: keeping n and ∆ fixed, a larger value of p implies a larger
4
value of Bκ2 in order for that assumption to hold. In turn, this leads to a larger bound in Theorem 2.
Furthermore, the dimension p appears in the probability of the event that WBSIP fails to localize
all the change points. We remark that, for the different problem of high-dimensional mean change
point detection, Wang and Samworth (2016) also obtained a localization rate independent of the
dimension: see Theorem 3 there. In Section 3 below we will prove that Assumption 3 is in fact
essentially necessary for any algorithm to produce a vanishing localization rate.
Proof. Since ǫn is the desired order of localization rate, by induction, it suffices to consider any
generic (s, e) ⊂ (0, T ) that satisfies
ηr−1 ≤ s ≤ ηr ≤ . . . ≤ ηr+q ≤ e ≤ ηr+q+1 ,

q ≥ −1,

max{min{ηr − s, s − ηr−1 }, min{ηr+q+1 − e, e − ηr+q }} ≤ ǫn ,

where q = −1 indicates that there is no change point contained in (s, e).
Note that under Assumption 3, ǫn ≤ ∆/4; it, therefore, has to be the case that for any change
point ηp ∈ (0, T ), either |ηp − s| ≤ ǫn or |ηp − s| ≥ ∆ − ǫn ≥ 3∆/4. This means that min{|ηp −
e|, |ηp − s|} ≤ ǫn indicates that ηp is a detected change point in the previous induction step, even
if ηp ∈ (s, e). We refer to ηp ∈ [s, e] as an undetected change point if min{ηp − s, ηp − e} ≥ 3∆/4.
In order to complete the induction step, it suffices to show that WBSIP((s, e), {(αm , βm )}M
m=1 , τ, δ)
(i) will not detect any new change point in (s, e) if all the change points in that interval have been
previous detected, and (ii) will find a point b in (s, e) (in fact, in (s + δ, e− δ)) such that |ηp − b| ≤ ǫn
if there exists at least one undetected change point in (s, e).
Let
n
M
{um }M
m=1 = P C({Xi }i=1 , {(αm , βm )}m=1 ).

n
n
Since the intervals {(αm , βm )}M
m=1 are generated independently from {Xi }i=1 ∪ {Wi }i=1 , the rest
of the argument is made under M defined in Equation (20) of Appendix A, which has no effects
on the distribution of {Xi }ni=1 ∪ {Wi }ni=1 .
p
Step 1. Let λ1 = B 2 p log(n). In this step, we are to show that, on the event A1 ({Wi }ni=1 , λ1 )
and for some c′1 > 0,
′
′
sup |u⊤
m (Σηk − Σηk−1 )um | ≥ c1 kΣηk − Σηk−1 kop = c1 κk

for every

k = 1, . . . , K + 1

(9)

1≤m≤M

On the event M, for any ηk ∈ (0, n), without loss of generality, there exists αk ∈ [ηk −3∆/4, ηk −∆/2]
and βk ∈ [ηk + ∆/2, ηk + 3∆/4]. Thus [αk , βk ] contains only one change point ηk . Using Lemma 20
in Appendix D and the inequality p log(n) ≤ ∆/8, we have that
√
e ηαk ,βk kop ≥ (1/2)kΣη − Ση kop ∆.
e αk ,βk kop = kΣ
max
kΣ
(10)
t
k
k−1
k
t=⌈αk +δ⌉,...,⌊βk −δ⌋

15

Let bk ∈ arg maxt=⌈αk +δ⌉,...,⌊βk −δ⌋ kSetαk ,βk kop , where Sets,e denote the covariance CUSUM statistics
of {Wi }ei=s+1 at evaluated t. Since kΣηk − Σηk−1 kop = κk , by definition,
e αk ,βk uk | ≥ |u⊤ Seαk ,βk uk | − λ1
|u⊤
k Σbk
k bk
=

≥

max

t=⌈αk +δ⌉,...,⌊βk −δ⌋

max

e αk ,βk kop − 2λ1
kΣ
t
√
− Σηk−1 kop ∆ − 2λ1
√
− Σηk−1 kop ∆

t=⌈αk +δ⌉,...,⌊βk −δ⌋

≥ (1/2)kΣηk

≥ (1/4)kΣηk

kSetαk ,βk kop − λ1

where the first and second inequalities hold on the event A1 ({Wi }ni=1 , λ1 ), the third inequality
follows from (10) and the last inequality from Assumption 3. Next, observe that
q
t−αk

t ≤ ηk ,
(βk −αk )(βk −t) (βk − ηk )(Σηk − Σηk−1 ),
αk ,βk
e
q
=
Σt
βk −t

(βk −αk )(t−αk ) (ηk − αk )(Σηk − Σηk−1 ), t ≥ ηk .
Using the above expression, for bk ≥ ηk , we have that
√
e αk ,βk uk
(1/4)kΣηk − Σηk−1 kop ∆ ≤ u⊤
k Σbk
s
βk − bk
(ηk − αk ) u⊤
=
k (Σηk − Σηk−1 )uk
(βk − αk )(bk − αk )
s
(βk − ηk )(ηk − αk ) ⊤
uk (Σηk − Σηk−1 )uk
≤
βk − αk
√
≤ 2∆ u⊤
k (Σηk − Σηk−1 )uk .
√
Therefore (9) holds with c′1 = 1/(2 2). The case of bk < ηk follows from very similar calculations.

Step 2. In this step, we will show that WBSIP((s, e), {(αm , βm )}M
m=1 , τ, δ) will consistently detect
or reject the existence of undetected change points within
(s,
e),
provided that (9) holds and on
p
2 log(n), and M, given in Equation (18)
,
λ
),
where
λ
=
B
the two events B1 ({Xi }ni=1 , {um }M
2
m=1 2
and Equation (20) in Appendix A, respectively.
Let am , bm and m∗ be defined as in WBSIP((s, e), {(αm , βm )}M
m=1 , τ, δ). Denote Yi (um ) =
s,e
s,e
⊤
2
⊤
e
e
(um Xi ) and fi (um ) = um Σi um . Let Yt (um ) and ft (um ) be defined as in (53) and (54) of
Appendix D.2 respectively.
Suppose there exists a change point ηp ∈ (s, e) such that min{ηp − s, e − ηp } ≥ 3∆/4. Let
δ ≤ 3∆/32. Then, on the event M, there exists an interval (αm , βm ) selected by WBSIP such that
αm ∈ [ηp − 3∆/4, ηp − ∆/2] and βm ∈ [ηp + ∆/2, ηp + 3∆/4].
Then [s′m , e′m ] = [αm , βm ] ∩ [s, e], and [sm , em ] = [s′m + δ, e′m − δ] (see details of the WBSIP
procedure in Algorithm 3). Moreover, we have that min{ηp −sm , em −ηp } ≥ (1/2)∆. Thus, [sm , em ]
contains at most one change point of the time series {fi (um )}ni=1 . A similar calculation as the one
shown in the proof of Lemma 20 gives that
√
max
|fetsm ,em (um )| ≥ (1/8) ∆|u⊤
m (Σηp − Σηp−1 )um |,
⌈sm +log(n)⌉≤t≤⌊em −log(n)⌋

16

where em − sm ≤ (3/2)∆ is used in the last inequality. Therefore
am =
≥

max

⌈sm +log(n)⌉≤t≤⌊em −log(n)⌋

max

⌈sm +log(n)⌉≤t≤⌊em −log(n)⌋

|Yetsm ,em (um )|

|fetsm ,em (um )| − λ2

√
≥ (1/8) ∆|u′m (Σηp − Σηp−1 )um | − λ2 ,

where the first inequality holds on the event B1 ({Xi }ni=1 , {um }M
m=1 , λ2 ). Thus for any undetected
change point ηp within (s, e), it holds that
am∗ =
≥

sup am

(11)

√
√
sup (1/8) ∆|u′m (Σp − Σp−1 )um | − λ2 ≥ (c′1 /8)κp ∆ − λ2

(12)

1≤m≤M
1≤m≤M

√
≥ (c′1 /16)κp ∆

(13)

where the second inequality follows from (9), and the last inequality from
p
√
λ2 = B 2 log(n) ≤ (c′1 /16)κ ∆,
√
by choosing the constant C in Assumption 3 to be at least 32.
Then, WBSIP((s, e), {(αm , βm )}M
m=1 , τ, δ) correctly accepts the existence of undetected change
points on the events (9), B1 ({Xi }ni=1 , {um }M
m=1 , λ2 , δ) and M.
Suppose there does not exists any undetected change points within (s, e), then for any (s′m , e′m ) =
(αm , βm ) ∩ (s, e), one of the following situations must hold.
(a) There is no change point within (s′m , e′m );
(b) there exists only one change point ηr within (sm , em ) and min{ηr − s′m , e′m − ηr } ≤ ǫn ; or
(c) there exist two change points ηr , ηr+1 within (sm , em ) and max{ηr − s′m , e′m − ηr+1 } ≤ ǫn .
Observe that if (a) holds, then, on the event B1 ({Xi }ni=1 , {um }M
m=1 , λ2 ) given in Equation (18)
′
′
in Appendix A, for (sm , em ) = (sm + δ, em − δ), we have
max

⌈sm +log(n)⌉≤t≤⌊em −log(n)⌋

|Yetsm ,em (um )| ≤

max

⌈sm +log(n)⌉≤t≤⌊em −log(n)⌋

|fetsm ,em (um )| + λ2 = 0 + λ2 .

If (b) or (c) holds, then since (sm , em ) = (s′m + δ, e′m − δ) and δ ≤ ǫn , it must be the case that
(sm , em ) does not contain any change points. This reduces to case (a). Therefore if (8) holds, then
WBSIP((s, e), {(αm , βm )}M
m=1 , τ, δ) will always correctly reject the existence of undetected change
points, on the event B1 ({Xi }ni=1 , {um }M
m=1 , λ2 ).
Step 3. Assume that there exists a change point ηp ∈ (s, e) such that min{ηp − s, ηp − e} ≥ 3∆/4.
Let am , bm and m∗ be defined as in WBSIP((s, e), {(αm , βm )}M
m=1 , τ ).
To complete the proof it suffices to show that, on the events B1 ({Xi }ni=1 , {um }M
m=1 , λ2 ) and
n
M
B2 ({Xi }i=1 , {um }m=1 , λ2 ) given in Equation (18) and Equation (19) respectively of Appendix A,
there exists a change point ηk ∈ [sm∗ , em∗ ] such that min{ηk −s, ηk −e} ≥ 3∆/4 and |bm∗ −ηk | ≤ ǫn .
Consider the univariate time series {Yi (um∗ )}ni=1 and {fi (um∗ )}ni=1 defined in (53) and (54) of
m∗
is a
Appendix D.2. Since the collection of the change points of the time series {fi (um∗ )}ei=s
m∗
17

em∗
subset of that of {ηk }K+1
k=0 ∩ [s, e], we may apply Corollary 13 to the time series {Yi (um∗ )}i=sm∗
em∗
and {fi (um∗ )}i=sm∗ . To that end, we will need to ensure that the assumptions of Corollary 13
are verified. Let δ′ = log(n) and λ = λ2 . Observe that (48) and (49) are straightforward consequences of Assumption 3, (46) and (47) follow from the definiiton of B1 ({Xi }ni=1 , {um }M
m=1 , λ2 ) and
n
M
B2 ({Xi }i=1 , {um }m=1 , λ2 ), and that (45) follows from (13).
Thus, all the conditions in Corollary 13 are met, and we therefore conclude that there exists a
m∗
change point ηk , which is also a change point of {fi (v)}ei=s
, satisfying
m∗

min{em∗ − ηk , ηk − sm∗ } > ∆/4

(14)

and
|bm∗ − ηk | ≤ max{C3 λ22 κ−2 , δ′ } ≤ ǫn ,

where the last inequality holds because λ22 κ−2 = B 4 log(n)κ−2 ≥ log(n), which is a consequence of
the inequality B 2 ≥ κ.
The proof is complete with the following two observations: i) The change points of {fi (um∗ )}ei=s
belong to (s, e) ∩ {ηk }K
k=1 ; and ii) Equation (14) and (sm∗ , em∗ ) ⊂ (s, e) imply that
min{e − ηk , ηk − s} > ∆/4 > ǫn .
As discussed in the argument before Step 1, this implies that ηk must be an undetected change
point of {Xi }ni=1 in the covariance structure.

3

Lower bounds

In this section, we provide lower bounds for the problem of chenge point estimation with high
dimensional covariance matrices.
In Theorem 2 we showed that if the distribution of {Xi }ni=1 satisfies Assumption 1 and additionally the condition that ∆ ≥ CB 4 κ−2 p log(n) for sufficiently large C as given in Assumption 3, then
the WBSIP algorithm can, with high probability, detect all the change points with a localization
rate of the order
B 4 log n
.
(15)
κ2 n
Assumption 3 might seem a bit arbitrary at first glance. However, we will show in the next result
that if the class of distribution of interest allows for a spacing parameter ∆ = cB 4 κ−2 p for some
sufficiently small constant c, then it is not possible to estimate the location of the change point at
a vanoishing rate. This result furher implies that the WBSIP algorithm is optimal in the sense of
requiring the minimax scaling for the problem parameters.
To that effect, we will consider the following class of data generating distribution.
Assumption 4. Let X1 , . . . , Xn ∈ Rp be independent Gaussian random vectors such that Xi ∼
Np (0, Σi ) with kΣi kop ≤ 2σ. Let {ηk }K+1
k=0 ⊂ {0, . . . , n} be a collection of change points, such that
η0 = 0 and ηK+1 = n and that
Σηk +1 = Σηk +2 = . . . = Σηk+1 , for any k = 1, . . . , K + 1.
Assume there exists parameters κ = κ(n) and ∆ = ∆(n) such that
inf

{ηk − ηk−1 } ≥ ∆ > 0,

k=1,...,K+1

kΣηk − Σηk−1 kop = κk ≥ κ > 0, for any k = 1, . . . , K + 1.
18

n
We will denote with Pκ,∆,σ
the joint distribution of the data when the above assumption is in
n
effect. Notice that Pκ,∆,σ satisfies Assumption 1 with B = 8σ. Below we prove that consistent
2

4
estimation of the locations of the change point requires ∆κ
p B to diverge, as in Assumption 3. We
recall that all the parameters of interest, ∆, κ, σ and p are allowed to change with n. The proof is
based on a construction used in Cai and Ma (2013) to obtain minimax lower bounds for a class of
hypothesis testing problems involving covariance matices.

Lemma 3. Let {Xi }ni=1 be a time series satisfies Assumption 4 with only one change point and let
n
Pκ,∆,σ
denote the corresponding joint distribution. Consider the class of distributions
P = P(n) =



n
Pκ,∆,σ


σ4 p
2
:
≤ ∆ ≤ n/3, κ ≤ σ /4 ,
33κ2

Then
inf sup EP (|η̂ − η|) ≥ n/6.
η̂ P ∈P

Remark 9 (Phase transiton for the localization rate in covariance change point detection problem). In light of Lemma 3 and Theorem 2, we conclude that, in the covariance change
point detection problem, the solution undergoes a phase transition, which are able to characterize
up to logarithmic factor (in n). Specifically,
• if ∆ ≥ CB 4 p log(n)/κ2 for a sufficiently large constant C > 0, then it is possible to estimate
the locations of the change points with a localization rate vanishing in n;
• on the other hand, if ∆ = CB 4 p/κ2 for a sufficiently small constant c > 0, then the localization
rate of any algorithm remains, at least in the worst case, bounded away from 0.
To the best of our knowledge, this phase transition effect is new and unique to our settings.
We conclude this section by showing that the localization rate that we have obtained for the
WBSIP algorithm, given above in (15), is, up to a logarithmic factor, minimax optimal.
n
o
n
: ∆ < n/2, κ ≤ σ 2 /4 . Then,
Lemma 4. Consider the class of distribution Q = Pκ,∆,σ
inf sup EP (|η̂ − η|) ≥ cσ 4 κ−2 .
η̂ P ∈Q

4

Discussion

In this paper, we tackle the problem of change point detection for a time series of length n of
independent p-dimensional random vectors with covariance matrices that are pirece-wise constant.
We allow the dimension, as well as other parameters quantifying the difficulty of the problem, to
change with n. We have devised two procedures based on existing algorithms for change point
detection – binary segmentation and wild binary segmentation – and show that the localization
rates they yield are consistent with those in the univariate time series mean change point detection
problems. In particular we demonstrate the algorithm WBSIP, which applies wild binary segmentation to carefully chosen univariate projections of the data, produces a localization rate that is,
up to a logarithmic factor, minimax optimal.
19

The model setting adopted in the paper allows for the dimension p to grow with n. However,
in order for the localization rates of any procedure to vanish with n it must be the case that p
is of smaller order than n. One possible future direction is to consider different high dimensional
settings whereby p is permitted to grow even faster than n, with additional structural assumptions
on the underlying covariance matrices. For instance, we may model the covariance matrices as
spiked matrices with sparse leading eigenvectors.
A common, undesirable feature of the WBSIP algorithms is the fact that, for given interval
(s, e), the search of the next change point is limited to points inside the interval that are δ away
from the endpoints, where δ is an inout parameter that is larger than the localization error. Such
restriction, which appears also in the algorithm of Wang and Samworth (2016) for mean change
point localization of high-dimensional time series, is made in order to prevent the algorithms from
returning spurious change points in the proximity of a true change point. The reason for such
phenomenon is subtle, and is ultimately due to the fact that the rate of decay of the expected value
of the covariance CUSUM statistics around the true change points is in general not uniform, as it
depends on the magnitude of the change. A possible solution to such an issue – which appears to
be unavoidable – would be to design adaptive algorithm yielding local rates, one for each change
point. We will pursue this line of research in future work.
One key assumption used in this paper is the the time series of comprised of independent
observations. This is of course a rather strong condition, which might not apply to many real-life
problems. In order to handle time dependence, a natural approach is introduce mixing conditions
(see, e.g. Wu, 2005) or assume the observations come from certain well-defined time series models.
Further extensions that will be worth pursuing include the cases in which model is mis-specified or
the observations are contaminated We leave these interesting extensions to future work but expect
that many of the results derived in this manuscript will provide the theoretical underpinning for
devising and studying more complicated algorithms.

Acknowledgments
We thank Dr Haeran Cho for her constructive comments.

References
Anderson, T. W. (2003). An Introduction to Multivariate Statistical Analysis. 3rd ed. Wiley
Series in Probability and Statistics, John Wiley & Sons.
Aston, J. A. D. and Kirch, C. (2014). Efficiency of change point tests in high dimensional
settings. arXiv preprint arXiv: 1409.1771.
Aue, A., Hömann, S., Horváth, L. and Reimherr, M. (2009). Break detection in the covariance
structure of multivariate nonlinear time series models. The Annals of Statistics, 37 4046–4087.
Avanesov, V. and Buzun, N. (2016). Change-point detection in high-dimensional covariance
structure. arXiv preprint arXiv:1610.03783.
Baranowski, R., Chen, Y. and Fryzlewicz, P. (2016). Narrowest-Over-Threshold detection
of multiple change-points and change-point-like feature. arXiv preprint arXiv: 1609.00293.
20

Barigozzi, M., Cho, H. and Fryzlewicz, P. (2016). Simultaneous multiple change-point and
factor analysis for high-dimensional time series. arXiv preprint arXiv: 1612.06928.
Berkes, I., Gabrys, R., Horváth, L. and Kokoszka, P. (2009). Detecting changes in the
mean of functional observations. Journal of the Royal Statistical Society: Series B (Statistical
Methodology), 71 927–946.
Berthet, Q. and Rigollet, P. (2013). Optimal detection of sparse principal components in high
dimension. The Annals of Statistics, 41 1780–1815.
Birke, M. and Dette, H. (2005). A note on testing the covariance matrix for large dimension.
Statistics and Probability Letters, 74 281–289.
Cai, T. T. and Ma, Z. (2013). Optimal hypothesis testing for high dimensional covariance
matrices. Bernoulli, 19 2359–2388.
Chan, H. P. and Walther, G. (2013). Detection with the scan and the average likelihood ratio.
Statistica Sinica, 1 409–428.
Cho, H. (2015). Change-point detection in panel data via double cusum statistic. Electronic
Journal of Statistics in press.
Cho, H. and Fryzlewicz, P. (2012). Multiscale and multilevel technique for consistent segmentation of nonstationary time series. Statistica Sinica, 22 207–229.
Cho, H. and Fryzlewicz, P. (2015). Multiple change-point detection for high-dimensional time
series via Sparsified Binary Segmentation. Journal of the Royal Statistical Society: Series B
(Statistical Methodology), 77 475–507.
Davies, P. L. and Kovac, A. (2001). Local extremes, runs, strings and multiresolution. Ann.
Statist., 29 1–65.
Davis, R. A., Lee, T. C. M. and Rodriguez-Yam, G. A. (2006). Structural break estimation for
nonstationary time series models. Journal of the American Statistical Association, 101 223–239.
Frick, K., Munk, A. and Sieling, H. (2014). Multiscale change point inference. Journal of the
Royal Statistical Society: Series B (Statistical Methodology), 76 495–580.
Fryzlewicz, P. (2014). Wild binary segmentation for multiple change-point detection. The Annals
of Statistics, 42 2243–2281.
Gombay, E., Horváth, L. and Hušková, M. (1996). Estimators and tests for change in variances. Statistics and Risk Modeling, 14 145–160.
Harchaoui, Z. and Lévy-Leduc, C. (2010). Multiple change-point estimation with a total
variation penalty. Journal of American Statistical Association, 105 1480–1493.
Harchaoui, Z. and Lvy-Leduc, C. (2010). Multiple change-point estimation with a total variation penalty. Journal of the American Statistical Association, 105 1480–1493.
Horváth, L. and Hušková, M. (2012). Change-point detection in panel data. Journal of Time
Series Analysis, 33 631–648.
21

Inclan, C. and Tiao, G. C. (1994). Use of cumulative sums of squares for retrospective detection
of changes of variance. Journal of the American Statistical Association, 89 913–923.
James, B., James, K. L. and Siegmund, D. (1987). Tests for a change-point. Biometrika, 74
71–83.
Jirak, M. (2015). Uniform change point tests in high dimension. The Annals of Statistics, 43
2451–2483.
Johnstone, I. M. (2001). On the distribution of the largest eigenvalue in principal components
analysis. The Annals of Statistics, 29 295–327.
Killick, R., Fearnhead, P. and Eckley, I. A. (2012). Optimal detection of changepoints with
a linear computational cost. Journal of the American Statistical Association, 107 1590–1598.
Korkas, K. K. and Fryzlewicz, P. (2017). Multiple change-point detection for non-stationary
time series using wild binary segmentation. Statistica Sinica, 27 287–311.
Lavielle, M. (1999). Detection of multiple changes in a sequence of dependent variables. Stochastic
Processes and their Applications, 83 79–102.
Lavielle, M. and Moulines, E. (2000). Least-squares estimation of an unknown number of shifts
in a time series. Journal of Time Series Analysis, 21 33–59.
Ledoit, O. and Wolf, M. (2002). Some hypothesis tests for the covariance in principal components analysis. The Annals of Statistics, 30 1081–1102.
Li, H., Guo, Q. and Munk, A. (2017). Multiscale change-point segmentation: Beyond step
functions. arXiv preprint arXiv: 1708.03942.
Lin, K., Sharpnack, J. L., Rinaldo, A. and Tibshirani, R. J. (2017). A sharp error analysis
for the fused lasso, with application to approximate changepoint screening. In Advances in Neural
Information Processing Systems 30 (I. Guyon, U. V. Luxburg, S. Bengio, H. Wallach, R. Fergus,
S. Vishwanathan and R. Garnett, eds.). 6887–6896.
Olshen, A. B., Venkatraman, E., Lucito, R. and Wigler, M. (2004). Circular binary
segmentation for the analysis of array-based dna copy number data. Biostatistics, 5 557–572.
Page, E. S. (1954). Continuous inspection schemes. Biometrika, 41 100–115.
Picard, D. (1985). Testing and estimating change-points in time series. Advances in Applied
Probability, 17 841–867.
Qian, J. and Jia, J. (2012).
Https://arxiv.org/abs/1211.5194.

On pattern recovery of the fused lasso.

Tech. rep.

Rojas, C. R. and Wahlberg, B. (2014). On change point detection using the fused lasso method.
Tech. rep. Https://arxiv.org/abs/1401.5408.
Roy, S. N. (1957). Some Aspects of Multivariate Analysis. New York: Wiley.

22

Siegmund, D. and Venkatraman, E. S. (1995). Using the generalized likelihood ratio statistic
for sequential detection of a change-point. The Annals of Statistics 255–271.
Tibshirani, R., Saunders, M., Rosset, S., Zhu, J. and Knight, K. (2005). Sparsity and
smoothness via the fused lasso. Journal of the Royal Statistical Society: Series B (Statistical
Methodology), 67 91–108.
Venkatraman, E. S. (1992). Consistency results in multiple change-point problems. Ph.D. thesis,
Stanford University.
Vershynin, R. (2010). Introduction to the non-asymptotic analysis of random matrices. arXiv
preprint arXiv:1011.3027.
Vostrikova, L. (1981). Detection of the disorder in multidimensional random-processes. Doklady
Akademii Nauk SSSR, 259 270–274.
Wald, A. (1945). Sequential tests of statistical hypotheses. The Annals of Mathematical Statistics,
16 117–186.
Wang, T. and Samworth, R. J. (2016). High-dimensional changepoint estimation via sparse
projection. arXiv preprint arXiv:1606.06246.
Wang, Y. (1995). Jump and sharp cusp detection by wavelets. Biometrika, 82 385–397.
Wu, W. B. (2005). Nonlinear system theory: Another look at dependence. Proceedings of the
National Academy of Sciences of the United States of America, 102 14150–14154.
Yao, Y.-C. and Au, S.-T. (1989). Least-squares estimation of a stop function. Sankhyā: The
Indian Journal of Statistics, Series A 370–381.
Yu, B. (1997). Festschrift for Lucien Le Cam, vol. 423, chap. Assouad, Fano, and Le Cam. Springer
Science & Business Media, 435.

Appendices
A

Probabilistic bounds

In this section, we give basic high-probability concentration bounds on the fluctuations of te covariance CUSUM statistics using the notions of sub-Gaussian and sub-Exponential random vectors.
We also state some properties of the randomly selected intervals {sm , em } in the WBS algorithm,
which hold with high probability.
We start by introducing the definitions of sub-Gaussian and sub-Exponential random variables
through Orlicz norms. See, e.g., Vershynin (2010) for more details.
Definition 2.

(i) A random variable X ∈ R is sub-Gaussian if
kXkψ2 := sup k−1/2 {E(|X|k )}1/k < ∞.
k≥1

23

A random vector X ∈ Rp is sub-Gaussian if
kXkψ2 := sup kv ⊤ Xkψ2 < ∞,
v∈S p−1

where S p−1 denote unit sphere in Eucliden norm in Rp .
(ii) A random variable Y ∈ R is sub-Exponential if
kY kψ1 = sup k−1 E(|Y |k )1/k < ∞.
k≥1

A random vector Y ∈ Rp is sub-Exponential if
kY kψ1 = sup kv ⊤ Y kψ1 < ∞.
v∈S p−1

We note that if X ∈ R is sub-Gaussian, then X 2 is sub-Exponential, due to the easily verifiable
fact that
kXk2ψ2 ≤ kX 2 kψ1 ≤ 2kXk2ψ2 .
(16)
Recall the sample and the population versions of the covariance CUSUM statistics given in
Definition 1. For λ > 0, we define the events, which depend on {Xi }ni=1 .


s,e
s,e
n
e
e
≤ λ, min{t − s, e − t} ≥ p log(n)
(17)
A1 ({Xi }i=1 , λ) =
sup
St − Σ t
op

0≤s<t<e≤n

and

A2 ({Xi }ni=1 , λ) =

(

Pe

⊤
i=s+1 (Xi Xi

sup

0≤s<e≤n

√

e−s

− Σi )

op

≤ λ,

)

e − s ≥ p log(n) .

p
Next, for an arbitary collection {vm }M
m=1 of deterministic unit vectors in R we define the events
)
(
s,e
n
M
⊤ es,e
e )vm ≤ λ, min{t − s, e − t} ≥ log(n)
B1 ({Xi }i=1 , {vm }m=1 , λ) =
sup
sup
vm (St − Σ
t
1≤m≤M 0≤s<t<e≤n

(18)

and
B2 ({Xi }ni=1 , {vm }M
m=1 , λ)

=

(

sup

Pe

⊤
⊤
i=s+1 vm (Xi Xi

sup

1≤m≤M 0≤s<e≤n

√

e−s

− Σi )vm

≤ λ,

e − s ≥ log(n) .

Lemma 5. Suppose {Xi }ni=1 ⊂ Rp are i.i.d sub-Gaussian centered random vectors such that
sup kXi kψ2 ≤ B.

1≤i≤n

24

)

(19)

There exists an absolute constant c > 0 such that,
p
P(A1 ({Xi }ni=1 , B 2 p log(n))) ≥ 1 − 2 × 9p n3 n−cp ,
p
P(A2 ({Xi }ni=1 , B 2 p log(n))) ≥ 1 − 2 × 9p n2 n−cp ,
p
2
log(n))) ≥ 1 − 2n3 M n−c ,
P(B1 ({Xi }ni=1 , {vm }M
,
B
m=1
p
2
P(B2 ({Xi }ni=1 , {vm }M
log(n))) ≥ 1 − 2n2 M n−c ,
m=1 , B
for any set {vm }M
m=1 of deterministic unit vectors.

Proof. We first tackle A1 . For any v ∈ Rp such that kvk2 = 1, we can write
e s,e)v =
v (Sets,e − Σ
t
⊤

where

e
X

i=s+1

e


X
⊤
2
⊤
2
ai Zi ,
ai (v Xi ) − E((v Xi ) ) =
i=s+1

q
e−t

(e−s)(t−s)
q
ai =
t−s

(e−s)(e−t)

s + 1 ≤ i ≤ t,
t + 1 ≤ i ≤ e,

and Zi = (v ⊤ Xi )2 − E[(v ⊤ Xi )2 ]. Since min{t − s, e − t} ≥ p log(n), we further have that
e
X

a2i = 1,

and

i=s+1

p
max |ai | ≤ 1/ p log(n).

s+1≤i≤e

Thus by Proposition 5.16 in Vershynin (2010), for any ǫ ≥ 0,
!
(
)!
p
e
X
ǫ2 ǫ p log(n)
ai Zi ≥ ǫ ≤ 2 exp −c min
P
,
,
K2
K
i=s+1

where c > 0 is an absolute constant and
K = max kZi kψ1 ≤ 2k(v ⊤ Xi )2 kψ1 ≤ 4kXi k2ψ2 ≤ 4B 2 .
i

Therefore,



p

e s,e v ≥ B 2 p log(n) ≤ 2n−cp .
P v ⊤ Sets,e − Σ
t

Let N1/4 be a minimal 1/4-net (with respect to the Eucidean norm) of the unit sphere in Rp . Then,
card(N1/4 ) ≤ 9p and, by a standard covering argument followed by a union bound, we arrive at the
inequality
p
P(A1 ({Xi }ni=1 , B 2 p log(n))) ≥ 1 − 2 × 9p n3 n−cp ,
for a universal constant c > 0. Following the same arguments we have,
p
P(A2 ({Xi }ni=1 , B 2 p log(n))) ≥ 1 − 2 × 9p n2 n−cp ,
p
2
P(B1 ({Xi }ni=1 , {vm }M
,
B
log(n))) ≥ 1 − 2n3 M n−c ,
m=1
and

for some c > 0.

2
P(B2 ({Xi }ni=1 , {vm }M
m=1 , B

p
log(n))) ≥ 1 − 2n2 M n−c ,

25

M
Let {sm }M
m=1 , {em }m=1 be two sequences independently selected at random in [s, e], and

M=

K
\

sm ∈ Sk , em ∈ Ek , for some m ∈ {1, . . . , M } ,

(20)

k=1

where Sk = [ηk − 3∆/4, ηk − ∆/2] and Ek = [ηk + ∆/2, ηk + 3∆/4], k = 1, . . . , K. In the lemma
below, we give a lower bound on the probability of M. Under the scaling assumed in our setting,
this bound approaches 1 as n grows.
Lemma 6. For the event M defined in (20), we have


∆2
n
.
P(M) ≥ 1 − exp log − M
∆
16n2
Proof. Since the number of change points are bounded by n/∆,
M
K Y
X


P Mc ≤
1 − P sm ∈ Sk , em ∈ Ek ≤ K(1 − ∆2 /(16n2 ))M ≤ n/∆(1 − ∆2 /(16n2 ))M .
k=1 m=1

B

Properties of the univariate CUSUM statistics

In this section, we derive some important and useful properties of the univariate CUSUM statistic.
Our results and proofs build upon the existing literature on univariate mean change point detection; see in partiucular, Venkatraman (1992) and Fryzlewicz (2014), whose notation will be used
throughout. It is imprtnat however to note that we have made several non-trivial modifications
of those arguments, and have made a special effort in keeping track of the changes in all the key
parameters. This careful treatment eventually allows us to achieve tight upper bounds for the
the localization rate implied by the WBSIP algorithm and which in turn have revealed a phase
transition in the problem parameters (see Section 3). In particular, the results of this section can
be used to sharpen existing analyses of the BS and WBS algorithms.

B.1

Results from Venkatraman (1992)

We start by introducing some notations for one dimensional change point detection and the corresponding CUSUM statistics. Let {Yi }ni=1 , {fi }ni=1 ⊂ R be two univariate,
sequences. We will make the following assumptions on.
Assumption 5 (Univariate mean change points). Let {ηk }K+1
k=0 ⊂ {0, . . . , n}, where η0 = 0 and
ηK+1 = n, and
fηk−1 +1 = fηk−1 +2 = . . . = fηk for all 1 ≤ k ≤ K + 1,

Assume

inf

{ηk − ηk−1 } ≥ ∆ = ∆(n) > 0,

k=1,...,K+1

|fηk+1 − fηk | := κk > 0, k = 1, . . . , K,
sup

k=1,...,K+1

|fηk | < B1 .
26

as

We also have the corresponding CUSUM statistics over any generic interval [s, e] ⊂ [1, T ] defined
Yets,e

=

fets,e =

s

s

t
X
e−t
Yi −
(e − s)(t − s)
i=s+1

e−t
(e − s)(t − s)

t
X

i=s+1

fi −

s

s

e
X
t−s
Yi ,
(e − s)(e − t)
i=t+1

e
X
t−s
fi .
(e − s)(e − t)
i=t+1

Throughout this Appendix B, all of our results are proven by regarding {Yi }Ti=1 and {fi }Ti=1 as two
deterministic sequences. We will frequently assume that fets,e is a good approximation of Yets,e in
ways that we will specify through appropriate assumptions.
Observe that the function fets,e is only well defined on [s, e] ∩ Z. Our first result, which is
taken from Venkatraman (1992), shows that the there exists a continuous realization of the discrete
function fets,e
Lemma 7. Suppose [s, e] ⊂ [1, T ] satisfies

ηr−1 ≤ s ≤ ηr ≤ . . . ≤ ηr+q ≤ e ≤ ηr+q+1 ,

q ≥ 0.

Then there exists a continuous function Fets,e : [s, e] → R such that Fers,e = fers,e for every r ∈ [s, e]∩ Z
with the following additional properties.
1. |Fets,e | is maximized at the change points within [s, e]. In other words,
arg max |Fets,e | ∩ {ηr , . . . , ηr+q } =
6 ∅.
s≤t≤e

2. If Fets,e > 0 for some t ∈ (s, e), then Fets,e is either monotonic or decreases and then increases
within each of the interval [s, ηr ], . . . , [ηr+q , e].

The proof of this lemma can be found in Lemma 2.2 and 2.3 of Venkatraman (1992). We remark
that if Fets,e ≤ 0 for all t ∈ (s, e), then it suffices to consider the time series {−fi }Ti=1 and a similar
result as in the second part of Lemma 7 still holds. Throughout the entire section, we always view
fets,e as a continuous function and frequently invoke Lemma 7 as a basic property of the CUSUM
statistics without further notice.
Our next lemma is an adaptation of a result first obtained by Venkatraman (1992), which quantifies characterizes how fast the CUSUM statistics decays around a good change point. An analogous
result, derive using different arguments, can be found in Proposition 21 in Wang and Samworth
(2016).

Lemma 8 (Venkatraman (1992) Lemma 2.6). Let [s, e] ⊂ [1, T ] be any generic interval. For some
c1 , c2 > 0 and λ > 0 such that
min{ηk − s, e − ηk } ≥ c1 ∆
fes,e ≥ c2 κ∆(e − s)−1/2 ,
ηk

(21)
(22)

suppose there exists sufficient small c3 such that

max |fets,e | − feηs,e
≤ 2λ ≤ c3 κ∆3 (e − s)−5/2
k

s≤t≤e

27

(23)

Then there exists an absolute constant c > 0 such that if the point d ∈ [s, e] is such that |d − ηk | ≤
c1 ∆/16, then
feηs,e
− feds,e > cfeηs,e
|ηk − d|∆(e − s)−2
k
k

Remark 10. If feηs,e
k < 0, and d ∈ [s, e] is such that |d − ηk | ≤ c1 /16, then by considering the
sequence {−fi }ni=1 , it holds that
(−feηs,e
) − (−feds,e ) > c(−feηs,e
)|ηk − d|∆(e − s)−2
k
k

Proof. Without loss of generality, assume that d ≥ ηk . Following the argument of Venkatraman
(1992) Lemma 2.6, it suffices to consider two cases: (1) ηk+1 > e, and (2) ηk+1 ≤ e.
Case 1. Let El be defined as in the case 1 in Venkatraman (1992) Lemma 2.6. There exists
es,e (which in the notation of Venkatraman
a c > 0 such that, for every d ∈ [ηk , ηk + c1 ∆/16], feηs,e
k − fd
(1992) is the term El ) can be written as
feηs,e
|d−ηk | √
k

e − ηk

p

ηk − s + (d − ηk )

p

e−s

(ηk − s + (d − ηk ))(e − ηk ) +

p

(ηk − s)(e − ηk − (d − ηk ))

Using the inequality (e − s) ≥ 2c1 ∆, the previous expression is lower bounded by

.

≥ c′ |d − ηk |feηs,e
∆(e − s)−2 .
k

Case 2. Let h = c1 ∆/8 and l = d − ηk ≤ h/2. Then, following closely the initial calculations for
case 2 of Lemma 2.6 of Venkatraman (1992), we obtain that
− feds,e ≥ E1l (1 + E2l ) + E3l ,
feηs,e
k

where
E1l = p

(ηk − s + l)(e − ηk − l)

E2l = p

and

E3l = −

p

feηs,e
k l(h − l)

(ηk − s + l)(e − ηk − l) +

p

(ηk − s)(e − ηk )

,

((e − ηk − h) − (ηk − s))((e − ηk − h) − (ηk − s) − l)
 p

p
p
(ηk − s + l)(e − ηk − l) + (ηk − s + h)(e − ηk − h)
(ηk − s)(e − ηk ) + (ηk − s + h)(e − ηk − h)

(feηs,e
− feηs,e
k )l
k +h
h

s

(ηk − s + h)(e − ηk − h)
.
(ηk − s + l)(e − ηk − l)

Since h = c′′ ∆ and l ≤ h/2,

Observe that

E1l ≥ c′′ feηs,e
|d − η|∆(e − s)−2 .
k

ηk − s ≤ ηk − s + l ≤ ηk − s + h ≤ 9(ηk − s)/8,

e − ηk ≥ e − ηk − l ≥ e − ηk − h ≥ 7(e − ηk )/8.
(24)
28

Thus
E2l
= p

((e − ηk − h) − (ηk − s))2 + l(h + ηk − s) − l(e − ηk )
 p

p
p
(ηk − s + l)(e − ηk − l) + (ηk − s + h)(e − ηk − h)
(ηk − s)(e − ηk ) + (ηk − s + h)(e − ηk − h)

−l(e − ηk )
(ηk − s + h)(e − ηk − h)
−l(e − ηk )
≥
≥ −1/2
(ηk − s)(7/8)(e − ηk )
≥

where (24) is used in the second inequality and the fact that l ≤ h/2 ≤ c1 ∆/16 ≤ (ηk − s)/16 is
used in the last inequality. For E3l , observe that
≤ 2λ
≤ max |fets,e| − feηs,e
≤ |feηs,e
| − feηs,e
feηs,e
− feηs,e
k
k
k
k +h
k +h
s≤t≤e

and that (21) and l/2 ≤ h = c1 ∆/8 imply that

ηk − s ≤ ηk − s + l ≤ ηk − s + h ≤ 9(ηk − s)/8 and e − ηk ≥ e − ηk − l ≥ e − ηk − h ≥ 7(e − ηk )/8.
Therefore,
2(d − ηk )λ
E3l ≥ −
c1 ∆/8

s

(9/8)(ηk − s)(e − ηk )
(ηk − s)(7/8)(e − ηk )

32(d − ηk )λ
c1 ∆
′′
≥ −(c /4)fes,e (d − ηk )∆(e − s)−2 ,

≥−

ηk

where the first inequality follows from (24) and the last inequality follows from (22) and (23), for
a sufficiently small c3 . Thus,
feηs,e
− feds,e ≥ E1l (1 + E2l ) + E3l ≥ (c′′ /4)feηs,e
|ηk − d|∆(e − s)−2 .
k
k
The following proposition is a direct consequence of Lemma 8 and essentially characterized the
localization error rate of the BS algorithm.
Proposition 9. Consider any generic interval (s, e) ⊂ (0, T ) such that
ηr−1 ≤ s ≤ ηr ≤ . . . ≤ ηr+q ≤ e ≤ ηr+q+1 ,

q ≥ 0.

(25)

Let b ∈ arg maxs≤t≤e |Yets,e |. Suppose for some c1 > 0 and κ > 0,

max{min{ηr − s, s − ηr−1 }, min{ηr+q+1 − e, e − ηr+q }} = ǫn ,

where
ǫn < min{(3c1 /8)2 κ2 ∆2 (e − s)−1 B1−2 , ∆/4},
29

(26)

and
|Yebs,e | ≥ c1 κ∆(e − s)−1/2 .

(27)

Assume also that there exists sufficient small c3 > 0 such that
sup |fets,e − Yets,e | = λ1 < min{c3 κ∆3 (e − s)−5/2 ,

s≤t≤e

(c1 /4)κ∆(e − s)−1/2 }.

(28)

Then there exists a change point ηk ∈ [s, e] and an absolute constant C1 > 0 such that

and

min{e − ηk , ηk − s} > (3c1 /8)2 κ2 ∆2 (e − s)−1 B1−2

(29)

|feηs,e
| ≥ |Yebs,e | − λ1 ≥ max |fets,e| − 2λ1 .
k

(30)

|ηk − b| ≤ C1 λ1 (e − s)5/2 ∆−2 κ−1 ,
s≤t≤e

Proof. Observe that from (28),

max |fets,e| ≤ max |Yets,e | + λ1 ≤ |Yebs,e | + λ1 ≤ |febs,e | + 2λ1 .

s<t<e

s<t<e

(31)

Suppose ηk ≤ b ≤ ηk+1 for some r − 1 ≤ k ≤ r + q. Observe that

|febs,e | ≥ |Yebs,e | − λ1 > (3c1 /4)κ∆(e − s)−1/2 > 0,

where the second inequality follows from (27) and (28). It suffices to consider the case in which
febs,e > 0, since, if febs,e < 0, then the same arguments can be applied to the time series {−fi }ni=1 .
From Lemma 7, fets,e is either monotonic or decreasing and then increasing on [ηk , ηk+1 ]. Thus
} ≥ febs,e.
, feηs,e
max{feηs,e
k+1
k

es,e
If fets,e is locally decreasing at b, then feηs,e
k ≥ fb . Therefore

≥ febs,e > (3c1 /4)κ∆(e − s)−1/2 .
feηs,e
k

(32)

Step 1. We first show that (32) implies (29). For the sake of contradiction, suppose that min{e −
ηk , ηk − s} ≤ ((3c1 /8)κ∆(e − s)−1/2 B1−1 )2 . Then
r
r
(e
−
η
)(η
−
s)
(e − ηk )(ηk − s)
k
k
s,e
feηk ≤
B1 +
B1
e−s
e−s
p
≤ 2 min{e − ηk , ηk − s}B1 ≤ (3c1 /4)κ∆(e − s)−1/2 .
This is a contradiction to (32). Therefore (29) holds for ηk .

Step 2. We now apply Lemma 8, since (22) and (23) hold in virtue of (32) and (28), respectively.
Thus, we will need to prove that
min{ηk − s, e − ηk } ≥ (3/4)∆.
30

For the sake of contradiction, assume that min{ηk − s, e − ηk } < (3/4)∆. Suppose ηk − s < (3/4)∆.
Since ηk − ηk−1 ≥ ∆, one has s − ηk−1 ≥ ∆/4. This also means that ηk is the first change
point within [s, e]. Therefore k = r in (25). By (26), min{ηr − s, s − ηr−1 } ≤ ǫ < ∆/4. Since
s − ηk−1 = s − ηp−1 ≥ ∆/4, it must be the case that
ηk − s = ηr − s ≤ ǫ ≤ (3c1 /8)2 κ2 ∆2 (e − s)−1 B1−2 .
This is a contradiction to (29). Therefore ηk − s ≥ (3/4)∆. The argument of e − ηk ≥ (3/4)∆ can
be made analogously.
Step 3. By Lemma 8, if there exists ad and a sufficient large constant C1 satisfying
d ∈ [ηk , ηk + C1 λ1 (e − s)5/2 ∆−2 κ−1 ],
then

|ηk − d|∆(e − s)−2 ≥ λ1 ,
− feds,e > cfeηs,e
feηs,e
k
k

where the last inequality follows from (32) and (28).
For the sake of contradiction, suppose b ≥ d. Then

− λ1 ≤ max |fets,e| − 2λ1 ,
febs,e ≤ feds,e < feηs,e
k
s<t<e

where the first inequality follows from Lemma 7 which ensures that fets,e is decreasing on [ηp , b] and
d ∈ [ηp , b]. This is a contradiction to (31). Thus b ∈ [ηk , ηk + C1 λ1 (e − s)5/2 ∆−2 κ−1 ].
(30) follows from (31) and (32).
The argument for the case when febs,e is locally increasing at b is similar and therefore we omit the
details.
√
Corollary 10. Let δ > 0 be such that 2 δB1 ≤ (3c1 /4)κ∆(e−s)−1/2 . Let b′ = arg max⌈s+δ⌉≤t≤⌊e−δ⌋ |Yets,e |.
Suppose all the assumption in Proposition 9 hold except that (27) and (28) are replaced by
and

−1/2
|Yebs,e
′ | ≥ c1 κ∆(e − s)

sup

⌈s+δ⌉,...,⌊e−δ⌋

(33)

|fets,e − Yets,e | = λ1 < min{c3 κ∆3 (e − s)−5/2 ,

(c1 /4)κ∆(e − s)−1/2 },

(34)

respectively. Then, all the conclusions of Proposition 9 still hold for b′ .
Proof. For any t ∈ [s + δ] ∪ [e − δ],
p
√
|fets,e | ≤ 2 min{e − t, t − s}B1 ≤ 2 δB1 ≤ (3c1 /4)κ∆(e − s)−1/2 .

Let ηk and b defined as in the proof of Proposition 9. Then by (32) ηk , b ∈ ⌈s + δ⌉, . . . , ⌊e − δ⌋.
Therefore b = b′ and all the conclusions of Proposition 9 still hold.

31

B.2

Results from Fryzlewicz (2014)

Below, we will derive some further properties of the CUSUM statistic using the ANOVA decompositiontype of arguments first introduced by Fryzlewicz (2014), which are particularly effective at separating the noise from the signal. Some of the results below are contained in Fryzlewicz (2014), but
others requires different, subtle arguments. For completeness, we include all the proofs.
For a pair (s, e) of positive integers with s < e, let Wds,e be the two dimensional linear subspace
of Re−s spanned by the vectors
u1 = (1, . . . , 1, 0, . . . , 0)⊤
| {z } | {z }
d−s

u2 = (0, . . . , 0, 1, . . . , 1)⊤ .
| {z } | {z }

e−d

d−s

e−d

For clarity, we will use h , i to denote the inner product of two vectors in the Euclidean space.
Lemma 11. For x = (xs+1 , . . . , xe )⊤ ∈ Re−s , let Pds,e (x) be the projection of x onto Wds,e .
1. The projection Pds,e (x) satisfies
Pds,e (x)

e
1 X
=
xi + hx, ψds,e iψds,e ,
e−s
i=s+1

where h·, ·i is the inner product in Euclidean space, and ψds,e = ((ψds,e )s , . . . , (ψds,e )e−s )⊤ with
q
e−d

i = s + 1, . . . , d,
(e−s)(d−s) ,
s,e
q
(ψd )i =
d−s
−
(e−s)(e−d) , i = d + 1, . . . , e,
i.e. the i-th entry of Pds,e (x) satisfies
(
Pds,e (x)i

2. Let x̄ =

1
e−s

Pe

i=s+1 xi .

=

1
d−s
1
e−d

Pd
xj ,
Pj=s+1
e
j=d+1 xj ,

i = s + 1, . . . , d,
i = d + 1, . . . , e.

Since hx̄, ψds,e i = 0,

kx − Pds,e (x)k2 = kx − x̄k2 − hx, ψds,e i2 .

Proof. The results hold following the fact that the projection matrix of subspace Wds,e is


1/(d − s) · · · 1/(d − s)
0
···
0


..
..
..
..
..
..


.
.
.
.
.
.




1/(d
−
s)
·
·
·
1/(d
−
s)
0
·
·
·
0
s,e
.
PW s,e = 

0
···
0
1/(e − d) · · · 1/(e − d) 
d




..
..
..
..
..
..


.
.
.
.
.
.
0
···
0
1/(e − d) · · · 1/(e − d)

32

(35)

For any pair d1 , d2 ∈ {s + 1, . . . , e} and f ∈ Re−s , the following two statements are equivalent:
hf, ψds,e
i2 ≤ hf, ψds,e
i2 ⇐⇒ kf − Pds,e
(f )k2 ≥ kf − Pds,e
(f )k2 .
1
2
1
2
Lemma 12. Assume Assumption 1. Let [s0 , e0 ] be an interval with e0 − s0 ≤ CR ∆ and containing
at lest one change point ηr such that
ηr−1 ≤ s0 ≤ ηr ≤ . . . ≤ ηr+q ≤ e0 ≤ ηr+q+1 ,

q ≥ 0.

Suppose that min{ηp′ − s0 , e0 − ηp′ } ≥ ∆/16 for some p′ and let κs,e
max = max{κp : min{ηp − s0 , e0 −
ηp } ≥ ∆/16}. Consider any generic [s, e] ⊂ [s0 , e0 ], satisfying
min{ηr − s0 , e0 − ηr } ≥ ∆/16

for all ηr ∈ [s, e].

Let b ∈ arg maxs<t<e |Yets,e |. For some c1 > 0, λ > 0 and δ > 0, suppose that
√
|Yebs,e | ≥ c1 κs,e
max ∆
sup |Ye s,e − fes,e | ≤ λ
s<t<e

t

(36)

t

and

sup √

s1 <t<e1

1
e1 − s1

e1
X

(Yt − ft ) ≤ λ

for every

t=s1 +1

If there exists a sufficient small c2 > 0 such that
√
λ ≤ c2 κs,e
max ∆ and

e1 − s1 ≥ δ.

δ ≤ c2 ∆,

(37)

(38)

then there exists a change point ηk ∈ (s, e) such that
min{e − ηk , ηk − s} > ∆/4

and

|ηk − b| ≤ min{C3 λ2 κ−2
k , δ}.

Proof. Without loss of generality, assume that febs,e > 0 and that fets,e is locally decreasing at b.
Observe that there has to be a change point ηk ∈ [s, b], or otherwise febs,e > 0 implies that fets,e is
decreasing, as a consequence of Lemma 18.
Thus if s ≤ ηk ≤ b ≤ e, then
√
√
√
s,e
s,e
≥ febs,e ≥ |Yebs,e | − λ ≥ c1 κs,e
feηs,e
(39)
max ∆ − c2 κmax ∆ ≥ (c1 /2)κmax ∆.
k

Observe that e − s ≤ e0 − s0 ≤ CR ∆ and that that (s, e) has to contain at least one change point
or otherwise |feηs,e
k | = 0 which contradicts (39).

Step 1. In this step, we are to show that min{ηk − s, e − ηk } ≥ min{1, c21 }∆/16.
Suppose ηk is the only change point in (s, e). So min{ηk − s, e − ηk } ≥ min{1, c21 }∆/16 must
hold or otherwise it follows from Lemma 17, we have
√
c1 √
c1
| < κk ∆ ≤ κs,e
|feηs,e
max ∆,
k
4
2

which contradicts (39).

33

Suppose (s, e) contains at least two change points. Then ηk − s ≤ min{1, c21 }∆/16 implies that
ηk is the first change point in [s, e]. Therefore
√
√
1 es,e
1
es,e| + c1 κr ∆
|
+
2κ
|
f
max
|
f
|
≤
|feηs,e
η
−
s
≤
r
r
η
t
k
4 k+1
4 s<t<e
2
√
3 e s,e
c1 s,e
1 e s,e
≤ |Yb | + λ + κmax ∆ ≤ |Yb | + λ < |Yebs,e | − λ
4
2
4
where the first inequality follows from Lemma 18, the fourth inequality follows from (36), and the
last inequality holds when c2 is sufficiently small. This contradicts (39).
Step 2. By Lemma 8 there exists d such that

√
−1
d ∈ [ηk , ηk + λ ∆(κs,e
max ) ]

es,e > 2λ. For the sake of contradiction, suppose b ≥ d. Then
and that feηs,e
k − fd

febs,e ≤ feds,e < feηs,e
− 2λ ≤ max |fets,e| − 2λ ≤ max |Yets,e | + λ − 2λ = |Yebs,e | − λ,
k
s<t<e

s<t<e

that fets,e is decreasing on [ηp , b] and
where the first inequality follows from Lemma 7, which ensures √
−1
d ∈ [ηp , b]. This is a contradiction to (39). Thus b ∈ [ηk , ηk + λ ∆(κs,e
max ) ].

Step 3. Let f s,e = (fs+1 , . . . , fe )⊤ ∈ R(e−s) and Y s,e = (Ys+1 , . . . , Ye )⊤ ∈ R(e−s) . By the definition
of b, it holds that
Y s,e − Pbs,e (Y s,e )

2

(Y s,e )
≤ Y s,e − Pηs,e
k

2

2

≤ Y s,e − Pηs,e
(f s,e) .
k

For the sake of contradiction, throughout the rest of this argument suppose that, for some sufficient
large constant C3 to be specified,
ηk + max{C3 λ2 κ−2
k , δ} < b.
(This will of course imply that ηk +
the bound

−2
max{C3 λ2 (κs,e
max ) , δ}

Y s,e − Pbs,e (Y s,e )

2

(40)

< b). We will show that this leads to
2

(f s,e ) ,
> Y s,e − Pηs,e
k

(41)

which is a contradiction.
derive (41) from (40), we note that min{e−ηk , ηk −s} ≥ min{1, c21 }∆/16
√ Tos,e
and that |b − ηk | ≤ λ ∆(κmax )−1 implies that
√
−1
min{e − b, b − s} ≥ min{1, c21 }∆/16 − λ ∆(κs,e
≥ min{1, c21 }∆/32,
(42)
max )
where the last inequality follows from (38) and holds for an appropriately small c2 > 0.
Equation (41) is in turn implied by
2hεs,e , Pb (Y s,e ) − Pηk (f (s,e) )i < kf s,e − Pb (f s,e )k2 − kf s,e − Pηk (f s,e )k2 ,

where εs,e = Y s,e − f s,e. By (35), the right hand side of (43) satisfied the relationships
kf s,e − Pb (f s,e )k2 − kf s,e − Pηk (f s,e )k2 = hf s,e, ψηk i2 − hf s,e, ψb i2
= (fes,e)2 − (fes,e )2
ηp

b

|
− febs,e)|feηs,e
≥ (feηs,e
k
k
s,e
2
≥ c|d − ηk |(fe ) ∆−1
′

≥ c |d −

34

ηk
2
ηk |(κs,e
max ) ,

(43)

where Lemma 8 and (39) are used in the second and third inequalities. The left hand side of (43)
can in turn be rewritten as
2hεs,e , Pb (X s,e ) − Pηk (f s,e)i = 2hεs,e , Pb (X s,e ) − Pb (f s,e)i + 2hεs,e , Pb (f s,e ) − Pηk (f s,e)i.

(44)

The second term on the right hand side of the previous display can be decomposed as


ηk
e
b
X
X
X
 εs,e (Pb (f s,e)i − Pηk (f s,e )i )
+
+
hεs,e , Pb (f s,e ) − Pηk (f s,e )i = 
i
i=s+1

i=ηk +1

i=b+1

= I + II + III.

In order to bound the terms I, II and III, observe that, since e − s ≤ e0 − s0 ≤ CR ∆, the interval
[s, e] must contain at most CR + 1 change points. Let
ηr′ −1 < s ≤ ηr′ ≤ . . . ≤ ηp′ +q′ < e ≤ ηp′ +q′ +1 .

Then p′ + q ′ + 1 − r ′ ≤ CR + 1.
Step 4. We can write
√

I = ηk − s

√

1
ηk − s

ηk
X

i=s+1

εs,e
i

!

ηk
b
X
1 X
1
fi −
fi
b−s
ηk − s
i=s+1

i=s+1

!

.

Thus,
Pk
P
Pk
ηk
b
X
(ηk − s)( ηi=s+1
fi + bi=ηk +1 fi ) − (b − s) ηi=s+1
fi
1 X
1
fi =
fi −
b−s
ηk − s
(b − s)(ηk − s)
i=s+1
i=s+1
P
Pk
Pk
fi + (ηk − s) bi=ηk +1 fi )
(ηk − b) ηi=s+1
(ηk − b) ηi=s+1
fi + (ηk − s)(b − ηk )fηk +1 )
=
=
(b − s)(ηk − s)
(b − s)(ηk − s)
=

ηk
X
b − ηk
b − ηk
1
fi + fηk+1 ≤
(CR + 1)κs,e
−
max
b−s
ηk − s
b−s
i=s+1

where Lemma 16 is used in the last inequality. It follows from Equation (37) that
√
√
|b − ηk |
4 2
s,e
(CR + 1)κmax ≤
|b − ηk |∆−1/2 λ(CR + 1)κs,e
|I| ≤ ηk − sλ
max ,
b−s
min{1, c1 }
where (42) is used in the last inequality.

Step 5. For the second term II, we have that



d
b
e
X
X
X
p
1
1
1
εs,e  
fi −
fi 
|II| = b − ηk  √
b−s
e − ηk
b − ηk i=η +1 i
i=s+1
i=ηk +1
k


b
e
X
X
p
1
1
fi − fηk +
fi − fηk+1 
≤ b − ηk λ  fηk − fηk+1 +
b−s
e − ηk
i=s+1
i=ηk +1
p
s,e
s,e
≤ b − ηk (κs,e
max + (CR + 1)κmax + (CR + 1)κmax ),
35

where the first inequality follows from (42) and (37), and the second inequality from Lemma 16.
Step 6. Finally, we have that
III =
Therefore,
|III| ≤

√

√

e−b

1
e−b

e
X

i=b+1

!
εs,e 
i

1
e − ηk

e
X

i=ηk +1

fi −

1
e−b

e
X

i=b+1



fi  .

√
b − ηk
4 2
s,e
(CR + 1)κmax ≤
|b − ηk |∆−1/2 λ(CR + 1)κs,e
e − bλ
max .
e−b
min{1, c1 }

Step 7. Using the first part of Lemma 11, the first term on the right hand side of (44) can be
bounded as
hεs,e , Pd (X s,e ) − Pd (f s,e )i ≤ λ2
. Thus (43) holds if

n
o
p
2
−1/2
s,e
s,e
2
b
−
η
|b − ηk |(κs,e
)
≥
C
max
|b
−
η
|∆
λκ
,
.
λκ
,
λ
k
k
max
max
max
√
Since λ ≤ c3 ∆κ, the first inequality holds. The second inequality follows from |b − ηk | ≥
−2
C3 λ2 (κk )−2 ≥ C3 λ2 (κs,e
max ) , as assumed in (40). This completes the proof.
Corollary 13. Let [s0 , e0 ] be a generic interval satisfying e0 − s0 ≤ CR ∆ and containing at lest
one change point ηr such that
ηr−1 ≤ s0 ≤ ηr ≤ . . . ≤ ηr+q ≤ e0 ≤ ηr+q+1 ,
κs,e
max

q ≥ 0.

= max{κp : min{ηp −s0 , e0 −ηp } ≥
Suppose min{ηp′ −s0 , e0 −ηp′ } ≥ ∆/16 for some p′ and, let
∆/16}. Consider a generic interval (s, e) ⊂ (s0 , e0 ), satisfying
min{ηp − s0 , e0 − ηp } ≥ ∆/16

for all ηp ∈ [s, e].

Let δ′ > 0 be some constant and b′ = arg max⌈s+δ′ ⌉,...,⌊e−δ′ ⌋ |Yets,e |. Suppose in addition that, for
some positive constants c1 and c2 ,
√
s,e
∆,
(45)
|Yebs,e
′ | ≥ c1 κmax
s,e
s,e
sup
|Ye − fe | ≤ λ,
(46)
t=⌈s+δ′ ⌉,...,⌊e−δ′ ⌋

t

e1
X
1
(Yt − ft )| ≤ λ,
|
e1 − s1 t=s +1
s1 ≤t≤e1
1
√
s,e
λ ≤ c2 κmax ∆,

sup

and

t

√

for every

δ′ ≤ c2 ∆.

e1 − s1 ≥ δ′ ,

(47)

(48)
(49)

Then there exists a change point ηk ∈ [s, e] such that

min{e − ηk , ηk − s} > ∆/4

′
|ηk − b′ | ≤ max{C3 λ2 κ−2
k , δ }.

Proof. By the same proof of Corollary 10, if b is defined as in Lemma 12, then b = b′ if c2 is
sufficiently small.
36

C

Proofs of the Results from Section 3

e u = σ 2 Ip + κuu′ . Observe that if κ ≤ σ 2 /4,
Proof of Lemma 3. For any vector u ∈ Rp , denote Σ
e u kop = σ 2 + κ ≤ 2σ 2 .
kΣ
Step 1.
p
n denote the joint distribution of independent random vectors {X }n
Let Pe0,u
i i=1 in R such that
i.i.d.
e u)
X1 , . . . , X∆ ∼ Np (0, Σ

i.i.d.

and X∆+1 , . . . , Xn ∼ Np (0, σ 2 I).

n denote the joint distribution of independent random vectors {X }n
p
Similarly, let Pe1,u
i i=1 in R with
i.i.d.

X1 , . . . , Xn−∆ ∼ Np (0, σ 2 I)

P

i.i.d.
e u ).
and Xn−∆+1 , . . . , Xn ∼ Np (0, Σ

n
√ en
u∈{±1}p / p Pi,u . Let η(Pi,u ) denote the location of the
n . Then since η(P
en ) = ∆ and η(Pen ) = n − ∆
to the distribution Pei,u
0,u
1,u
n
n
e
e
|η(P0,u ) − η(P1,u )| ≥ n/3. By Le cam’s lemma (see, e.g. Yu, 1997),

Let Pein =

1
2p

inf

sup

η̂ P ∈P n
κ
e ,∆,σ

change point associated
√
for any u ∈ {±1}p / p,

EP (|η̂ − η|) ≥ (n/3)(1 − dT V (Pe0n , Pe1n )),

where T V (Pe0n , Pe1n ) = 21 kPe0n − Pe1n k1 .
n
Let Σu = Ip + κ
euu′ , where κ
e = κ/σ 2 . Observe that by assumption κ
e ≤ 1/4. Denote with P0,u
n the joint distributions of independent samples {X }n
p
and P1,u
i i=1 in R where
i.i.d.

i.i.d.

X1 , . . . , X∆ ∼ Np (0, Σu )
and

i.i.d.

X1 , . . . , Xn−∆ ∼ Np (0, I)

and X∆+1 , . . . , Xn ∼ Np (0, I)
and

ii.d.

Xn−∆+1 , . . . , Xn ∼ Np (0, Σu ),

respectively. Since total variation distance is invariant under rescaling of the covariance, then
kPe0n − Pe1n k1 = kP0n − P1n k1 . Therefore
inf

sup

η̂ P ∈P n
κ
e ,∆,σ

EP (|η̂ − η|) ≥ (n/3)(1 − dT V (P0n , P1n ))

(50)

Step 2.
Let x = (x1 , . . . , x∆ ), y = (x∆+1 , . . . , xn−∆ ) and z = (xn−∆+1 , . . . , xn ). Let
• f0 (x) denote the joint distribution of X1 , . . . , X∆ ∼ Np (0, I) and fu (x) denote the joint
i.i.d.

distribution of X1 , . . . , X∆ ∼ Np (0, Σu );

i.i.d.

• g0 (x) denote the joint distribution of X∆+1 , . . . , Xn−∆ ∼ Np (0, I);
i.i.d.

• h0 (x) denote the joint distribution of Xn−∆+1 , . . . , Xn ∼ Np (0, I) and hu (x) denote the
i.i.d.

joint distribution of Xn−∆+1 , . . . , Xn ∼ Np (0, Σu ).

37

Then,
kP0n − P1n k1
Z Z Z
1
=
2p
=

=

≤

Z

g0 (y)dy

Z Z

1
2p

Z Z

1
2p

X

√
u∈{±1}p / p

Z Z

1
2p

fu (x)g0 (y)h0 (z) −
X

√
u∈{±1}p / p

1
2p

fu (x)h0 (z) −

X

√
u∈{±1}p / p

1
2p

f0 (x)g0 (y)hu (z) dxdydz

X

√
u∈{±1}p / p

X

fu (x)h0 (z) −

X

fu (x)h0 (z) − f0 (x)h0 (z) + f0 (x)g0 (z) −

√
u∈{±1}p / p

√
u∈{±1}p / p

X

f0 (x)hu (z) dxdz

1
2p

√
u∈{±1}p / p

f0 (x)hu (z) dxdz
1
2p

X

√
u∈{±1}p / p

=2kP0∆ − P1∆ k1 ,

i.i.d.

where P0∆ is the joint distribution of X1 , . . . , X∆ ∼ Np (0, I) and P1∆ =

f0 (x)hu (z) dxdz

1
2p

i.i.d.

P

∆
√
u∈{±1}p / p P1,u ,

∆ is the joint distribution of X , . . . , X
where P1,u
1
∆ ∼ Np (0, Σu ). Thus (50) becomes

inf

sup

η̂ P ∈P n
κ
e ,∆,σ

EP (|η̂ − η|) ≥ (n/3)(1 − kP0∆ − P1∆ k1 ).

(51)

Step 3.
To bound 2kP0∆ − P1∆ k1 , let P0 = Np (0, Ip ) and Pu = Np (0, Σu ). It is easy to see that
χ

2

(P1∆ , P0∆ )

For any u, v ∈ {±1}p ,
EP ∆
0



= EP ∆
0



dPu∆ dPv∆
dP0∆ dP0∆

2
1
dP1∆
−1 = p
∆
4
dP0



=



EP0



X

EP ∆
0

u,v∈{±1}p

dPu dPv
dP0 dP0

∆



dPu∆ dPv∆
dP0∆ dP0∆



− 1.

= (1 − (e
κu′ v)2 )−∆/2 ,

where the last equality follows from Lemma 14. Denote U and V to be the p dimensional
Rademacher variables with U being independent of V and εp = (1′ V /p)2 . Thus
χ2 (P1∆ , P0∆ ) =

1
4p

X

u,v∈{±1}p

(1 − (e
κu′ v)2 )−∆/2 − 1


∆/2 
κU ′ V /p)2
−1
= EU,V 1 − (e


∆/2
= EV 1 − (e
−1
κ1′ V /p)2

≤ E exp(εp κ
e2 ∆) − 1,
38

where √
the inequality follows from inequality (1 − t)−∆/2 ≤ exp(∆t) for any t ≤ 1/2 and that
κ
e ≤ 1/ 2. The Höffding’s equality, applied to Rademacher variables, gives
P (εp ≥ λ) ≤ 2e−2pλ .

(52)

Thus
Z

∞


P exp εp κ
e2 ∆ ≥ u du
0
Z ∞

≤1+
P ε ≥ log(u)/(e
κ2 ∆) du


Z1 ∞
log(u)2p
2 exp −
≤1+
du
κ
e2 ∆
1
2
,
=1−
1 − κe2p
2∆


E exp(εp κ
e ∆) =
2

where the second inequality follows from (52), and the last equality holds if κe2p
2 ∆ > 1. Thus if
q
4
2p
= 2pσ
, then, using the well-known fact that kP0∆ − P1∆ k1 ≤ 2 χ2 (P1∆ , P0∆ ), we obtain
∆ = 33e
κ2
33κ2
the bounds
s
q
2
∆
∆
= 1/2.
kP0 − P1 k1 ≤ 2 χ2 (P1∆ , P0∆ ) ≤ 2
2p
−1
κ
e2 ∆
This and (51) give
inf sup EP (|η̂ − η|) ≥ n/6.
η̂ P ∈P

Proof of Lemma 4. Let

 n
Q = Pκ,∆
: ∆ < n/2 .

Using the same argument as in Step 1 of the proof of Lemma 3, it suffices to take
κ = κ/σ 2 ≤ 1/4,
e

and consider Σu = Ip + κ
euu′ .
Let u to be any unit vector in Rp and δ a positive number. Let P0n denote the joint distribution
of independent samples {Xi }ni=1 in Rp where
iid.

iid.

X1 , . . . , X∆−1 ∼ Np (0, I),

X∆ , . . . , Xn ∼ Np (0, Σu )

and P1n the joint distribution of independent samples {Xi }ni=1 in Rp where
iid.

X1 , . . . , X∆+δ ∼ Np (0, I),

iid.

X∆+δ+1 , . . . , Xn ∼ Np (0, Σu ).

n . Then since η(P n ) = ∆ and
Let η(Pin ) denote the location of the change point of distribution Pi,u
0
η(P1n ) = ∆ + δ, thus |η(P0n ) − η(P1n )| ≥ δ. By Le cam’s lemma (Yu, 1997),

1
inf sup EP (|η̂ − η|) ≥ δ(1 − dT V (P0n , P1n )) = δ(1 − kP0n − P1n k1 )
η̂ P ∈Q
2
39

Since P0n and P1n are identically distributed on X1 , . . . X∆ − 1 and on X∆+δ , . . . , Xn , kP0n − P1n k1 =
kP0δ − P1δ k1 , where P0δ is the joint distribution of X1 , Xδ ∼ Np (0, I) and P1δ is the joint distribution
of X1 , Xδ ∼ Np (0, Σu ). By Lemma 14,
χ2 (P1δ , P0δ ) = (1 − κ
e2 )−δ/2 − 1 ≤ 4e
κ2 δ.

if κ2 δ/2 ≤ 1/4 and δ/2 ≥ 2. Thus by taking δ = e
κ−2 /4 ≥ 4 we obtain that


q
1
δ
δ
2
χ (P1 , P0 ) ≥ e
inf sup EP (|η̂ − η|) ≥ δ 1 −
κ−2 /32 = σ 4 κ−2 /32.
η̂ P ∈Q
2
Lemma 14. Let P0 = Np (0, Ip ) and Pu = Np (0, Ip + κuu′ ). Then


−1/2
dPu dPv
EP0
= 1 − (κu′ v)2
.
dP0 dP0
Proof. See, e.g., lemma 5.1 in Berthet and Rigollet (2013).

Lemma 15. For t ≥ 2, x ≥ 0, if tx ≤ 1/4, (1 − x)−t − 1 ≤ 4tx
Proof. There exists s ∈ [0, x] such that
(1 − x)−t − 1 = tx + t(t + 1)x2 (1 − s)−t−2 ≤ tx + 4t2 x2 (1 − x)−t ≤ 2tx + 4t2 x2 ((1 − x)−t − 1).
Thus
(1 − x)−t − 1 ≤

D

2tx
≤ 4tx
1 − 4t2 x2

Properties of the covariance CUSUM statistic

D.1

Properties of 1d CUSUM statistics

Lemma 16. Suppose [s, e] ⊂ [1, T ] such that e − s ≤ CR ∆, and that
ηr−1 ≤ s ≤ ηr ≤ . . . ≤ ηr+q ≤ e ≤ ηr+q+1 ,

q ≥ 0.

Denote
κs,e
max = max{ηp − ηp−1 : r ≤ p ≤ r + q}.
Then for any r − 1 ≤ p ≤ r + q,
e

1 X
fi − fηp ≤ CR κs,e
max .
e−s
i=s

40

Proof. Since e − s ≤ CR ∆, the interval [s, e] contains at most CR + 1 change points. Observe that
e

1 X
fi − fηp
e−s
i=s

ηr+1
ηr
X
1 X
=
(fηr−1 − fηp ) +
(fηr − fηp ) + . . . +
e−s
i=s

i=ηr +1

e
X

ηr+1
ηr
X
1 X
s,e
≤
|p − r|κmax +
|p − r − 1|κs,e
max + . . . +
e−s

≤

1
e−s

i=s
e
X

(fηr+q − fηp )

i=ηr+q +1

i=ηr +1

e
X

i=ηr+q +1

|p − r − q − 1|κs,e
max

(CR + 1)κs,e
max ,

i=s

where |p1 − p2 | ≤ CR + 1 for any ηp1 , ηp2 ∈ [s, e] is used in the last inequality.
Lemma 17. If ηp is the only change point in [s, e], then
r
q
(ηp − s)(e − ηp )
s,e
e
κp ≤ min{ηp − s, e − ηp }κp
|fηp | =
e−s
Lemma 18. Let [s, e] contains two or more change points such that
ηr−1 ≤ s ≤ ηr ≤ . . . ≤ ηr+q ≤ e ≤ ηr+q+1 ,

q ≥ 1.

If
ηr − s ≤ c21 ∆
then

√
|feηs,e
| ≤ c1 |feηs,e
| + 2κr ηr − s.
r
r+1

This can be useful in testing when there are exactly two change points with
ηr − s ≤ λ2 κ−2
r ,

e − ηr+1 ≤ λ2 κ−2
r+1 .

It is also useful to show ηr − s ≥ ∆/4 for some absolute constant c when
| ≥ max |fets,e | − 2λ.
|feηs,e
r
s≤t≤e

Proof. Consider the sequence {gt }et=s+1 be such that
(
fηr+1 if s + 1 ≤ t ≤ ηr ,
gt =
ft if ηr + 1 ≤ t ≤ e.
For any t ≥ ηr ,
feηs,e
− geηs,e
=
r
r

s

√
(e − s) − t
(ηr − s)(fηr+1 − fηr ) ≤ ηr − sκr
(e − s)(t − s)
41

Thus
√
| + ηr − sκr
| ≤ |e
gηs,e
|feηs,e
r
r
s
√
(ηr − s)(e − ηr+1 ) s,e
≤
|e
gηr+1 | + ηr − sκr
(ηr+1 − s)(e − ηr )
r
√
c21 ∆ s,e
|e
gηr+1 | + ηr − sκr
≤
∆
√
e
≤ c1 |fηs,e
| + 2 ηr − sκr .
r+1

where the first inequality follows from the observation that the first change point of gt in [s, e] is
at ηr+1 .

D.2

Properties of the covariance CUSUM statistics

All of our consistency results heavily rely on the properties of population quantity of the CUSUM
statistic. In the covariance change point detection problem, however, it is not trivial to analyze the
e s,e
properties of the function t 7→ Σ
in the multiple change point case. For example, it is difficult
t
op
e s,e
to determine the regions of monotonicity of Σ
as a function of t as is done in Venkatraman
t

op

(1992, Lemma 2.2). As a remedy, we introduce the concept of shadow vector, which is defined as
a maximizer of the operator norm of the CUSUM statistics in all the following results. In this
way, we turn the covariance change point detection problem into a mean change point detection
problem.
For any v ∈ Rp with kvk = 1, let Yi (v) = (v ⊤ Xi )2 and fi (v) = v ⊤ Σi v, i = 1, . . . , n. Note
that both {Yi (v)} and {fi (v)} are univariate sequences, we hence have the corresponding CUSUM
statistics defined below
s
s
t
e
X
X
e
−
t
t−s
s,e
e
Yi (v) −
Yi (v),
(53)
Yt (v) =
(e − s)(t − s)
(e − s)(e − t)
i=s+1
i=t+1
s
s
t
e
X
X
e
−
t
t−s
s,e
fi (v) −
fi (v).
(54)
fet (v) =
(e − s)(t − s)
(e − s)(e − t)
i=s+1

i=t+1

The key rationale of the CUSUM based BS algorithm or any variants thereof being a powerful
tool selecting the change points is that the population version of the CUSUM statistic achieving
its maxima at the true change points. In Lemma 19, we show the same holds for the covariance
CUSUM statistic.

e s,e
Lemma 19. Assume (s, e) ∩ {ηk }K
k=1 6= ∅ and Assumption 1. The quantity Σt
of t achieves its maxima at the true change points, i.e.
arg max

t=s+1,...,e−1

e s,e
Σ
t

op

op

as a function

∩ {ηk }K
k=1 6= ∅.

Proof. For the sake of contradiction suppose that there exists t∗ ∈ (s, e) \ {ηk }K
k=1 such that
t∗ ∈ arg max

t=s+1,...,e−1

42

e s,e
Σ
t

op

,

and

e s,e
Σ
t∗

op

>

e s,e
Σ
ηk

max

k: ηk ∈(s,e)

op

.

n
′
n
e s,e
Let v ∈ arg maxkuk=1 |u′ Σ
t∗ u|, and consider the sequence {fi (v)}i=1 = {v Σi v}i=1 . By the above
display, we have

e s,e
fets,e
Σ
∗ (v) =
t∗

>

op

max

k: ηk ∈(s,e)

e s,e
Σ
ηk

op

≥

max

k: ηk ∈(s,e)

feηk (v) ,

(55)

where fets,e (v) is defined in (54). It follows from Lemma 2.2 of Venkatraman (1992), the quantities
fets,e (v) is maximized at the change points of the time series {ft (v)}et=s+1 . Note that the change
points of the sequence {ft (v)}et=s+1 is a subset of {ηk }K
k=1 , this contradicts (55).
Lemma 19 shows that the population version of the covariance CUSUM statistic is maximized
at the true change points in terms of the operator norm. In Lemma 20 below, we give the lower
bound of the maxima thereof. One can interpret it as the signal strength.
Lemma 20. Under Assumption 1, let 0 ≤ s < ηk < e ≤ n be any interval satisfying
min{ηk − s, e − ηk } ≥ c1 ∆.
Then for any 0 < δ < (c1 /8)∆,
max

t=⌈s+δ⌉,...,⌊e−δ⌋

Proof. Let

e s,e kop ≥ (c1 /2)κ∆(e − s)−1/2 .
kΣ
t


v ∈ arg max u⊤ Σηk − Σηk+1 u .
kuk=1

Therefore

we have

e s,e
Σ
t
op

≥

fets,e(v)

. Since


fηk (v) − fηk+1 (v) = v ⊤ Σηk − Σηk+1 v = Σηk − Σηk+1
max

t=⌈s+δ⌉,...,⌊e−δ⌋

e s,e kop ≥
kΣ
t

max

t=⌈s+δ⌉,...,⌊e−δ⌋

op

≥ κ,

|fets,e(u)| ≥ (c1 /2)κ∆(e − s)−1/2 ,

where the second inequality follows from the same arguments in Venkatraman (1992, Lemma 2.4)
by regarding fets,e(u) as the CUSUM statistic for a univariate time series.

In Lemma 21, we show that we can conform the problem of change point detection on covariance
to the one on mean via shadow vectors. This translation allows us to convert a high dimensional
covariance problem to a univariate mean problem.
Lemma 21. Under the same assumptions as Lemma 20, let
b∈

and denote a shadow vector by

arg max
t=⌈s+δ⌉,...,⌊e−δ⌋

Sets,e

op

v ∈ arg max u⊤ Sebs,e u .

,

kuk=1

Then

b∈

arg max
t=⌈s+δ⌉,...,⌊e−δ⌋

43

Yets,e (v) .

(56)

Proof. It suffices to show that
Yebs,e (v) = Sebs,e

op

≥

max

t=⌈s+δ⌉,...,⌊e−δ⌋

Sets,e

op

≥

max

t=⌈s+δ⌉,...,⌊e−δ⌋

Yets,e (v) .

For the shadow vector v defined in (56), it is tempting to argue that we can estimate v and
then use {Yi (v)}ni=1 as the new sequence to derive an upper for the rate of localization. This ideal
approach, however, suffers from two non-trivial obstacles.
• We don’t have any guarantee that the estimation of v is consistent. This is because estimating
the first principle direction of any sample matrix in general requires a non-vanish gap between
the first and the second eigenvalues of the corresponding population matrix. Since b depends
on the data, without more involved additional assumption, it is difficult to show Sebs,e converges
to a population quantity as n → ∞.
• Suppose v given in (56) has a well defined population quantity, say E(v). The estimation of
v depends on {Xi }ni=1 and therefore E(Yi (v)) 6= fi (E(v)).

On one hand our knowledge of the population version of the shadow vector v is very limited; on the
other hand, in Lemma 22 below we show that without estimating the population shadow vector,
the CUSUM statistics and its corresponding population quantity are close enough and that the
maximum of the CUSUM statistic is detectable.
Lemma 22. Under the same assumptions as Lemma 20, assume
max

t=⌈s+δ⌉,...,⌊e−δ⌋

e s,e
Sets,e − Σ
t

op

≤ λ,

where λ1 > 0. For the shadow vector v defined in (56), then
max

t=⌈s+δ⌉,...,⌊e−δ⌋

Yets,e (v) ≥ (c1 /2)κ∆(e − s)−1/2 − λ.

Proof. Observe that for all compatible t and v,


e s,e
e s,e v ≤ Ses,e − Σ
Yets,e (v) − fets,e (v) = v ⊤ Sets,e − Σ
t
t
t

op

≤ λ.

In addition, for the shadow vector v, we have
max

t=⌈s+δ⌉,...,⌊e−δ⌋

≥

max

t=⌈s+δ⌉,...,⌊e−δ⌋

Yets,e (v) =
e s,e
Σ
t

op

max

t=⌈s+δ⌉,...,⌊e−δ⌋

Sets,e

op

− λ1 ≥ (c1 /2)κ∆(e − s)−1/2 − λ,

where the last inequality follows from Lemma 20.
We finally prove a simple property of covariance CUSUM statistics.
p
e s,e kop ≤ 2 min{e − t, t − s}B 2 for any t ∈ (s, e).
Lemma 23. kΣ
t
Proof. Observe that

r
p
(e − t)(s − t) 2
s,e
e
kΣt kop ≤ 2
B ≤ 2 min{e − t, t − s}B 2 .
e−s
44

