Confirmatory Bayesian Online Change Point Detection in the Covariance Structure of Gaussian Processes

TL;DR

This paper introduces a Bayesian online change point detection method that confirms covariance structure changes in Gaussian Processes, improving detection accuracy and prediction in nonstationary time series.

Contribution

It proposes a statistically justified confirmatory approach to Bayesian change point detection in GPs, enhancing robustness and accuracy without precise hyperparameter tuning.

Findings

Successfully detects covariance structure changes in synthetic and real data.

Outperforms existing methods in prediction error and log likelihood.

Provides theoretical thresholds for change detection tests.

Abstract

In the analysis of sequential data, the detection of abrupt changes is important in predicting future changes. In this paper, we propose statistical hypothesis tests for detecting covariance structure changes in locally smooth time series modeled by Gaussian Processes (GPs). We provide theoretically justified thresholds for the tests, and use them to improve Bayesian Online Change Point Detection (BOCPD) by confirming statistically significant changes and non-changes. Our Confirmatory BOCPD (CBOCPD) algorithm finds multiple structural breaks in GPs even when hyperparameters are not tuned precisely. We also provide conditions under which CBOCPD provides the lower prediction error compared to BOCPD. Experimental results on synthetic and real-world datasets show that our new tests correctly detect changes in the covariance structure in GPs. The proposed algorithm also outperforms existing methods for the prediction of nonstationarity in terms of both regression error and log likelihood.

Figures

Tables

Table 1. Table 1 : Comparison of BOCPD, CS-BOCPD, and CBOCPD with NLL and MSE on synthetic and Gazebo robot simulation datasets.

	LEN-CHANGE		VAR-CHANGE		Gazebo:Env1		Gazebo:Env2		Gazebo:Env3
Method	NLL	MSE	NLL	MSE	NLL	MSE	NLL	MSE	NLL	MSE
BOCPD	0.99$\pm$0.38	0.47$\pm$0.32	1.13$\pm$0.61	0.55$\pm$0.58	2.07$\pm$0.51	0.14$\pm$0.05	2.24$\pm$0.48	0.57$\pm$0.26	0.28$\pm$0.12	0.11$\pm$0.03
CS-BOCPD	0.95$\pm$0.32	0.43$\pm$0.30	1.05$\pm$0.42	0.48$\pm$0.40	-0.11$\pm$0.33	0.12$\pm$0.04	1.82$\pm$1.51	0.55$\pm$0.27	0.25$\pm$0.16	0.10$\pm$0.02
CBOCPD	0.79$\pm$0.33	0.34$\pm$0.24	0.89$\pm$0.36	0.41$\pm$0.34	-0.31$\pm$0.34	0.11$\pm$0.04	0.69$\pm$0.36	0.45$\pm$0.19	-0.99$\pm$0.47	0.10$\pm$0.04

Table 2. Table 2 : Results of predictive performance on Nile data, Well Log Data and Snow Data. The results are provided with 95% of confidence interval and the p-value of the null hypothesis that a method is equivalent to the best performing method according to NLL, using a one sided t-test.

Method	NLL	p-value	MSE	p-value
Nile Data (200 training points, 463 test points)
ARGP	1.07$\pm$0.64	0.0001	5.06$\pm$0.86	0.0005
ARGP-BOCPD	0.78$\pm$0.72	0.0001	4.94$\pm$0.87	0.0017
GPTS	0.86$\pm$0.64	0.0001	4.78$\pm$0.81	0.0100
BOCPD	0.57$\pm$0.77	0.0014	4.73$\pm$0.82	0.0115
CBOCPD	0.00$\pm$0.80	N/A	4.32$\pm$0.74	N/A
Well Log Data (500 training points, 3050 test points)
ARGP	7.20$\pm$0.60	0.0001	17.3$\pm$2.6	0.0001
ARGP-BOCPD	0.00$\pm$0.30	N/A	4.68$\pm$0.46	N/A
GPTS	3.73$\pm$0.42	0.0001	8.27$\pm$0.61	0.0001
BOCPD	4.35$\pm$0.31	0.0001	19.2$\pm$1.3	0.0001
CBOCPD	0.30$\pm$0.27	0.0010	4.92$\pm$0.44	0.2124
Snow Data (500 training points, 13380 test points)
ARGP	17.48$\pm$0.82	0.0001	14.82$\pm$0.57	0.0001
ARGP-BOCPD	0.06$\pm$0.39	0.0001	9.65$\pm$0.39	0.0001
GPTS	16.60$\pm$0.22	0.0001	8.76$\pm$0.36	0.0001
BOCPD	0.00$\pm$0.39	N/A	9.43$\pm$0.38	N/A
CBOCPD	1.92$\pm$0.37	0.0001	6.34$\pm$0.27	0.0001

Equations

$${\phi }_n(\mathfrak{T})=\mathbb{P}(\mathfrak{T}=1|{\mathbb{H}}_0)+\underset{t\in {\mathcal{C}}_n}{\max }\mathbb{P}(\mathfrak{T}=0|{\mathbb{H}}_{1,t}).$$

$$X^T{\Sigma }^{-1}X-\underset{t\in {\mathcal{C}}_n}{\min }\underset{b\ne 0}{\min }\left[{\left(X-\frac{b}{2}{\zeta }_t\right)}^T{\Sigma }^{-1}\left(X-\frac{b}{2}{\zeta }_t\right)\right],$$

$${\mathfrak{T}}_{GLRT}=\mathbb{I}\left(\underset{t\in {\mathcal{C}}_n}{\max }{\left|\frac{{\zeta }_t^T{\Sigma }^{-1}X}{\sqrt{{\zeta }_t^T{\Sigma }^{-1}{\zeta }_t}}\right|}^2\ge {\mathfrak{R}}_{n,\delta }\right).$$

$${\mathfrak{R}}_{n,\delta }=1+2\left[log\left(\frac{2n}{\delta }\right)+\sqrt{log\left(\frac{2n}{\delta }\right)}\right].$$

$${\mathbb{H}}_{1,t}:\mathrm{Cov}(X_i,X_j)=\{\begin{array}{ll}K(i,j),& i,j<t\\ K^{\prime }(i,j),& i,j\ge t\\ K^{\prime \prime }(i,j),& \text{otherwise}\end{array}$$

$$\underset{t\in {\mathcal{C}}_n}{\max }\left[X^T(\Sigma {)}^{-1}X-X^T({\Sigma }_t^{{}^{\prime }}{)}^{-1}X+\ln \left(\frac{|\Sigma |}{|{\Sigma }_t^{{}^{\prime }}|}\right)\right].$$

$$2\mathfrak{L}-\ln \left(\frac{|\Sigma |}{|{\Sigma }_{t^{\ast }}^{{}^{\prime }}|}\right)=\sum _{i=1}^n{v}_i-\sum _{i=1}^n{\lambda }_i{u}_i$$

$$2\mathfrak{L}-\ln \left(\frac{|\Sigma |}{|{\Sigma }_{t^{\ast }}^{{}^{\prime }}|}\right)=\sum _{i=1}^n{\lambda }_i{u}_i-\sum _{i=1}^n{v}_i$$

$$\Sigma =\left(\begin{array}{cc}{K}_{aa}& K_{ab}\\ K_{ba}& K_{bb}\end{array}\right),{\Sigma }_t^{{}^{\prime }}=\left(\begin{array}{cc}{K}_{aa}& 0\\ 0& K_{bb}^{\prime }\end{array}\right).$$

$${\phi }_n(\mathfrak{T})=\mathbb{P}(2\mathfrak{L}\ge {\mathfrak{R}}_{\delta }|{\mathbb{H}}_0)+\underset{t\in {\mathcal{C}}_n}{\max }\mathbb{P}(2\mathfrak{L}\le {\mathfrak{R}}_{\delta }|{\mathbb{H}}_{1,t})\le \delta .$$

$$\mathbb{P}(x_{t+1}|x_{1:t})=\sum _{r_t}\mathbb{P}(x_{t+1}|x_t^{(r)})\mathbb{P}(r_t|x_{1:t}).$$

$$\sum _{r_{t-1}}\mathbb{P}(r_t|r_{t-1},x_{t-1}^{(r)})\mathbb{P}(x_t|r_{t-1},x_{t-1}^{(r)})\mathbb{P}(r_{t-1},x_{1:t-1}).$$

$$\mathbb{P}(r_t|r_{t-1})=\{\begin{array}{ll}H(r_{t-1}+1)& r_t=0\\ 1-H(r_{t-1}+1)& r_t=r_{t-1}+1\\ 0& \text{otherwise.}\end{array}$$

$$\mathbb{P}(r_t=0|r_{t-1},x_{t-1}^{(r)})=\{\begin{array}{ll}1-\delta ,& {\tau }^{\ast }=t\text{ and }{\mathfrak{T}}_{GLRT}^{\ast }=1\\ \delta ,& {\mathfrak{T}}_{GLRT}^{\ast }=0\\ H_{const},& \text{otherwise.}\end{array}$$

$$\forall i\in [0,t-1],-{ϵ}_U\le \mathbb{E}[x_t|\varnothing ]-\mathbb{E}[x_t|x_{i:t-1}]\le -{ϵ}_L$$

$$\forall i\in [0,t-1],{ϵ}_L\le \mathbb{E}[x_t|\varnothing ]-\mathbb{E}[x_t|x_{i:t-1}]\le {ϵ}_U$$

$$\frac{{ϵ}_U}{{ϵ}_L}\le {\alpha }_0\left(1+\frac{(1-{\delta }_0^{\text{II}})(1-{\delta }_1^{\text{II}})}{{\delta }_0^{\text{II}}{\delta }_1^{\text{II}}}\right)$$

$$\mathbb{E}[|{\mu }_1-{\mathbb{E}}_{BO}[x_t|x_{1:t-1}]|]\ge \mathbb{E}[|{\mu }_1-{\mathbb{E}}_{CBO}[x_t|$$

$$\forall i\in [0,t-1],-{ϵ}_U\le \mathbb{E}[x_t|x_{1:t-1}]-\mathbb{E}[x_t|x_{i:t-1}]\le -{ϵ}_L,$$

$$\forall i\in [0,t-1],{ϵ}_L\le \mathbb{E}[x_t|x_{1:t-1}]-\mathbb{E}[x_t|x_{i:t-1}]\le {ϵ}_U$$

$$\frac{{ϵ}_U}{{ϵ}_L}\le \frac{{\alpha }_{t-1}((1-{\delta }_0^I)(1-{\delta }_1^I)+{\delta }_0^I{\delta }_1^I)}{{\beta }_{t-1}(1-{\delta }_0^I)(1-{\delta }_1^I)+{\delta }_0^I{\delta }_1^I}$$

$$\mathbb{E}[|{\mu }_2-{\mathbb{E}}_{BO}[x_t|x_{1:t-1}]|]\ge \mathbb{E}[|{\mu }_2-{\mathbb{E}}_{CBO}[x_t|$$

$$\mathbb{P}(X\ge ϵ)\le \exp (-\frac{{ϵ}^2}{2{\sigma }^2}),\text{and }\mathbb{P}(X\le -ϵ)\le \exp (-\frac{{ϵ}^2}{2{\sigma }^2}).$$

$${\lambda }_1(\Sigma )\ge {\lambda }_2(\Sigma )\ge \cdots \ge {\lambda }_n(\Sigma ),\forall k\in [1,n].$$

$$C_{t,(i,j)}=\{\begin{array}{ll}K(i,j),& i,j<t\\ K^{\prime }(i,j),& i,j\ge t\\ 0,& \text{otherwise}\end{array}$$

$$X^T(\Sigma {)}^{-1}X-X^T({\Sigma }_{t^{\ast }}^{{}^{\prime }}{)}^{-1}X=X^T{\Sigma }^{-\frac{1}{2}}{\Sigma }^{\frac{1}{2}}{\Sigma }^{-1}{\Sigma }^{\frac{1}{2}}{\Sigma }^{-\frac{1}{2}}X-X^T{\Sigma }^{-\frac{1}{2}}{\Sigma }^{\frac{1}{2}}({\Sigma }_{t^{\ast }}^{{}^{\prime }}{)}^{-1}{\Sigma }^{\frac{1}{2}}{\Sigma }^{-\frac{1}{2}}X$$

$${X^{\prime }}^T{X}^{\prime }-{X^{\prime }}^T{\Sigma }^{\frac{1}{2}}({\Sigma }_{t^{\ast }}^{{}^{\prime }}{)}^{-1}{\Sigma }^{\frac{1}{2}}{X}^{\prime }$$

$${X^{\prime }}^T{X}^{\prime }-{X^{\prime \prime }}^T\Lambda X^{\prime \prime }=\sum _{i=1}^n{v}_i-\sum _{i=1}^n{\lambda }_i{u}_i$$

$$X^T(\Sigma {)}^{-1}X-X^T({\Sigma }_{t^{\ast }}^{{}^{\prime }}{)}^{-1}X=X^T\Sigma {{}^{\prime }}_{t\ast }^{-\frac{1}{2}}\Sigma {{}^{\prime }}_{t\ast }^{\frac{1}{2}}{\Sigma }^{-1}\Sigma {{}^{\prime }}_{t\ast }^{\frac{1}{2}}\Sigma {{}^{\prime }}_{t\ast }^{-\frac{1}{2}}X-X^T\Sigma {{}^{\prime }}_{t\ast }^{-\frac{1}{2}}\Sigma {{}^{\prime }}_{t\ast }^{\frac{1}{2}}({\Sigma }_{t^{\ast }}^{{}^{\prime }}{)}^{-1}\Sigma {{}^{\prime }}_{t\ast }^{\frac{1}{2}}\Sigma {{}^{\prime }}_{t\ast }^{-\frac{1}{2}}X$$

$${X^{\prime }}^T\Sigma {{}^{\prime }}_{t\ast }^{\frac{1}{2}}{\Sigma }^{-1}\Sigma {{}^{\prime }}_{t\ast }^{\frac{1}{2}}{X}^{\prime }-{X^{\prime }}^T{X}^{\prime }argmaxtZt$$

Full text

1 Introduction

In sequential data, the change point detection (CPD) problem, i.e., analysis of the data with the aim of detecting abrupt changes, is an important component in improving the prediction of future events. When the underlying distribution does not change (stationary), predicting future events becomes tractable. However, the stationarity assumption does not hold in practice.

A change point (CP) is a specific sequential position at which the underlying distribution changes. CPs play critical roles in numerous real-world applications, including climate modeling Manogaran and Lopez (2018), speech recognition Panda and Nayak (2016), image analysis Tewkesbury et al. (2015), and human activity recognition Cleland et al. (2014). In econometrics, structural breaks, which essentially apply CPD to regression models, have been studied for decades in an attempt to identify structural stability in the forecasting of time series Chow (1960); Ewing and Malik (2016). Trend filtering determines CPs by assuming piecewise linearity in the sequential data Kim et al. (2009). CPD also plays an important role in the domain adaptation problem, where it is known as a covariate shift Sugiyama et al. (2008), as the distribution of the test data changes from that of the training data.

Most existing CPD methods are based on either statistical (Bayesian) inference or hypothesis tests. Statistical inference methods compute the probability of the occurrence of a change. Bayesian CPD algorithms Barry and Hartigan (1993); Xuan and Murphy (2007) identify CPs using Bayesian framework. BOCPD algorithm Adams and MacKay (2007); Garnett et al. (2010); Saatçi et al. (2010); Kim and Choi (2015) detects CPs sequentially considering the correlated interval, the so-called run length, between CPs. Such probabilistic methods, however, do not provide a statistical error bound leading to a lack of reliability and are highly sensitive to selected hyperparameters.

Hypothesis-test-based approaches determine the existence of changes based on a statistical test, where the error probability is naturally determined during the computation. Representative hypothesis test techniques include kernel methods such as two-sample tests based on the maximum mean discrepancy Gretton et al. (2012); Li et al. (2015), kernel Fisher discriminant ratio Eric et al. (2008), and likelihood ratio tests such as the cumulative sum (CUSUM) test Chernoff and Zacks (1964); Gombay et al. (1996).

While the conventional GPs only deal with stationary (globally smooth) functions, GPs with CPs can mimic locally smooth functions Saatçi et al. (2010), allowing them to represent many real-world time series data. There have been a number of studies of CPD in GPs using hypothesis tests. For instance, CPD in GPs using p-value test on the likelihood given the null hypothesis of stationary time series have been investigated Isupova (2017). If the null distribution is not legitimate, this approach cannot decisively say there is a change even if the null hypothesis is rejected. Another work has studied to detect the mean changes in GPs through likelihood ratio test Keshavarz et al. (2018).

In this paper, we propose new statistical likelihood ratio tests that detect changes in the covariance structure of GPs and build a theoretically justified online detection algorithm, Confirmatory BOCPD, which detects CPs with a reasonable time delay. We also present sufficient conditions under which CBOCPD provides the lower prediction error compared to BOCPD. Moreover, our algorithm adjusts the parameter of BOCPD to avoid false alarms and missed detections when the results of hypothesis tests are sound. When the results are not sound, our algorithm takes advantages of Bayesian inference from BOCPD.

2 Background

2.1 GP Models and Assumptions

A GP is a random process produced by a collection of random variables such that any finite set of those variables follows a multivariate Gaussian distribution. A GP is completely specified by the mean function $\mu(\cdot)=\mathbb{E}[f(\cdot)]$ and the kernel function $K(\cdot,\cdot)=\text{Cov}(f(\cdot),f(\cdot))$. The kernel function describes distinctive characteristics of the sequential data, such as variance, length scales, and periodicity. For example, in the case of the Radial Basis Function kernel $K(t,t^{\prime})=\sigma^{2}\exp(-\frac{(t-t^{\prime})^{2}}{2l^{2}})$, the length scale hyperparameter $l$ controls the smoothness of the sequential data.

In this paper, we assume that observations are corrupted by the white Gaussian noise with the variance $\sigma_{no}^{2}$. In modeling sequential data using GPs, the index $t$ is considered as the input, and the corresponding output is modeled as $x_{t}\sim N(f(t),\sigma_{no}^{2})$. Given the GP hyperparameters $\theta_{m}$ and $\sigma_{no}$, the log marginal likelihood of the GP over the observed samples $\mathbf{x}$ is specified by the mean vector $\mathbf{\mu}$ and the covariance matrix $\Sigma$ for the multivariate Gaussian distribution.222$\log p(\mathbf{x}|\theta_{m},\sigma_{no})=-\frac{1}{2}(\mathbf{x}{-}\mathbf{\mu})^{T}(\Sigma+\sigma_{no}^{2}I)^{-1}(\mathbf{x}{-}\mathbf{\mu})-\frac{1}{2}\log|\Sigma+\sigma_{no}^{2}I|-\frac{n}{2}\log 2\pi$, where $n$ is the number of observed samples.

2.2 Optimal CPD of the Mean

Here, we briefly review some results on the detection of a single CP in the mean function of a GP Keshavarz et al. (2018). The goal is to set a formal statistical test that can identify a change in the mean function and to quantify the confidence of a change.

We denote the $n$ observed sequential data as $X=\{\it{X}_{t}\}_{t=\text{1}}^{n}$. Let $t\in\mathcal{C}_{n}\subseteq\{1,...,n\}$ represents the point at which a sudden change occurs. Given sequential points, two hypotheses are used for the likelihood ratio test. One is the null hypothesis, $\mathbb{H}_{0}$, which insists there is no CP. The other is the alternative hypothesis, $\mathbb{H}_{1}$, which declares there is at least one CP. The hypothesis testing problem is constructed as follows. Given the two hypotheses, the likelihood ratio is defined as $2\mathfrak{L}=2(\text{sup}_{\theta\in\Theta_{1}}\ell(\theta_{1})-\text{sup}_{\theta\in\Theta_{0}}\ell(\theta_{0}))$ where $\ell$ is the log likelihood function and $\Theta_{0}$ and $\Theta_{1}$ are the parameter spaces of $\mathbb{H}_{0}$ and $\mathbb{H}_{1}$ respectively. The generalized likelihood ratio test (GLRT) is formulated as $\mathfrak{T}_{GLRT}=\mathbb{I}\left(2\mathfrak{L}\geq\mathfrak{R}_{n,\delta}\right)$ with a proper threshold $\mathfrak{R}_{n,\delta}$ and the indicator function $\mathbb{I}(\cdot)$. Here, $\delta$ is the upper bound of the corresponding conditional detection error probability $\varphi_{n}(\mathfrak{T})$, which is defined as

[TABLE]

We reject $\mathbb{H}_{0}$ when $\mathfrak{T}_{GLRT}{=}1$, otherwise, we fail to reject $\mathbb{H}_{0}$. In the mean change detection problem in a GP, we assume that the samples are generated by a GP with zero mean under the null hypothesis, which can be stated as $\mathbb{H}_{0}:\mathbb{E}X=\boldsymbol{0}$. Under the associative alternative hypothesis at time $t$, we assume that there is a change of size $b$ in the mean as $\mathbb{H}_{1,t}:\exists\text{ }b\neq 0,\text{ }\mathbb{E}X=\frac{b}{2}\zeta_{t}$, where $\zeta_{t}\in\mathbb{R}^{n}$ is given by $\zeta_{t}(k):=sign(k-t)$ for any $t\in\mathcal{C}_{n}$. Combining all possible cases, the alternative hypothesis can be written as $\mathbb{H}_{1}:\bigcup_{t\in\mathcal{C}_{n}}\mathbb{H}_{1,t}$, which states that there exists at least one CP with jump size $b$. Further, we can write $2\mathfrak{L}$ as

[TABLE]

where $\Sigma$ is the covariance matrix of $X$. Maximizing Equation (2) with respect to $b$ and plugging in the test, we obtain the following test.

[TABLE]

With a suitable threshold $\mathfrak{R}_{n,\delta}$, the error probability is bounded as $\varphi_{n}(\mathfrak{T}_{GLRT})\leq\delta$ under the sufficient condition on $b$ Keshavarz et al. (2018), and $\mathfrak{R}_{n,\delta}$ could be chosen to be

[TABLE]

3 CPD of the Covariance

This section presents our new hypothesis tests to detect CPs in the covariance structure of a GP.

3.1 Motivational Examples

Figure 1 shows how CPD of a covariance structure could affect the quality of a GP regression. The left plot in Figure 1 shows samples from a GP with an intended CP in the middle. The middle plot shows samples from a GP model after the hyperparameters have been learnt using the whole datasets. The right plot shows samples from a GP model whose covariance structure breaks and the hyperparameters have been learnt separately. Figure 1 suggests that fitting nonstationary data to a time-invariant GP results in an imprecise model. GP regression with a structural break in the covariance structure is more expressive and better suited to the analysis of nonstationary data.

3.2 Likelihood Ratio Test

To construct a test for detecting changes in the covariance structure, we define the null hypothesis as $\mathbb{H}_{0}:\mathrm{Cov}(X_{i},X_{j})=K(i,j)$ and the alternative hypothesis as $\mathbb{H}_{1}=\bigcup_{t\in\mathcal{C}_{n}}\mathbb{H}_{1,t}$, with

[TABLE]

where $K$, $K^{\prime}$ and $K^{\prime\prime}$ are the kernel functions. Let $\Sigma$ and $\Sigma_{t}^{\prime}$ denote the covariance matrices for $\mathbb{H}_{0}$ and $\mathbb{H}_{1,t}$, respectively. The likelihood ratio $2\mathfrak{L}$ is written as

[TABLE]

Theorem 3.1**.**

Let $t^{*}=\text{argmin}_{t}\left[X^{T}(\Sigma_{t}^{{}^{\prime}})^{-1}X+\ln|\Sigma_{t}^{{}^{\prime}}|\right]$. Under the null hypothesis, i.e., $X\sim N(\boldsymbol{0},\Sigma)$,

[TABLE]

where $\lambda_{1},...,\lambda_{n}$ are the eigenvalues of $\Sigma^{\frac{1}{2}}(\Sigma_{t^{*}}^{{}^{\prime}})^{-1}\Sigma^{\frac{1}{2}}$ and $u_{i},v_{i}\sim\chi_{1}^{2}$, a chi-squared distribution with degree 1.

Theorem 3.2**.**

Let $t^{*}=\text{argmin}_{t}\left[X^{T}(\Sigma_{t}^{{}^{\prime}})^{-1}X+\ln|\Sigma_{t}^{{}^{\prime}}|\right]$. Under the alternative hypothesis, i.e., $X\sim N(\boldsymbol{0},{\Sigma_{t^{*}}^{{}^{\prime}}})$,

[TABLE]

where $\lambda_{1},...,\lambda_{n}$ are the eigenvalues of $\Sigma{{}^{\prime}}_{t^{*}}^{\frac{1}{2}}\Sigma^{-1}\Sigma{{}^{\prime}}_{t^{*}}^{\frac{1}{2}}$ and $u_{i},v_{i}\sim\chi_{1}^{2}$, a chi-squared distribution with degree 1.

Theorems 3.1 and 3.2 show that the difference of two positive semi-definite quadratic terms can be expressed as a subtraction between a chi-square random variable with $n$ degrees of freedom and a linear combination of independent chi-square random variables having one degree of freedom each. See Appendix333Appendix is available at https://arxiv.org/abs/1905.13168 A for the proofs of Lemmas and Theorems in Sections 3 and 4.

3.3 Tests for the Covariance Structural Break

For the case where the covariance structure breaks into two different kernels, $\mathbb{H}_{0}$ and $\mathbb{H}_{1}$ are similarly defined as in Section 3.2 except that $K^{\prime\prime}(i,j)=0$. The corresponding covariance matrices can then be written as

[TABLE]

Here, $K_{rc}$ for $r,c\in\{a,b\}$ indicates the covariance matrix between $X_{r}$ and $X_{c}$ with kernel $K$ where $X_{a}:=X_{1:t}$ and $X_{b}:=X_{t+1:n}$. We now define a likelihood ratio test as $\mathfrak{T}_{GLRT}=\mathbb{I}\left(2\mathfrak{L}\geq\mathfrak{R}_{\delta}\right)$. For further Lemmas and Theorems, we define constant $C_{t}$ as follows.

Definition 1**.**

For the covariance matrices $\Sigma$ and $\Sigma_{t}^{{}^{\prime}}$, $C_{t}=\frac{1}{\lambda_{\text{min}}(\Sigma)}+\frac{1}{\lambda_{\text{min}}(\Sigma_{t}^{\prime})}$ where $\lambda_{\text{min}}(M)$ represents the smallest eigenvalue of matrix $M$.

Lemma 3.1**.**

$C_{t}\leq C_{0}$* for all $t\in[1,n]$ when $C_{0}=\frac{1}{\lambda_{\text{min}}(\Sigma)}+\frac{1}{\lambda_{\text{min}}(\Sigma)\wedge\lambda_{\text{min}}(\Sigma_{t=n}^{\prime})}$ where $\wedge$ stands for the minimum operator.*

Lemma 3.2**.**

Let $Q_{t}$ be $X^{T}((\Sigma)^{-1}-(\Sigma_{t}^{\prime})^{-1})X$. When $X_{t}$ is bounded with $X_{t}\in[-V,V]$ for all $t\in[1,n]$, $Q_{t}-\mathbb{E}[Q_{t}]$ is $\frac{C_{t}V^{2}n}{2}$-subgaussian.

Lemma 3.2 shows that, for bounded $X_{t}$, the tail probability of the centered quadratic term in the likelihood ratio decays approximately as fast as the Gaussian with zero mean and the same variance as the centered quadratic term.

Lemma 3.3**.**

The probability that the $\mathfrak{T}_{GLRT}$ is correct under the null hypothesis (absence of CP) is at least $1-\delta/2$, i.e., $\mathbb{P}(2\mathfrak{L}\geq\mathfrak{R}_{n,\delta,\mathbb{H}_{0}}|\mathbb{H}_{0})\leq\delta/2,$ for $\mathfrak{R}_{n,\delta,\mathbb{H}_{0}}=\max_{t}\left(n-Tr(\Sigma(\Sigma_{t}^{{}^{\prime}})^{-1})+\ln\left(\frac{|\Sigma|}{|\Sigma_{t}^{{}^{\prime}}|}\right)\right)+C_{0}V^{2}n\;\sqrt[]{0.5\ln(2/\delta)}$.

Lemma 3.4**.**

The probability that the $\mathfrak{T}_{GLRT}$ is correct under the alternative hypothesis (existence of a CP) is at least $1-\delta/2$, i.e., $\mathbb{P}(2\mathfrak{L}\leq\mathfrak{R}_{n,\delta,\mathbb{H}_{1}}|\mathbb{H}_{1})\leq\delta/2,$ for $\mathfrak{R}_{n,\delta,\mathbb{H}_{1}}=\min_{t}\left(Tr(\Sigma_{t}^{{}^{\prime}}(\Sigma)^{-1})-n+\ln\left(\frac{|\Sigma|}{|\Sigma_{t}^{{}^{\prime}}|}\right)\right)-C_{0}V^{2}n\;\sqrt[]{0.5\ln(2/\delta)}$.

Using Lemma 3.2 and the concentration inequality, Lemmas 3.3 and 3.4 show that we can control the type I or type II errors to be below $\delta/2$.

Theorem 3.3**.**

For $\mathfrak{R}_{n,\delta,\mathbb{H}_{0}}$, $\mathfrak{R}_{n,\delta,\mathbb{H}_{1}}$ in Lemmas 3.3 and 3.4, when $\mathfrak{R}_{n,\delta,\mathbb{H}_{1}}\geq\mathfrak{R}_{n,\delta,\mathbb{H}_{0}}$ and $\mathfrak{R}_{n,\delta,\mathbb{H}_{0}}\leq\mathfrak{R}_{\delta}\leq\mathfrak{R}_{n,\delta,\mathbb{H}_{1}}$, the conditional detection error probability is bounded as

[TABLE]

Proof.

It follows directly from Lemmas 3.3 and 3.4. ∎

Using Theorem 3.3, we guarantee that the likelihood ratio test for a general covariance kernel change is statistically correct for any error bound $\delta$ under specified conditions. If we set the threshold to be greater than or equal to the upper epsilon bound of the null distribution, $\mathfrak{R}_{n,\delta,\mathbb{H}_{0}}$, we can guarantee a bounded type I error. If we set the threshold to be less than or equal to the lower epsilon bound of the alternative distribution, $\mathfrak{R}_{n,\delta,\mathbb{H}_{1}}$, we can guarantee a bounded type II error.

The inequalities for $\mathfrak{R}_{n,\delta,\mathbb{H}_{0}}$ and $\mathfrak{R}_{n,\delta,\mathbb{H}_{1}}$ can have three possible cases. If $\mathfrak{R}_{n,\delta,\mathbb{H}_{0}}>\mathfrak{R}_{n,\delta,\mathbb{H}_{1}}$, there is no threshold guaranteeing both type I and type II errors. If $\mathfrak{R}_{n,\delta,\mathbb{H}_{0}}=\mathfrak{R}_{n,\delta,\mathbb{H}_{1}}$, there is only one threshold that can guarantee both type I and type II errors. If $\mathfrak{R}_{n,\delta,\mathbb{H}_{0}}<\mathfrak{R}_{n,\delta,\mathbb{H}_{1}}$, the thresholds that can guarantee both type I and type II errors are indicated by the shaded area in Figure 3 in Appendix A.

4 Confirmatory Bayesian Online CPD

This section presents our algorithmic contribution which improves the performance of Bayesian online CPD (BOCPD) with statistical hypothesis tests.

4.1 Bayesian Online CPD

Many real-world applications require CPs to be detected in an online manner. BOCPD Adams and MacKay (2007) uses Bayesian inference to update the distribution of future data based on the distribution of the run length. We assume that, if we divide data with CPs, the data in different partitions will be independent whereas the data in the same partition will be autocorrelated.

The formulation of BOCPD is as follows. Let $r_{t}$ be the run length at time $t$ and $x_{t}$ be the data observation at time $t$. We denote a set of consecutive discrete observations between time $a$ and $b$ inclusive as $x_{a:b}$. Let ${x_{t}}^{(r_{t})}$ denote the data observations since the most recent CP. We simply use ${x_{t}}^{(r)}$ instead of ${x_{t}}^{(r_{t})}$ when there is no ambiguity. Then, we calculate the distribution of the next data $x_{t+1}$ given the data up to time $t$. Through marginalization of $r_{t}$, this can be written as

[TABLE]

The joint distribution of the run length and data $\scriptstyle\mathbb{P}(r_{t},x_{1:t})$ is written as,

[TABLE]

We can see that $\scriptstyle\mathbb{P}(r_{t},x_{1:t})$ has a recursive form with respect to time. The first term in Equation (6) can be substituted by $\scriptstyle\mathbb{P}(r_{t}|r_{t-1})$ under the assumption that the current run length only depends on the previous run length. Here, the conditional prior444$\scriptstyle H(\tau)$ is the hazard function with $\scriptstyle H(\tau)=\frac{\mathbb{P}_{gap}(g=\tau)}{\sum_{t=\tau}^{\infty}\mathbb{P}_{gap}(g=t)}.$ of $\scriptstyle r_{t}$ is given by

[TABLE]

$\scriptstyle\mathbb{P}_{gap}(g)$ is the a priori probability distribution over the interval between CPs. When $\scriptstyle\mathbb{P}_{gap}(g)$ is a geometric distribution with timescale $\scriptstyle\lambda$, the hazard function becomes constant as $\scriptstyle H(\tau)=1/\lambda$. The second term in Equation (6) can be calculated using a GP. Thus, we can calculate the overall distribution $\scriptstyle\mathbb{P}(x_{t+1}|x_{1:t})$ by recursive message passing scheme with $\scriptstyle\mathbb{P}(r_{t},x_{1:t})$.

The BOCPD framework efficiently determines changes in GP modeling. The change is modeled by considering all the possible run lengths. One of the main assumptions of this framework is that parameters such as the kernel parameters or the hazard function do not change. However, this online CPD framework is vulnerable to the parameters, because inappropriate parameters degrade predictive performance on non-stationary data leading imprecise run length distribution as shown in Figure 2.

4.2 Confirmatory BOCPD

Algorithm 1 presents a theoretically justified online change detection algorithm, CBOCPD. The main idea of CBOCPD is to overcome the limitations of the assumption that the run length is independent of the data. However, we claim that the first term in Equation (6) can be directly calculated by the likelihood ratio test in Section 3.3. The first two lines initialize the parameters. In lines 3–13, Equation (7) is altered as

[TABLE]

Here, we have two likelihood ratio tests, $\scriptstyle\mathfrak{T}_{GLRT}^{0}=\mathbb{I}\left(2\mathfrak{L}\geq\tilde{\mathfrak{R}}_{\delta,\mathbb{H}_{0}}\right)$ and $\scriptstyle\mathfrak{T}_{GLRT}^{1}=\mathbb{I}\left(2\mathfrak{L}\geq\tilde{\mathfrak{R}}_{\delta,\mathbb{H}_{1}}\right)$, with thresholds in line 4. Empirical thresholds are used since the theoretically calculated thresholds in Lemmas 3.3 and 3.4 are not tight enough to use in practice. We further define $\scriptstyle\mathfrak{T}_{GLRT}^{*}=1$ if $\scriptstyle\mathfrak{T}_{GLRT}^{1}=1$, and $\scriptstyle\mathfrak{T}_{GLRT}^{0}=1$ and $\scriptstyle\mathfrak{T}_{GLRT}^{*}=0$ if $\scriptstyle\mathfrak{T}_{GLRT}^{1}=0$ and $\scriptstyle\mathfrak{T}_{GLRT}^{0}=0$. The likelihood ratio tests are applied to the window around $\scriptstyle t$ with $\scriptstyle W=x_{t-m:t+m}$. $\scriptstyle\tau^{*}$ is the time point in the window that maximizes the likelihood ratio, $\scriptstyle\tau^{*}=\text{argmax}_{\tau\in\mathcal{C}_{W}}2\mathfrak{L}_{\tau}$. Here $\scriptstyle\mathcal{C}_{W}\subseteq\{t-m,...,t+m\}$ is a set of change point candidates for the window. If both likelihood ratio tests at time $\scriptstyle\tau^{*}$ pass and $\scriptstyle\tau^{*}$ coincides with $\scriptstyle t$, we decide that $\scriptstyle t$ is a CP and set $\scriptstyle\mathbb{P}(r_{t}=0|r_{t-1},x_{t-1}^{(r)})=1-\delta$, which enhances the probability of change in the BOCPD framework. In contrast, if neither test passes, we strongly believe there is no change and reduce the probability of change in the BOCPD framework. This is why we name this algorithm Confirmatory BOCPD. We add $\scriptstyle\tau^{*}=t$ to avoid situations where the same time point maximizes the likelihood ratio in several consecutive windows, resulting in duplicate CPs. Lines 14–19 follow the BOCPD framework Adams and MacKay (2007).

4.3 Theoretical Analysis of CBOCPD

We present sufficient conditions under which CBOCPD provides the lower prediction error compared to BOCPD. Here, the prediction error is defined as the expected absolute difference between the predictive mean given the true run length and the predictive mean under BOCPD and CBOCPD at a (detected) CP $\scriptstyle t$, as a good prediction at a CP is a key factor in overall performance. For simplicity, we denote the expected value of $\scriptstyle x_{t}$ under BOCPD and CBOCPD as $\scriptstyle\mathbb{E}_{BO}[x_{t}|x_{1:t-1}]$ and $\scriptstyle\mathbb{E}_{CBO}[x_{t}|x_{1:t-1}]$, respectively. Further we define $\scriptstyle\alpha_{i}=1-\mathbb{P}_{BO}(r_{t-1}=i|x_{1:t-1})$ under BOCPD and $\scriptstyle\beta_{i}=1-\mathbb{P}_{CBO}(r_{t-1}=i|x_{1:t-1})$ under CBOCPD.

4.3.1 Non-stationary Case

We first investigate the case where a change exists.

Theorem 4.1**.**

Consider BOCPD (Sec. 4.1) and CBOCPD (Sec. 4.2) algorithms where CBOCPD uses two statistical tests $\scriptstyle\mathfrak{T}_{GLRT}^{0}$ and $\scriptstyle\mathfrak{T}_{GLRT}^{1}$ with type II error of $\scriptstyle\delta_{0}^{\text{II}}$ and $\scriptstyle\delta_{1}^{\text{II}}$ respectively. When there exists a CP at $\scriptstyle t$ with mean value of $\scriptstyle\mu_{1}$ which satisfies either

[TABLE]

or

[TABLE]

for $\scriptstyle 0<\epsilon_{L}<\epsilon_{U}$ where $\scriptstyle\mathbb{E}[x_{t}|\emptyset]$ indicates the expected value of $\scriptstyle x_{t}$ with prior distribution, if

[TABLE]

with $\scriptstyle\alpha_{0}$ indicating the probability of the run length not to be zero, then the absolute error of expected value of CBOCPD at $\scriptstyle t$ is less than or equal to the one of BOCPD as

[TABLE]

4.3.2 Stationary Case

Next, we investigate the case where a change does not exist.

Theorem 4.2**.**

Consider BOCPD (Sec. 4.1) and CBOCPD (Sec. 4.2) algorithms where CBOCPD uses two statistical tests $\scriptstyle\mathfrak{T}_{GLRT}^{0}$ and $\scriptstyle\mathfrak{T}_{GLRT}^{1}$ with type I error of $\scriptstyle\delta_{0}^{I}$ and $\scriptstyle\delta_{1}^{I}$ respectively. When there exists a confirmed non-CP at $\scriptstyle t$ with mean value of $\scriptstyle\mu_{2}$ which satisfies either

[TABLE]

or

[TABLE]

for $\scriptstyle 0<\epsilon_{L}<\epsilon_{U}$, if

[TABLE]

with $\scriptstyle\alpha_{t-1}$ and $\scriptstyle\beta_{t-1}$ indicating the probability of the run length not to be $\scriptstyle t-1$ under BOCPD and CBOCPD respectively, then the absolute error of the expected value of CBOCPD at $\scriptstyle t$ is less than or equal to the one of BOCPD as

[TABLE]

5 Experimental Evaluations

5.1 Synthetic Data

In this experiment, we investigate two synthetic datasets generated by GPs with changes in the length scale and the variance of a Radial Basis Function kernel, respectively. Observations are obtained by adding the white Gaussian noise with the variance $\scriptstyle\sigma_{no}^{2}=0.1$. For both datasets, two CPs are drawn uniformly from time intervals $\scriptstyle(75,125)$ and $\scriptstyle(275,325)$ with end time $\scriptstyle T=400$.

For the first experiment, the length scale $\scriptstyle l$ switches from $\scriptstyle 3$ to $\scriptstyle 20$ and from $\scriptstyle 20$ to $\scriptstyle 1$ with the variance $\scriptstyle\sigma^{2}=1$. For the second experiment, the variance $\scriptstyle\sigma^{2}$ switches from $\scriptstyle 1$ to $\scriptstyle 4$ and from $\scriptstyle 4$ to $\scriptstyle 0.3$ with the length scale $\scriptstyle l=3$. We compare the proposed CBOCPD with BOCPD and CS (CUSUM)-BOCPD. For BOCPD, we use a modified version of the stationary GP-BOCPD Saatçi et al. (2010) by fixing the timescale $\scriptstyle\lambda$. CS-BOCPD uses CUSUM instead of the our proposed likelihood ratio test in Algorithm 1. For all algorithms, the timescale of $\scriptstyle\lambda=200$ is used.

Figure 2(a) and Figure 2(b) show that CBOCPD identifies the length scale change in the data with the help of statistical tests, whereas BOCPD captures the change too less or too many times. The first two columns in Table 1 present the quantitative results from 100 runs of each algorithm after training on the first 100 points and testing on the remaining 300 points. CBOCPD outperforms BOCPD in terms of both negative log-likelihood (NLL) and mean squared error (MSE). The differences in performances become larger when there is a stronger correlation among the samples.

5.2 Gazebo Robot Simulation Data

We also conducted experiments on the Gazebo robot simulator to detect changes in the environment of a robot. We gathered data by moving the robot through a changing environment (i.e., varying ground properties). We used the Pioneer3AT robot with kinetic ROS. There are three environments. In the first environment (Env1), the ground changes at the midpoint from ‘Plane ground’ to ‘Bumpy ground 1’. In the second environment (Env2), the robot moves from ‘Bumpy ground 1’ to ‘Bumpy ground 2’ where ‘Bumpy ground 2’ is more coarser than ‘Bumpy ground 1’. Finally, in the third environment (Env3), the robot moves from ‘Bumpy ground 2’ back to the ‘Plane ground’. For simplicity, we used only z-directional position data as the inputs. The last three columns in Table 1 show that CBOCPD outperforms other methods in all environments. Interestingly, the results show that CUSUM test does not help BOCPD much when the variance decreases as in Env3 but proposed likelihood ratio test improves BOCPD in all the cases.

5.3 Real World Data

Further, we compare our proposed algorithm with autoregressive GP (ARGP), autoregressive GP with BOCPD (ARGP-BOCPD), GP time series with BOCPD (BOCPD) and GP time series model (GPTS). For baseline methods we used the publicly available source codes555https://sites.google.com/site/wwwturnercomputingcom/software/ThesisCodeAndData.zip.

We first consider Nile data which records the lowest annual water levels on the Nile river during AD 622-1284 measured at the island of Roda. There is domain knowledge suggesting a CP in the year 715 due to the construction of a ‘nilometer’ Popper (1951); Balek (1977). Secondly, we used Well Log data which records 4050 nuclear magnetic resonance measurements while drilling a well. We also considered Snowfall data which records daily snowfall data in Whistler, BC, Canada. In this dataset each method is evaluated by the ability to predict snowfall of next day using 35 years of test data.

Table 2 shows the predictive performance of our proposed algorithm compared to other GP based BOCPD methods. In Nile data, we see that combining BOCPD method with GP improves the performance. CBOCPD further improves the performance by leveraging the proposed statistical hypothesis tests and outperforms all other algorithms. In Well Log data, the nonlinear temporal correlations within each regime give a slight advantage to underlying predictive model of ARGP. However, CBOCPD still shows the competitive result. In Snowfall data, CBOCPD well detects the difference in noise levels and achieves the best performance in MSE.

6 Conclusion

This paper presents a novel framework for detecting changes in the covariance structure of the GP. We propose a semi-online CPD algorithm, Confirmatory BOCPD, which is an improved version of BOCPD with statistical hypothesis tests. Experiments using synthetic and real-world datasets demonstrate that CBOCPD outperforms conventional BOCPDs.

Acknowledgments

This work was supported by the Industrial Convergence Core Technology Development Program funded by the Ministry of Trade, Industry & Energy (MOTIE), Korea (No. 10063172, Development of robot intelligence technology for mobility with learning capability toward robust and seamless indoor and outdoor autonomous navigation) and Institute for Information & communications Technology Planning & Evaluation (IITP) grant funded by the Ministry of Science and ICT (MSIT), Korea (No. 2017-0-01779, A machine learning and statistical inference framework for explainable artificial intelligence).

Appendix A Proofs for Sections 3 and 4

Definition 2** (Subgaussianity).**

A random variable $\scriptstyle X$ is $\scriptstyle\sigma$-subgaussian if for all $\scriptstyle\lambda\in\mathbb{R}$ it holds that $\scriptstyle\mathbb{E}[\exp(\lambda X)]\leq\exp(\lambda^{2}\sigma^{2}/2)$.

Lemma A.1**.**

If $\scriptstyle X$ is a bounded and centered random variable, with $\scriptstyle X\in[a,b]$, then $\scriptstyle X$ is $\scriptstyle\frac{b-a}{2}$-subgaussian.

Please refer to Ch. 5.2 in Lattimore and Szepesvári (2018).

Theorem A.1** (Chernoff Bound).**

If $\scriptstyle X$ is $\scriptstyle\sigma$-subgaussian, then for any $\scriptstyle\epsilon\geq 0$,

[TABLE]

Theorem A.2** (Cauchy Interlace Theorem).**

Let $\scriptstyle A$ be a Hermitian matrix of order $\scriptstyle n$, and let $\scriptstyle B$ be a principal submatrix of $\scriptstyle A$ of order $\scriptstyle n-1$, and denote the $\scriptstyle k$th largest eigenvalues of a general $\scriptstyle n\times n$ matrix by $\scriptstyle\lambda_{k}(\Sigma)$ so that

[TABLE]

If $\scriptstyle\lambda_{n}(A)\leq\lambda_{n-1}(A)\leq\dots\leq\lambda_{2}(A)\leq\lambda_{1}(A)$ and $\scriptstyle\lambda_{n-1}(B)\leq\lambda_{n-2}(B)\leq\dots\leq\lambda_{2}(B)\leq\lambda_{1}(B)$, then $\scriptstyle\lambda_{n}(A)\leq\lambda_{n-1}(B)\leq\lambda_{n-1}(A)\leq\lambda_{n-2}(B)\leq\cdots\leq\lambda_{2}(A)\leq\lambda_{1}(B)\leq\lambda_{1}(A)$.

Corollary A.1**.**

Let $\scriptstyle A,B,C_{t}\in\mathbb{R}^{n\times n}$ be positive semi-definite matrix defined element wise as $\scriptstyle A_{i,j}=K(i,j)$, $\scriptstyle B_{i,j}=K^{\prime}(i,j)$ and $\scriptstyle C_{t}$ is defined as

[TABLE]

for all $\scriptstyle t\in[1,n]$. Then, $\scriptstyle\forall 1\leq k\leq n,(\lambda_{n}(A)\wedge\lambda_{n}(B))\leq\lambda_{k}(C_{t})\leq(\lambda_{1}(A)\vee\lambda_{1}(B))$, where $\scriptstyle\wedge$ and $\scriptstyle\vee$ stand for minimum and maximum operators.

Proof of Lemma 3.1.

The proof follows from Corollary A.1. ∎

Proof of Theorem 3.1.

The quadratic forms in $\scriptstyle 2\mathfrak{L}$ can be reorganized as

[TABLE]

Letting $\scriptstyle X^{\prime}=\Sigma^{-\frac{1}{2}}X$, one has $\scriptstyle X^{\prime}\sim N(\boldsymbol{0},I)$ since $\scriptstyle X\sim N(\boldsymbol{0},\Sigma)$, and then it gives:

[TABLE]

The matrix $\scriptstyle\Sigma^{\frac{1}{2}}(\Sigma_{t^{*}}^{{}^{\prime}})^{-1}\Sigma^{\frac{1}{2}}$ is symmetric as both matrices $\scriptstyle\Sigma$ and $\scriptstyle\Sigma_{t^{*}}^{{}^{\prime}}$ are symmetric. Performing eigenvalue decomposition gives $\scriptstyle\Sigma^{\frac{1}{2}}(\Sigma_{t^{*}}^{{}^{\prime}})^{-1}\Sigma^{\frac{1}{2}}=Q\Lambda Q^{T}$ where $\scriptstyle Q$ is an orthogonal matrix and diagonal elements of $\scriptstyle\Lambda$ is non-negative since $\scriptstyle\Sigma$ and $\scriptstyle\Sigma_{t^{*}}^{{}^{\prime}}$ are positive semi-definite matrices. Letting $\scriptstyle X^{\prime\prime}=Q^{T}X^{\prime}$, one has $\scriptstyle X^{\prime\prime}\sim N(\boldsymbol{0},I)$ since $\scriptstyle Q$ is an orthogonal matrix, and then the equation is further transformed as:

[TABLE]

where $\scriptstyle\lambda_{1},...,\lambda_{n}$ are the eigenvalues of $\scriptstyle\Lambda$ and $\scriptstyle u_{i},v_{i}\sim\chi_{1}^{2}$, a chi-squared distribution with degree 1. The last equality holds as Q is orthonormal.

∎

Proof of Theorem 3.2.

The proof proceeds in a similar manner as Theorem 3.1 except for that $\scriptstyle X\sim N(\boldsymbol{0},{\Sigma_{t^{*}}^{{}^{\prime}}})$. The quadratic forms in $\scriptstyle 2\mathfrak{L}$ can be reorganized as

[TABLE]

Letting $\scriptstyle X^{\prime}=\Sigma{{}^{\prime}}_{t*}^{-\frac{1}{2}}X$, one has $\scriptstyle X^{\prime}\sim N(\boldsymbol{0},I)$ since $\scriptstyle X\sim N(\boldsymbol{0},\Sigma{{}^{\prime}}_{t*})$, and then it gives:

[TABLE]

The matrix $\scriptstyle\Sigma{{}^{\prime}}_{t*}^{\frac{1}{2}}\Sigma^{-1}\Sigma{{}^{\prime}}_{t*}^{\frac{1}{2}}$ is symmetric as both matrices $\scriptstyle\Sigma$ and $\scriptstyle\Sigma_{t^{*}}^{{}^{\prime}}$ are symmetric. Performing eigenvalue decomposition gives $\scriptstyle\Sigma{{}^{\prime}}_{t*}^{\frac{1}{2}}\Sigma^{-1}\Sigma{{}^{\prime}}_{t*}^{\frac{1}{2}}=Q\Lambda Q^{T}$ where $\scriptstyle Q$ is an orthogonal matrix and diagonal elements of $\scriptstyle\Lambda$ is non-negative since $\scriptstyle\Sigma$ and $\scriptstyle\Sigma_{t^{*}}^{{}^{\prime}}$ are positive semi-definite matrices. Letting $\scriptstyle X^{\prime\prime}=Q^{T}X^{\prime}$, one has $\scriptstyle X^{\prime\prime}\sim N(\boldsymbol{0},I)$ since $\scriptstyle Q$ is an orthogonal matrix, and then the equation is further transformed as:

[TABLE]

where $\scriptstyle\lambda_{1},...,\lambda_{n}$ are the eigenvalues of $\scriptstyle\Lambda$ and $\scriptstyle u_{i},v_{i}\sim\chi_{1}^{2}$, a chi-squared distribution with degree 1.

∎

Proof of Lemma 3.2.

We will first show the bound of a quadratic term, $\scriptstyle X^{T}(\Sigma^{-1}-(\Sigma_{t}^{\prime})^{-1})X$. For a positive semi-definite matrix $\scriptstyle M$, its eigenvectors are orthonormal and all of its eigenvalues are non-negative. By performing eigenvalue decomposition

[TABLE]

where $\scriptstyle\lambda_{k}$ is the $\scriptstyle k$th diagonal element of $\scriptstyle\Lambda$ and $\scriptstyle q_{k}$ is the $\scriptstyle k$th column of $\scriptstyle Q$. Since eigenvectors are orthonormal

[TABLE]

By the restriction that $\scriptstyle X_{t}$ is bounded with $\scriptstyle X_{t}\in[-V,V]$ for all $\scriptstyle t$,

[TABLE]

where $\scriptstyle\lambda_{\text{min}}(M)$ and $\scriptstyle\lambda_{\text{max}}(M)$ represent the smallest and biggest eigenvalues of $\scriptstyle M$ respectively.

As $\scriptstyle\Sigma$ is symmetric, its inverse $\scriptstyle(\Sigma)^{-1}$ is also symmetric and positive semi-definite. Substituting $\scriptstyle M$ by $\scriptstyle\Sigma^{-1}$ gives

[TABLE]

because eigenvalues of the inverse matrix are inverses of eigenvalues of an original matrix.

Similarly for $\scriptstyle(\Sigma_{t}^{\prime})^{-1}$,

[TABLE]

Then the difference of the quadratic term of our interest is bounded as

[TABLE]

with $\scriptstyle C_{t}=\frac{1}{\lambda_{\text{min}}(\Sigma)}+\frac{1}{\lambda_{\text{min}}(\Sigma_{t}^{\prime})}$. Thus, we can conclude proof by Lemma A.1. ∎

Proof of Lemma 3.3.

For brevity, let $\scriptstyle Z_{t}=X^{T}(\Sigma^{-1}-{\Sigma^{\prime}_{t}}^{-1})X+\ln\left(\frac{|\Sigma|}{|\Sigma_{t}^{{}^{\prime}}|}\right),\;t=1,...,n$. Lemma 3.2 implies that $\scriptstyle Z_{t}-\mathbb{E}[Z_{t}]$ is $\scriptstyle\frac{C_{t}V^{2}n}{2}$-subgaussian. Under the null hypothesis the expectation of $\scriptstyle Z_{t}$ is defined as

[TABLE]

Theorem A.1 implies that $\scriptstyle\mathbb{P}(Z_{t}\geq n-Tr(\Sigma(\Sigma_{t}^{{}^{\prime}})^{-1})+\ln\left(\frac{|\Sigma|}{|\Sigma_{t}^{{}^{\prime}}|}\right)+C_{t}V^{2}n\;\sqrt[]{0.5\ln(2n/\delta)})\leq\frac{\delta}{2n}$. By Lemma 3.1, $\scriptstyle C_{t}\leq C_{0}$ for all $\scriptstyle t\in[1,n]$. Thus,

[TABLE]

As the above inequality is satisfied for all $\scriptstyle t\in[1,n]$, it is also satisfied for $\scriptstyle t*=\text{argmax}_{t}Z_{t}$. Thus we can conclude as follows.

[TABLE]

∎

Proof of Lemma 3.4.

For brevity, let $\scriptstyle Z_{t}=X^{T}(\Sigma^{-1}-{\Sigma^{\prime}_{t}}^{-1})X+\ln\left(\frac{|\Sigma|}{|\Sigma_{t}^{{}^{\prime}}|}\right),\;t=1,...,n$. Lemma 3.2 implies that $\scriptstyle Z_{t}-\mathbb{E}[Z_{t}]$ is $\scriptstyle\frac{C_{t}V^{2}n}{2}$-subgaussian. Under the alternative hypothesis the expectation of $\scriptstyle Z_{t}$ is defined as

[TABLE]

Theorem A.1 implies that $\scriptstyle\mathbb{P}(Z_{t}\leq Tr(\Sigma_{t}^{{}^{\prime}}(\Sigma)^{-1})-n+\ln\left(\frac{|\Sigma|}{|\Sigma_{t}^{{}^{\prime}}|}\right)-C_{t}V^{2}n\;\sqrt[]{0.5\ln(2/\delta)})\leq\frac{\delta}{2}$. By Lemma 3.1, $\scriptstyle C_{t}\leq C_{0}$ for all $\scriptstyle t\in[1,n]$. Thus,

[TABLE]

As the above inequality is satisfied for all $\scriptstyle t\in[1,n]$, it is also satisfied for $\scriptstyle t*=\text{argmax}_{t}Z_{t}$. Thus we can conclude as follows.

[TABLE]

∎

Proof of Theorem 4.1.

Let’s define the gain of CBOCPD over BOCPD as

[TABLE]

In the case $\scriptstyle\mathfrak{T}_{GLRT}^{*}=1$, the gain is written as

[TABLE]

The inequality in the last line comes from the assumption $\scriptstyle\epsilon_{L}\leq|\mathbb{E}[x_{t}|\emptyset]-\mathbb{E}[x_{t}|x_{i:t-1}]|\leq\epsilon_{U}$.

In the case $\scriptstyle\mathfrak{T}_{GLRT}^{*}=0$, the loss of CBOCPD is written as

[TABLE]

The equation in the last line comes from the fact that $\scriptstyle\Sigma_{r_{t-1}}^{t-1}\mathbb{P}_{CBO}(r_{t-1}|x_{1:t-1})=1$ under the CBOCPD when non-change is detected. Then, the gain is bounded as

[TABLE]

As $\scriptstyle\mathbb{P}(\mathfrak{T}_{GLRT}^{*}=1)=(1-\delta_{0}^{\text{II}})(1-\delta_{1}^{\text{II}})$ and $\scriptstyle\mathbb{P}(\mathfrak{T}_{GLRT}^{*}=0)=\delta_{0}^{\text{II}}\delta_{1}^{\text{II}}$ in non-stationary case, the expected gain is bounded from below as

[TABLE]

where the last inequality follows the assumption. Thus we can conclude that the expected gain is non-negative.

∎

Proof of Theorem 4.2.

Let’s define the gain of CBOCPD over BOCPD as

[TABLE]

The loss of BOCPD is written as

[TABLE]

In the case $\scriptstyle\mathfrak{T}_{GLRT}^{*}=1$, the loss of CBOCPD is written as

[TABLE]

Then, the gain is bounded as

[TABLE]

In the case $\scriptstyle\mathfrak{T}_{GLRT}^{*}=0$, the loss of CBOCPD is written as

[TABLE]

Then, the gain is bounded as

[TABLE]

As $\scriptstyle\mathbb{P}(\mathfrak{T}_{GLRT}^{*}=1)=\delta_{0}^{\text{I}}\delta_{1}^{\text{I}}$ and $\scriptstyle\mathbb{P}(\mathfrak{T}_{GLRT}^{*}=0)=(1-\delta_{0}^{\text{I}})(1-\delta_{1}^{\text{I}})$ in stationary case, the expected gain is bounded from below as

[TABLE]

where the last inequality follows the assumption. Thus we can conclude that the expected gain is non-negative.

∎