Optimal change point detection in Gaussian processes
Hossein Keshavarz1 , Clayton Scott2,1 and XuanLong Nguyen1,2

arXiv:1506.01338v2 [math.ST] 8 Apr 2017

1

2

Department of Statistics, University of Michigan
Department of Electrical Engineering and Computer Science, University of Michigan
Email: {hksh,clayscot,xuanlong}@umich.edu
April 11, 2017
Abstract
We study the problem of detecting a change in the mean of one-dimensional Gaussian process data.
This problem is investigated in the setting of increasing domain (customarily employed in time series
analysis) and in the setting of fixed domain (typically arising in spatial data analysis). We propose a
detection method based on the generalized likelihood ratio test (GLRT), and show that our method
achieves nearly asymptotically optimal rate in the minimax sense, in both settings. The salient feature
of the proposed method is that it exploits in an efficient way the data dependence captured by the
Gaussian process covariance structure. When the covariance is not known, we propose the plug-in
GLRT method and derive conditions under which the method remains asymptotically near optimal. By
contrast, the standard CUSUM method, which does not account for the covariance structure, is shown to
be asymptotically optimal only in the increasing domain. Our algorithms and accompanying theory are
applicable to a wide variety of covariance structures, including the Matern class, the powered exponential
class, and others. The plug-in GLRT method is shown to perform well for maximum likelihood estimators
with a dense covariance matrix.

1

Introduction

Change point detection is the problem of detecting an abrupt change or changes arising in a sequence of
observed samples. A common problem of this type involves detecting shifts in the mean of a temporal
process. This problem has found a variety of applications in many ﬁelds, including audio analysis [20], EEG
segmentation [36], structural health monitoring [28, 42] and environment sciences [35, 56]. Despite advances
in the development of algorithms [33,36,39,48] and asymptotic theory [6,38,49,54] for a number of contexts,
such studies are mainly conﬁned to the setting of (conditionally) independently distributed data. Existing
works on optimal detection of shifts in the mean in temporal data with statistically dependent observations
are far less common.
Incorporating dependence structures into the modelling of random processes is a natural approach. In
fact, this has been considered in detecting changes of remotely collected data [11, 19]. For instance, Chandola et al. [11] proposed a Gaussian process based algorithm to identify changes in Normalized Diﬀerence
Vegetation Index (NDVI) time series for a particular location in California. Despite such statistical modelling considerations, to our knowledge, most researchers have not exploited the dependence structures of
the underlying temporal process, e.g., its covariance function and spectral density, in designing a (minimax)
optimal detection method.
In this paper, we focus on the detection of a single change in the mean of a Gaussian process data
sequence. Our main contribution is to show that neglecting the dependence structures in the data samples
leads to suboptimal detection procedures, particularly in the presence of strong correlation among samples.

1

Moreover, it is possible to exploit the underlying dependence structures to design asymptotically optimal
detection algorithms.
Consider a simpliﬁed setting in which we let G be a Gaussian process in a domain D ⊆ R and Dn :=
n
n
{tk }k=1 ⊂ D denote the ﬁnite index set of sampling points. Denote the observed samples by X = {Xk }k=1
in which Xk = G (tk ) for k = 1, . . . , n. Moreover, let t ∈ Cn,α ⊆ {1, . . . , n} (the parameter α is a positive
scalar which will be introduced in Section 2.3) and b > 0 represents the point of sudden change and the
jump/shift value, respectively. Namely, there is µ ∈ R (which will be assumed to be 0 for now) such that




b
b
EXk = µ −
1 (k ≤ t) + µ +
1 (k > t) , k ∈ {1, . . . , n} .
(1.1)
2
2
To analyze the performance of a detection procedure as sample size n grows to inﬁnity, one is confronted with
two fundamentally diﬀerent theoretical frameworks: increasing domain asymptotics and fixed domain (infill)
asymptotics ([46], Chapter 13). The former arises naturally in time series analysis, which is distinguished
by the constraint that the distance between consecutive sampling time points are bounded away from zero.
The simplest instance of the sampling scenario in this regime arises when the diameter of Dn is of order
n and min |ti+1 − ti | > δ for some strictly positive, ﬁxed scalar δ. In our notation
S the index set for the
Gaussian process represents the sampling time points. Typically we set D = R and ∞
n=1 Dn = N or Z. See,
[5,25,26,34,47,59] for examples of works studying change point detection via increasing domain asymptotics.
Fixed domain asymptotics, one the other hand, is a more suitable setting when the index set of sampling
points D is bounded, so that the observations get denser in D as n increases. Particularly for D ⊂ R,
we have that min |ti+k − ti | = O (k/n) for positive integers i, k with i, (i + k) ∈ {1, . . . , n}, and it can be
extended to multidimensional domains in a straightforward way. This is the case for spatially distributed
data [53], where the domain of the index set is typically of one, two or three dimensions. But there are also
other examples: a particularly useful approach to change point detection in non-stationary processes widely
adopted in speech processing and ﬁnance can also be cast in this framework. In fact, this approach involves
piecewise locally stationary (PLS) processes, which can be interpreted as processes that are approximately
piecewise stationary, due to the gradual and smooth change of the spectrum [1,14]. Abrupt change detection
in PLS processes has been considered in, for instance, [1, 35]. Since the future samples of a PLS process
may not carry any information about the current state of the process, the increasing domain setting cannot
suitably capture the dependence structure of the observed data sequence. As a result, the theory of change
detection and inference in PLS processes is also based upon ﬁxed domain asymptotics.
Our goal in this work is to derive a computationally eﬃcient detection procedure that eﬀectively accounts
for the underlying dependence structures of the observed sequential data. For the analysis of such a procedure
we shall adopt both aforementioned asymptotic frameworks. Obviously, speciﬁc applications ultimately
determine which one among the two is more suitable. There are also several motivations for this dual
treatment. First, we are not aware of any established detection algorithm and accompanying asymptotic
theory under the ﬁxed-domain setting. Second, although there is a vast literature on the change point
detection in the increasing domain scenario, existing work focus mostly on independent (or conditionally
independent) data sequences, so our analysis in this setting still carries some notable novelty. Third, there is
a fundamental diﬀerence in the behavior of the same detection algorithm when applied to the two asymptotic
settings. This point is worth highlighting: we will show that in the ﬁxed domain setting, ignoring the
underlying dependence structure may result in suboptimal detection performance, but this is not the case for
the increasing domain setting. Finally, we note that the ﬁxed domain minimax detection problem considered
in this paper also serves as a useful starting point in the study of optimal detection of discontinuities in
Gaussian spatial processes, as considered recently by [19, 50].
The amount of correlation between the nearby samples is one of the diﬀerences among the two asymptotic
frameworks that will plays a crucial role in our results. In the ﬁxed domain regime, regardless of how large
the sample size is, if |j − i| is of order nβ for some β ∈ (0, 1), the correlation among Xi and Xj is still so
close to one. Roughly speaking, the eﬀective sample size is much smaller than n. However, in the increasing
domain setting, even for long range dependent processes the correlation among samples at time points i, j
is small when |j − i| is large. Accordingly, the asymptotic behavior of a statistical method applied to the
2

ﬁxed-domain regime is expected to be distinct from that of the increasing domain regime. Moreover, new
techniques are needed to address the statistical dependence intrinsic in the former.
From here on, for ﬁxed domain asymptotic results, we assume that G is a one dimensional Gaussian
n
process in D = [0, 1] with n regularly spaced samples, i.e. Dn = {k/n}k=1 . Under the increasing domain
asymptotics framework, we let D = [0, ∞) and Dn = {1, . . . , n}. The remaining notations are the same for
both.
Previous works. One of the earliest attempts to study shift in mean detection was perhaps that of
Chernoﬀ et al. [12]. More general settings of this problem have been studied in subsequent works, e.g.,
[15, 41, 59]. For instance it is assumed in [59] that the sequence of Xk ’s are independent Gaussian variables.
They proposed a detection method based on generalized likelihood ratio test (GLRT), also known as the
cumulative sum (CUSUM) test, and given by
)
(
r

n
t
X
1
1X
t (n − t)
TCUSUM = 1 max
(1.2)
Xk −
Xk ≥ Rn .
t∈Cn,α
n
n−t
t
k=t+1

k=1

CUSUM compares the maximum of a test statistic over Cn,α with a critical values Rn . Non-asymptotic
upper bounds on the error probabilities of this simple test were obtained by the authors under the Gaussian
and i.i.d. assumptions. Interestingly, due to test’s simplicity, when such assumptions do not hold, one can
still apply the same test statistic to the data sequence. Most subsequent works appeared to follow this
direction, in addition to adhering to the increasing domain asymptotic framework, e.g. [5, 25, 26, 47]. We
wish to mention Rencova ([47], chapter 4), who studied the same CUSUM test as [59], but working with
the assumption that X be a strong mixing time series. Kokoszka [34] also analyzed the CUSUM test, but
working with a diﬀerent dependent observation model P
with sub-squared growth of the variance of partial
δ
sum, i.e. there is δ ∈ (0, 2) such that for any k < m, var m
j=k Xj . (m − k + 1) . Horváth et al. [25,26] and
Antoch [5] studied the performance
of the CUSUM test for the detection of a sudden change in the mean in
P∞
∞
linear processes, i.e. Xt = j=0 wj ǫt−j , in which {ǫt }t=−∞ are i.i.d. and zero mean random variables and
∞
the weights {wj }j=0 satisfy some properties such as absolute or square summability.
At ﬁrst sight, it may seem puzzling how the CUSUM test continues to admit nearly optimal detection
performance even as its test statistics apparently ignore the dependence among data samples.
S∞ A high-level
explanation can be made regarding this phenomenon: when G is a Gaussian process with n=1 Dn = N and
the covariance function cov (Xs , Xt ) → 0 as |t − s| grows to inﬁnity, the percentage of pairs (Xs , Xt )ns,t=1
whose covariance is non-negligible tends to zero as n → ∞. As a result, there is very little gain in accounting for the dependence structures underlying the sequence, and so the CUSUM statistic provides a good
approximation of the likelihood ratio test for large n, leading to the asymptotic optimality of TCUSUM in
the increasing domain setting. This is of course not the case for the ﬁxed domain setting. Indeed, one of
the contributions of this paper is to show that the CUSUM test is suboptimal when applied to the ﬁxed
domain setting of the detection problem. Moreover, to achieve mimimax optimal detection performance we
shall develop a new test statistic that account for the underlying dependence structures in the data sequence.
We also note in passing that in comparison to the increasing domain analysis typically encountered in the
literature, the theoretical analysis for the ﬁxed domain setting is considerably more challenging, as one needs
to take into account the statistical dependence in the data sequence in a more fundamental way.
CUSUM test also applied to one dimensional processes with highly correlated samples, after a proper
standardization. For instance, Horváth et al. [26] used a diﬀerent normalizing factor for applying CUSUM to
one dimensional Gaussian time series with long range dependence. However apart from standardizing factor,
they do not directly incorporate the correlation structures of the data in the formulation of the test statistic.
Note that, the proposed test in this paper achieve consistency under weaker condition on minimal detectable
jump. Furthermore, Lai [37] adjusts CUSUM test for detecting abrupt changes in the mean of a stochastic
process. However we are not aware of the analysis of change point problem, for the case of unknown before
and after distribution parameters, in the inﬁll regimes and its comparison with the increasing domain setup.
There are also some notable works (see e.g., [52]) on estimating the volatility parameter of non-stationary
time series using the change point models.
3

Overview of main results. In this paper we study the change point detection problem that arises in
Gaussian processes in both settings of increasing and ﬁxed domains. We show when it is important to
account for the dependence structure in the data sequence, and analyze a number of detection algorithms
based upon a generalized likelihood ratio test. More speciﬁcally, our contributions are as follows.
1. Given an n-sample drawn from a one dimensional Gaussian process data with a known covariance
structure, we propose a generalized likelihood ratio test for detecting a sudden shift in the mean. This
method requires the knowledge of the dependence structure (via the covariance matrix), and will be
shown to achieve asymptotically optimal detection performance in both ﬁxed and increasing domain
settings. Our theory holds for a variety of covariance structures, such as the Matern class, powered
exponential class, and several others speciﬁed in terms of the covariance kernel’s spectral density. The
smoothness parameter for the Gaussian process (which determines how fast the corresponding spectral
density decays) plays a central role in characterizing the minimax optimal detection performance —
but this is shown to be the case only in the ﬁxed domain setting, not in the increasing domain setting.
2. We provide an upper bound guarantee for the CUSUM detection method. This result conﬁrms that
the CUSUM is asymptotically optimal in the increasing domain setting, but it also suggests that the
CUSUM is suboptimal in the ﬁxed domain setting. The suboptimality is conﬁrmed in our simulation
study, which demonstrates a wide gap between the CUSUM and our method. This result makes sense,
in light of the minimax result described in the preceding paragraph.
3. In practice, the covariance structure is not known a priori, and often has to be estimated from the
data. To address this scenario, we propose a plug-in GLRT method, and analyze its performance. In
particular, we derive detection performance bounds which also account for the quality of a particular
covariance estimation method (such as MLE based methods with dense or tapered covariances) used
in the plug-in GLRT. Most surprisingly, we show that as long as a consistent covariance estimate is
employed (the notion of consistency will be deﬁned in Section 4), regardless of its estimation rate, the
plug-in GLRT achieves asymptotically optimal detection performance. Moreover, in some situations a
plug-in GLRT with an inconsistent covariance estimate is shown to perform almost as well as the case
of known covariance.
In addition to the above contributions, our proof methods contain several useful techniques. Our proofs
integrate four major technical tools: we exploit properties of the mutually absolutely Gaussian measures,
the decorrelation of samples drawn from Gaussian processes in ﬁxed domain, the non-asymptotic analysis
of the inverse of large Toeplitz matrices, in addition to the classical theory of minimax detection. We
also develop novel techniques for analyzing diﬀerent norms of the decorrelation matrix of X in either of
the asymptotic regimes. The Appendices provide several beneﬁcial and easy-to-reference technicalities for
theoretical problems in the area of Gaussian random ﬁelds and time series that may be of independent
interest.
Structure of the paper. Section 2 precisely formulates detection of a shift in the mean, focusing on
one dimensional Gaussian process data, and then introduces our proposed detection algorithm as hypothesis
testing in the case of fully known spectral density of the underlying process. We also adapt our proposed
detection technique, which will be referred as the plug-in test, to the much more realistic case of unknown
spectral density. Note that Sections 3, 5 and 6 are divided into two subsections focusing on theoretical results
in the ﬁxed domain and the increasing domain setting, respectively. Section 3 presents suﬃcient conditions
on shift value b and spectral density to detect the existence of shift in mean with high probability. Section 4
is devoted to analyze the performance of plug-in test by imposing some suﬃcient conditions on the estimated
spectral density. Section 5 serves as a comprehensive study of the CUSUM test. The minimax optimality
of the proposed algorithms will be discussed in section 6. Section 7 is devoted to the numerical experiments
and assessing the proposed algorithms using simulation studies. Section 8 contain concluding remarks with a
concise discussion of future directions. Appendix A contains the proofs of the main results and Appendix B

4

presents and proves some auxiliary results used in Appendix A. Lastly, Appendix C develops non-asymptotic
results on the inverse of large Toeplitz matrices which are useful in the study of CUSUM test.
Notation. ∧ and ∨ stand for minimum and maximum operators and the indicator function is represented
by 1 (·). For any m ∈ N, Im , 0m and 1m respectively denote the m by m identity matrix, all zeros column
vector of length
P m, and all ones column vector of length m. For two matrices of the same size M1 and M2 ,
hM1 , M2 i := i,j (M1 )ij (M2 )ij denotes their usual inner product. For any symmetric matrix M , λmin (M )
represents the smallest eigenvalue of M . We will use the following matrix norms on M ∈ Rm×n . For any 1 ≤
1/p
P
p
p < ∞, kM kℓp :=
|M
|
stands for the element-wise ℓp norm of M , while kM kℓ∞ := maxi,j |Mij |
ij
i,j
R
p
p
represents the sup norm of M . For a function f : D 7→ R and p > 0, kf kp := D |f (u)| du. The special
case of p = ∞ is deﬁned by kf k∞ := supu∈D |f (u)|. For any f ∈ L1 (R), fˆ represents its Fourier transform
deﬁned by
Z∞
fˆ (ω) =
f (t) e−jωt dt, ∀ ω ∈ R,
−∞

2

where j = −1 denotes the imaginary unit. Moreover, for a symmetric function f ∈ L∞ ([−π, π]), {fm }m∈Z
denotes the set of the its Fourier coeﬃcients. Assuming that {fm }m∈Z is absolutely summable, TN (f )
represents the inﬁnite Toeplitz matrix generated by the Fourier coeﬃcient of f , i.e.


f0 f1 f2 · · ·


 f1 f0 f1 . . . 


TN (f ) := 
(1.3)
.
 f2 f1 f0 . . . 


.. . .
.
.
. .. ..
.

h
in
Moreover, Tn (f ) = (TN (f ))ij

i,j=1

denotes a n × n truncated Toeplitz matrix generated by f . For two
f (t)
t→t0 g(t)

functions f and g on R, we write f (t) ≍ g (t) as t → t0 , if C1 ≤ lim

≤ C2 for some strictly

positive bounded scalars C1 ≤ C2 . In particular, we write f (t) ∼ g (t) as t → t0 to indicate the case that
C1 = C2 = 1. Furthermore, for sequences an and bn , we write bn = Ω (an ) when bn is bounded below by
an asymptotically, i.e. lim |bn /an | ≥ C for some positive C. In case an and bn are random, bn = oP (an )
n→∞

means that bn /an converges in probability to zero as n → ∞. For two sets Ω1 , Ω2 ⊂ Rm , dist (Ω1 , Ω2 ) :=
inf ωi ∈Ωi , i=1,2 kω1 − ω2 kℓ2 stands for their mutual distance with respect to Euclidean distance. Lastly, Γ (·)
denotes the gamma function.

2

Problem formulation and detection algorithms

In this section we present a formulation of the shift-in-mean detection problem associated with Gaussian
process data, and then describe detection algorithms that account for the underlying process modelling
n
assumptions. Recall that G is a Gaussian process in D ⊂ R and {Xk = G (tk )}k=1 represents the set of n
n
samples of G at times Dn = {tk }k=1 ⊂ D. At this point we proceed to split the problem formulation into
two subsections according to the two diﬀerent asymptotic settings.

2.1

Gaussian process model in fixed domain setting

In the ﬁxed domain setting, G − EG is assumed to be a mean-zero stationary Gaussian process on a bounded
domain D = [0, 1] and regularly sampled at tk = k/n for k = 1, . . . , n. Let the symmetric real functions K :

5

R 7→ R and K̂ : R 7→ R respectively denote the covariance function and spectral density of G. Accordingly,
n
Σn := cov ({Xk }k=1 ) is a symmetric Toeplitz matrix given by
n
 
r−s
n
.
(2.1)
Σn = {cov Xr Xs }r,s=1 = K
n
r,s=1
In the next section, we impose some regularity conditions on K.

2.2

Gaussian process model in increasing domain setting

In the increasing domain setting, G − EG is a mean-zero Gaussian process in D = [0, ∞) and so {Xk }k
endowed with a Toeplitz covariance function. This is a common setting for time series data, where the
observed samples are indexed by time points in D. It is customary to assume that the covariance of the
observed samples decreases as the temporal distance increases. Deﬁne cov (X1 , Xk ) = fk for any k, in
∞
which {fm }m=0 be a absolutely summable sequence with f0 = 1. Due to the stationarity assumption,
∞
n
ΣN := cov ({Xk }k=1 ) is an inﬁnite symmetric Toeplitz matrix. We view {Xk }k=1 as the observed part of an
∞
n
inﬁnite stationary time series, {Xk }k=1 . Accordingly, the covariance matrix of {Xk }k=1 , denoted by Σn , is
a symmetric (truncated) Toeplitz matrix.
It is a known fact that (Chapter 4, [22]) there is a symmetric and almost surely (with respect to Lebesgue
measure) positive function, f : [−π, π] 7→ R such that ΣN = TN (f ). Thus Σn = Tn (f ). For studying the
asymptotic properties of the change detection algorithm, certain regularity conditions are required on f .

2.3

Detection procedure based on generalized likelihood ratio test

Now we proceed to formulate the detection of the existence of a sudden change in the mean of a one
dimensional Gaussian process as a composite hypothesis testing problem. As noted above, we are dealing
with two diﬀerent settings of the domains. That is, we assume that G satisﬁes either of the two conditions:
(a) Fixed domain setting: G is restricted to D = [0, 1] where its spectral density admits Assumption 3.1
n
and Dn = {k/n}k=1 .
(b) Increasing domain setting: The domain of G is D = R and the samples are taken at Dn = {1, . . . , n}.
Moreover, Σn = Tn (f ) for some f fulﬁlling Assumption 3.2.
Although the domain settings are diﬀerent, the detection procedure that we propose based on a generalized likelihood ratio test will be the same. The composite hypothesis testing problem is set out as follows.
Under the null hypothesis, all the random variables have zero mean, i.e. EX = 0n . To specify the alternative
hypothesis H1 we ﬁrst introduce a few additional notations. Let t ∈ Cn,α denote the occurrence time of the
single change point. The set Cn,α ⊆ {1, . . . , n} contains plausible occurrence time of the change, and we
assume there is α ∈ (0, 1/2) such that Cn,α = {t : t ∧ (n − t) > αn}. Another important parameter b denotes
the amount of shift in the mean before and after the change point. Thus, for a ﬁxed t ∈ Cn , the associated
alternative hypothesis to t can be stated as,
H1,t : ∃ b 6= 0, EX =

b
ζt ,
2

(2.2)

where ζt ∈ Rn is given by ζt (k) := sign (k − t) for any t ∈ Cn,α . Since t is not known a priori, the alternative
hypothesis is speciﬁed by taking the union of H1,t . Thus, the composite hypothesis testing problem is given
by
[
b
(2.3)
H0 : EX = 0n , v.s. H1 =
H1,t , i.e. ∃ t ∈ Cn,α , b 6= 0, s.t. EX = ζt .
2
t∈Cn,α

Next, we propose a test statistic which is constructed by the generalized likelihood ratio (GLR). Note that
the GLR is an explicit function of the joint density of samples and so the Gaussian process assumption is
essential to its calculation.
6

Proposition 2.1. Assuming that Σn is known, there exists Rn,δ > 0 for which the GLRT is given by


TGLRT

−1

2

ζ ⊤ (Σn ) X

= I  max qt
t∈Cn,α
−1
ζt⊤ (Σn ) ζt


≥ Rn,δ  .

(2.4)

The threshold value Rn,δ depends only on n and some parameter δ determining the false alarm rate.
The precise form of Rn,δ will be presented in subsequent sections. We also note that setting µ = 0 in (1.1)
results in a substantially simpliﬁed expression of the GLR, which eases the exposition of our analysis of the
computational and theoretical properties of the proposed test. The general form of the GLRT, when µ is
unknown, is presented as Proposition A.1 in Appendix A.
Unlike the CUSUM test, cf. Eq. (1.2), the covariance function of G is explicitly taken into account in the
GLRT. As a result, it will be shown in the sequel that the proposed detection method is optimal, while the
same cannot be said for the CUSUM test, speciﬁcally in the setting of ﬁxed domain asymptotics. In practice,
however, the covariance is not known and needs to be estimated. To address such scenarios, we propose to
approximate the likelihood ratio by plugging in a positive deﬁnite estimate of the covariance matrix, which
will be indicated by Σ̃n , in Eq. (2.4). The plug-in detection technique will be called plug-in GLRT. Here
is the formulation of the plug-in GLRT, while the choice of Σ̃n and the accompanying theory will be given
later in Section 4.
Definition 2.1. Let Σ̃n be a positive deﬁnite estimate of Σn . The plug-in GLRT is given by

T̃GLRT



 −1
X
Σ̃n

r
= I
max
t∈Cn,α
 −1
ζt
ζt⊤ Σ̃n
ζt⊤

2




≥ R̃n,δ 
,

(2.5)

for some strictly positive threshold value R̃n,δ .

3

Detection rate of GLRT: known Σn

This section is devoted to a treatment of the detection rate of the GLRT, given that Σn is known. Section
3.1 addresses this problem in the ﬁxed domain, while Section 3.2 considers the increasing domain. We ﬁrst
deﬁne the risk measure that we will consider in the subsequent sections.
Definition 3.1. For any change detection algorithm T ∈ {0, 1}, the conditional detection error probability
(CDEP) of T which is denoted by ϕn (T ), is deﬁned as
ϕn (T ) = P (T = 1 | H0 ) + max P (T = 0 | H1,t ) .
t∈Cn,α

Remark 3.1. In words, ϕn is the sum of the false positive rate and the maximal miss detection rate taken
over the set of possible change point locations Cn,α . Clearly, CDEP hinges on the choices of Cn,α – the value
of ϕn increases as Cn,α becomes a larger proper subset of {1, . . . , n}. CDEP as a risk measure has been
considered in the literature for detecting abnormal clusters in a network (see e.g., [4, 10]). It also provides
an upper bound on the Bayesian risk measure. We refer the reader to [2] for a prudent comparison of CDEP
and the Bayesian risk measure.
In the following results, we seek for suﬃcient conditions on the shift value b under which the CDEP of
GLRT is bounded above by some small δ ∈ (0, 1).

7

3.1

Detection rate in fixed domain setting

The results in this sub-section are guaranteed for two common classes of covariance function, one of which
admits polynomially decaying spectral density, and the other is Gaussian covariance function. Throughout
this section we assume that G is a one dimensional Gaussian process restricted to D = [0, 1] whose covariance
function and spectral density are denoted by K and K̂, respectively. We ﬁrst focus on the case of polynomially
decaying K̂.
Assumption 3.1. K is an integrable positive deﬁnite covariance function. Moreover, there exist ν ∈ (0, ∞)
and CK (depending on K) for which K̂ satisﬁes the following condition:
CK := sup K̂ (ω) 1 + ω 2
ω∈R

ν+ 12

< ∞.

(3.1)

We shall always choose the largest possible ν that satisﬁes (3.1). It is simple to see that Assumption 3.1
holds if and only if K̂ is bounded at the origin and K̂ (ω) ≍ ω −(2ν+1) as ω tends to inﬁnity. It is well-known
that the tail behavior of K̂ is closely linked to the smoothness of K at the origin (e.g., Section 2.8, [53]).
The following are a few examples of common covariance functions that will be studied in this paper.
(a) Matern: This class is widely used in geostatistics, and has a fairly simple explicit form of spectral
density.
√
−(ν+1/2)

4πΓ (ν + 1/2) 2 −2ν 1
2
K̂ (ω) =
+ω
σ ρ
,
(3.2)
Γ (ν)
ρ2
in which ρ, ν, σ ∈ (0, ∞). Regardless of the choice of ρ and σ, condition (3.1) holds for Matern spectral
density with parameter ν.
(b) Powered exponential : Another versatile class of covariance functions is
!
β
r
K (r) = σ 2 exp −
ρ

(3.3)

for some β ∈ (0, 2) and ρ, σ ∈ (0, ∞). Although the spectral density does not have a closed form in
terms of simple functions, Lemma B.3 shows that K̂ admits Assumption 3.1 with ν = β/2.
(c) Rational spectral densities: Rational spectral densities form a general class admitting Assumption 3.1.
For any K̂ in this class, there are two polynomials, Qn and Qd , with real coeﬃcients, unit leading
coeﬃcients and p := deg (Qd ) − deg (Qn ) ∈ N, such that
K̂ (ω) = λ

|Qn (jω)|2
|Qd (jω)|

2.

(3.4)

Moreover, we assume that Qd has no root on the imaginary axis and λ is a strictly positive scalar. Since
K (0) < ∞ and K̂ (ω) ≍ ω −2p as ω → ∞, Assumption 3.1 holds with ν = p − 1/2.
(d) Triangular : For T, σ ∈ (0, ∞), the covariance function and spectral density are given by


 ρω  2
r
ρσ 2
2
K (r) = σ 1 −
sinc
, K̂ (ω) =
.
ρ +
2
2
Triangular covariance is less favorable than the aforestated cases due to the oscillatory behaviour of K̂
(p. 31, [53]). One can easily show that this covariance fulﬁls Assumption 3.1 with ν = 1/2.

8

Theorem 3.1. Let δ ∈ (0, 1). Suppose that G is a one dimensional Gaussian process restricted to D = [0, 1]
whose associated spectral density K̂ admits Assumption 3.1 for some ν and CK . G is regularly sampled on
i/n, i = 1, . . . , n. There exist Rn,δ > 0 (depending only on n and δ), n0 := n0 (K) and a positive universal
constant C such that if n ≥ n0 and
s




n (1 − 2α)
1
−ν
log
,
(3.5)
|b| ≥ Cn
CK 1 +
ν
δ
we have
ϕn (TGLRT ) ≤ δ.
See Appendix A.2 for the proof of Theorem 3.1. We now make several comments regarding the roles of
various quantities embedded in Theorem 3.1.
(a) Rn,δ in Theorem 3.1 can be chosen as
" 
 s 
#
2n (1 − 2α)
2n (1 − 2α)
Rn,δ = 1 + 2 log
+ log
.
δ
δ

(3.6)

We guarantee that CDEP is less than or equal δ by controlling false alarm and miss detection probabilities below δ/2. Our trick provides an evidence to choose Rn,δ . Notice that under null hypothesis, the
test statistic in Eq. (2.4) has the same distribution as the supremum of a χ21 process over Cn,α , which
is represented by {Ψ (t) : t ∈ Cn,α }. Strictly speaking for controlling the false alarm probability below
δ/2, Rn,δ needs to be chosen in such a way that


δ
P max Ψ (t) ≥ Rn,δ ≤ .
t∈Cn,α
2
The standard χ21 tail inequality in [7] implies that if Rn,δ is chosen based upon Eq. (3.6), then Ψ (t) is
below δ/ {2n (1 − 2α)} for any t ∈ Cn,α . Thus, the union bound inequality yields
!
δ |Cn,α |
δ
P sup Ψ (t) ≥ Rn,δ ≤ |Cn,α | max P (Ψ (t) ≥ Rn,δ ) ≤
= .
t∈C
2n
(1
−
2α)
2
n,α
t∈Cn,α
√
(b) The minimal detectable shift is proportional to CK , as deﬁned in (3.1). Note that CK is determined by
both low frequency and √tail behaviour of spectral density via ν. (E.g., for Matern covariance functions

given by (3.2), CK = 4πΓ(ν+1/2)
σ 2 1 ∨ ρ−2ν ). It is easily veriﬁable by (3.1) that CK is linearly
Γ(ν)
p
proportional to K (0), meaning that CK also captures the notion of the standard deviation of the
observations. Thus, Theorem 3.1 implicitly expresses that change detection is more challenging for
Gaussian processes with larger variance.
(c) Sample size n has two opposing eﬀects on the detection rate. On the one hand, n (1 − 2α) appearing in
the logarithmic function, is closely connected to the size of alternative hypothesis which is determined
by |Cn,α | = n (1 − 2α). On the other hand, the term n−ν indicates the possibility of small shift detection
as more observations are available.
We will see in Section 3.2 that parameters δ, variance of observations and sample size have almost
analogous roles in the increasing domain change detection. The main diﬀerence between the two asymptotic
settings is the role of the decay rate of K̂ in ﬁxed domain, which is encapsulated as ν. Note that ν is closely
related to the smoothness of G with larger values of ν corresponds to smoother Gaussian process in the mean
squared sense (cf. [53], Chapter 2). For smooth Gaussian processes, G (t0 ) can be interpolated using the
observations in the vicinity of t0 with small estimation error. This leads to a simpler shift-in-mean detection
for smoother processes. More precisely, as n → ∞ the lower bound on detectable b, (3.5), vanishes more
rapidly for larger ν.
9

Remark 3.2. We describe the rate of the minimal detectable jump for some speciﬁc commonly used classes
of spectral densities, all of which satisfy Assumption 3.1.
(a) Matern: For the Matern class with
p parameters (σ, ρ, ν) satisfying Assumption 3.1, the smallest detectable shift in mean is |b| ≍ n−ν log (n (1 − 2α) /δ).

(b) Powered
exponential : For the
p
 power exponential class, any jump size of magnitude at least |b| =
−β
Ω
n log (n (1 − 2α) /δ) is detectable. On the contrary to Matern class, obtaining a closed form
for CK is quite diﬃcult for powered exponential class.
−2p

(c) Rational spectral densities: It has been discussed
previously that K̂ (ω)
as ω → ∞ and

 ≍ |ω|
p
−(p−1/2)
ν = p − 1/2, revealing that each |b| = Ω n
log (n (1 − 2α) /δ) can be detected with high
probability.
p

n−1 log (n (1 − 2α) /δ)
(d) Triangular : Since K̂ satisﬁes Assumption 3.1 with ν = 1/2, so any |b| = Ω
is detectable. Although it needs a great algebraic eﬀort to ﬁnd CK , it can be shown easily that
CK ≤ σ (2/ρ + ρ/2).
We conclude this section by a comprehensive explanation the role of α in Theorem 3.1. The dependence
on α in Eq. (3.5) is logarithmic, which encoding how the size of Cn,α aﬀecting the detection rate.
Remark 3.3. The minor role of α in Eq. (3.5) may seem a bit surprising. Strictly speaking the asymptotic
behavior of the smallest detectable jump remain the same, regardless of how small α has been chosen (even
if α tends ro zero). It is also notable to mention that we did not use the assumption that α is a ﬁxed and
strictly positive scalar in our proof. This puzzling aspect can be resolved by a deeper look at the formulation
of the hypothesis testing problem (2.3). For algebraic convenience, we assume that the mean of G ﬂuctuates
around µ = 0 in Eq. (2.3). The fact that we assumed µ is known is the main reason that α parameter in the
detraction rate of GLRT, as we do not need to estimate µ from the data. That is why in this particular case α
even be chosen as small as O (1/n). We want to emphasize that the generic form of GLRT test for unknown
µ are presented in Proposition A.1. We believe that in the analysis of the extended version of GLRT, the
constant C in Eq. (3.5) depends on α (without changing the dependence on n and other parameters).
Gaussian covariance function. The Gaussian covariance function is given by
"
#
"
 2 #
√
(ρω)2
1 r
2
2
, K̂ (ω) = ρσ 2π exp −
.
K (r) = σ exp −
2 ρ
2

(3.7)

It is widely used in practice for modeling of smooth Gaussian processes, e.g. in [40]. Regarding this choice
of covariance, we have the following result:
Theorem 3.2. Let G be a Gaussian process on [0, 1] which is observed at i/n, i = 1, . . . , n, whose covariance
function is given by Eq. (3.7). Choose δ ∈ (0, 1). There are Rn,δ > 0, n0 := n0 (ρ), C0 := C0 (ρ) > 0 and a
universal constant C > 0, such that if n ≥ n0 and
s


h
i
n (1 − 2α)
,
(3.8)
|b| ≥ C exp − n log (C0 n) log
δ
then
ϕn (TGLRT ) ≤ δ.
The details of the proof is given in Appendix A.2. Because of the super-exponential decay of Gaussian
spectral density, Assumption 3.1 is actually satisﬁed for any ν > 0. This result shows that it is possible to
detect exponentially small jump size b as n increases. Moreover, the result of Theorem 3.2 is also compatible
with that of Theorem 3.1, in quantifying precisely the assertion that the smoother the Gaussian process is,
the easier it is to detect the presence of the shift in the mean.
10

3.2

Detection rate in increasing domain setting

Turning now to the increasing domain setting, recall that G is assumed to follow the setting described in
Section 2.2. We also assume that Σn = Tn (f ) for some function f satisfying the following conditions.
Assumption 3.2. f : [−π, π] 7→ R is a real symmetric function such that
(a) There are two positive universal scalars, 0 < mf ≤ Mf < ∞ such that
mf :=

inf

ω∈[−π,π]

f (ω) ≤ Mf :=

sup

f (ω) .

ω∈[−π,π]

(b) There exist positive universal constants c and λ such that
−(1+λ)

|fk | ≤ c (1 + k)

.

(3.9)

Note that the ﬁrst condition regarding the inﬁmum of f is necessary to have a positive deﬁnite inﬁnite
covariance matrix, i.e., ν ⊤ ΣN ν > 0 for any non-zero ν ∈ RN . Moreover, the polynomial decay of fk ’s as
stated in (3.9) is a suﬃcient condition to ensure that f can be equivalently expressed by its Fourier series.
Such condition is common in the non-asymptotic analysis of Toeplitz matrices (see, e.g., [22]).
Theorem 3.3. Let δ ∈ (0, 1) and suppose that Σn = Tn (f ) in which f admits Assumption 3.2 for some
positive scalars c and λ. There exist n0 ∈ N, C > 0 (depending only on c and λ) and Rn,δ > 0 such that for
any n ≥ n0 , if
s


n (1 − 2α)
|b| ≥ C f (0) n−1 log
,
(3.10)
δ
then
ϕn (TGLRT ) ≤ δ.
See Appendix A.2 for the proof of Theorem 3.3. Some comments are in order. First, the threshold Rn,δ
in Theorem 3.3 is chosen in exactly the same way as in the ﬁxed domain setting, as given by Eq. (3.6).
Second, in contrast to the ﬁxed domain setting, the dependence structure for G no longer plays the central
role in the characterization of the detection performance. In particular, f (0) is the only factor in (3.10)
that captures the correlation in the samples, but this scalar quantity evidently has an insigniﬁcant eﬀect:
the asymptotic behaviour of GLRT remains the same (up to some constant factor) for diﬀerent Gaussian
processes satisfying Assumption 3.2. A related observation that arises by comparing between (3.5) and (3.10)
is that the correlation structure of observations, which is encapsulated into ν or f (0), and the quantities
encoding the marginal density information such as n have been completely decoupled in the rate of GLRT in
increasing domain. An examination of the proof reveals that the decoupling eﬀect in the increasing domain
setting arises due to the short-range correlation assumption (cov (Xr , Xs ) → 0 polynomially in |r − s|). It
follows that as n increases the correlation for most pairs of observed sample become negligible.
Remark 3.4. Although for algebraic convenience, throughout the paper we assume that α is a ﬁxed and
positive constant, it is not necessary for our proof technique of Theorem 3.3. We basically can generalize
Theorem 3.3 to the case that α tends to zero as n → ∞. A deeper look at our proof (see A.2) and the
auxiliary results in Appendix C reveals that Eq. (3.10) can be replaced by a complicated from
s 

n (1 − 2α) p
|b| ≥ C log
f (0) − ξn ,
δ
In which ξn is a vanishing sequence of n, which can be easily eliminated from the detection rate by adjusting
the universal constant C. We refer the reader to Remark 3.3 for a detailed discussion on the role of α in the
case that µ is unknown in the formulation of the shift in mean problem ((1.1)).
11

4

Detection rate of plug-in GLRT

As we have shown in Proposition 2.1, full knowledge of Σn is central to computing the generalized likelihood
ratio. In practice, the spectral density and covariance function of G are not known a priori, and so we take
a plug-in approach, by approximating the GLRT by estimating the covariance estimate Σn (see Deﬁnition
2.1). This section serves to investigate various ways of constructing plug-in GLRT and assessing its detection
performance. We focus only on the ﬁxed domain setting, because the dependence structure underlying the
hypothesis plays an important role in determining the detection error rate, as shown in Section 3.2.
We ﬁrst assume that G is a Matern Gaussian process on D = [′, ∞] with unknown parameters η = (σ, ρ) ∈
n
Ω (see Eq. (3.2)), and is regularly observed on {i/n}i=1 . We use η̃m = (σ̃m , ρ̃m ) to indicate the estimated
parameters using m regularly spaced samples in D. We also assume that Cn,α ⊆ {k : αn ≤ k ≤ (1 − α) n}.
Namely, the Gaussian process is under control for a certain number of observations. The controlled samples
before the sudden change, XB := {Xk : k ≤ αn}, will be used to estimate η. The parameter estimation
stage is typically called the burn-in period in the literature.
It is known in the Gaussian processes literature that η is not consistently estimable in the ﬁxed domain
setting when the number of the observations in D grows to inﬁnity — see [60, 61] for further details. Zhang
[61] showed that neither σ or ρ are consistently estimable but the quantity σρ−ν can be consistently estimated
using MLE. The profound reason behind the inconsistency is the existence of a class of mutually absolutely
continuous models for G which are almost surely impossible to discern by observing one realization of G.
Strictly speaking, the induced measures corresponding to two Matern Gaussian processes with parameters η
and η ′ are absolutely continuous with respect to each other, whenever σρ−ν = σ ′ ρ′−ν . Furthermore Zhang
[61] showed that if one ﬁxes ρ at an arbitrary value, then the maximum likelihood estimator for σρ−ν is
consistent. We shall show that despite the inconsistency in estimating η, plug-in GLRT exhibits an analogous
performance as GLRT with fully known covariance function whenever the estimate of η is consistent up to
the equivalence class.
It has been discussed in that ﬁxing ρ̃m at large values has triﬂing impact on predictive performance. Note
that due to the complicated dependence of Matern covariance function to ρ, estimating ρ is a computationally
challenging task, particularly for large data sets. So we can accelerate the whole detection procedure without
estimating ρ. Plug-in GLRT change detector is a two stage algorithm as follows:
• Estimation step:
1. Fix ρ̃m at the largest possible element in Ω. Namely, ρ̃m is a deterministic quantity given by ρ̃m =
sup {ρ : (σ, ρ) ∈ Ω}.

2. Estimate σρ−ν given the controlled samples XB , using any consistent procedure such as maximum
likelihood (MLE) [61], weighted local Whittle likelihood [58], averaging quadratic variation [3]. We
P

use the term consistent to refer to the cases that |σρ−ν − σ̃m ρ̃−ν
m | → 0 as m grows to inﬁnity.
n

:= ⌊αn⌋).
3. Construct the approximated covariance matrix of X, as Σ̃n = K r−s
n , η̃m r,s=1 (here m

• Detection step:

1. Applying the GLRT by plugging Σ̃n in place of Σn into (2.4), as described in Deﬁnition 2.1.
Now we turn to state the main result of this section regarding the rate of the plug-in GLRT.
Theorem 4.1. Let δ ∈ (0, 1). Let G be Gaussian process whose associated spectral density K̂ has Matern
form with unknown parameters (σ, ρ) ∈ Ω. Given regular samples of one realization of G, there are ﬁnite
scalar C, n0 ∈ N, a non-negative sequence lim τm = 0, and threshold level Rn,δ > 0 such that for any
m→∞
n ≥ n0 ,


ϕn T̃GLRT ≤ δ + 2τm ,
whenever

12

|b| ≥ Cn

−ν

s





1
n (1 − 2α)
CK 1 +
log
.
ν
δ

(4.1)

See Appendix A.3 for the proof of Theorem 4.1.
Remark 4.1. The threshold value for the plug-in GLRT is chosen exactly same as in Theorem 3.1.
"

 s 
!#
2n (1 − 2α)
2n (1 − 2α)
Rn,δ = 1 + 2 log
.
+ log
δ
δ

(4.2)

Since {τk }k∈N is a vanishing sequence and m = ⌊αn⌋ is an increasing function of n, (δ + 2τm ) lies in
the vicinity of δ for large n. It is worthwhile to mention that C appearing in Theorem 4.1 is larger than
the previously introduced scalar in Theorem 3.1, which can be viewed as the cost of mis-specifying ρ. The
most interesting aspect of Theorem 4.1 is perhaps that if some consistent estimate of σρ−ν is available then
regardless of its rate, the plug-in GLRT has asymptotically the same rate as the GLRT with fully known
covariance function.
We conclude this section by studying the performance of plug-in GLRT when both variance and range
parameter are consistently estimable. In this case the plug-in GLRT test is constructed by replacing the
estimated parameters in the GLRT test statistic. Suppose that G has a powered exponential covariance
function, introduced in Eq. (3.3). Anderes [3] proposed a consistent estimate of covariance parameters using
empirical average of the quadratic variation of G. According to Theorem 5 of [3], unlike the Matern class,
P

both σ0 and ρ0 are consistently estimable when β ∈ (0, 1/2). Namely, |ρ − ρ̃m | ∨ |σ − σ̃m | → 0, for the
method introduced in [3]. The following result, which has a similar ﬂavor as Theorem 4.1, determines the
detection rate of plug-in GLRT for one dimensional powered exponential Gaussian processes.
Theorem 4.2. Let δ ∈ (0, 1). Let G be Gaussian with powered exponential covariance function with
unknown parameters (σ, ρ) ∈ Ω and known β ∈ (0, 1/2). Given regular samples of one realization of G, there
are ﬁnite scalar n0 ∈ N and C (which depends on the covariance parameters β, σ and ρ), a non-negative
sequence lim τm = 0, such that for any n ≥ n0 ,
m→∞



ϕn T̃GLRT ≤ δ + 2τm ,

whenever

|b| ≥ C
and

"

Rn,δ = 1 + 2 log



s

n−β log



2n (1 − 2α)
δ


n (1 − 2α)
.
δ



+

s

log



2n (1 − 2α)
δ

!#

.

Theorem 4.2 states that given a consistent estimate of η = (σ0 , ρ0 ), the plug-in GLRT procedure has
the same asymptotic behavior as GLRT with fully known parameters (see part (b) of Remark 3.2 for the
detection rate of GLRT with known σ0 and ρ0 ).

5

Detection rate of CUSUM

In this section we revisit the classical CUSUM test and present several results regarding its detection rate in
both asymptotic settings. These results should be contrasted with our earlier theorems on the performance
of the proposed exact and plug-in GLRT tests, and highlight the need for accounting for the dependence
13

structures underlying the data, especially in the ﬁxed domain setting. In the following, Theorem 5.1 introduces suﬃcient condition on |b| under which CUSUM can distinguish null and alternative hypotheses with
high probability. Theorem 5.2 studies the performance of CUSUM in the increasing domain setting.
Theorem 5.1. Let G be a Gaussian process in [0, 1] satisfying kKk1 < ∞ and K̂ ′

∞

< ∞. Moreover let

δ ∈ (0, 1), α ∈ (0, 1/2) and Cn,α = [αn, (1 − α) n]. Given n samples of one realization of G at i/n, i = 1, . . . , n,
there are Rn,δ > 0, and n0 := n0 (δ, α) such that if n ≥ n0 and
v


u
u log 2n(1−2α)
t
δ
|b| ≥ 4
,
(5.1)
α (1 − α)
then,

ϕn (TCUSUM ) ≤ δ.
We refer the reader to Appendix 4.1 for the proof of above result. The risk of ﬁxed domain-CUSUM
has been controlled from above under mild conditions on K, which holds true for all considered examples
of covariance functions in this paper. It is indeed obvious form the following inequality that K satisﬁes the
assumptions in Theorem 5.1 if a (r) := rK (r) is absolutely integrable:
K̂

′

∞

= sup
ω∈R

Z∞

a (r) e

−jωr

−∞

dr ≤

Z∞

−∞

|rK (r)| dr.

The main feature of the above theorem is the suﬃcient condition that the jump size increases (at the order of
log n at least) in order to have an upper bound guarantee on the detection error. Although we do not have a
proof that this suﬃcient condition is also necessary, our result suggests that the CUSUM test is inconsistent
in the ﬁxed domain setting: the detection error may not vanish as data sample size increases, when the jump
size is a constant. This statement is in fact veriﬁed by simulations. By contrast, we have shown earlier that
using the GLRT based approach, we can guarantee the detection error to vanish as long as the jump size is
either constant or (better yet) bounded from below by a suitable vanishing term.
Remark 5.1. Let us give a qualitative argument for the inconsistency of the CUSUM test in the ﬁxed
domain setting. Suppose that b tends to zero as n → ∞. Deﬁne
"
#
r
n
t
X
1
1X
t (n − t)
Xk −
Xk .
Ut :=
n
n−t
t
k=t+1

k=1

The expected value of Ut is zero, under the null hypothesis and for any t. Regardless of the existence of a
shift in the mean, the standard deviation of Ut remains the same. A careful look at the proof
of Theorem
√
5.1 reveals that the smallest value of the standard deviation of Ut over t ∈ Cn,α is order n. Moreover, if
there
mean occuring at the change point t̄ ∈ Cn,α , then the expected value of Ut̄ is given by
p is a shift in the √
b t̄ (n − t̄) /n = O (b n) (Recall that αn ≤ t̄ ≤ (1 − α) n). Generally speaking, as the mean of Ut under
the null hypothesis, denoted by E (Ut | H0 ), is zero for any t ∈ Cn,α , the CUSUM test cannot distinguish
between the null and the alternative (even for large sample size), since
E (Ut̄ | H1 )
p
=O
var (Ut̄ )

 √ 
b n
√
= O (b) → 0,
n

as n ր ∞.

Here, E (Ut̄ | H1 ) represents the expected value of Ut under the alternative. This suggests that regardless
of the sample size, the CUSUM test cannot detect the existence of a small shift in the mean in the ﬁxed
domain setting.

14

Remark 5.2. The threshold value of the CUSUM test in Theorem 5.1 is given by
v
s 
u


!
u
2n (1 − 2α)
2n (1 − 2α)
t
Rn,δ = n 1 + 2 log
+ 2 log
.
δ
δ

This threshold has diﬀerent form of dependence to n than that of the threshold of GLRT in Eq. (3.6), since
unlike GLRT the CUSUM test do not reduce the correlation among the samples. In order to remove the
gap between the threshold
of GLRT and CUSUM in the ﬁxed domain setting, we further normalize Ut by
√
considering Un⋆ = Un / n. So equivalently CUSUM test in this regime can be written as
!
2
R
n,δ
2
⋆
:=
.
TCUSUM = 1 max |Un⋆ | > Rn,δ
t∈Cn,α
n
⋆
Here Rn,δ
is exactly same as the critical value of GLRT test.

Now we aim to study CUSUM test in the increasing domain asymptotic regime. Recall that for this
scenario, G is a Gaussian process with covariance matrix Σn = Tn (f )
Theorem 5.2. Let δ ∈ (0, 1), ϑ > 0 and Cn,α = [αn, (1 − α) n]. Assume that f satisﬁes Assumption 3.2 for
some c and λ. There are n0 = n0 (f, ϑ) and universal constant C (λ, c) > 0, such that if n ≥ n0 and
s


(1 + ϑ) f (0)
n (1 − 2α)
|b| ≥ C
,
(5.2)
log
nα (1 − α)
δ
then
ϕn (TCUSUM ) ≤ δ.
Th reader can ﬁnd the detailed proof of Theorem 5.2 in Appendix A.4. By comparing between Theorems
3.3 and 5.2, it is clear that the CUSUM test exhibits a similar detection performance as the GLRT test in
the increasing domain setting. In fact, we will show in the next section that both tests achieve minimax
optimality in that setting. However, according to the numerical studies, using GLRT slightly improve the
detection performance comparing to CUSUM, especially in the presence of strong long range dependence.
Thus, in practice, one can aﬀord to ignore the covariance structure of G in the increasing domain setting,
and due to its simplicity, CUSUM is to be preferred.
Remark 5.3. The parameter α, which characterizes the prior knowledge on the location of possible abrupt
change, plays a similar role in both Theorems 5.1 and 5.2. The fact that α is a ﬁxed, strictly positive
p quantity,
means that the detectable jump for CUSUSM in the increasing domain setting is of order n−1 log n.
However, such restriction on α is not critical for our proof approach and the CDEP of CUSUM is still less
than δ, if condition (5.2) holds for the case of α → 0 with the sample size growing. p
For instance even if
α ≍ nβ−1 for some β ∈ (0, 1), CUSUSM is still consistent with the detection rate |b| & n−β log n.

6

Minimax lower bound on detection rate

In this section, we establish minimax lower bounds on the detectable jump in the mean of G. Theorem 6.1
shows that the obtained rate for GLRT and plug-in GLRT (Theorems 3.1 and 4.1) are nearly optimal (up
to some logarithmic term in n) for rational spectral densities. Section 6.2 demonstrates the near minimax
optimality of the GLRT and CUSUM algorithms in increasing domain setting. Before jumping to the main
result of this section, let us rigorously introduce the notion of near minimax optimality.
Definition 6.1. Given n samples, let T ∈ {0, 1} be a shift in mean detection algorithm whose CDEP
is denoted by ϕn (T ). T is said to be near minimax optimal (up to some logarithmic term in n) in the
∞
asymptotic sense, if for any δ ∈ (0, 2), there are two vanishing sequences {h1,n }∞
n=1 and {h2,n }n=1 relying on
n, δ and spectral density, such that
15

1. ϕn (T ) ≤ δ whenever |b| ≥ h1,n .
2. When |b| ≤ h2,n , there is no algorithm whose CDEP is strictly less than δ.


3. There exists a positive, bounded scalar β, for which h1,n /h2,n = O logβ n as n → ∞.

Put simply, when the sample size is large enough, no algorithm is considerably superior to a near optimal
T , regardless of how complicated its formulation might be.

6.1

Lower bound in the fixed domain

We begin this section by recalling that in the ﬁxed domain regime, G is a Gaussian process in [0, 1] which is
n
observed at {i/n}i=1 . We formally introduce the class of spectral densities that we consider in this section.
While the following conditions on K̂ are more restrictive than Assumption 3.1, it is still provides a rich class
of commonly used spectral densities.
Assumption 6.1. There are constants p ∈ N and β ∈ (1/2, ∞) such that
1. lim K̂ (ω) |ω|
ω→∞

2. lim sup
ω→∞



2p

2p
′ :=
exists and CK
lim K̂ (ω) |ω| ∈ (0, ∞).

K̂(ω)|ω|2p
′
CK

ω→∞



− 1 ω β < ∞.

Generally speaking, Assumption 6.1 contains the class of spectral densities K̂ (ω) for which there is some
−2p
p ∈ N such that K̂ (ω) ≍ |ω|
as ω tends to inﬁnity. Note that the second condition in Assumption
6.1 is of theoretical purposes and does not have a simple qualitative interpretation. It can be observed
that Assumption 6.1 excludes any K̂ (ω) satisfying Assumption 3.1 with (ν + 1/2) ∈
/ N. For instance,
Assumption6.1 does not hold for Matern covariance functions with (ν + 1/2) ∈
/ N.
Remark 6.1. Here, we name a salient class of spectral densities satisfying Assumption 6.1.
• Simple calculations show that any rational spectral density K̂ (See (3.4)) admits Assumption 6.1 with
′
CK
= λ, β = 1 and p = deg (Qd ) − deg (Qn ) ∈ N. Moreover, K̂ satisﬁes Assumption 3.1 with
ν = p − 1/2. We discussed in Remark that Matern covariance function with p := (ν + 1/2) ∈ N
has indeed a rational spectral density. These particular instances of Matern covariance, which are
commonly used in machine learning and geostatistics, are of the form K (r) = Q (|r|) e−d|r| , where
Q (·) is a polynomial of degree p − 1.
Theorem 6.1. Let
and assume that Assumption 6.1 holds for K. Consider the problem (2.3) in
2)
 δ ∈ (0,
n
. There are positive scalars C̄K and n0 := n0 (K) such that if n ≥ n0 and
which cov (X) = K r−s
n
r,s=1
|b| ≤ C̄K n

−p+1/2

s

log




1
,
δ (2 − δ)

(6.1)

then for any test T ,
ϕn (T ) ≥ δ.
See Appendix A.5 for the proof of Theorem 6.1.
Remark 6.2. Comparing the detection rate of GLRT (see Theorem 3.1) and plug-in GLRT (recall from
Theorem 4.1), with the rate described in Eq. (6.1) establishes near minimax optimality of the GLRT with
known covariance structure and plug-in GLRT for the class of spectral densities considered in Remark
√ 6.1,
in the asymptotic sense. Strictly speaking, under the ﬁxed domain setting, there is a gap of order log n
between (6.1) and the detection rate of formerly studied GLRT based algorithms. Although we do not have
a proof to establish the near minimax optimality of the GLRT and plug-in GLRT for spectral densities
satisfying Assumption 3.1, our conjecture is that Theorem 6.1 can be extended to this broader class.
16

6.2

Lower bound in the increasing domain setting

Turning to the increasing domain setting, we give a condition on jump size |b| according to which no algorithm
in the increasing domain can properly detect the existence of a shift in the mean. Unlike Section 6.1, there
is no distinction between the assumptions used to obtain the minimax lower bound and Assumption 3.2.
Theorem 6.2. Let δ ∈ (0, 2), ϑ > 0 and Cn,α = [αn, (1 − α) n]. Suppose that Σn = Tn (f ) in which f
satisﬁes Assumption 3.2. There are n0 := n0 (f, ϑ) and a universal constant C > 0 such that if n ≥ n0 and
v


u
u (1 + ϑ) f (0) log
1
t
δ(2−δ)
|b| ≤ C
,
αn

then for any test T ,

ϕn (T ) ≥ δ.
The direct comparison between the detection rate of both CUSUM (in Theorem 5.2) and GLRT (see
Theorem 3.3) test with the above result indicates the minimax optimality (up to some order log n term) of
both os these procedures in the increasing domain setting.

7

Simulation study

To illustrate the performance of the proposed shift-in-mean detection algorithms, we conduct a set of controlled simulation studies for verifying the results in Sections 3, 4 and 5. Our goals are two-fold:
(a) comparing the performance of the GLRT based algorithms with the standard CUSUM test.
(b) assessing the sensitivity of algorithm (2.5) to the parameters of the covariance function and tapering of
Σn .
In all the numerical studies in this section we ﬁx n = 500 and α = 0.1.
The area under the receiver operating characteristic (ROC) curve, which will be referred as AUC, is a
standard way for assessing the performance of a test. The ROC curve plots the power against the false alarm
probability. Since the ROC curve is conﬁned in the unit square, the AUC ranges in [0,1]. The ROC curve of
a test based on pure random guessing is the diagonal line between origin and (1, 1) and so the AUC of any
realistic test is at least 0.5.
The subsequent ﬁgures in this section exhibit empirical AUC versus b. For a ﬁxed value of b, covariance
function K and a detection algorithm T , we apply the following method to compute the AUC of T :
1. Set T1 = 500 and T2 = 50.
2. For k = 1 to T2 repeat independently
(a) For ℓ = 1 to T1 repeat independently
i. Choose p ∈ {0, 1} with equal probability which denotes null or alternative hypotheses. Thus,
approximately T1 /2 = 250 experiments correspond to both null and alternative.
ii. If p = 0, generate zero mean X ∈ Rn according to covariance function K. That is, X are sampled
from a Gaussian process with no abrupt shift in mean. Otherwise, choose t ∈ [αn, (1 − α) n] =
{50, 51, · · · , 450} uniformly at random (recall that t represents the location of the mean shift)
and generate X ∈ Rn according to H1,t .
iii. Compute T score.
(b) Numerically obtain the ROC curve of T based upon T1 experiments in part i.
(c) Given the ROC curve, compute AU Ck using trapezoidal integration method.
17

3. Compute the average AUC by AU C =

1
T2

T2
P

AU Ck .

k=1

The ﬁrst simulation study aims to compare CUSUM and GLRT based algorithms in the ﬁxed domain
regime and assess the role of smoothness and other parameters of K in the performance of the GLRT. For this
n
experiment G is a Gaussian process in [0, 1] which is observed at regularly spaced samples, Dn = {k/n}k=1 ,
i.e., Xk = G (k/n) , k = 1, . . . , n . The covariance function of G is assumed to has Matern form with
parameters (σ0 , ρ0 , ν). Strictly speaking,


i−l
cov (Xi , Xl ) = σ02 Kν
, i, l = 1, . . . , n,
nρ0
√
Z
−(ν+1/2)
4πΓ (ν + 1/2) ∞ −jωx
dω, ∀x ≥ 0.
1 + ω2
e
Kν (x) =
Γ (ν)
−∞

We consider three diﬀerent scenarios on ν, 0.5, 1, and 1.5. We also set σ0 = 1 and ρ0 = 1/2. As customary
in the literature, we assume that ν is known and so ν will not be estimated. For conducting the plug-in
GLRT procedure, both parameters (σ0 , ρ0 ) are estimated using full MLE. Due to the low dimensionality of
unknown parameters, the most eﬀective way to estimate (σ0 , ρ0 ) is to apply brute force grid search over a
pre-speciﬁed set P. Here, we choose P = {0.2, 0.4, · · · , 2} × {1/4, 1/3.9, . . . , 1/0.1}. The ﬁnal results of this
numerical study is exhibited in Figure 1. We observe the following:

• GLRT and plug-in GLRT have a signiﬁcantly better detection performance than CUSUM. This performance improvement is more pronounced for smoother covariance function (larger ν). In particular,
the CUSUM is completely impractical for detection of an small change when ν = 1 or 1.5.
• In each panel of Figure 1, the GLRT has a slightly larger AUC than that of plug-in GLRT. Thus,
Figure 1 veriﬁes the existence of a small gap between the smallest detectable jump of GLRT and plugin GLRT. Note that this fact can also be observed by comparing (3.5) and (4.1). In short, having full
knowledge of covariance parameters slightly improves the detection performance and so our proposed
algorithm is robust to the estimation error of the unknown parameters of K.
• Comparing the range of b in each panel of Figure 1 discloses that more rapid decay of the spectral
density can decrease the smallest detectable jump. This observation substantiates the role of ν in the
theory established in Sections 3 and 4.
Next, we compare the performance of the GLRT with known parameters and the CUSUM in the increasing
domain setting. Recalling from Theorems 3.3 and 5.2, these two methods have analogous asymptotic rates.
In the left panel of Figure 2, we choose an exponentially decaying covariance function


|i − l|
, i, l = 1, . . . , n,
cov (Xi , Xl ) = f|i−l| = σ02 exp −
ρ0
in which σ0 = 1 and ρ0 = 2. That is Σn has exponentially decaying oﬀ-diagonal entries. However, in the
right panel, the chosen covariance function has a polynomially decaying tail given by
cov (Xi , Xj ) = f|i−j| =

σ02


−(1+λ)
|i − l|
1+
,
ρ0

with σ0 = 1, ρ0 = 2 and λ = 0.5. In this case, Σn has heavier oﬀ-diagonal terms. Note that Assumption
3.2 is satisﬁed in either of the two cases. It is evident from Figure 2 that the GLRT exhibits a slightly
better performance than the CUSUM, and the gap between the two AUC curves is more visible in the case
of polynomially decaying covariance function. Thus, we still recommend the use of GLRT in the presence of
strong correlation among samples in applications described by the increasing domain regime.

18

1.0

1.0
0.9

CUSUM
GLRT
PGLRT

0.8

AU C

AU C

0.9

0.7
0.6
0.5

CUSUM
GLRT
PGLRT

0.8
0.7
0.6

0.1

0.2

0.3
b

0.4

0.5

0.6

0.5

0.01

0.02

0.03

0.04

0.05

b

1.0

AU C

0.9
0.8

CUSUM
GLRT
PGLRT

0.7
0.6
0.5

0.005 0.010 0.015 0.020 0.025 0.030
b

Figure 1: The above ﬁgures assess the performance of diﬀerent detection algorithms when G is one dimensional Matern Gaussian process, with parameters (ν, σ0 , ρ0 ), and regularly sampled in [0, 1]. From left to
right then from top to bottom, (ν, σ0 , ρ0 ) = (0.5, 1, 0.5) , (1, 1, 0.5) , (1.5, 1, 0.5). In each panel horizontal axis
displays jump value b and the three curves (dashed black, solid blue and green) respectively exhibit AUC of
GLRT with known covariance structure, plug in GLRT (PGLRT) using full MLE and CUSUM.

19

1.0

0.9

0.9

0.8

0.8

AU C

AU C

1.0

0.7
0.6
0.5
0.0

0.6

CUSUM
GLRT
0.2

0.4
b

0.6

0.7

0.8

0.5
0.0

CUSUM
GLRT
0.2

0.4
b

0.6

0.8

Figure 2: The above ﬁgure assesses the performance of increasing domain detection algorithms. In each panel
the horizontal axis displays the jump value b and the two curves (dashed black and solid blue) respectively
exhibit the AUC of the GLRT with known covariance structure and CUSUM. In the right panel, we choose
−(1+λ)
cov (Xi , Xl ) = σ02 (1 + |i − l| /ρ0 )
in which (σ0 , ρ0 ) = (1, 2) and λ = 0.5. For the left panel, the
covariance function is given by cov (Xi , Xl ) = σ02 exp (− |i − l| /ρ0 ) where (σ0 , ρ0 ) = (1, 2).

8

Discussion

As indicated in the Introduction, the comprehensive analysis of the detection of shift-in-mean of a Gaussian
process in the ﬁxed domain regime has remained relatively unexplored. However, the considered model in
(1.1) is only one of several plausible scenarios which should be subject to more thorough investigation. We
note that the probabilistic model of G can be extended in some possible ways for future research.
(a) Here we deal with a single abrupt change in EG. However, we believe that our techniques can be
extended to rigorously formulate the minimax optimal rate of the GLRT and plug-in GLRT for detecting
multiple shifts in EG.
(b) Recently, Ivanoﬀ et al. [30] studied the problem of change-set detection in two dimensional Poisson
processes. Speciﬁcally, G is a Poisson process in R2 and there are scalars µ0 6= µ1 and Ω ⊂ R2 such
that the intensity of G can be formulated by
EG (s) = µ0 1s∈Ω + µ1 1s∈Ω
/ ,

∈ R2 .

(8.1)

The objective is to detect Ω as well as possible based on observation of G. However, a comprehensive study of minimax optimal change-set detection methods for multi-dimensional Gaussian processes
remains unavailable. Note that the ﬁxed domain setting is the natural way to study asymptotic behaviour of algorithms regarding spatial processes. So, this paper can provide valuable intuition about
the minimax rate of change-set detection in Gaussian spatial processes. We expect that aside from the
sample size and smoothness of the covariance functions, the geometric properties of change-sets will
have a crucial role in the design and analysis of detection algorithms.

Appendices
Appendix A contains the proofs of the main results in Sections 2-6. Appendix B states and proves the
technical results required in the proofs of the main results. In addition, Appendix C establishes properties
20

of the inverses of large Toeplitz matrices. Such results are useful for the proofs of Theorems 3.3, 5.2 and 6.2,
and may also be of independent interest.

A
A.1

Proofs
Proofs for Section 2

Proof of Proposition 2.1. In the following L stands for the generalized negative log-likelihood ratio.
"
⊤

#
b
b
−1
−1
⊤
X − ζt
.
(Σn )
2L = X (Σn ) X − min min X − ζt
t∈Cn,α b6=0
2
2

(A.1)

Note that the objective function in (A.1) is quadratic in terms of b. The explicit form of 2L can be obtained
with a bit of algebraic derivations. The algebra has been skipped to save space; we arrive at
−1

ζ ⊤ (Σn )
2L = max max − t
t∈Cn,α b6=0
4

ζt

b2 + bζt⊤ (Σn )

−1

X

!

2

−1

ζ ⊤ (Σn ) X
= max qt
t∈Cn,α
−1
ζt⊤ (Σn ) ζt

.

So, there is a threshold value, Rn,δ > 0, for which the GLRT is given by (2.4).

The following result expressing the form of GLRT in the generic case of unknown µ can be proved in an
analogous way as Proposition 2.1.
Proposition A.1. There is Rn,δ > 0 for which the GLRT is given by


2
hY , ζt − B1 (t) 1n i
p
≥ Rn,δ  ,
TGLRT = I  max
t∈Cn,α
B2 (t)
where Y = (Σn )

−1

X and

B1 (t) =

A.2

(A.2)

−1
ζt⊤ (Σn )
⊤ (Σ )−1
n
n

1



−1

ζt⊤ (Σn )

1n

2

1n
, B2 (t) = ζt⊤ (Σn )−1 ζt −
.
−1
1n
1⊤
1n
n (Σn )

Proofs for Section 3

Proof of Theorem 3.1. Let p = ⌈ν + 1/2⌉, P = {1, . . . , p} and θn = exp (−1/n). Construct a banded
triangular matrix An ∈ Rn×n by the following procedure.
 
p
j
An [k, k − j] =
(−θn ) , j ∈ {0, . . . , p} , k ∈ {p + 1, . . . , n} ,
j
(An )P,P = n−2ν Ip .

⊤

ζ (Σ )
It is relatively simple to verify that An is invertible. In addition, for brevity let Zt = √t ⊤ n
n

−1

X

ζt (Σn )−1 ζt

for

any t ∈ Cn,α , in which ζt has been deﬁned in (2.2). Lastly, deﬁne Un,t := An ζt ∈ R , W := An X and
Dn := cov (W ).
Easy calculations show that under the null hypothesis {Zt }t∈Cn,α is a set of standard Gaussian random

variables and so, by Lemma B.1, we have P maxt∈Cn,α Zt2 ≥ Rn,δ ≤ δ/2. That is, the false alarm probability
21


is less than δ/2. Moreover if the alternative hypothesis H1,t̃ (for some t̃ ∈ Cn,α ) holds then Zt2
non-central

χ21

Zt̃2

random variables and the non-centrality parameter of
 |b| q ⊤
−1
E Zt̃ | H1,t̃ =
ζt̃ (Σn ) ζt̃ .
2

t∈Cn,α

are

is given by

Applying Lemma B.1 (σ0 = σk = 1 for any k) demonstrates that ϕn (T2 ) ≤ δ, whenever
s  
q
q
4n
−1
−1
ζt⊤ (Σn ) ζt ≥ 8 log
.
|b| ζt̃⊤ (Σn ) ζt̃ ≥ |b| min
t∈Cn,α
δ

(A.3)

Thus, in order to get a suﬃcient condition on detectable b, it suﬃces to ﬁnd a tight uniform lower bound on
ζt⊤ (Σn )−1 ζt for t ∈ Cn,α .
−1
⊤
The identity Σ−1
An can be shown using the linearity of covariance operator and nonn = An (Dn )
singularity of A. Choose t ∈ Cn,α in an arbitrary way. As a result of this alternative representation of
−1
−1
⊤
⊤
(Dn ) Un,t . Applying Kantorovich inequality (cf. Appendix B) and the
Σ−1
ζt = Un,t
n , we have ζt (Σn )
triangle inequality yields
#2
"
kUn,t k4ℓ2
kUn,t k2ℓ2
1
−1
−1
⊤
⊤
.
(A.4)
ζt (Σn ) ζt = Un,t (Dn ) Un,t ≥ ⊤
≥
kUn,t kℓ1
kDn kℓ∞
Un,t Dn Un,t
Now, we show that

kUn,t k2ℓ
2
kUn,t kℓ

1

2

kUn,t kℓ2

≥ 31 , for large enough n. Indeed, after some algebra, we can get

t+p
X

≥

2
Un,t
(k) =

p
X

k=1

k=t+1





− (1 − θn )p + 2
2

k−1  
X
p
j=0

j

2

j
(−θn ) 

2



p 
p
k−1
X
X
X p
p−1
2 (p − 1)
k−1
j

=2
(−1)
=2
≥ 2p ,
2
(−1)
k
−
1
p
−
1
j
j=0

(a)

≥

k=1

(A.5)

k=1

where inequality (a) follows from the fact that for large enough n, θn is arbitrarily close to 1. To get an
upper bound on kUn,t kℓ1 ,
kUn,t kℓ1

=

p
X

k=1

=

p
X

|Un,t (k)| +
n

−2ν

k=1

+

t
X

k=p+1

t
X

k=p+1

≤ pn−2ν

|Un,t (k)| +
p

(1 − θn ) +

n
X

k=t+p+1

n
X

k=t+p+1

|Un,t (k)| +
p

(1 − θn ) +

t+p
X

k=t+1

p
X

k=1

|Un,t (k)|
p

− (1 − θn ) + 2

k−1  
X
p
j=0

j

(−θn )

j

p k−1
p k−1
X
X p
X
X p
(b)
j
j
1−p
(−θn )
(−θn ) ≤ 2 + 2
+n
+2
j
j
j=0
j=0
k=1

k=1


p 
p k−1
X
X p
X
p−1
k−1
j
(−1)
= 2 + 2p+1 ≤ 3 2p .
(−1) = 2 + 4
≤ 2+4
k
−
1
j
j=0
k=1

(A.6)

k=1

Note that inequality (b) is valid when pn−2ν +n1−p ≤ 2, which obviously holds for suﬃciently large n = O (1).
The remaining inequalities and identities in (A.6) can be easily veriﬁed via basic properties of the binomial
coeﬃcients. Combining (A.5) and (A.6) yields the desired goal. Now, inequality (A.4) can be rewritten as

−1
1
−1
⊤
ζt (Σn ) ζt ≥
= 9 max var (Wk )
.
(A.7)
1≤k≤n
9 kDn kℓ∞
22

In the ﬁnal phase of the proof, we achieve a tight upper bound on max1≤k≤n var (Wk ). It is obvious from
the formulation of An and the stationarity of X − EX that max1≤k≤n var (Wk ) = n−2ν ∨ var (Wp+1 ). So,
the goal is reduced to give an upper bound on the variance of Wp+1 .
!
2

Z
p  
p  
X
X
p
p
−jrω
(c) 1
r
r
(−θn ) Xp+1−r =
var (Wp+1 ) = var
dω
(−θn ) exp
K̂ (ω)
r
2π
n
r
r=0
r=0
R

=

r


p
− (1 + jω)
− exp
n
r

2

1
2π

K̂ (ω)

1
2π

Z

(d) C
p
K
K̂ (ω) 1 + θn2 − 2θn cos (ω/n) dω ≤
2π

r=0

R

=

p 
X

Z

R



dω =

1
2π

Z

K̂ (ω) 1 − e

−(1+jω)
n

2p

dω

R

Z 
1 + θn2 − 2θn cos
R

(1 + ω 2 )

ω
n

ν+1/2

p

dω,

(A.8)

where, identity (c) is implied by Bochner theorem (cf. [53], Chapter 2) and (d) is immediate consequence of
Assumption 3.1. Notice that

 ω 
hω
 ω i2

 ω 
ω 
1
1
2
sinc
.
≤ (1 − θn ) + 2θn 1 − cos
≤ 2 + 2 1 − cos
= 2+
1 + θn2 − 2θn cos
n
n
n
n
n
n
2n

Let ξ = p − (ν + 1/2) < 1. Henceforth, for any R > 0,
n
n
op
2 op
ω
1
ω
Z 1/n2 + ω sinc ω 2
Z
sinc 2n
2 +
n
n
2
2πn2ν
dω =
dω
var (Wp+1 ) ≤ n2ν
CK
(1 + ω 2 )ν+1/2
(1/n2 + ω 2 )ν+1/2
R
R
n
n
op

 op
ω 2
2
Z
ZR 1/n2 + ω sinc ω 2
1/n
+
ω
sinc
2
2
dω +
dω
=
ν+1/2
ν+1/2
2
2
2
2
(1/n + ω )
(1/n + ω )
−R
|ω|≥R
n

 op
ω 2
2
Z
ZR
1/n
+
ω
sinc
(e)
ξ
2
≤
1/n2 + ω 2 dω +
dω
ν+1/2
2
2
(1/n + ω )
−R

(f )

≤

ZR

−R

|ω|≥R

1/n2 + ω 2

ξ

dω + 5p

Z

|ω|

−(2ν+1)

(g)

dω ≤ 3R3 + 5p

R−2ν
.
ν

(A.9)

|ω|≥R

Inequality (e) follows form the fact that supω∈R |sinc (ω/2)| ≤ 1. In order to justify (f ), observe that
|ω sinc (ω/2)| ≤ 2 for any ω ∈ R. Thus, for large enough n and |ω| ≥ R, we get
n

2 op
1/n2 + ω sinc ω2
p
−(2ν+1)
−(2ν+1)
≤ |ω|
1/n2 + 4 ≤ 5p |ω|
.
ν+1/2
(1/n2 + ω 2 )
ξ
Note that there is some n0 := n0 (R, ν) such that supω∈R 1/n2 + ω 2 ≤ 3/2R2 for all n > n0 . This
immediately entails inequality (g).
Finally, minimizing the obtained upper bound in (A.9) over R > 0, we get


1
−2ν
(A.10)
1+
var (Wp+1 ) ≤ CCK n
ν
for some universal constant C > 0. Thus, there is another strictly
positive universal constant, C ′ , for which

1
−2ν
′
−2ν
max1≤k≤n var (Wk ) = n
∨ var (Wp+1 ) ≤ C CK n
1 + ν . So, (A.7) implies that
−1

ζt⊤ (Σn )

ζt &

23

CK

n2ν
.
1 + ν1

(A.11)

The combination of (A.3) and (A.11) completes our proof.
Proof of Theorem 3.2. The proof proceeds in a similar manner as that
 of2the
 preceding theorem, in the sense
ρ
that it is required to show that inequality (A.3) holds. Let θn = exp − n2 . An represents the inverse of the
Cholesky factorization of Σn . For any k ≤ j and q ∈ [0, 1], G (k, j; q) denotes the following rational function.
j
Y

G (k, j; q) =

ℓ=j−k+1

1−q

" k
Y

ℓ

ℓ=1

1−q

#−1

ℓ

,


and G (k, j; 1) = kj . G (k, j; q) is usually referred to Gaussian binomial coeﬃcients in the combinatorics
literature. Finally, let Un,t := An ζt . Similar to (A.3), the aim is to obtain a universal lower bound on
−1
−1
2
ζt⊤ (Σn ) ζt for t ∈ Cn,α . Observe that, ζt⊤ (Σn ) ζt = kUn,t kℓ2 .
In order to achieve a tight lower bound on kUn,t kℓ2 , it is pivotal to study the non-asymptotic behaviour
of the entries of An . According to Proposition 1 of [40], the entries of An are given by
 p (j−k) G (k − 1, j − 1; θ )
n
s
1{j≥k} .
(An )jk = − θn
j−1
Q
(1 − θnℓ )
ℓ=1

Since

ℓρ2
n2

1−e−x
x
xց0

tends to 0 as n gets large for any ℓ ∈ {0, . . . , n} and lim
"j−1
Y
ℓ=1

1−

θnℓ

#−1


= 1, we get

"j−1 
 2(j−1)
#−1

Y
n
1
ℓρ2
.
≍
=
1 − exp − 2
n
(j − 1)! ρ

(A.12)

ℓ=1

Direct calculations show that G (k − 1, j − 1; θn ) ≍
j, k ∈ {1, . . . , n}. Thus,

j−1
k−1



for any θn in a small neighborhood of 1 and



 p (j−k)
(j−k) j − 1
.
− θn
G (k − 1, j − 1; θn ) ≍ (−1)
k−1

(A.13)

2

The asymptotic identities (A.12) and (A.13) come in handy to analyze kUn kℓ2 :
kUn k2ℓ2

≥

n
X

(Un )2j

=

n
X

j=t+1

j=t+1

"

j
X

k=t+1

(An )jk −

t
X

k=1

(An )jk

#2

 2(j−1) " X


#2
 X
j
t
n
(j−k) j − 1
(j−k) j − 1
≍
−
(−1)
(−1)
ρ
k−1
k−1
j=t+1
k=1
k=t+1
"

#2
 2(j−1) X


j
t
n
X
X
n
1
(j−k) j − 1
(j−k) j − 1
−2
(−1)
(−1)
=
k−1
(j − 1)! ρ
k−1
j=t+1
k=1
k=1


2
 2(j−1) 
n
n
j−1 2  2(j−1)
X
X
1
n
n
j j−1
t
=
0 − 2 (−1)
≍
.
t
(j − 1)! ρ
(j − 1)! ρ
j=t+1
j=t+1
n
X

1
(j − 1)!

Thus, there are universal constants C, C ′ > 0 and C0 depending on α and ρ such that
2
kUn kℓ2

≥C

n
X

j=t+1


j−1 2
t


 2(j−1)
n−1 2  2 (n−1) (a)
 t 1  en2 (n−1) (b)
n
n
n
′ n
t
√
≥C
≥ C
≥ (C0 n) .
(j − 1)! ρ
(n − 1)! ρ2
t
n nρ2
24

Note that inequality (a) can be shown using Stirling’s formula and (b) is obvious implication of the fact that
t ≤ (1 − α) n (Recall Cn,α from Section 2.3). In summary, we have that
q
−1
n/2
|b| ζn⊤ (Σn ) ζn & |b| (C0 n)
.
We conclude the proof by appealing to Lemma B.1.

Proof of Theorem 3.3. The proof is also similar to that of Theorem 3.1. By applying Lemma B.1, we need
to show (A.3) holds. That is,
s 

q
4n (1 − 2α)
−1
.
(A.14)
|b| min
ζt⊤ (Σn ) ζt ≥ 8 log
t∈Cn,α
δ
−1

So, the suﬃcient detectability condition will be obtained by ﬁnding a tight uniform lower bound on ζt⊤ (Σn ) ζt
on Cn,α . For any t ∈ Cn,α , deﬁne at , a′t ∈ {0, 1}n by at (i) = 1{i≤t} and a′t = 1n − at . Observe that at , a′t
have non-overlapping support and ζt = at − a′t . Thus,
−1

ζt⊤ (Σn )

ζt = a⊤
t (Σn )

−1

at + a′⊤
t (Σn )

−1

−1

a′t − 2a⊤
t (Σn )

a′t ≥ a⊤
t (Σn )

−1

−1

at − 2a⊤
t (Σn )

a′t

The last inequality leads to the following key result
min

t∈Cn,α

−1

ζt⊤ (Σn )
n

ζt

≥
(b)

≥

min

t∈Cn,α

−1

a⊤
t (Σn )
n

at

− 2 max

t∈Cn,α

a⊤
t (Σn )
n

−1

a′t

(a)

≥ min

t∈Cn,α

(c)
1
1
− ξn′ − ξn ≥
,
f (0)
2f (0)

−1

a⊤
t (Σn )
n

at

− ξn
(A.15)

′ ∞
in which {ξn }∞
n=1 and {ξn }n=1 are appropriately chosen non-negative vanishing sequences, based on the
developed results in Appendix C. In (A.15), inequality (a) and the explicit form ξn are obtained from
Lemma C.1. One can ﬁnd a closed form of ξn′ and veriﬁes inequality (b) using Corollary C.1. Furthermore,
(c) holds whenever n is greater than some n0 , which depends on c and λ. The proof of Theorem 3.3 is
completed by combining (A.14) and (A.15).

A.3

Proofs for Section 4

Proof of Theorem 4.1. Recall dK,∞ from (??) and deﬁne the event Am by Am := [dK,∞ (η̃m , Ωη ) ≤ ρm ].
−1
ζ ⊤ (Σ̃n ) X
Moreover, we use Z̃t to denote qt ⊤
. Notice that η̃m , Σ̃n and are measurable functions on the
−1
ζt (Σ̃n ) ζt
sample space. Furthermore for any u ∈ Am , η̃m (u) and Σ̃n (u) represent the value of these measurable
functions at u. Lastly, for u ∈ Am deﬁne the random variable Z̃t (u) by

−1
ζt⊤ Σ̃n (u)
X
Z̃t (u) = r 
−1 .
⊤
ζt
ζt Σ̃n (u)

As the range of ϕn is [0, 2] (See Deﬁnition 3.1), we have by Assumption ?? that



 

ϕn T̃GLRT ≤ 2τm + E ϕn T̃GLRT | Am .

(A.16)

As lim supm→∞ ρm < 1, there are γ ∈ (0, 1) and m0 ∈ N such that ρm ≤ γ for any m ≥ m0 . Choose
γ < γ0 < 1 in an arbitrary fashion. We aim to obtain a suﬃcient condition on b to control the second
term in the right hand side of (A.16) below δ. Notice that throughout the proof we assume that m ≥ m0 .
Conditioning on the occurrence of Am , the following statement trivially holds for m.
25

∃ η̄m ∈ Ωη s.t.

K̂ (·, η̃m )
K̂ (·, η̄m )

−1

∞

≤

ρm
< 1.
γ0

(A.17)

Now choose
Deﬁne the covariance matrix of observations associated to η̄m (u) by
n
 u ∈ Am arbitrarily.
,
η̄
(u)
.
Similar
to the proof of Theorem 3.1, it is necessary to study the two
Σ̄n (u) = K r−s
m
n
r,s=1

following quantities: 1. variance of Z̃ (u) and 2. expected value of Z̃ (u) under the alternative hypothesis,
to control the false alarm and miss detection probabilities. Notice that

σt (u) :=
(b)

≤



−1
−1

−1

−1
ζt⊤ Σ̃n (u)
Σn Σ̃n (u)
ζt (a) ζt⊤ Σ̃n (u)
Σ̄n (u) Σ̃n (u)
ζt
var Z̃ (u) =
≤ B
−1
−1


ζt
ζt
ζt⊤ Σ̃n (u)
ζt⊤ Σ̃n (u)

−1
ρm
B 1−
.
(A.18)
γ0

Lemma B.4 ensures the existence of some scalar B ∈ (1, ∞) for which inequality (a) in (A.18) holds (since η
and η̄m (u) belong to the same equivalence class). Moreover, (b) can be easily deduced from the combination
−1

. Now, using
of (A.17) and the second inequality in Lemma B.2. Namely, σt (u) ≤ σ0 := B 1 − ργm0
Lemma B.1, we get


δ
2
P max Z̃t ≥ Rn,m | Am ≤ .
(A.19)
1≤t≤n
2
h

 q
i
2n
+
. So, we have controlled
log
Note that Lemma B.2 suggests to take Rn,δ,ρm = σ0 1 + 2 log 2n
δ
δ
type 1 error from above in (A.19). Now we turn to control the type II error from above. Assume that there
is a sudden change in the mean of G at t̄ ∈ Cn,α . According to Lemma B.1, type II error is less than δ/2
whenever for any u ∈ Am
s  
s  
r 
−1
2n
2n
|b|
4B
 log
ζt̄ ≥ 4σ0 log
ζt̄⊤ Σ̃n (u)
=
.
(A.20)
ρm
2
δ
δ
1−
γ0

Applying Lemma B.2 and then Lemma B.4, one can easily show that



−1
1 − ργm0 ⊤
−1
−1
ρm
⊤
⊤
ζt̄ Σ̃n (u)
ζt̄ ≥ 1 −
ζt̄ Σ̄n (u)
ζt̄ ≥
ζt̄ Σ̄n (u)
ζt̄ .
γ0
B

(A.21)

The combination of last inequality and (A.20) along one line of algebra leads to the following suﬃcient
condition to control type II error holds if
B
1 − ργm0

q
−1
|b| min ζt⊤ (Σn ) ζt ≥ 8
t

!3/2 s

log




2n
,
δ

We conclude the proof by invoking (A.11) to obtain (A.22).

A.4

Proofs for Section 5

Proof of Theorem 5.1. Choose t ∈ Cn,α and deﬁne
Ut⋆

:=

r

t (n − t)
n2

n
n
X
1
1X
Xk −
Xk
n−t
t
k=t+1

26

k=1

!

.

(A.22)

Moreover set
Rn,δ

v
s 
u


!
u
2n
(1
−
2α)
2n (1 − 2α)
t
= n 1 + 2 log
+ 2 log
.
δ
δ

Note that under the null hypothesis, Ut⋆ is a zero mean random variable and
lim var (Ut⋆ )
n→∞

=

t (n − t)
lim
n→∞
n2

=

Z∞

(a)

lim

n→∞
−∞

Z∞

(a)

=

−∞

Z∞

−∞

K̂ (ω)
2π

K̂ (ω)
2π

s

2

n
n
X
K̂ (ω)
1
1X
exp (−jkω/n) −
exp (−jkω/n) dω
2π n − t
t

s

k=t+1

k=1

n
X
exp (−jkω/n)
β
−
1−β
n

β
1−β

k=t+1

Z1

e

−jωu

du −

β

s

1−β
β

Zβ

e

s

n
1 − β X exp (−jkω/n)
β
n

dω

k=1

2

−jωu

2

du dω =

Z∞

−∞

0

K̂ (ω) Gβ (ω)
dω,
2π

where
2 
2





ω
(1 − β) ω
(1 − β) ω
βω
Gβ (ω) := (1 − β) sinc
sinc
sin2
.
+ β sinc
+4β (1 − β) sinc
2
2
2
2
(A.23)
The identity (a) is implied by Bochner Theorem and (b) follows from the dominated convergence theorem.
It is easy to see that kGβ k∞ ≤ 1 and so lim var (Un⋆ ) ≤ 1 by the triangle inequality. Moreover, Lemma




βω
2

n→∞

B.8 shows that the achieved upper bound on σn2 = var (Un⋆ ) is tight up to some constant whenever K̂ has a
uniformly bounded derivative. Namely, there is a universal constant c ∈ (0, 1) such that c ≤ lim var (Un⋆ ) ≤ 1
n→∞

⋆
2
for any β ∈ (0, 1). Let Rn,δ
= Rn,δ
/n. Thus

P (T = 1 | H0 ) = P



max |Ut | ≥ Rn,δ | H0

t∈Cn,α



=P



max

t∈Cn,α

2
|Ut⋆ |

≥

⋆
Rn,δ



| H0 .

(A.24)

2

For any t ∈ Cn,α , |Ut⋆ | is a (non-normalized) χ21 random variable, as σn2 ≤ 1. Moreover |Cn,α | = n (1 − 2α).
So the part (a) of Lemma B.1 says that


δ
2
⋆
P max |Ut⋆ | ≥ Rn,δ
| H0 ≤ .
t∈Cn,α
2
Now we turn to control the miss detection probability. Without loss of generality assume that b > 0. Choose
an arbitrary t ∈ Cn,α . A line of algebra shows that
p
E (Ut⋆ | H1,t ) ≥ b α (1 − α).
(A.25)
p
Eq. (5.1) on b implies that E (Ut⋆ | H1,t ) ≥ 4 log (2n (1 − 2α) /δ). In other words, given a sudden jump at
2
t, |Us⋆ | , s ∈ Cn,α are non-central χ21 random variables satisfying the conditions of the part (b) of Lemma
B.1. Hence


δ
2
⋆
(A.26)
P (T = 0 | H1,t ) = P max |Us⋆ | ≤ Rn,δ
| H1,t ≤ .
s∈Cn,α
2
Proof of Theorem 5.2. We continue to use the samenotation as the proof of Theorem 5.1. Note that there
are three appropriately chosen vanishing sequences a1n n∈N , a2n n∈N and a3n n∈N such the
27

var (Ut )

=

t
var
n (n − t)

n
X

Xk

!

Xk

!

k=t+1
n
X

t
X
n−t
+
var
Xk
nt
k=1

t
X

(a)

t
var
n (n − t)

(b)

t
n−t
f (0) +
f (0) + a1n + a2n = f (0) + a3n ,
n
n

≤

=

k=t+1

+

n−t
var
nt

Xk

k=1

!
!

t
n
X
X
2
− cov
Xk ,
Xk
n
k=1

k=t+1

!

+ a1n
(A.27)

in which inequality (a) follows from Lemma C.2 and (b) is implied by identity (2.1.3) in [62]. Thus, there
is n0 ∈ N (depending on f and ϑ) such that for any n ≥ n0 , maxt∈Cn,α var (Ut ) ≤ (1 + ϑ) f (0). The rest of
proof will be omitted because of the analogy to (A.24)-(A.26) in the proof of Theorem 5.1.

A.5

Proofs for Section 6

Proof of Theorem 6.1. We follow the standard method for bounding the Bayes risk from below. Observe
that
inf ϕn (T ) =
T

≥

1 − sup inf [P (T = 0 | H0 ) − P (T = 0 | H1,t )]
T

t∈Cn,α

(a)

1 − inf sup |P (T = 0 | H0 ) − P (T = 0 | H1,t )| ≥ 1 − inf H (P0 , P1,t ) ,
t∈Cn,α T

t∈Cn,α

2

where (a) follows from inequality 2.27 in [55]. So, it suﬃces to show that inf t∈Cn,α H 2 (P0 , P1,t ) ≤ (1 − δ) . A
few lines of straightforward algebra on the explicit form of Hellinger distance of Gaussian measures indicates
that inf T ϕn (T ) ≥ δ, whenever


1
−1
2
⊤
.
(A.28)
b inf ζt (Σn ) ζt ≤ 32 log
t∈Cn,α
δ (2 − δ)
Henceforth, it is enough to obtain a tight upper bound on inf t∈Cn,α ζt⊤ (Σn )
2p−1

−1

ζt .

′
CK
Γ(p−1/2)
√

Let σ = 1 and choose d > 0 by d
=
. Furthermore, let F̂d,p,σ : R 7→ R denote the
4πΓ(p)
Matern spectral density parametrized by p and d and σ as (3.2). Note that d is well deﬁned due to the
ﬁrst condition in Assumption 6.1. Deﬁne ξt ∈ Rn by ξt (k) = 1k>t and let ξt′ = ξt − ζt for any t ∈ Cn,α .
Moreover, let θn = exp (−d/n) and St = {t + 1, . . . , n}. Finally, deﬁne the covariance matrix Ψn ∈ Rn×n by
n
Ψn = [Fd,p,σ ((r − s) /n)]r,s=1 . Observe that
−1

ζt⊤ (Σn )

ζt



−1
−1
−1
= 2 ξt⊤ (Σn ) ξt + ξt′⊤ (Σn ) ξt′ − 1⊤
1n
n (Σn )


≤ 4 ξt⊤ (Σn )−1 ξt ∨ ξt′⊤ (Σn )−1 ξt′ .

(A.29)

−1
We aim to prove that there is a constant C := C (p) > 0 for which ξt⊤ (Σn ) ξt ≤ Cn2p−1 . The same upper
−1
bound can be obtained for ξt′⊤ (Σn ) ξt′ in an analogous manner.
We ﬁrst show that
!
K̂
− 1 ∈ L2 (R) .
(A.30)
F̂d,p,σ

Let M represent the ﬁnite lim sup in the second condition of Assumption 6.1. Without loss of generality, we

28

can assume that β < 2 in Assumption 6.1. Using a few lines of algebra along with this condition, we get
!
!
!


2 p
2p
2p
d
K̂
(ω)
K̂
(ω)
ω
K̂
(ω)
ω
1+ 2
−1
lim sup ω β
− 1 ≤ lim sup ω β
−1
= lim sup ω β
′
′
CK
ω
CK
ω→∞
ω→∞
ω→∞
F̂d,p (ω)
"

p
#
2p
d2
β K̂ (ω) ω
1+ 2
+ lim sup ω
−1
′
CK
ω
ω→∞
(a)

=

M+

2pd2
2p−2+β (b)
lim sup K̂ (ω) |ω|
= M.
′
CK ω→∞
p

Notice that, identity (a) follows from Assumption 6.1 and ﬁrst order Taylor expansion of (1 + x) for inﬁnitesimal x > 0. Moreover, (b) follows from the combination of β < 2 and the ﬁrst condition in Assumption
6.1. Namely, there is R > 0 such that
K̂ (ω)
F̂d,p (ω)

−1 ≤

2M
β

|ω|

,

∀ |ω| ≥ R,

which substantiates (A.30) as β > 1/2.
It is known (4.31, Chapter III, [29]) that there is a function φ ∈ L2 (R) with bounded support such

that F̂d,p,σ (ω) ≍ φ̂ (ω)

2

as |ω| → ∞. Theorem 4 of Skorokhod [51] implies that the associated zero mean

Gaussian measures to spectral densities K̂ and F̂d,p,σ are equivalent. Based upon Lemma B.4, there exists
a constant H ∈ (0, ∞) such that
ξ ⊤ (Σn )−1 ξt
1
≤ H.
≤ lim t⊤
n→∞ ξ (Ψ )−1 ξ
H
n
t
t
−1

So, it suﬃces to show that ξt⊤ (Ψn ) ξt ≤ C ′ n2p−1 for some appropriately chosen C ′ > 0 depending on H
and C.
Letting ν = p − 1/2 and recalling An , W and Dn form the proof of Theorem 3.1, we have
ξt⊤ (Ψn )

−1

(b)

⊤

ξt = (An ξt ) Dn−1 (An ξt ) ≤

2

kAn ξt kℓ2
.
λmin (Dn (St , St ))

(A.31)

Note that inequality (b) is inferred from supp (An ξt ) = St . Applying a similar technique as (A.5), we get
= (n − t − p) (1 − θn )p +

kAn ξt k2ℓ2

p
X

k=1




k−1
X
j=0

2

2

 
p
k−1
X
X
p
p

(−θn )j  ≤ n (1 − θn )p + 2
(−1)j 
j
j
j=0
k=1



2p − 2
p
p−1
p
≤ 2
+ n (1 − θn ) ≤ αp n−(p−1) + 2 (2e)
≤ (2e) .
p−1
−1

2p

(A.32)

−1

So, ξt⊤ (Ψn ) ξt ≤ (2e) [λmin (Dn (St , St ))] .
Next, we control the smallest eigenvalue of Dn (St , St ) from the below. We ﬁrst control the diagonal
entries from below. Note that all the diagonal entries of Dn (St , St ) are the same and given by (cf. (A.8))
Q =

=

−p
Z
p

p
d2 + ω 2
(c)
F̂d,p (ω) 
2
1 + θn − 2θn cos (ω/n) dω ∝
1 + θn2 − 2θn cos (ω/n) dω
2π
d−2ν
R
R
#p
Z "
Z
2
2
(d) n−2ν
(1 − θn ) + 4θn sin (dω/2)
n−2ν
dω
≥
[sinc (dω/2)]2p dω = Cd′ n−2ν ,
(A.33)
1/n2 + ω 2
2
Z

R

R

29

where (c) is obtained from (3.2) and the inequality (d) follows from the fact that for any γ ∈ (0, 1) (here we
1
put γ = 2− p ), there is n0 (γ) such that for any n ≥ n0 (γ),
(1 − θn )2 + 4θn sin2 (dω/2)
2
≥ γd2 [sinc (dω/2)] .
1/n2 + ω 2
The proof of the last inequality will be skipped due to its simplicity.
Now, let Ξ := Dn (St , St ) /Q. The combination of (A.31), (A.32) and (A.33) shows that
ξt⊤ (Ψn )−1 ξt ≤

C0 n−2ν
C0′ n2ν
C ′ n2p−1
⇒ ξt⊤ (Σn )−1 ξt ≤
= 0
,
λmin (Ξ)
λmin (Ξ)
λmin (Ξ)

for some constants, C0 (p) and C0′ depending on C0 , H and K. It can be shown using identity 1.2 of [8]
that there is some integrable function g : [−π, π] 7→ R with mg := essinf (g) > 0 such that Ξ is a p−banded
correlation matrix, i.e. Ξ (r, s) = 0 for |r − s| ≥ p, and Ξ = Tn (f ). In remains to note that Lemma 6 of [22]
implies that λmin (Ξ) > mg for any n, which concludes the proof.
Proof of Theorem 6.2. The proof is similar to the proof of Theorem 6.1, with some minor diﬀerence in the
detail. Applying the classical technique of bounding ϕn (T ) from below in terms of Hellinger distance, we
need to verify (A.28). Recalling the formulation of ξt and ξt′ from the proof of Theorem 6.1, we get the
following inequality for any t ∈ Cn,α ,
−1

ζt⊤ (Σn )

 (b) 4n
(a) 
(c) 4n (1 + ϑ)
−1
−1
ζt ≤ 4 ξt⊤ (Σn ) ξt ∨ ξt′⊤ (Σn ) ξt′ ≤
+ o (n) ≤
.
f (0)
f (0)

(A.34)

Notice that (a) is exactly the same as inequality (A.29). Moreover, (b) follows from Corollary C.1 and lastly,
there is n0 ∈ N, which depends on f and ϑ, for which inequality (c) holds. Using (A.34), one can easily
verify inequality (A.28) and concluding the proof.

B

Auxiliary results

This section contains several technical results needed in Appendix A.
2
Lemma B.1. Let σ0 ≥ 1 and n ≥ 2. Let Z ∈Rn be a Gaussian
random
q
 vector with EZ = µ and var Zk ≤ σ0


. For any δ ∈ (0, 1) and any n ∈ N, the
+ log 2n
for any 1 ≤ k ≤ n. Moreover, let Rn = 1 + 2 log 2n
δ
δ
following results hold.


1. If µ = 0, then P max Z2j ≥ σ02 Rn ≤ δ2 .
1≤j≤n

2. If max |µj | ≥ 4σ0
1≤j≤n

q
log

2n
δ


, then P




max Z2j ≤ σ02 Rn ≤ 2δ .

1≤j≤n



Zj
σj

2

are standard χ21 random variables, for

δ
δ
any j = 1, . . . , n. Lemma 8.1 in [7] implies that P Z2j ≥ σj2 Rn ≤ 2n
. Thus, P Z2j ≥ σ02 Rn ≤ 2n
due to
σj ≤ σ0 . We conclude the proof of the ﬁrst part by a union bound argument. Now, we turn to prove the
second part. Deﬁne k := arg max |µj |. It is easy to verify that Rn ≤ 4 log 2n
δ . Observe that
Proof. For brevity, let σj = var Zj , j = 1, . . . , n. Notice that



1≤j≤n

P



max Z2j ≤ σ02 Rn

1≤j≤n



"

Z2
≤ P k2 ≤
σk



σ0
σk

2

30

Rn

#

"

Z2
≤ P k2 ≤ 4
σk



σ0
σk

2

log



2n
δ

#

.

Z2k
2
σk

is a non-central χ21 random variables with non-centrality parameter Bk := µσkk . The lower
q

bound condition on |µk | implies that Bk ≥ 4 σσk0 log 2n
δ . We ﬁnish the proof by the following inequality,
s
"
 # (b)
 2
 # (a) " 2
Z2k
2
δ
Zk
σ0
2n
2
2
P 2 ≤4
≤ P 2 ≤ 1 + Bk − 2 (1 + 2Bk ) log
≤ .
log
σk
σk
δ
σk
δ
2
q
 2


In order to demonstrate inequality (a), we need to show that 1+Bk2 −2 (1 + 2Bk2 ) log 2δ ≥ 4 σσk0 log 2n
δ
which can be shown by obvious inequality σ0 /σk ≥ 1 and a few lines of algebra. Inequality (b) can be inferred
from Lemma 8.1 of [7].
Moreover,

Proposition B.1 (Kantorovich inequality, (p. 452, [24])). Let Σ ∈ Rn×n be a non-singular covariance

matrix and let V ∈ Rn be a non-zero vector. Then, V ⊤ Σ−1 V ≥

kV k4ℓ
2
V ⊤ ΣV

.

Lemma B.2. Let K0 and
function
with spectral densities K̂0 and K̂1 , respectively.
n K1 be two covariance


r−s n
:=
Deﬁne, Σ0 := K0 r−s
K
and
Σ
.
Suppose that there exists ρ ∈ (0, 1) such that
1
1
n
n
r,s=1
r,s=1
K1
K0

−1

∞

≤ ρ. Then,

(a) Σ0 − Σ1 

ρ
1−ρ Σ1 .

−1
(b) Σ−1
1 (Σ0 − Σ1 ) Σ1 

ρ
−1
1−ρ Σ1 .
K1
K0

Proof. Trivial calculations on

−1

∞

≤ ρ shows that for any ω ∈ R,

−ρ
ρ
K̂1 (ω) ≤ K̂0 (ω) − K̂1 (ω) ≤
K̂1 (ω) .
1+ρ
1−ρ

(B.1)

Choose v ∈ Rn arbitrarily. The basic properties of spectral density and inequality (B.1) imply that
⊤

2πv (Σ0 − Σ1 ) v

=

Z 
R

=

2
n
 X
jωℓ/n
K̂0 (ω) − K̂1 (ω)
vℓ e
dω ≤
ℓ=1

ρ
(1 + ρ)

Z

K̂1 (ω)

R

2πρ ⊤
v Σ1 v.
1−ρ

n
X

2

vℓ e

jωℓ/n

dω

ℓ=1

ρ
1−ρ Σ1 .

The second inequality is an obvious implication of the ﬁrst inequality.


δ
. Then,
Lemma B.3. Let δ ∈ (0, 2), d ∈ (0, ∞) and deﬁne K : R 7→ R by K (r) = σ 2 exp − dr
Thus, Σ0 − Σ1 

lim K̂ (ω) |ω|

1+δ

ω→∞

σ 2 δΓ (δ) sin
= Cδ (d, σ) :=
πdδ

πδ
2



.

Proof. It is obvious that Cδ (d, σ) = σ 2 Cδ (d, 1), so without loss of generality assume that σ = 1. It is trivial
δ
that K (r) is of index δ as |r| → 0, i.e. lim 1−K(rλ)
1−K(r) = λ for any λ > 0. Thus, applying the Tauberian
Theorem (p. 35, [53]) leads to
lim [1 − K (1/ω)]

ω→∞

|r|→0

−1

Z∞

K̂ (u) du =

ω

31

Γ (δ) sin
π

πδ
2



=

Cδ (d, 1) dδ
.
δ

(B.2)

−1

−δ

Moreover, the ﬁrst order Taylor expansion of e−x is at 0, implies that [1 − K (1/ω)] (dω)
Thus, (B.2) can be rewritten by last limiting identity and applying L’Hospital’s rule.
−δ

Cδ (d, 1) = lim δd
ω→∞

δ

(dω) δω

δ

Z∞

K̂ (u) du = lim δω
ω→∞

ω

δ

Z∞
ω

→ 1 as ω → 0.

1+δ

K̂ (u) du = lim K̂ (ω) |ω|
ω→∞

.

The following Lemma is probably well-known in the literature of Gaussian processes (e.g. the identity 2
of [53] (p.112) is analogous but not exactly same as the part (a) of Lemma B.4). Because of the absence of
direct references, we include and prove the following result in this section.
Lemma B.4. Let Gi , i = 1, 2 be two zero mean stationary Gaussian process in [0, 1] associated to covariance
functions
For

 any n ∈ N, deﬁne two positive deﬁnite covariance matrices by
 Ki , i= 1, 2, respectively.
r−s
:=
K
and
Ψ
. If G1 and G2 induce equivalent measures on the Hilbert space of
Σn := K1 r−s
2
n
n
n
L2 ([0, 1]), then there exists an scalar B ∈ [1, ∞) for which
1.

1
B

≤ lim inf v6=0n

2.

1
B

≤ lim inf v6=0n

n→∞

n→∞

v ⊤ Σn v
v ⊤ Ψn v

≤ lim supv6=0n
n→∞

v ⊤ Σ−1
n v
v ⊤ Ψ−1
n v

v ⊤ Σn v
v ⊤ Ψn v

≤ lim supv6=0n
n→∞

≤ B.

v ⊤ Σ−1
n v
v ⊤ Ψ−1
n v

≤ B.

Proof. We use Pi , i = 1, 2 to denote the probability measures with respect to Gi , i = 1, 2, respectively.
Abusing the notation, X ∈ Rn represents the random vector generated by sampling Gaussian process at
⊤
Σn v
≤ B1 .
{k/n}nk=1 for any n ∈ N. We prove the existence of a ﬁnite scalar B1 for which lim supv6=0n vv⊤ Ψ
nv
n→∞

Assume toward contradiction that lim supv6=0n
vectors {vn ∈

∞
Rn }n=1

n→∞

v ⊤ Σn v
v ⊤ Ψn v

tends to inﬁnity. So, there is a sequence of non-zero

such that

vn⊤ Σn vn
= ∞.
⊤
n→∞ vn Ψn vn
h
i
p
Consider the measurable event En = |hvn , Xi| ≥ vn⊤ Σn vn . Simple calculations shows that
lim sup

P1 (En ) = Q (1) ,

P2 (En ) = Q

in which Q (·) stands for the Q-function, i.e. Q (r) =
(B.4) leads to lim sup
n→∞

P1 (En )
P2 (En )

R

√1
2π

s

vn⊤ Σn vn
vn⊤ Ψn vn

exp −x2 /2



!

,

(B.3)

(B.4)

1 (|x| ≥ r) dx. Combining (B.3) and

= ∞ which contradicts the absolute continuity of P1 with respect to P2 . One

cam show using the same technique that there is B2 ∈ (1, ∞) such that
v ⊤ Σn v
1
≤ lim inf ⊤
.
n→∞ v6=0n v Ψn v
B2
We conclude the proof by choosing B = B1 ∨ B2 . Now, we turn to substantiate the second claim. Pick a
non-zero vector v ∈ Rn . According to Lemma B.5, there is an suitably chosen n−dimensional vector u (The
inner product of u and v is necessarily 1) such that
v ⊤ Σn−1 v
u⊤ Ψn u (a)
u⊤ Ψn u
⊤ −1
⊤
≤ B.
≤
=
v
Σ
vu
Ψ
u
=
n
n
maxhω,vi=1 ω ⊤ Σn ω
u ⊤ Σn u
v ⊤ Ψ−1
n v
Note that the inequality (a) is obtained from the ﬁrst part of this Lemma. Taking supremum over all
non-zero v ∈ Rn and n ∈ N terminates the proof.
32

Lemma B.5. Let Σ ∈ Rn×n be a non-singular covariance matrix and let ω ∈ Rn be a non-zero vector.
Then,
−1
ω ⊤ Σ−1 ω
= min v ⊤ Σv.
(B.5)
hv,ωi=1

Proof. Since the optimization problem in (B.5) is a convex program with continuously diﬀerentiable objective
function and constraint, so its minimal value can be obtained solving the KKT equations. That is, there are
λ̂ ≥ 0 and v̂ such that
2Σv̂ − λ̂ω = 0, λ̂ (hv̂, ωi − 1) = 0.
Solving the above set of equations yields, v̂ =
into the right hand side of (B.5).

Σ−1 ω
.
ω ⊤ Σ−1 ω

The desired result will be established by replacing v̂

2

Lemma B.6. Let Ω = {η = (d, σ) : η ∈ Ω} ⊂ (0, ∞) be a compact set such that
dist (Ω, {(x, y) : x = 0 or y = 0}) > 0.

(B.6)

The conditions of Proposition ?? are satisﬁed for the following scenarios.
1. K (·, η) is the powered exponential covariance function with known β ∈ (0, 2) and η ∈ Ω.
2. K̂ (·, η) is the Matern spectral density with known ν, given by (3.2).
Proof. We ﬁrst substantiate part 1. It is easy to verify the ﬁrst condition in Proposition ?? for powered
exponential covariance. For proving the continuity condition, let ηm be a convergent sequence in Ω to η. So
Z
≤ |K (r, ηm ) − K (r, η)| dr → 0, as m → ∞
K̂ (·, ηm ) − K̂ (·, η)
∞

R

due to the dominated convergence Theorem. Finally, Lemma B.3 and a few lines of algebra imply that
q
2
σ −2 + (β/d) < ∞.
sup k∇ log CK,Ω (η)kℓ2 = 2 sup
η∈Ω

(σ,d)∈Ω

Notice that the last inequality is a consequence of (B.6). Now we turn to the proof of second part. The
veriﬁcation of the ﬁrst two conditions
in Proposition ?? is analogous to the powered exponential covariance
√
2
ρ−2ν for Matern covariance. Thus,
function. Note that CK,Ω (η) = σ 4πΓ(ν+1/2)
Γ(ν)
sup k∇ log CK,Ω (η)kℓ2

η∈Ω

q
= 2 sup
σ −2 + (ν/d)2 < ∞.
(σ,d)∈Ω

Lemma B.7. Suppose that K̂ (ω, η) satisﬁes the conditions in the part 2 of Lemma B.6 and the taper
function fˆtap (ω) admits Assumption ??. Let K̂tap (·, η) be the convolution of K̂ (ω, η) and fˆtap . Then,
g (η) =

K̂tap (·, η)
K̂ (·, η)

−1

,
2

is a continuous function of η in Ω.
Proof. Let dmax = sup {d : (σ, d) ∈ Ω} and dmin = inf {d : (σ, d) ∈ Ω}. Note that 0 < dmin ≤ dmax < ∞ as
Ω is a compact set in the interior of the upper right half-plane. deﬁne h : R × Ω 7→ (0, ∞) by
h (ω, η) :=

K̂tap (ω, η)
K̂ (ω, η)

−1=

Z∞

−∞

33

K̂ (u, η)
fˆtap (ω − u)
du − 1.
K̂ (ω, η)

We ﬁrst show that for any ﬁxed ω ∈ R, h (ω, η) is a continuous function of η. Deﬁne two real valued function
 2
ν+1/2
d
+ω 2
K̂(u,η)
hd,ω (u) = dmax
. Note that hd,ω is integrable and for any u ∈ R,
and p (u, η) := fˆtap (ω − u) K̂(ω,η)
2
2
+u
min
p (u, ·) is a continuous function of η. Moreover,
sup
η=(σ,d)∈Ω

p (u, η) = fˆtap (ω − u) sup

η∈Ω



d2 + ω 2
d2 + u2

ν+1/2

≤ sup

η∈Ω



d2 + ω 2
d2 + u2

ν+1/2

≤ hd,ω (u) .

∞

Choose η ∈ Ω and a convergent sequence {ηn }n=1 ⊂ Ω to η, in an arbitrary fashion. The dominated
convergence theorem shows that
lim h (ω, ηn ) = lim

n→∞

Z∞

n→∞
−∞

p (u, ηn ) du = h (ω, η) ,

which conﬁrms the continuity of h with respect to η. Now, we demonstrate the continuity of g. Notice that
g (η) = kh (·, η)k2 . One can show using the exact same technique as Lemma 4 of [16] (with a minor algebraic
−(1+r)
diﬀerence) that there is a function q : R 7→ R given by q (ω) = c1 ∨ c2 |ω|
, where c1 , c2 and r are
appropriately chosen ﬁnite positive scalars, uniformly dominating h (ω, η) over Ω. Namely,
sup h (ω, η) ≤ q (ω) .

η∈Ω

Thus,


lim g (ηn ) − g (η) ≤ lim kh (·, η) − h (·, ηn )k2 = 

n→∞

n→∞

Z∞

−∞



1/2

h2 (ω, η) − h2 (ω, ηn ) dω 


(a)

= 0,

terminating the proof. Notice that identity (a) follows from the continuity of h with respect to η and the
dominated convergence theorem.
Lemma B.8. Let K be a covariance function such that K̂ ′
Then, there is a universal constant c > 0 such that
inf

Z∞

β∈(0,1)
−∞

∞

< ∞ and deﬁne Gβ : R 7→ [0, 1] by (A.23).

K̂ (ω) Gβ (ω) dω ≥ c.

Proof. Observe that for any ω ∈ R, Gβ (ω) is a quadratic function of β in the compact interval [0, 1] and
lim kGβn − Gβ k∞ = 0 for any convergent sequence βn → β. This property implies that

n→∞

inf

Z∞

β∈(0,1)
−∞



Z∞
Z∞
1
K̂ (ω) G0.5 (ω) dω  .
K̂ (ω) Gβ (ω) dω ∧
inf
K̂ (ω) Gβ (ω) dω ≥
2 β∈(0,1), |β−1/2|≥r

(B.7)

−∞

−∞

2

for some suﬃciently small r > 0. Observe that, Gβ (0) = (1 − 2β) > 0 for β 6= 1/2. The diﬀerentiability of
Gβ and K̂ (ω) implies the existence of a non-degenerate open interval Iβ centered at 0 such that,
2

(1 − 2β) K̂ (0)
inf K̂ (ω) Gβ (ω) ≥
⇒
ω∈Iβ
2

Z∞

−∞

34

2

K̂ (ω) Gβ (ω)
(1 − 2β) K̂ (0)
dω ≥
|Iβ | .
2π
4π

2

Notice that inf |β−1/2|≥r (1 − 2β) |Iβ | > 0. So, we just need to show that the corresponding term to β = 1/2
2
in the right hand side of (B.7) is strictly positive. For β = 1/2, Gβ (ω) = [sinc (ω/4) sin (ω/2)] and so
Z∞

−∞

K̂ (ω) Gβ (ω)
dω ≥
2π

Z2π

−2π

2

(b) 2
K̂ (ω) [sinc (ω/4) sin (ω/2)]
dω ≥ 3
2π
π

Z2π

(c)

K̂ (ω) sin2 (ω/2) dω > 0.

−2π

Note that (b) is a consequence of monotonicity of sinc (·) in the interval (0, π/2) and inequality (c) follows
from the combination of K̂ ′ (0) < ∞ and K̂ (0) > 0.

C

Non-asymptotic behaviour of the inverse of Toeplitz matrices

In this section, we investigate some non-asymptotic properties of the inverse of Toeplitz matrices with
polynomially decaying oﬀ-diagonal entries. The developed results plays a crucial role in the analysis of
GLRT in increasing domain. We ﬁrst introduced some simplifying notation. For any symmetric and periodic
function f ∈ L∞ (R) with period 2π, deﬁne


f0 f1 f2 · · ·


 f1 f0 f1 . . . 


TN (f ) := 
.
 f2 f1 f0 . . . 


.. . .
..
..
.
.
.
.

in which {fm }m∈Z denotes the set of Fourier coeﬃcients of f and f−m = fm for any m ∈ Z. Moreover, let
Tn (f ) = [TN (f )]nr,s=1 . Finally, for any Sn , Sn′ ⊆ {1, . . . , n} and any Ln ∈ Rn×n deﬁne

1⊤ Ln1Sn′
,
τ (Ln , Sn , Sn′ ) := pSn
|Sn | |Sn′ |

in which 1Sn ∈ Rn represents the indicator vector of Sn in {1, . . . , n}. For the sake of brevity, we use
the shorthand notation τ (Ln , Sn ) if Sn = Sn′ . Throughout this section we assume that there exists some
α ∈ (0, 1/2) such that Cn,α = [αn, (1 − α) n] and Ln is a Toeplitz covariance matrix where its generator
satisﬁes Assumption 3.2 for some c, λ, mf and Mf . We also use LN to indicate the inﬁnite sized version of
Ln .
Lemma C.1. Let Sn = {1, . . . , t} for some t ∈ Cn,α and Snc = {1, . . . , n} \ Sn . There is a bounded constant
C ≥ 0 depending on f and α such that
c
τ L−1
n , Sn , Sn



C
≤ p
n−(λ∧1)
α (1 − α)

1{λ6=1} + 1{λ=1} log n .


Proof. Without loss of generality assume that t ≤ n − t. The proof is based upon the well known fact that
for a positive deﬁnite matrix with polynomially decaying oﬀ-diagonal entries, the corresponding elements of
its inverse shrink with same rate (See e.g. [31]). For any d ∈ {1, . . . , n − 1}, let
ζ (d, t) := {(r, s) : r ∈ Sn , j ∈
/ Sn , |r − s| = d} .
A simple counting argument leads to



d,
t,
|ζ (d, t)| =

n − d,


d ∈ {1, . . . , t}

d ∈ {t, . . . , n − t}
.

d ∈ {n − t + 1, . . . , n − 1}
35

(C.1)


c
Using triangle inequality and rearranging the diﬀerent components in τ L−1
n , Sn , Sn , we get
c
τ L−1
n , Sn , Sn



n−1
X
1
≤p
Ln−1 (r, s) .
|ζ (d, t)| max
(r,s)∈ζ(d,t)
t (n − t) d=1

(C.2)

According to Lemma A1 of [23], there is a C (f ) for which the entries of L−1
n admits the following inequality,
C̃ (f )

L−1
n (r, s) ≤

(1 + |r − s|)

1+λ

(C.3)

.

We terminate the proof by substituting inequality (C.3) and identity (C.1) into inequality (C.2). We will
skip the algebraic details due to lack of space.
!
n−1
n−t
t
X
X
X

C̃
n−d
t
d
−1
c
p
≤
τ L n , Sn , Sn
+
+
1+λ
1+λ
1+λ
t (n − t) d=1 (1 + d)
d=n−t+1 (1 + d)
d=t+1 (1 + d)
(a)

≤

C
p
n−(λ∧1)
α (1 − α)

Note that inequality (a) follows from the fact that

1{λ6=1} + 1{λ=1} log n .


p
p
t (n − t) ≥ n α (1 − α).

The proof of the following result is omitted due to its analogy to the proof of Lemma C.1.
Lemma C.2. With the same conditions and notation as Lemma C.1, there is a non-negative constant
C ′ := C ′ (f ) such that
C′
|τ (Ln , Sn , Snc )| ≤ p
n−(λ∧1)
α (1 − α)

1{λ6=1} + 1{λ=1} log n .


Proposition C.1. Deﬁne the inﬁnite column vector νn ∈ RN by νn (r) = 1{r≤n} . Then, there is a constant
C := C (f, α) such that


1
1 ⊤ −1
≤ C n−(λ−ǫ) 1{λ∈(0,1]} + n−1 1{λ>1} ,
νn L N νn −
n
f (0)
in which ǫ ∈ (0, λ) is chosen in an arbitrary way.

−1
Proof. We obtain an upper bound on νn⊤ LN
νn /n − 1/f (0) and the lower bound can be achieved using
akin techniques. For brevity, let g = 1/f . According to Proposition 1.12 of [9], L−1
N satisﬁes the following
identity,
−1
L−1
N = TN (g) + LN HN (f ) HN (g) ,

in which HN (f ) is the generated Hankel matrix by f

f1
 f2

HN (f ) =  f3

..
.

as
f2
f3
f4
..
.

36

f3
f4
f5
..
.

···
···
···
..
.





.


(C.4)

HN (g) can also be deﬁned in a similar way. A simple algebraic manipulation on identity C.4 leads to
1 ⊤ −1
ν L νn
n n N

ν ⊤ L−1 HN (f ) HN (g) νn
1 ⊤
νn TN (g) νn + n N
n
n
n X
n
X
1 ⊤
1
ν TN (g) νn +
hL−1 HN (f ) eN (r) , HN (g) eN (k)i
n n
n r=1 s=1 N

=
=
≤

L−1
1 ⊤
N
νn TN (g) νn +
n
n

2→2

≤

1
1 ⊤
ν TN (g) νn +
n n
nmf

n
X

(a)

n
X
r=1

r=1

kHN (g) eN (r)kℓ2

n
X
s=1

n
X

kHN (g) eN (r)kℓ2 ∨

s=1

kHN (f ) eN (s)kℓ2

kHN (f ) eN (s)kℓ2

!2

.
(C.5)

Note that (a) is direct consequence of the combination of generalized Cauchy-Schwartz inequality and simple
properties of operator norm; The last inequality of C.5 can be obtained by inequality (1.14) of [9] on
n
P
the operator norm of Toeplitz matrices. As next step, we control
kHN (g) eN (r)kℓ2 from above. It is
r=1

known that if mf > 0 and f satisﬁes condition 3.9, then g does as well. So, there bounded are constants
c (λ) , c′ (λ) , c′′ (λ) > 0 such that
v
v
uZ∞
u∞
n
n
n
n
u
X uX
X
X
X
u
t
kHN (g) eN (r)kℓ2 =
j −(1/2+λ)
|gk |2 ≤ c
t x−(2+2λ) dx ≤ c′
r=1

j=1

r=1

k=r

r=1

r




c′′ n1/2−λ 1{λ∈(0,1/2)} + log n1{λ=1/2} + 1{λ>1/2} .

≤

(C.6)

Henceforth, there is some constant strictly positive c̃ (depending on λ) such that


1
1 ⊤ −1
νn LN νn ≤ νn⊤ TN (g) νn + c̃ n−2λ 1{λ∈(0,1/2)} + n−1 log2 n1{λ=1/2} + n−1 1{λ>1/2} .
n
n

Lastly, we obtain an upper bound on

1 ⊤
n νn TN

(g) νn . νn⊤ TN (g) νn can be viewed as variance of Sn =

(C.7)
n
P

Xk

k=1

where {Xk }k∈N is a stationary process with spectral density g. Thus, identity 2.1.2 of [62] shows that
π

1 ⊤
ν TN (g) νn
n n

=

+

(b)

≤

1
πn
2
πn

π

Z2 

−π
2

π
2

1
+
f (0)

2

g (0)
g (2ω) dω =
πn

sin (nω)
sin (ω)

2

1
2
[g (2ω) − g (0)] dω =
+
f (0) πn

−π
2

sin (nω)
sin (ω)

2

dω
π

Z 

−π
2

Z2 

sin (nω)
sin (ω)

2

P

m∈Z

|fm | |m|

πnmf f (0)

(a)

λ′

Z2 

−π
2

sin (nω)
sin (ω)

2

[g (2ω) − g (0)] dω

π

Z2 

−π
2

sin (nω)
sin (ω)

2

′

|ω|λ dω.

(C.8)

Identity (a) follows from the following results which can be proved by applying Parseval’s identity on the
triangular pulse centred at 0.
π
2
Z2 
1
sin (nω)
dω = 1
πn
sin (ω)
−π
2

37

Based upon the following inequality, g inherits the Holder property from f .
|g (2ω) − g (0)| =

|f (2ω) − f (0)|
f (2ω) − f (0)
≤
.
f (0) f (2ω)
f (0) mf

Hence, (b) is implied by this inheritance property and Proposition 3.2.12 of [21] regarding that f is a (λ′ )Holder function for any 0 ≤ λ′ < λ.
For any λ ∈ (0, 1], choose ǫ = ǫ (λ) > 0 such that (λ − ǫ) > 0. Let λ′ = (λ − ǫ) 1{λ∈(0,1]} + 1+λ
2 1{λ>1} .
P
λ′
′
′
|fm | |m| is bounded, so,
Obviously λ < λ and λ > 1 for any λ > 1. Notice that
m∈Z

π

π

1
n

Z2 

−π
2

sin (nω)
sin (ω)

2

λ′

|ω| dω

(a0 )

≍

(a1 )

≤

Z2
0

nπ

sin2 (nω)
2

(nω)

′

λ

′

|ω| dnω = n−λ

Z2

sin2 (u)
du
u2−λ′

0


ĉ n−λ+ǫ 1{λ∈(0,1]} + n−1 1{λ>1} .

Inequality (a0 ) follows from the fact that π2 ≤ sinω ω ≤ 1 for any ω ∈ 0, π2
integration techniques. Combination of (C.7)-(C.9) and concludes the proof.



(C.9)

and (a1 ) is given by simple

Proposition C.2. Let Sn = {1, . . . , n}. Under the same notation and assumptions as Proposition C.1, we
have



1 ⊤ −1
−1
−(λ−ǫ)
1
+
n
1
.
≤
C
n
νn LN νn − τ L−1
,
S
n
{λ∈(0,1]}
{λ>1}
N
n

Proof. Using Widom’s theorem (Theorem 2.14, [9]) and one line of straightforward algebra, there is a matrix
Dn ∈ Rn×n such that

1 ⊤ −1
νn LN νn − τ L−1
N , Sn
n

=
≤


1 ⊤
νn Vn L−1
N − TN (g) Vn νn + τ (Dn , Sn )
n

1
kDn k2→2 + νn⊤ Vn L−1
N − TN (g) Vn νn .
n

Vn is a bounded operator which is deﬁned by Vn (ν) = (νn , · ·· , ν1 , 0, 0, · · · ) for any ν ∈ ℓ2 . Based upon
Theorem 2.15 of [9] (ﬁrst equation, p. 44), kDn k2→2 = o n−λ . Moreover, the special form of Vn , gives


1
1
−1
(Vn νn )⊤ LN
− TN (g) Vn νn = νn⊤ L−1
N − TN (g) νn .
n
n

Ultimately, combining the same tricks as inequality (a) in (C.5), identity (C.4) and inequality (C.6) yields



1 ⊤
−(λ−ǫ)
′
−1
n
−
T
(g)
V
ν
νn Vn L−1
≤
C
1
+
n
1
,
N
n
n
{λ∈(0,1]}
{λ>1}
N
n

for some constant C ′ (λ) > 0. We end the proof by using triangle inequality.
Corollary C.1. There is a constant C > 0 depending on f such that

τ Σ−1
n , {1, . . . , n} −

1
f (0)



≤ C ′ n−(λ−ǫ) 1{λ∈(0,1]} + n−1 1{λ>1} .

(C.10)

Acknowledgment. This research is partially supported by NSF grant ACI-1047871. Additionally, CS
is partially supported by NSF grants 1422157, 1217880, and 0953135, and LN by NSF CAREER award
DMS-1351362, NSF CNS-1409303, and NSF CCF-1115769.

38

References
[1] S. Adak, “Time-dependent spectral analysis of non-stationary time series”, Journal of the American
Statistical Association 93, no. 444 (1998): 1488-1501.
[2] L. Addario-Berry, N. Broutin, L. Devroye, and G. Lugosi, “On combinatorial testing problems”, The
Annals of Statistics, 2010;38(5):3063-92.
[3] E. Anderes, “On the consistent separation of scale and variance for Gaussian random ﬁelds”, The Annals
of Statistics (2010): 870-893.
[4] E. Arias-Castro, E.J. Candès, and A. Durand, “Detection of an anomalous cluster in a network”, The
Annals of Statistics, 2011 Feb 1:278-304.
[5] J. Antoch, M. Hušková and Z. Prášková, “Eﬀect of dependence on statistics for determination of change”,
Journal of Statistical Planning and Inference 60, No. 2 (1997): 291-310.
[6] P. Bertrand, P. Raphael, M. Fhima and A. Guillin, “Oﬀ-line detection of multiple change points by the
ﬁltered derivative with p-value method”, Sequential Analysis 30, no. 2 (2011): 172-207.
[7] L. Birgé, “An alternative point of view on Lepski’s method”, Lecture Notes-Monograph Series, 2001, p.
113-133.
[8] D. Bolin and F. Lindgren, “Spatial models generated by nested stochastic partial diﬀerential equations,
with an application to global ozone mapping”, The Annals of Applied Statistics 5, no. 1 (2011): 523-550.
[9] A. Böttcher and B. Silbermann, “Introduction to large truncated Toeplitz matrices”, Springer, 1999.
[10] C. Butucea and Y.I. Ingster, “Detection of a sparse submatrix of a high-dimensional noisy matrix”,
Bernoulli, 2013;19(5B):2652-88.
[11] V. Chandola and R.R. Vatsavai, “A Gaussian process based online change detection algorithm for
monitoring periodic time series." In SDM, pp. 95-106. 2011.
[12] H. Chernoﬀ, S. Zacks, “Estimating the current mean of a normal distribution which is subjected to
changes in time”, The Annals of Mathematical Statistics, 1964 Sep 1;35(3):999-1018.
[13] M. Csörgő, L. Horváth, “Nonparametric methods for changepoint problems”, Handbook of statistics,
1988 Dec 31;7:403-25.
[14] R. Dahlhaus, “Fitting time series models to nonstationary processes”, The Annals of Statistics 25, no.
1 (1997): 1-37.
[15] J. Deshayes, D. Picard, Oﬀ-line statistical analysis of change-point models using non parametric and
likelihood methods”, In Detection of Abrupt Changes in Signals and Dynamical Systems 1985 (pp.
103-168), Springer Berlin Heidelberg.
[16] J. Du, H. Zhang and V.S. Mandrekar, “Fixed-domain asymptotic properties of tapered maximum likelihood estimators”, The Annals of Statistics 37, no. 6A (2009): 3330-3361.
[17] H. Exton, “q-Hypergeometric functions and applications”, Horwood, 1983.
[18] R. Furrer, M.G. Genton and D. Nychka, “Covariance tapering for interpolation of large spatial datasets”,
Journal of Computational and Graphical Statistics 15, no. 3 (2006).
[19] E. Gabriel, D. Allard and J.N. Bacro, “Estimating and testing zones of abrupt change for spatial data”,
Statistics and Computing 21, no. 1 (2011): 107-120.

39

[20] O. Gillet, S. Essid and G. Richard, “On the correlation of automatic audio and visual segmentations
of music videos”, Circuits and Systems for Video Technology, IEEE Transactions on 17, no. 3 (2007):
347-355.
[21] L. Grafakos, “Classical Fourier analysis”, Vol. 86, Springer, 2008.
[22] R.M. Gray, “Toeplitz and circulant matrices: A review”, now publishers inc, 2006.
[23] P. Hall and J. Jin, “Innovated higher criticism for detecting sparse signals in correlated noise”, The
Annals of Statistics 38, No. 3 (2010), 1686-1732.
[24] R.A. Horn and C.R. Johnson, "Matrix analysis", Cambridge university press, 2012.
[25] L. Horváth, “Detection of changes in linear sequences”, Annals of the Institute of Statistical Mathematics
49, No. 2 (1997): 271-283.
[26] L. Horváth and M. Hušková, “Change point detection in panel data”, Journal of Time Series Analysis
33, no. 4 (2012): 631-648.
[27] L. Horváth, P. Kokoszka, “The eﬀect of long-range dependence on change-point estimators”, Journal of
Statistical Planning and Inference, 1997 Oct 30;64(1):57-81.
[28] X. Hu, H. Qiu and N. Iyer, “Multivariate change detection for time series data in aircraft engine fault
diagnostics”, In Systems, Man and Cybernetics, ISIC. IEEE International Conference on, pp. 2484-2489.
IEEE, 2007.
[29] I.A. Ibragimov and Y.A. Rozanov, “Gaussian random processes”, Springer, 1978.
[30] B.G. Ivanoﬀ and E. Merzbach, “Optimal detection of a change-set in a spatial Poisson process”, The
Annals of Applied Probability (2010): 640-659.
[31] S. Jaﬀard, “Propri’et’es des matrices bien localis’ees pr’es de leur diagonale et quelques applications”,
Ann. Inst. Henri Poincar’e, 7, 1990, pp. 461–476
[32] C.G. Kaufman, M.J. Schervish and D.W. Nychka, “Covariance tapering for likelihood-based estimation
in large spatial data sets”, Journal of the American Statistical Association 103, no. 484 (2008): 15451555.
[33] Y. Kawahara and M. Sugiyama, “Change-point detection in time-series data by direct density-ratio
estimation”, In Proceedings of SIAM International Conference on Data Mining (SDM), pp. 389-400.
2009.
[34] P. Kokoszka and R. Leipus, “Change-point in the mean of dependent observations”, Statistics and
probability letters 40, No. 4 (1998): 385-393.
[35] M. Last and R. Shumway, “Detecting abrupt changes in a piecewise locally stationary time series”,
Journal of multivariate analysis 99, no. 2 (2008): 191-214.
[36] M. Lavielle, “Using penalized contrasts for the change-point problem”, Signal Processing 85, no. 8 (2005):
1501-1510.
[37] T.L. Lai, “Information bounds and quick detection of parameter changes in stochastic systems”, IEEE
Transactions on Information Theory, 1998 Nov;44(7):2917-29.
[38] C. Levy-leduc, “Catching change-points with lasso”, In Advances in Neural Information Processing
Systems(NIPS), pp. 617-624. 2007.

40

[39] J. Liu, L. Yuan, and J. Ye, “An eﬃcient algorithm for a class of fused lasso problems”, In Proceedings of
the 16th ACM SIGKDD international conference on Knowledge discovery and data mining, pp. 323-332.
ACM, 2010.
[40] W.L. Loh and T.K. Lam, “Estimating structured correlation matrices in smooth Gaussian random ﬁeld
models”, The Annals of Statistics 28, no. 3 (2000): 880-904.
[41] I.B. MacNeill, “Test for change of parameter at unknown time and distributions of some related functionals on Brownian motion”, The Annals of Statistics 2, 950-962, 1974.
[42] H.Y. Noh, R. Rajagopal and A.S. Kiremidjian, “Damage diagnosis algorithm using a sequential change
point detection method with an unknown distribution for damage”, In SPIE Smart Structures and
Materials+ Nondestructive Evaluation and Health Monitoring, pp. 834507-834507. International Society
for Optics and Photonics, 2012.
[43] M.P. Paciﬁco, C. Genovese, I. Verdinelli and L. Wasserman, “False discovery control for random ﬁelds”,
Journal of the American Statistical Association 99, no. 468 (2004): 1002-1014.
[44] M.B. Priestley, “Evolutionary spectra and non-stationary processes”, Journal of the Royal Statistical
Society, Series B (Methodological) (1965): 204-237.
[45] J. Qian and L. Su, “Shrinkage estimation of regression models with multiple structural change”, Working
paper, Department of Economics, Singapore Management University, 2013.
[46] T. S. Rao, S. S. Rao and C. R. Rao, “Handbook of statistics: time series analysis, methods and applications”, Vol. 30. Elsevier, 2012.
[47] M. Rencova, “Change point detection in temperature series”, Doctoral dissertation, Czech Technical
University, 2009.
[48] G. Rigaill, “Pruned dynamic programming for optimal multiple change-point detection”, arXiv preprint
arXiv:1004.0887 (2010).
[49] X. Shao and X. Zhang, “Testing for change points in time series”, Journal of the American Statistical
Association 105, no. 491 (2010).
[50] X. Shen, H.C. Huang and N. Cressie, “Nonparametric hypothesis testing for a spatial signal”, Journal
of the American Statistical Association 97, no. 460 (2002): 1122-1140.
[51] A.V. Skorokhod and M.I. Yadrenko, “On absolute continuity of measures corresponding to homogeneous
Gaussian ﬁelds”, Theory of Probability and Its Applications 18, no. 1 (1973): 27-40.
[52] V. Spokoiny, “Multiscale local change point detection with applications to value-at-risk”, The Annals of
Statistics, 2009 Jun 1:1405-36.
[53] M.L. Stein, “Interpolation of spatial data: some theory for kriging”, Springer, 1999.
[54] A.G. Tartakovsky, B.L. Rozovskii, R.B. Blazek and H. Kim, “A novel approach to detection of intrusions
in computer networks via adaptive sequential and batch-sequential change-point detection methods”,
Signal Processing, IEEE Transactions on 54, no. 9 (2006): 3372-3382.
[55] A.B. Tsybakov, “Introduction to non-parametric estimation”, Vol. 11. New York, Springer, 2009.
[56] J. Verbesselt, R. Hyndman, G. Newnham and D. Culvenor, “Detecting trend and seasonal changes in
satellite image time series”, Remote sensing of Environment 114, no. 1 (2010): 106-115.
[57] H. Wendland, “Error estimates for interpolation by compactly supported radial basis functions of minimal degree”, Journal of approximation theory 93, no. 2 (1998): 258-272.
41

[58] W.Y. Wu, C.Y. Lim and Y. Xiao, “Tail estimation of the spectral density for a stationary Gaussian
random ﬁeld”, Journal of Multivariate Analysis 116 (2013): 74-91.
[59] Y.C. Yao and R.A. Davis, “The asymptotic behaviour of the likelihood ratio statistic for testing a shift in
mean in a sequence of independent normal variates”, The Indian Journal of Statistics, Series A (1986):
339-353.
[60] Z. Ying, “Asymptotic properties of a maximum likelihood estimator with data from a Gaussian process”,
Journal of Multivariate Analysis 36, no. 2 (1991): 280-296.
[61] H. Zhang, “Inconsistent estimation and asymptotically equal interpolations in model-based geostatistics”,
Journal of the American Statistical Association 99, no. 465 (2004): 250-261.
[62] L. Zhengyan and L. Chuanrong, “Limit theory for mixing dependent random variables”, Vol. 378,
Springer, 1996.

42

