Sequential Change-Point Detection via Online Convex Optimization

Abstract

Sequential change-point detection when the distribution parameters are unknown is a fundamental problem in statistics and machine learning. When the post-change parameters are unknown, we consider a set of detection procedures based on sequential likelihood ratios with non-anticipating estimators constructed using online convex optimization algorithms such as online mirror descent, which provides a more versatile approach to tackling complex situations where recursive maximum likelihood estimators cannot be found. When the underlying distributions belong to a exponential family and the estimators satisfy the logarithm regret property, we show that this approach is nearly second-order asymptotically optimal. This means that the upper bound for the false alarm rate of the algorithm (measured by the average-run-length) meets the lower bound asymptotically up to a log-log factor when the threshold tends to infinity. Our proof is achieved by making a connection between sequential change-point and online convex optimization and leveraging the logarithmic regret bound property of online mirror descent algorithm. Numerical and real data examples validate our theory.

1. Introduction

Sequential analysis is a classic topic in statistics concerning online inference from a sequence of observations. The goal is to make statistical inference as quickly as possible, while controlling the false-alarm rate. An important sequential analysis problem commonly studied is sequential change-point detection [1]. It arises from various applications including online anomaly detection, statistical quality control, biosurveillance, financial arbitrage detection and network security monitoring (see, e.g., [2,3,4]).

We are interested in the sequential change-point detection problem with known pre-change parameters but unknown post-change parameters. Specifically, given a sequence of samples $X_1$, $X_2$, …, we assume that they are independent and identically distributed (i.i.d.) with certain distribution $f_{\theta }$ parameterized by $\theta$, and the values of $\theta$ are different before and after some unknown time called the change-point. We further assume that the parameters before the change-point are known. This is reasonable since usually it is relatively easy to obtain the reference data for the normal state, so that the parameters in the normal state can be estimated with good accuracy. After the change-point, however, the values of the parameters switch to some unknown values, which represent anomalies or novelties that need to be discovered.

1.1. Motivation: Dilemma of CUSUM and Generalized Likelihood Ratio (GLR) Statistics

Consider change-point detection with unknown post-change parameters. A commonly used change-point detection method is the so-called CUSUM procedure [4] that can be derived from likelihood ratios. Assume that before the change, the samples $X_i$ follow a distribution $f_{{\theta }_0}$ and after the change the samples $X_i$ follow another distribution $f_{{\theta }_1}$. CUSUM procedure has a recursive structure: initialized with $W_0=0$, the likelihood-ratio statistic can be computed according to $W_{t+1}=max\{W_t+log(f_{{\theta }_1}\left(X_{t+1}\right)/f_{{\theta }_0}\left(X_{t+1}\right)),0\}$, and a change-point is detected whenever $W_t$ exceeds a pre-specified threshold. Due to the recursive structure, CUSUM is memory and computation efficient since it does not need to store the historical data and only needs to record the value of $W_t$. The performance of CUSUM depends on the choice of the post-change parameter ${\theta }_1$; in particular, there must be a well-defined notion of “distance” between ${\theta }_0$ and ${\theta }_1$. However, the choice of ${\theta }_1$ is somewhat subjective. Even if in practice a reasonable choice of ${\theta }_1$ is the “smallest” change-of-interest, in the multi-dimensional setting, it is hard to define what the “smallest” change would mean. Moreover, when the assumed parameter ${\theta }_1$ deviates significantly from the true parameter value, CUSUM may suffer a severe performance degradation [5].

An alternative approach is the Generalized Likelihood Ratio (GLR) statistic based procedure [6]. The GLR statistic finds the maximum likelihood estimate (MLE) of the post-change parameter and plugs it back into the likelihood ratio to form the detection statistic. To be more precise, for each hypothetical change-point location k, the corresponding post-change samples are $\{X_{k+1},\dots ,X_t\}$. Using these samples, one can form the MLE denoted as ${\widehat{\theta }}_{k+1,t}$. Without knowing whether the change occurs and where it occurs beforehand when forming the GLR statistic, we have to maximize k over all possible change locations. The GLR statistic is given by ${max}_{k<t}{\sum }_{i=k+1}^{t}log(f_{{\widehat{\theta }}_{k,t}}\left(X_i\right)/f_{{\theta }_0}\left(X_t\right))$, and a change is announced whenever it exceeds a pre-specified threshold. The GLR statistic is more robust than CUSUM [7], and it is particularly useful when the post-change parameter may vary from one situation to another. In simple cases, the MLE ${\widehat{\theta }}_{k+1,t}$ may have closed-form expressions and may be evaluated recursively. For instance, when the post-change distribution is Gaussian with mean $\theta$ [8], ${\widehat{\theta }}_{k+1,t}=\left({\sum }_{i=k+1}^t{X}_i\right)/(t-k)$, and ${\widehat{\theta }}_{k+1,t+1}=(t-k)/(t-k+1)·{\widehat{\theta }}_{k+1,t}+X_{t+1}/(t-k+1)$. However, in more complex situations, in general MLE ${\widehat{\theta }}_{k+1,t}$ does not have recursive form and cannot be evaluated using simple summary statistics. One such instance is given in Section 1.2. Another instance is when there is a constraint on the MLE such as sparsity. In these cases, one has to store historical data and recompute the MLE ${\widehat{\theta }}_{k,t}$ whenever there is new data, which is not memory efficient nor computational efficient. For these cases, as a remedy, the window-limited GLR is usually considered, where only the past w samples are stored and the maximization is restricted to be over $k\in (t-w,t]$. However, even with the window-limited GLR, one still has to recompute ${\widehat{\theta }}_{k,t}$ using historical data whenever the new data are added.

Besides CUSUM or GLR, various online change-point detection procedures using one-sample updates have been considered, which replace with the MLE with a simple recursive estimator. The one-sample update estimate takes the form of ${\widehat{\theta }}_{k,t}=h(X_t,{\widehat{\theta }}_{k,t-1})$ for some function h that uses only the most recent data and the previous estimate. Then the estimates are plugged into the likelihood ratio statistic to perform detection. Online convex optimization algorithms (such as online mirror descent) are natural approach to construct these estimators (see, e.g., [9,10]). Such a scheme provides a more versatile approach to developing a detecting procedure for complex situations, where the exact MLE does not have a recursive form or even a closed-form expression. The one-sample update enjoys efficient computation, as information from the new data can be incorporated via low computational cost update. It is also memory efficient since the update only needs the most recent sample. The one sample update estimators may not correspond to the exact MLE, but they tend to result in good detection performance. However, in general there is no performance guarantees for such an approach. This is the question we aim to address in this paper.

1.2. Application Scenario: Social Network Change-Point Detection

The widespread use of social networks (such as Twitter) leads to a large amount of user-generated data generated continuously. One important aspect is to detect change-points in streaming social network data. These change-points may represent the collective anticipation of response to external events or system “shocks” [11]. Detecting such changes can provide a better understanding of patterns of social life. In social networks, a common form of the data is discrete events over continuous time. As a simplification, each event contains a time label and a user label in the network. In our prior work [12], we model discrete events using network point processes, which capture the influence between users through an influence matrix. We then cast the problem as detecting changes in an influence matrix, assuming that the influence matrix in the normal state (before the change) can be estimated from the reference data. After the change, the influence matrix is unknown (since it represents an anomaly) and has to be estimated online. Due to computational burden and memory constraint, since the scale of the network tends to be large, we do not want to store the entire historical data and rather compute the statistic in real-time. A simulated example to illustrate this case is shown later in Section 4.4.

1.3. Contributions

This paper has two main contributions. First, we present a general approach based on online convex optimization (OCO) for constructing the estimator for the one-sided sequential hypothesis test and the sequential change-point detection, in the non-anticipative approach of [8] if the MLE cannot be computed in a convenient recursive form.

Second, we provide a proof of the near second-order asymptotic optimality of this approach when a “logarithmic regret property” is satisfied and when the distributions are from an exponential family. The nearly second-order asymptotic optimality [4] means that the upper bound for performance matches the lower bound up to a log-log factor as the false-alarm rate tends to zero. Inspired by the existing connection between sequential analysis and online convex optimization in [13,14], we prove the near optimality leveraging the logarithmic regret property of online mirror descent (OMD) and the lower bound established in statistical sequential change-point literature [4,15]. More precisely, we provide a general upper bound for the one-sided sequential hypothesis test and change-point detection procedures with the one-sample update schemes. The upper bound explicitly captures the impact of estimation on detection by an estimation algorithm dependent factor. This factor shows up as an additional term in the upper bound for the expected detection delay, and it corresponds to the regret incurred by the one-sample update estimators. This establishes an interesting linkage between sequential change-point detection and online convex optimization. Although both fields, sequential change-point detection and online convex optimization, study sequential data, the precise connection between them is not clear, partly because the performance metrics are different: the former is concerned with the tradeoff between average run length and detection delay, whereas the latter focuses on bounding the cumulative loss incurred by the sequence of estimators through a regret bound [14,16]. Synthetic examples validate the performances of one sample update schemes. Here we focus on OMD estimators, but the results can be generalized to other OCO schemes such as the online gradient descent.

1.4. Literature and Related Work

Sequential change-point detection is a classic subject with an extensive literature. Much success has been achieved when the pre-change and post-change distributions are exactly specified. For example, the CUSUM procedure [17] with first-order asymptotic optimality [18] and exact optimality [19] in the minimax sense, and the Shiryayev-Roberts (SR) procedure [20] derived based on Bayesian principle that also enjoys various optimality properties. Both CUSUM and SR procedures rely on likelihood ratios between the specified pre-change and post-change distributions.

There are two main approaches in dealing with the unknown post-change parameters. The first one is a GLR approach [7,21,22,23,24], and the second is a mixture approach [15,25]. The GLR statistic enjoys certain optimality properties, but it can not be computed recursively in many cases [23]. To address the infinite memory issue, [7,21] studied the window-limited GLR procedure. The main advantage of the mixture approach is that it allows an easy evaluation of a threshold that guarantees the desired false alarm constraint. A disadvantage of this approach is that sometimes there may not be a natural way of selecting the weight function, in particular when there is no conjugate prior. This motivated a third approach to this problem, which was proposed first by Robbins and Siegmund in the context of hypothesis testing, and then Lorden and Pollak [8] in the sequential change detection problem. This approach replaces the unknown parameter with some non-anticipating estimator, which can be easier to find even if there is no conjugate prior, as in the Gamma example considered in [8,25]. This work developed a modified SR procedure by introducing a prior distribution to the unknown parameters. While the non-anticipating estimator approach [8,24] enjoys recursive and thus efficient computation for the likelihood ratio based detection statistics, their approaches to constructing recursive estimators (based on MLE or method-of-moments) cannot be easily extended to more complex cases (for instance, multi-dimensional parameters with constraints). Here, we consider a general and convenient approach for constructing non-anticipating estimators based on online convex optimization which is particularly useful for these complex cases. Our work provides an alternative proof for the nearly second-order asymptotic optimality by building a connection to online convex optimization and leveraging the regret bound type of results [14]. For one-dimensional Gaussian mean shift without any constraint, we replicate the second-order asymptotic optimality, namely, Theorem 3.3 in [24]. Recent work [26] also treats the problem when the pre-change distribution has unknown parameters.

Another related problem is sequential joint estimation and detection, but the goal is different in that one aims to achieve both good detection and good estimation performance, whereas in our setting, estimation is only needed for computing the detection statistics. These works include [27] and [28], which study the joint detection and estimation problem of a specific form that arises from many applications such as spectrum sensing [29], image observations [30], and MIMO radar [31]: a linear scalar observation model with Gaussian noise, and under the alternative hypothesis there is an unknown multiplicative parameter. The paper of [27] demonstrates that solving the joint problem by treating detection and estimation separately with the corresponding optimal procedure does not yield an overall optimum performance, and provides an elegant closed-form optimal detector. Later on [28] generalizes the results. There are also other approaches solving the joint detection-estimation problem using multiple hypotheses testing [30,32] and Bayesian formulations [33].

Related work using online convex optimization for anomaly detection includes [9], which develops an efficient detector for the exponential family using online mirror descent and proves a logarithmic regret bound, and [10], which dynamically adjusts the detection threshold to allow feedbacks about whether decision outcome. However, these works consider a different setting that the change is a transient outlier instead of a persistent change, as assumed by the classic statistical change-point detection literature. When there is persistent change, it is important to accumulate “evidence” by pooling the post-change samples (our work considers the persistent change).

Extensive work has been done for parameter estimation in the online-setting. This includes online density estimation over the exponential family by regret minimization [9,10,16], sequential prediction of individual sequence with the logarithm loss [13,34], online prediction for time series [35], and sequential NML (SNML) prediction [34] which achieves the optimal regret bound. Our problem is different from the above, in that estimation is not the end goal; one only performs parameter estimation to plug them back into the likelihood function for detection. Moreover, a subtle but important difference of our work is that the loss function for online detecting estimation is $-f_{{\widehat{\theta }}_i}\left(X_i\right)$, whereas our loss function is $-f_{{\widehat{\theta }}_{i-1}}\left(X_i\right)$ in order to retain the martingale property, which is essential to establish the nearly second-order asymptotic optimality.

2. Preliminaries

Assume a sequence of i.i.d. random variables $X_1,X_2,\dots$ with a probability density function of a parametric form $f_{\theta }$. The parameter $\theta$ may be unknown. Consider two related problems: one-sided sequential hypothesis test and sequential change-point detection. The detection statistic relies on a sequence estimator $\{{\widehat{\theta }}_t\}$ constructed using online mirror descent. The OMD uses simple one-sample update: the update from ${\widehat{\theta }}_{t-1}$ to ${\widehat{\theta }}_t$ only uses the current sample $X_t$. This is the main difference from the traditional generalized likelihood ratio (GLR) statistic [7], where each ${\widehat{\theta }}_t$ is estimated using historical samples. In the following, we present detailed descriptions for two problems. We will consider exponential family distributions and present our non-anticipating estimator based on the one-sample estimate.

2.1. One-Sided Sequential Hypothesis Test

First, we consider a one-sided sequential hypothesis test where the goal is only to reject the null hypothesis. This is a special case of the change-detection problem where the change-point can be either 0 or ∞ (meaning it never occurs). Studying this special case will given us an important intermediate step towards solving the sequential change-detection problem.

Consider the null hypothesis ${\mathsf{H}}_0:\theta ={\theta }_0$ versus the alternative ${\mathsf{H}}_1:\theta \ne {\theta }_0$. Hence, the parameter under the alternative distribution is unknown. The classic approach to solve this problem is the one-sided sequential probablity-ratio test (SPRT) [36]: at each time, given samples $\{X_1,X_2,\dots ,X_t\}$, the decision is either to reject ${\mathsf{H}}_0$ or taking more samples if the rejection decision can not be made confidently. Here, we introduce a modified one-sided SPRT with a sequence of non-anticipating plug-in estimators:

$${\widehat{\theta }}_t:={\widehat{\theta }}_t(X_1,\dots ,X_t),t=1,2,\dots .$$

(1)

Define the test statistic at time t as

$${\mathsf{\Lambda }}_t=\prod _{i=1}^t\frac{f_{{\widehat{\theta }}_{i-1}}\left(X_i\right)}{f_{{\theta }_0}\left(X_i\right)},i\ge 1.$$

(2)

The test statistic has a simple recursive implementation:

$${\mathsf{\Lambda }}_t={\mathsf{\Lambda }}_{t-1}·\frac{f_{{\widehat{\theta }}_{t-1}}\left(X_t\right)}{f_{{\theta }_0}(X_t}.$$

Define a sequence of $\sigma$-algebras ${\left\{{\mathcal{F}}_t\right\}}_{t\ge 1}$ where ${\mathcal{F}}_t=\sigma (X_1,\dots ,X_t)$. The test statistic has the martingale property due to its non-anticipating nature: $\mathbb{E}[{\mathsf{\Lambda }}_t\mid {\mathcal{F}}_{t-1}]={\mathsf{\Lambda }}_{t-1}$, where the expectation is taken when $X_1,\dots$ are i.i.d. random variables drawn from $f_{{\theta }_0}$. The decision rule is a stopping time

$$\tau \left(b\right)=min\{t\ge 1:log{\mathsf{\Lambda }}_t\ge b\},$$

(3)

where $b>0$ is a pre-specified threshold. We reject the null hypothesis whenever the statistic exceeds the threshold. The goal is to reject the null hypothesis using as few samples as possible under the false-alarm rate (or Type-I error) constraint.

2.2. Sequential Change-Point Detection

Now we consider the sequential change-point detection problem. A change may occur at an unknown time $\nu$ which alters the underlying distribution of the data. One would like to detect such a change as quickly as possible. Formally, change-point detection can be cast into the following hypothesis test:

$$\begin{array}{cc}\hfill {\mathsf{H}}_0:& X_1,X_2,\dots \stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}.}{\sim }{f}_{{\theta }_0},\hfill \\ \hfill {\mathsf{H}}_1:& X_1,\dots ,X_{\nu }\stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}.}{\sim }{f}_{{\theta }_0},X_{\nu +1},X_{\nu +2},\dots \stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}.}{\sim }{f}_{\theta },\hfill \end{array}$$

(4)

Here we assume an unknown $\theta$ to represent the anomaly. The goal is to detect the change as quickly as possible after it occurs under the false-alarm rate constraint. We will consider likelihood ratio based detection procedures adapted from two types of existing ones, which we call the adaptive CUSUM (ACM), and the adaptive SRRS (ASR) procedures.

For change-point detection, the post-change parameter is estimated using post-change samples. This means that, for each putative change-point location before the current time $k<t$, the post-change samples are $\{X_k,\dots ,X_t\}$; with a slight abuse of notation, the post-change parameter is estimated as

$${\widehat{\theta }}_{k,i}={\widehat{\theta }}_{k,i}(X_k,\dots ,X_i),i\ge k.$$

(5)

Therefore, for $k=1$, ${\widehat{\theta }}_{k,i}$ becomes ${\widehat{\theta }}_i$ defined in (2) for the one-sided SPRT. Initialize with ${\widehat{\theta }}_{k,k-1}={\theta }_0$. The likelihood ratio at time t for a hypothetical change-point location k is given by

$${\mathsf{\Lambda }}_{k,t}=\prod _{i=k}^t\frac{f_{{\widehat{\theta }}_{k,i-1}}\left(X_i\right)}{f_{{\theta }_0}\left(X_i\right)},$$

(6)

where ${\mathsf{\Lambda }}_{k,t}$ can be computed recursively similar to (2).

Since we do not know the change-point location $\nu$, from the maximum likelihood principle, we take the maximum of the statistics over all possible values of k. This gives the ACM procedure:

$$T_{\mathrm{ACM}}\left(b_1\right)=inf\left\{t\ge 1:\underset{1\le k\le t}{max}log{\mathsf{\Lambda }}_{k,t}>b_1\right\},$$

(7)

where $b_1$ is a pre-specified threshold. Similarly, by replacing the maximization over k in (7) with summation, we obtain the following ASR procedure [8], which can be interpreted as a Bayesian statistic similar to the Shiryaev-Roberts procedure.

$$T_{\mathrm{ASR}}\left(b_2\right)=inf\left\{t\ge 1:log\left(\sum _{k=1}^t{\mathsf{\Lambda }}_{k,t}\right)>b_2\right\},$$

(8)

where $b_2$ is a pre-specified threshold. The computations of ${\mathsf{\Lambda }}_{k,t}$ and estimator $\{{\widehat{\theta }}_t\}$, $\{{\widehat{\theta }}_{k,t}\}$ are discussed later in Section 2.4. For a fixed k, the comparison between our methods and GLR is illustrated in Figure 1.

Remark 1.

In practice, to prevent the memory and computation complexity from blowing up as time t goes to infinity, we can use window-limited version of the detection procedures in (7) and (8). The window-limited versions are obtained by replacing ${max}_{1\le k\le t}$ with ${max}_{t-w\le k\le t}$ in (7) and by replacing ${\sum }_{k=1}^t$ with ${\sum }_{k=t-w}^t$ in (8). Here w is a prescribed window size. Even if we do not provide theoretical analysis to the window-limited versions, we refer the readers to [7] for the choice of w the window-limited GLR procedures.

2.3. Exponential Family

In this paper, we focus on $f_{\theta }$ being the exponential family for the following reasons: (i) exponential family [10] represents a very rich class of parametric and even many nonparametric statistical models [37]; (ii) the negative log-likelihood function for exponential family $-log{f}_{\theta }\left(x\right)$ is convex, and this allows us to perform online convex optimization. Some useful properties of the exponential family are briefly summarized below, and full proofs can be found in [10,38].

Consider an observation space $\mathcal{X}$ equipped with a sigma algebra $\mathcal{B}$ and a sigma finite measure H on $(\mathcal{X},\mathcal{B})$. Assume the number of parameters is d. Let $x^{⊺}$ denote the transpose of a vector or matrix. Let $\varphi :\mathcal{X}\to {\mathbb{R}}^d$ be an H-measurable function $\varphi \left(x\right)={({\varphi }_1\left(x\right),\dots ,{\varphi }_d\left(x\right))}^{⊺}$. Here $\varphi \left(x\right)$ corresponds to the sufficient statistic for $\theta$. Let $\mathsf{\Theta }$ denote the parameter space in ${\mathbb{R}}^d$. Let $\{{\mathcal{P}}_{\theta },\theta \in \mathsf{\Theta }\}$ be a set of probability distributions with respect to the measure H. Then, $\{{\mathcal{P}}_{\theta },\theta \in \mathsf{\Theta }\}$ is said to be a multivariate exponential family with natural parameter $\theta$, if the probability density function of each $f_{\theta }\in {\mathcal{P}}_{\theta }$ with respect to H can be expressed as $f_{\theta }\left(x\right)=exp\{{\theta }^{⊺}\varphi \left(x\right)-\mathsf{\Phi }\left(\theta \right)\}.$ In the definition, the so-called log-partition function is given by

$$\mathsf{\Phi }\left(\theta \right):=log{\int }_{\mathcal{X}}exp\left({\theta }^{⊺}\varphi \left(x\right)\right)dH\left(x\right).$$

To make sure $f_{\theta }\left(x\right)$ a well-defined probability density, we consider the following two sets for parameters:

$$\mathsf{\Theta }=\{\theta \in {\mathbb{R}}^d:log{\int }_{\mathcal{X}}exp\left({\theta }^{⊺}\varphi \left(x\right)\right)dH\left(x\right)<+\infty \},$$

and

$${\mathsf{\Theta }}_{\sigma }=\{\theta \in \mathsf{\Theta }:{\nabla }^2\mathsf{\Phi }\left(\theta \right)⪰\sigma I_{d\times d}\}.$$

Note that $-log{f}_{\theta }\left(x\right)$ is $\sigma$-strongly convex over ${\mathsf{\Theta }}_{\sigma }$. Its gradient corresponds to $\nabla \mathsf{\Phi }\left(\theta \right)={\mathbb{E}}_{\theta }\left[\varphi \left(X\right)\right]$, and the Hessian ${\nabla }^2\mathsf{\Phi }\left(\theta \right)$ corresponds to the covariance matrix of the vector $\varphi \left(X\right)$. Therefore, ${\nabla }^2$$\mathsf{\Phi }\left(\theta \right)$ is positive semidefinite and $\mathsf{\Phi }\left(\theta \right)$ is convex. Moreover, $\mathsf{\Phi }$ is a Legendre function, which means that it is strongly convex, continuous differentiable and essentially smooth [38]. The Legendre-Fenchel dual ${\mathsf{\Phi }}^{\ast }$ is defined as

$${\mathsf{\Phi }}^{\ast }\left(z\right)=\underset{u\in \mathsf{\Theta }}{sup}\{u^{⊺}z-\mathsf{\Phi }\left(u\right)\}.$$

The mappings $\nabla {\mathsf{\Phi }}^{\ast }$ is an inverse mapping of $\nabla \mathsf{\Phi }$ [39]. Moreover, if $\mathsf{\Phi }$ is a strongly convex function, then $\nabla {\mathsf{\Phi }}^{\ast }={\left(\nabla \mathsf{\Phi }\right)}^{-1}$.

A general measure of proximity used in the OMD is the so-called Bregman divergence $B_F$, which is a nonnegative function induced by a Legendre function F (see, e.g., [10,38]) defined as

$$B_F(u,v):=F\left(u\right)-F\left(v\right)-\langle \nabla F\left(v\right),u-v\rangle .$$

(9)

For exponential family, a natural choice of the Bregman divergence is the Kullback-Leibler (KL) divergence. Define ${\mathbb{E}}_{\theta }$ as the expectation when X is a random variable with density $f_{\theta }$ and $I({\theta }_1,{\theta }_2)$ as the KL divergence between two distributions with densities $f_{{\theta }_1}$ and $f_{{\theta }_2}$ for any ${\theta }_1,{\theta }_2\in \mathsf{\Theta }$. Then

$$I({\theta }_1,{\theta }_2)={\mathbb{E}}_{{\theta }_1}\left[log(f_{{\theta }_1}\left(X\right)/f_{{\theta }_2}\left(X\right))\right].$$

(10)

It can be shown that, for exponential family, $I({\theta }_1,{\theta }_2)=\mathsf{\Phi }\left({\theta }_2\right)-\mathsf{\Phi }\left({\theta }_1\right)-{({\theta }_2-{\theta }_1)}^{⊺}\nabla \mathsf{\Phi }\left({\theta }_1\right).$ Using the definition (9), this means that $B_{\mathsf{\Phi }}$

$$B_{\mathsf{\Phi }}({\theta }_1,{\theta }_2):=I({\theta }_2,{\theta }_1)$$

(11)

is a Bregman divergence. This property is useful to constructing mirror descent estimator for the exponential family [39,40].

2.4. Online Convex Optimization (OCO) Algorithms for Non-Anticipating Estimators

Online convex optimization (OCO) algorithms [14] can be interpreted as a player who makes sequential decisions. At the time of each decision, the  outcomes are unknown to the player. After committing to a decision, the decision maker suffers a loss that can be adversarially chosen. An OCO algorithm makes decisions, which, based on the observed outcomes, minimizes the regret that is the difference between the total loss that has incurred relative to that of the best fixed decision in hindsight. To design non-anticipating estimators, we consider OCO algorithms with likelihood-based regret functions. We iteratively estimate the parameters at the time when one new observation becomes available based on the maximum likelihood principle, and hence the loss incurred corresponds to the negative log-likelihood of the new sample evaluated at the estimator ${\ell }_t\left(\theta \right):=-log{f}_{\theta }\left(X_t\right)$, which corresponds to the log-loss in [13]. Given samples $X_1,\dots ,X_t$, the regret for a sequence of estimators ${\{{\widehat{\theta }}_i\}}_{i=1}^t$ generated by a likelihood-based OCO algorithm $\mathsf{a}$ is defined as

$${\mathcal{R}}_t^{\mathsf{a}}=\sum _{i=1}^t\{-log{f}_{{\widehat{\theta }}_{i-1}}\left(X_i\right)\}-\underset{\tilde{\theta }\in \mathsf{\Theta }}{inf}\sum _{i=1}^t\{-log{f}_{\tilde{\theta }}\left(X_i\right)\}.$$

(12)

Below we omit the superscript $\mathsf{a}$ occasionally for notational simplicity.

In this paper, we consider a generic OCO procedure called the online mirror descent algorithms (OMD) [14,41]. Next, we discuss how to construct the non-anticipating estimators ${\{{\widehat{\theta }}_t\}}_{t\ge 1}$ in (1), and $\{{\widehat{\theta }}_{k,t}\},k=1,2,\dots ,t-1$ in (5) using OMD. The main idea of OMD is the following. At each time step, the estimator ${\widehat{\theta }}_{t-1}$ is updated using the new sample $X_t$, by balancing the tendency to stay close to the previous estimate against the tendency to move in the direction of the greatest local decrease of the loss function. For the loss function defined above, a sequence of OMD estimator is constructed by

$${\widehat{\theta }}_t=\underset{u\in \mathsf{\Gamma }}{argmin}[u^{⊺}\nabla {\ell }_t({\widehat{\theta }}_{t-1})+\frac{1}{{\eta }_i}{B}_{\mathsf{\Phi }}(u,{\widehat{\theta }}_{t-1})],$$

(13)

where $B_{\mathsf{\Phi }}$ is defined in (11). Here $\mathsf{\Gamma }\subset {\mathsf{\Theta }}_{\sigma }$ is a closed convex set, which is problem-specific and encourages certain parameter structure such as sparsity.

Remark 2.

Similar to (13), for any fixed k, we can compute ${\{{\widehat{\theta }}_{k,t}\}}_{t\ge 1}$ via OMD for sequential change-point detection. The only difference is that ${\{{\widehat{\theta }}_{k,t}\}}_{t\ge 1}$ is computed if we use $X_k$ as our first sample and then apply the recursive update (13) on $X_{k+1},\dots$. For ${\widehat{\theta }}_t$, we use $X_1$ as our first sample.

There is an equivalent form of OMD, presented as the original formulation in [40]. The equivalent form is sometimes easier to use for algorithm development, and it consists of four steps: (1) compute the dual variable: ${\widehat{\mu }}_{t-1}=\nabla \mathsf{\Phi }({\widehat{\theta }}_{t-1})$; (2) perform the dual update: ${\widehat{\mu }}_t={\widehat{\mu }}_{t-1}-{\eta }_t\nabla {\ell }_t({\widehat{\theta }}_{t-1})$; (3) compute the primal variable: ${\tilde{\theta }}_t={\left(\nabla \mathsf{\Phi }\right)}^{\ast }\left({\widehat{\mu }}_t\right)$; (4) perform the projected primal update: ${\widehat{\theta }}_t={argmin}_{u\in \mathsf{\Gamma }}{B}_{\mathsf{\Phi }}(u,{\tilde{\theta }}_t)$. The equivalence between the above form for OMD and the nonlinear projected subgradient approach in (13) is proved in [39]. We adopt this approach when deriving our algorithm and follow the same strategy as [9]. Algorithm 1 summarizes the steps [42].

For strongly convex loss function, the regret of many OCO algorithms, including the OMD, has the property that ${\mathcal{R}}_n\le Clogn$ for some constant C (depend on $f_{\theta }$ and ${\mathsf{\Theta }}_{\sigma }$) and any positive integer n [10,43]. Note that for exponential family, the loss function is the negative log-likelihood function, which is strongly convex over ${\mathsf{\Theta }}_{\sigma }$. Hence, we can have the logarithmic regret property.

Algorithm 1 Online mirror-descent for non-anticipating estimators.

Require: Exponential family specifications $\varphi \left(x\right),\mathsf{\Phi }\left(x\right)$ and $f_{\theta }\left(x\right)$; initial parameter value ${\theta }_0$; sequence of data $X_1,\dots ,X_t,\dots$; a closed, convex set for parameter $\mathsf{\Gamma }\subset {\mathsf{\Theta }}_{\sigma }$; a decreasing sequence ${\left\{{\eta }_t\right\}}_{t\ge 1}$ of strictly positive step-sizes.

1: ${\widehat{\theta }}_0={\theta }_0,{\mathsf{\Lambda }}_0=1$. {Initialization}

2: for all $t=1,2,\dots ,$ do

3: Acquire a new observation $X_t$

4: Compute loss ${\ell }_t({\widehat{\theta }}_{t-1})\triangleq -log{f}_{{\widehat{\theta }}_{t-1}}\left(X_t\right)=\mathsf{\Phi }({\widehat{\theta }}_{t-1})-{\widehat{\theta }}_{t-1}^{⊺}\varphi \left(X_t\right)$

5: Compute likelihood ratio ${\mathsf{\Lambda }}_t={\mathsf{\Lambda }}_{t-1}{\dot{f}}_{{\widehat{\theta }}_{t-1}}\left(X_t\right)/f_{{\theta }_0}\left(X_t\right)$

6: ${\widehat{\mu }}_{t-1}=\nabla \mathsf{\Phi }({\widehat{\theta }}_{t-1})$, ${\widehat{\mu }}_t={\widehat{\mu }}_{t-1}-{\eta }_t({\widehat{\mu }}_{t-1}-\varphi \left(X_t\right))$ {Dual update}

7: ${\tilde{\theta }}_t={\left(\nabla \mathsf{\Phi }\right)}^{\ast }\left({\widehat{\mu }}_t\right)$

8: ${\widehat{\theta }}_t={argmin}_{u\in \mathsf{\Gamma }}{B}_{\mathsf{\Phi }}(u,{\tilde{\theta }}_t)$ {Projected primal update}

9: end for

10: return ${\{{\widehat{\theta }}_t\}}_{t\ge 1}$ and ${\left\{{\mathsf{\Lambda }}_t\right\}}_{t\ge 1}$.

3. Nearly Second-Order Asymptotic Optimality of One-Sample Update Schemes

Below we prove the nearly second-order asymptotic optimality of the one-sample update schemes. More precisely, the nearly second-order asymptotic optimality means that the algorithm obtains the lower performance bound asymptotically up to a log-log factor in the false-alarm rate, as the false-alarm rate tends to zero (in many cases the log-log factor is a small number).

We first introduce some necessary notations. Denote ${\mathbb{P}}_{\theta ,\nu }$ and ${\mathbb{E}}_{\theta ,\nu }$ as the probability measure and the expectation when the change occurs at time $\nu$ and the post-change parameter is $\theta$, i.e., when $X_1,\dots ,X_{\nu }$ are i.i.d. random variables with density $f_{{\theta }_0}$ and $X_{\nu +1},X_{\nu +2},\dots$ are i.i.d. random variables with density $f_{\theta }$. Moreover, let ${\mathbb{P}}_{\infty }$ and ${\mathbb{E}}_{\infty }$ denote the probability measure when there is no change, i.e., $X_1,X_2,\dots$ are i.i.d. random variables with density $f_{{\theta }_0}$. Finally, let ${\mathcal{F}}_t$ denote the $\sigma$-algebra generated by $X_1,\dots ,X_t$ for $t\ge 1$.

3.1. “One-Sided” Sequential Hypothesis Test

Recall that the decision rule for sequential hypothesis test is a stopping time $\tau \left(b\right)$ defined in (3). The two standard performance metrics are the false-alarm rate, denoted as ${\mathbb{P}}_{\infty }(\tau \left(b\right)<\infty )$, and the expected detection delay (i.e., the expected number of samples needed to reject the null), denoted as ${\mathbb{E}}_{\theta ,0}\left[\tau \left(b\right)\right]$. A meaningful test should have both small ${\mathbb{P}}_{\infty }(\tau \left(b\right)<\infty )$ and small ${\mathbb{E}}_{\theta ,0}\left[\tau \left(b\right)\right]$. Usually, one adjusts the threshold b to control the false-alarm rate to be below a certain level.

Our main result is the following. As has been observed by [23], there is a loss in the statistical efficiency by using one-sample update estimators relative to the GLR approach using the entire samples $X_1,\dots ,X_t$ in the past. The theorem below shows that this loss corresponds to the expected regret given in (12).

Theorem 1 (Upper bound for OCO based SPRT).

Let ${\{{\widehat{\theta }}_t\}}_{t\ge 1}$ be a sequence of non-anticipating estimators generated by an OCO algorithm $\mathsf{a}$. As $b\to \infty$,

$${\mathbb{E}}_{\theta ,0}\left[\tau \left(b\right)\right]\le \frac{b}{I(\theta ,{\theta }_0)}+\frac{{\mathbb{E}}_{\theta ,0}\left[{\mathcal{R}}_{\tau \left(b\right)}^{\mathsf{a}}\right]}{I(\theta ,{\theta }_0)}+O\left(1\right)$$

(14)

Here $O\left(1\right)$ is a term upper-bounded by an absolute constant as $b\to \infty$.

The main idea of the proof is to decompose the statistic defining $\tau \left(b\right)$, $log\mathsf{\Lambda }\left(t\right)$, into a few terms that form martingales, and then invoke the Wald’s Theorem for the stopped process.

Remark 3.

The inequality (14) is valid for a sequence of non-anticipating estimators generated by an OCO algorithm. Moreover, (14) gives an explicit connection between the expected detection delay for the one-sided sequential hypothesis testing (left-hand side of (14)) and the regret for the OCO (the second term on the right-hand side of (14)). This illustrates clearly the impact of estimation on detection by an estimation algorithm dependent factor.

Note that in the statement of the Theorem 1, the stopping time $\tau \left(b\right)$ appears on the right-hand side of the inequality (14). For OMD, the expected sample size is usually small. By comparing with specific regret bound ${\mathcal{R}}_{\tau \left(b\right)}$, we can bound ${\mathbb{E}}_{\theta ,0}\left[\tau \left(b\right)\right]$ as discussed in Section 4. The most important case is that when the estimation algorithm has a logarithmic expected regret. For the exponential family, as shown in Section 3.3, Algorithm 1 can achieve ${\mathbb{E}}_{\theta ,0}\left[{\mathcal{R}}_n\right]\le Clogn$ for any positive integer n. To obtain a more specific order of the upper bound for ${\mathbb{E}}_{\theta ,0}\left[{\tau }_b\right]$ when b grows, we establish an upper bound for ${\mathbb{E}}_{\theta ,0}\left[{\tau }_b\right]$ as a function of b, to obtain the following Corollary 1.

Corollary 1.

Let ${\{{\widehat{\theta }}_t\}}_{t\ge 1}$ be a sequence of non-anticipating estimators generated by an OCO algorithm $\mathsf{a}$. Assume that ${\mathbb{E}}_{\theta ,0}\left[{\mathcal{R}}_n^{\mathsf{a}}\right]\le Clogn$ for any positive integer n and some constant $C>0$, we have

$${\mathbb{E}}_{\theta ,0}\left[\tau \left(b\right)\right]\le \frac{b}{I(\theta ,{\theta }_0)}+\frac{Clogb}{I(\theta ,{\theta }_0)}(1+o\left(1\right)).$$

(15)

Here $o\left(1\right)$ is a vanishing term as $b\to \infty$.

Corollary 1 shows that other than the well known first-order approximation $b/I(\theta ,{\theta }_0)$ [8,18], the expected detection delay ${\mathbb{E}}_{\theta ,0}\left[\tau \left(b\right)\right]$ is bounded by an additional term that is on the order of $log\left(b\right)$ if the estimation algorithm has a logarithmic regret. This $logb$ term plays an important role in establishing the optimality properties later. To show the optimality properties for the detection procedures, we first select a set of detection procedures with false-alarm rates lower than a prescribed value, and then prove that among all the procedures in the set, the expected detection delays of our proposed procedures are the smallest. Thus, we can choose a threshold b to uniformly control the false-alarm rate of $\tau \left(b\right)$.

Lemma 1 (false-alarm rate of $\tau \left(b\right)$).

Let ${\{{\widehat{\theta }}_t\}}_{t\ge 1}$ be any sequence of non-anticipating estimators. For any $b>0$, ${\mathbb{P}}_{\infty }(\tau \left(b\right)<\infty )\le exp(-b)$.

Lemma 1 shows that as b increases the false-alarm rate of $\tau \left(b\right)$ decays exponentially fast. We can set $b=log(1/\alpha )$ to make the false-alarm rate of $\tau \left(b\right)$ less than some $\alpha >0$. Next, leveraging an existing lower bound for general SPRT presented in Section 5.5.1.1 in [4], we establish the nearly second-order asymptotic optimality of OMD based SPRT as follows:

Corollary 2 (Nearly second-order optimality of OCO based SPRT).

Let ${\{{\widehat{\theta }}_t\}}_{t\ge 1}$ be a sequence of non-anticipating estimators generated by an OCO algorithm $\mathsf{a}$. Assume that ${\mathbb{E}}_{\theta ,0}\left[{\mathcal{R}}_n^{\mathsf{a}}\right]\le Clogn$ for any positive integer n and some constant $C>0$. Define a set $C\left(\alpha \right)=\{T:{\mathbb{P}}_{\infty }(T<\infty )\le \alpha \}$. For $b=log(1/\alpha )$, due to Lemma 1, $\tau \left(b\right)\in C\left(\alpha \right)$. For such a choice, $\tau \left(b\right)$ is nearly second-order asymptotic optimal in the sense that for any $\theta \in {\mathsf{\Theta }}_{\sigma }-\left\{{\theta }_0\right\}$, as $\alpha \to 0$,

$${\mathbb{E}}_{\theta ,0}\left[\tau \left(b\right)\right]-\underset{T\in C\left(\alpha \right)}{inf}{\mathbb{E}}_{\theta ,0}\left[T\right]=O(log(log(1/\alpha ))).$$

(16)

The result means that, compared with any procedure (including the optimal procedure) calibrated to have a false-alarm rate less than $\alpha$, our procedure incurs an at most $log(log(1/\alpha \left)\right)$ increase in the expected detection delay, which is usually a small number. For instance, even for a conservative case when we set $\alpha ={10}^{-5}$ to control the false-alarm rate, the number is $log(log(1/\alpha \left)\right)=2.44$.

3.2. Sequential Change-Point Detection

Now we proceed the proof by leveraging the close connection [18] between the sequential change-point detection and the one-sided hypothesis test. For sequential change-point detection, the two commonly used performance metrics [4] are the average run length (ARL), denoted by ${\mathbb{E}}_{\infty }\left[T\right]$; and the maximal conditional average delay to detection (CADD), denoted by ${sup}_{\nu \ge 0}{\mathbb{E}}_{\theta ,\nu }[T-\nu \mid T>\nu ]$. ARL is the expected number of samples between two successive false alarms, and CADD is the expected number of samples needed to detect the change after it occurs. A good procedure should have a large ARL and a small CADD. Similar to the one-sided hypothesis test, one usually choose the threshold large enough so that ARL is larger than a pre-specified level.

Similar to Theorem 1, we provide an upper bound for the CADD of our ASR and ACM procedures.

Theorem 2.

Consider the change-point detection procedure $T_{\mathrm{ACM}}\left(b_1\right)$ in (7) and $T_{\mathrm{ASR}}\left(b_2\right)$ in (8). For any fixed k, let ${\left\{{\widehat{\theta }}_{k,t}\right\}}_{t\ge 1}$ be a sequence of non-anticipating estimators generated by an OCO algorithm $\mathsf{a}$. Let $b_1=b_2=b$, as $b\to \infty$ we have that

$$\begin{array}{cc}& \underset{\nu \ge 0}{sup}{\mathbb{E}}_{\theta ,\nu }[T_{\mathrm{ASR}}\left(b\right)-\nu \mid T_{\mathrm{ASR}}\left(b\right)>\nu ]\le \underset{\nu \ge 0}{sup}{\mathbb{E}}_{\theta ,\nu }[T_{\mathrm{ACM}}\left(b\right)-\nu \mid T_{\mathrm{ACM}}\left(b\right)>\nu ]\hfill \\ \hfill \le & {\left(I(\theta ,{\theta }_0)\right)}^{-1}\left(b+{\mathbb{E}}_{\theta ,0}\left[{\mathcal{R}}_{\tau \left(b\right)}^{\mathsf{a}}\right]+O\left(1\right)\right).\hfill \end{array}$$

(17)

To prove Theorem 2, we relate the ASR and ACM procedures to the one-sided hypothesis test and use the fact that when the measure ${\mathbb{P}}_{\infty }$ is known, ${sup}_{\nu \ge 0}{\mathbb{E}}_{\theta ,\nu }[T-\nu \mid T>\nu ]$ is attained at $\nu =0$ for both the ASR and the ACM procedures. Above, we may apply a similar argument as in Corollary 1 to remove the dependence on $\tau \left(b\right)$ on the right-hand-side of the inequality. We establish the following lower bound for the ARL of the detection procedures, which is needed for proving Corollary 3:

Lemma 2 (ARL).

Consider the change-point detection procedure $T_{\mathrm{ACM}}\left(b_1\right)$ in (7) and $T_{\mathrm{ASR}}\left(b_2\right)$ in (8). For any fixed k, let ${\left\{{\widehat{\theta }}_{k,t}\right\}}_{t\ge 1}$ be any sequence of non-anticipating estimators. Let $b_1=b_2=b$, given a prescribed lower bound $\gamma >0$ for the ARL, we have

$${\mathbb{E}}_{\infty }\left[T_{\mathrm{ACM}}\left(b\right)\right]\ge {\mathbb{E}}_{\infty }\left[T_{\mathrm{ASR}}\left(b\right)\right]\ge \gamma ,$$

provided that $b\ge log\gamma$.

Lemma 2 shows that given a required lower bound $\gamma$ for ARL, we can choose $b=log\gamma$ to make the ARL be greater than $\gamma$. This is consistent with earlier works [8,25] which show that the smallest threshold b such that ${\mathbb{E}}_{\infty }\left[T_{ACM}\left(b\right)\right]\ge \gamma$ is approximate $log\gamma$. However, the bound in Lamma 2 is not tight, since in practice we can set $b=\rho log\gamma$ for some $\rho \in (0,1)$ to ensure that ARL is greater than $\gamma$.

Combing the upper bound in Theorem 2 with an existing lower bound for the CADD of SRRS procedure in [15], we obtain the following optimality properties.

Corollary 3 (Nearly second-order asymptotic optimality of ACM and ASR).

Consider the change-point detection procedure $T_{\mathrm{ACM}}\left(b_1\right)$ in (7) and $T_{\mathrm{ASR}}\left(b_2\right)$ in (8). For any fixed k, let ${\{{\widehat{\theta }}_{k,t}\}}_{t\ge 1}$ be a sequence of non-anticipating estimators generated by an OCO algorithm $\mathsf{a}$. Assume that ${\mathbb{E}}_{\theta ,0}\left[{\mathcal{R}}_n^{\mathsf{a}}\right]\le Clogn$ for any positive integer n and some constant $C>0$. Let $b_1=b_2=b$. Define $S\left(\gamma \right)=\{T:{\mathbb{E}}_{\infty }\left[T\right]\ge \gamma \}$. For $b=log\gamma$, due to Lemma 2, both $T_{\mathrm{ASR}}\left(b\right)$ and $T_{\mathrm{ACM}}\left(b\right)$ belong to $S\left(\gamma \right)$. For such b, both $T_{\mathrm{ASR}}\left(b\right)$ and $T_{\mathrm{ACM}}\left(b\right)$ are nearly second-order asymptotic optimal in the sense that for any $\theta \in \mathsf{\Theta }-\left\{{\theta }_0\right\}$

$$\begin{array}{c}\underset{\nu \ge 1}{sup}{\mathbb{E}}_{\theta ,\nu }[T_{\mathrm{ASR}}\left(b\right)-\nu +1\mid T_{\mathrm{ASR}}\left(b\right)\ge \nu ]\hfill \\ -\underset{T\left(b\right)\in S\left(\gamma \right)}{inf}\underset{\nu \ge 1}{sup}{\mathbb{E}}_{\theta ,\nu }[T\left(b\right)-\nu +1\mid T\left(b\right)\ge \nu ]=O(loglog\gamma ).\hfill \end{array}$$

(18)

A similar expression holds for $T_{\mathrm{ACM}}\left(b\right)$.

The result means that, compared with any procedure (including the optimal procedure) calibrated to have a fixed ARL larger than $\gamma$, our procedure incurs an at most $log(log\gamma )$ increase in the CADD. Comparing (18) with (16), we note that the ARL $\gamma$ plays the same role as $1/\alpha$ because $1/\gamma$ is roughly the false-alarm rate for sequential change-point detection [18].

3.3. Example: Regret Bound for Specific Cases

In this subsection, we show that the regret bound ${\mathcal{R}}_t$ can be expressed as a weighted sum of Bregman divergences between two consecutive estimators. This form of ${\mathcal{R}}_t$ is useful to show the logarithmic regret for OMD. The following result comes as a modification of [16].

Theorem 3.

Assume that $X_1,X_2,\dots$ are i.i.d. random variables with density function $f_{\theta }\left(x\right)$. Let ${\eta }_i=1/i$ in Algorithm 1. Assume that ${\{{\widehat{\theta }}_i\}}_{i\ge 1},{\left\{{\widehat{\mu }}_i\right\}}_{i\ge 1}$ are obtained using Algorithm 1 and ${\widehat{\theta }}_i={\tilde{\theta }}_i$ (defined in step 7 and 8 of Algorithm 1) for any $i\ge 1$. Then for any ${\theta }_0\in \mathsf{\Theta }$ and $t\ge 1$,

$${\mathcal{R}}_t=\sum _{i=1}^{t}i·B_{{\mathsf{\Phi }}^{\ast }}({\widehat{\mu }}_i,{\widehat{\mu }}_{i-1})=\frac{1}{2}\sum _{i=1}^{t}i·{({\widehat{\mu }}_i-{\widehat{\mu }}_{i-1})}^{⊺}\left[{\nabla }^2{\mathsf{\Phi }}^{\ast }\left({\tilde{\mu }}_i\right)\right]({\widehat{\mu }}_i-{\widehat{\mu }}_{i-1}),$$

where ${\tilde{\mu }}_i=\lambda {\widehat{\mu }}_i+(1-\lambda ){\widehat{\mu }}_{i-1}$, for some $\lambda \in (0,1)$.

Next, we use Theorem 3 on a concrete example. The multivariate normal distribution, denoted by $\mathcal{N}(\theta ,I_d)$, is parametrized by an unknown mean parameter $\theta$ and a known covariance matrix $I_d$ ($I_d$ is a $d\times d$ identity matrix). Following the notations in Section 2.3, we know that $\varphi \left(x\right)=x$, $dH\left(x\right)=(1/\sqrt{|2\pi I_d|})·exp\left(-x^{⊺}x/2\right)$, $\mathsf{\Theta }={\mathsf{\Theta }}_{\sigma }={\mathbb{R}}^d$ for any $\sigma <2$, $\mathsf{\Phi }\left(\theta \right)=(1/2){\theta }^{⊺}\theta$, $\mu =\theta$ and ${\mathsf{\Phi }}^{\ast }\left(\mu \right)=(1/2){\mu }^{⊺}\mu$, where $|·|$ denotes the determinant of a matrix, and H is a probability measure under which the sample follows $\mathcal{N}(0,I_d)$). When the covariance matrix is known to be some $\Sigma \ne I_d$, one can “whiten” the vectors by multiplying ${\Sigma }^{-1/2}$ to obtain the situation here.

Corollary 4 (Upper bound for the expected regret, Gaussian).

Assume $X_1,X_2,\dots$ are i.i.d. following $\mathcal{N}(\theta ,I_d)$ with some $\theta \in {\mathbb{R}}^d$. Assume that ${\{{\widehat{\theta }}_i\}}_{i\ge 1},{\left\{{\widehat{\mu }}_i\right\}}_{i\ge 1}$ are obtained using Algorithm 1 with ${\eta }_i=1/i$ and $\mathsf{\Gamma }={\mathbb{R}}^d$. For any $t>0$, we have that for some constant $C_1>0$ that depends on θ,

$${\mathbb{E}}_{\theta ,0}\left[{\mathcal{R}}_t\right]\le C_{1}dlogt/2.$$

The following calculations justify Corollary 4, which also serve as an example of how to use regret bound. First, the assumption ${\widehat{\theta }}_t={\tilde{\theta }}_t$ in Theorem 3 is satisfied for the following reasons. Consider $\mathsf{\Gamma }={\mathbb{R}}^d$ is the full space. According to Algorithm 1, using the non-negativity of the Bregman divergence, we have ${\widehat{\theta }}_t={argmin}_{u\in \mathsf{\Gamma }}{B}_{\mathsf{\Phi }}(u,{\tilde{\theta }}_t)={\tilde{\theta }}_t.$ Then the regret bound can be written as

$$\begin{array}{cc}\hfill {\mathcal{R}}_t=& \frac{1}{2}{({\widehat{\mu }}_1-{\widehat{\mu }}_0)}^{⊺}({\widehat{\mu }}_1-{\widehat{\mu }}_0)+\frac{1}{2}\sum _{i=2}^t[i·{({\widehat{\mu }}_i-{\widehat{\mu }}_{i-1})}^{⊺}({\widehat{\mu }}_i-{\widehat{\mu }}_{i-1})]\hfill \\ \hfill =& \frac{1}{2}{(X_1-{\theta }_0)}^{⊺}(X_1-{\theta }_0)+\frac{1}{2}\sum _{i=2}^t{({\widehat{\mu }}_i-{\widehat{\mu }}_{i-1})}^{⊺}(\varphi \left(X_i\right)-{\widehat{\mu }}_{i-1}).\hfill \end{array}$$

Since the step-size ${\eta }_i=1/i$, the second term in the above equation can be written as:

$$\begin{array}{cc}& \frac{1}{2}\sum _{i=2}^t{({\widehat{\mu }}_i-{\widehat{\mu }}_{i-1})}^{⊺}(\varphi \left(X_i\right)-{\widehat{\mu }}_{i-1})\hfill \\ \hfill =& \frac{1}{2}\sum _{i=2}^t{({\widehat{\mu }}_i-{\widehat{\mu }}_{i-1})}^{⊺}(\varphi \left(X_i\right)+{\widehat{\mu }}_i)-\sum _{i=2}^t\frac{1}{2}{({\widehat{\mu }}_i-{\widehat{\mu }}_{i-1})}^{⊺}({\widehat{\mu }}_{i-1}+{\widehat{\mu }}_i)\hfill \\ \hfill =& \sum _{i=2}^t\frac{1}{2(i-1)}{(\varphi \left(X_i\right)-{\widehat{\mu }}_i)}^{⊺}(\varphi \left(X_i\right)+{\widehat{\mu }}_i)+\sum _{i=2}^t\frac{1}{2}({∥{\widehat{\mu }}_{i-1}∥}^2-{∥{\widehat{\mu }}_i∥}^2)\hfill \\ \hfill =& \sum _{i=2}^t\frac{1}{2(i-1)}{∥X_i∥}^2-\sum _{i=2}^t\frac{1}{2(i-1)}{∥{\widehat{\mu }}_i∥}^2+\frac{1}{2}{∥{\widehat{\mu }}_1∥}^2-\frac{1}{2}{∥{\widehat{\mu }}_t∥}^2.\hfill \end{array}$$

Combining above, we have

$${\mathbb{E}}_{\theta ,0}\left[{\mathcal{R}}_t\right]\le \frac{1}{2}{\mathbb{E}}_{\theta ,0}\left[{(X_1-{\theta }_0)}^{⊺}(X_1-{\theta }_0)\right]+\frac{1}{2}\sum _{i=2}^t\frac{1}{i-1}{\mathbb{E}}_{\theta ,0}[{∥X_i∥}^2]+\frac{1}{2}{\mathbb{E}}_{\theta ,0}[{∥X_1∥}^2].$$

Finally, since ${\mathbb{E}}_{\theta ,0}[{∥X_i∥}^2]=d(1+{\theta }^2)$ for any $i\ge 1$, we obtain the desired result. Thus, with i.i.d. multivariate normal samples, the expected regret grows logarithmically with the number of samples.

Using the similar calculations, we can also bound the expected regret in the general case. As shown in the proof above for Corollary 4, the dominating term for ${\mathcal{R}}_t$ can be rewritten as

$$\sum _{i=2}^t\frac{1}{2(i-1)}{(\varphi \left(X_i\right)-{\widehat{\mu }}_i)}^{⊺}\left[{\nabla }^2{\mathsf{\Phi }}^{\ast }\left({\tilde{\mu }}_i\right)\right](\varphi \left(X_i\right)+{\widehat{\mu }}_i),$$

where ${\tilde{\mu }}_i$ is a convex combination of ${\widehat{\mu }}_{i-1}$ and ${\widehat{\mu }}_i$. For an arbitrary distribution, the term ${(\varphi \left(X_i\right)-{\widehat{\mu }}_i)}^{⊺}\left[{\nabla }^2{\mathsf{\Phi }}^{\ast }\left({\tilde{\mu }}_i\right)\right](\varphi \left(X_i\right)+{\widehat{\mu }}_i)$ can be viewed as a local normal distribution with the changing curvature ${\nabla }^2{\mathsf{\Phi }}^{\ast }\left({\tilde{\mu }}_i\right)$. Thus, it is possible to prove case-by-case the $O(logt)$-style bounds by making more assumptions about the distributions. Recall the notation ${\mathsf{\Theta }}_{\sigma }$ in Section 2.3 such that $-log{f}_{\theta }\left(x\right)$ is $\sigma$-strongly convex over ${\mathsf{\Theta }}_{\sigma }$. Let ${\parallel ·\parallel }_2$ denote the ${\ell }_2$ norm. Moreover, we assume that the true parameter belongs to a set $\mathsf{\Gamma }$ that is a closed and convex subset of ${\mathsf{\Theta }}_{\sigma }$ such that ${sup}_{\theta \in \mathsf{\Gamma }}{\parallel \nabla \mathsf{\Phi }\left(\theta \right)\parallel }_2\le M$ for some constant M. Thus, one can show that $-log{f}_{\theta }\left(x\right)$ is not only $\sigma$-strongly convex but also M-strongly smooth over $\mathsf{\Gamma }$. Theorem 3 in [10] shows that for all $\theta \in \mathsf{\Gamma }$ and $n\ge 1$, consider that ${\{{\widehat{\theta }}_i\}}_{i\ge 1}$ is obtained by OMD, then

$${\mathbb{E}}_{\theta ,0}\left[{\mathcal{R}}_n\right]\le \frac{{\mathbb{E}}_{\theta ,0}\left[{\left(\frac{1}{2}{max}_{1\le i\le n}{\parallel X_i\parallel }_2+\frac{1}{2}M\right)}^2\right]}{\sigma }·(logn+1).$$

Therefore, for any bounded distributions within the exponential family, we achieve a logarithmic regret. This logarithmic regret is valid for Bernoulli distribution, Beta distribution and some truncated versions of classic distributions (e.g., truncated Gaussian distribution, truncated Gamma distribution and truncated Geometric distribution analyzed in [44]).

4. Numerical Examples

In this section, we present some synthetic examples to demonstrate the good performance of our methods. We will focus on ACM and ASR for sequential change-point detection. In the following, we consider the window-limited versions (see Remark 1) of ACM and ASR with window size $w=100$. Recall that when the measure ${\mathbb{P}}_{\infty }$ is known, ${sup}_{\nu \ge 0}{\mathbb{E}}_{\theta ,\nu }[T-\nu \mid T>\nu ]$ is attained at $\nu =0$ for both ASR and ACM procedures (a proof can be found in the proof of Theorem 2). Therefore, in the following experiments we define the expected detection delay (EDD) as ${\mathbb{E}}_{\theta ,0}\left[T\right]$ for a stopping time T. To compare the performance between different detection procedures, we determine the threshold for each detection procedure by Monte-Carlo simulations such that the ARL for each procedure is about $10,000$. Below, we denote ${∥·∥}_2$, ${∥·∥}_1$ and ${∥·∥}_0$ as the ${\ell }_2$ norm, ${\ell }_1$ norm and ${\ell }_0$ norm defined as the number of non-zero entries, respectively. The following experiments are all run on the same Macbook Air with an Intel i7 Core CPU.

4.1. Detecting Sparse Mean-Shift of Multivariate Normal Distribution

We consider detect the sparse mean shift for multivariate normal distribution. Specifically, we assume that the pre-change distribution is $\mathcal{N}(0,I_d)$ and the post-change distribution is $\mathcal{N}(\theta ,I_d)$ for some unknown $\theta \in \{\theta \in {\mathbb{R}}^d:\parallel \theta {\parallel }_0\le s\}$, where s is called the sparsity of the mean shift. Sparse mean shift detection is of particular interest in sensor networks [45,46]. For this Gaussian case, the Bregman divergence is given by $B_{\mathsf{\Phi }}({\theta }_1,{\theta }_2)=I({\theta }_2,{\theta }_1)={\parallel {\theta }_1-{\theta }_2\parallel }_2^2/2.$ Therefore, the projection onto $\mathsf{\Gamma }$ in Algorithm 1 is a Euclidean projection onto a convex set, which in many cases can be implemented efficiently. As a frequently used convex relaxation of the ${\ell }_0$-norm ball, we set $\mathsf{\Gamma }=\{\theta :\parallel \theta {\parallel }_1\le s\}$ (it is known that imposing an ${\ell }_1$ constraint leads to sparse solution; see, e.g., [47]). Then, the projection onto ${\ell }_1$ ball can be computed very efficiently via a simple soft-thresholding technique [48].

Two benchmark procedures are the CUSUM and the GLR. For the CUSUM procedure, we specify a nominal post-change mean, which is an all-one vector. If knowing the post-change mean is sparse, we can also use the shrinkage estimator presented in [49], which performs hard or soft thresholding of the estimated post-change mean parameter. Our procedures are $T_{ASR}\left(b\right)$ and $T_{ACM}\left(b\right)$ with $\mathsf{\Gamma }={\mathbb{R}}^d$ and $\mathsf{\Gamma }=\{\theta :\parallel \theta {\parallel }_1\le 5\}$. In the following experiments, we run $10,000$ Monte Carlo trials to obtain each simulated EDD.

In the experiments, we set $d=20$. The post-change distributions are $\mathcal{N}(\theta ,I_d)$, where $100p\%$ entry of $\theta$ is 1 and others are 0, and the location of nonzero entries are random.

Table 1. Comparison of the EDDs in detecting the sparse mean shift of multivariate Gaussian distribution. Below, “CUSUM”: CUSUM procedure with pre-specified all-one vector as post-change parameter; “Shrinkage”: component-wise shrinkage estimator in [49]; “GLR”: GLR procedure; “ASR”: $T_{\mathrm{ASR}}\left(b\right)$ with $\mathsf{\Gamma }={\mathbb{R}}^d$; “ACM”: $T_{\mathrm{ACM}}\left(b\right)$ with $\mathsf{\Gamma }={\mathbb{R}}^d$; “ASR-L1”: $T_{\mathrm{ASR}}\left(b\right)$ with $\mathsf{\Gamma }=\{\theta :\parallel \theta {\parallel }_1\le 5\}$; “ACM-L1”: $T_{\mathrm{ACM}}\left(b\right)$ with $\mathsf{\Gamma }=\{\theta :\parallel \theta {\parallel }_1\le 5\}$. p is the proportion of non-zero entries in $\theta$. We run $10,000$ Monte Carlo trials to obtain each value. For each value, the standard deviation is less than one half of the value.

	$\mathit{p}=0.1$	$\mathit{p}=0.2$	$\mathit{p}=0.3$	$\mathit{p}=0.4$	$\mathit{p}=0.5$	$\mathit{p}=0.6$
CUSUM	188.60	146.45	64.30	18.97	7.18	3.77
Shrinkage	17.19	9.25	6.38	4.96	4.07	3.55
GLR	19.10	10.09	7.00	5.49	4.50	3.86
ASR	45.22	19.55	12.62	8.90	7.02	5.90
ACM	45.60	19.93	12.50	9.00	7.03	5.87
ASR-ℓ1	45.81	19.94	12.45	8.92	6.97	5.89
ACM-ℓ1	19.24	10.17	7.51	6.11	5.41	4.92

shows the EDDs versus the proportion p. Note that our procedures incur little performance loss compared with the GLR procedure and the CUSUM procedure. Notably, $T_{ACM}\left(b\right)$ with $\mathsf{\Gamma }=\{\theta :\parallel \theta {\parallel }_1\le 5\}$ performs almost the same as the GLR procedure and much better than the CUSUM procedure when p is small. This shows the advantage of projection when the true parameter is sparse.

4.2. Detecting the Scale Change in Gamma Distribution

We consider an example that detects the scale change in Gamma distributions. Assume that we observe a sequence $X_1,X_2\dots$ of samples drawn from $Gamma(\alpha ,\beta )$ for some $\alpha ,\beta >0$, with the probability density function given by $f_{\alpha ,\beta }\left(x\right)=exp(-x\beta )x^{\alpha -1}{\beta }^{\alpha }/\tilde{\mathsf{\Gamma }}\left(\alpha \right)$ (to avoid confusion with the $\mathsf{\Gamma }$ parameter in Algorithm 1 we use $\tilde{\mathsf{\Gamma }}(·)$ to denote the Gamma function). The parameter ${\alpha }^{-1}$ is called the dispersion parameter that scales the loss and the divergences. For simplicity, we fix $\alpha =1$ just like we fix the variance in the Gaussian case. The specifications in the Algorthm 1 are as follows: $\theta =-\beta$, $\mathsf{\Theta }=(-\infty ,0)$, $\varphi \left(x\right)=x$, $dH\left(x\right)=1$, $\mathsf{\Phi }\left(\theta \right)=-log(-\theta )$, $\mu =-1/\theta$ and ${\mathsf{\Phi }}^{\ast }\left(\mu \right)=-1-log\mu$. Assume that the pre-change distribution is $Gamma(1,1)$ and the post-change distribution is $Gamma(1,\beta )$ for some unknown $\beta >0$. We compare our algorithms with CUSUM, GLR and non-ancitipating estimator based on the method of moment (MOM) estimator in [8]. For the CUSUM procedure, we specify the post-change $\beta$ to be 2. The results are shown in

Table 2. Comparison of the EDDs in detecting the scale change in Gamma distribution. Below, “CUSUM”: CUSUM procedure with pre-specified post-change parameter $\beta =2$; “MOM”: Method of Moments estimator method; “GLR”: GLR procedure; “ASR”: $T_{\mathrm{ASR}}\left(b\right)$ with $\mathsf{\Gamma }=(-\infty ,0)$; “ACM”: $T_{\mathrm{ACM}}\left(b\right)$ with $\mathsf{\Gamma }=(-\infty ,0)$. We run $10,000$ Monte Carlo trials to obtain each value. For each value, the standard deviation is less than one half of the value.

	$\mathit{\beta }=0.1$	$\mathit{\beta }=0.5$	$\mathit{\beta }=2$	$\mathit{\beta }=5$	$\mathit{\beta }=10$
CUSUM	NaN	481.2	33.75	14.37	12.04
MOM	3.41	32.87	40.86	11.42	7.21
GLR	2.40	23.79	33.29	9.07	5.67
ASR	3.95	32.34	45.18	13.45	8.55
ACM	3.70	31.80	47.20	12.42	7.87

. CUSUM fails to detect the change when $\beta =0.1$, which is far away from the pre-specified post-change parameter $\beta =2$. We can see that performance loss of the proposed ACM method compared with GLR and MOM is very small.

4.3. Communication-Rate Change Detection with Erdos-Rényi Model

Next, we consider a problem to detect the communication-rate change in a network, which is a model for social network data. Suppose we observe communication between nodes in a network over time, represented as a sequence of (symmetric) adjacency matrices of the network. At time t, if node i and node j communicates, then the adjacency matrix has 1 on the $ij$th and $ji$th entries (thus it forms an undirected graph). The nodes that do not communicate have 0 on the corresponding entries. We model such communication patterns using the Erdos-Renyi random graph model. Each edge has a fixed probability of being present or absent, independently of the other edges. Under the null hypothesis, each edge is a Bernoulli random variable that takes values 1 with known probability p and value 0 with probability $1-p$. Under the alternative hypothesis, there exists an unknown time $\kappa$, after which a small subset of edges occur with an unknown and different probability $p^{\prime }\ne p$.

In the experiments, we set $N=20$ and $d=190$. For the pre-change parameters, we set $p_i=0.2$ for all $i=1,\dots ,d$. For the post-change parameters, we randomly select n out of the 190 edges, denoted by $\mathcal{E}$, and set $p_i=0.8$ for $i\in \mathcal{E}$ and $p_i=0.2$ for $i\notin \mathcal{E}$. As said before, let the change happen at time $\nu =0$ (since the upper bound for EDD is achieved at $\nu =0$ as argued in the proof of Theorem 2). To implement CUSUM, we specify the post-change parameters $p_i=0.8$ for all $i=1,\dots ,d$.

The results are shown in

Table 3. Comparison of the EDDs in detecting the changes of the communication-rates in a network. Below, “CUSUM”: CUSUM procedure with pre-specified post-change parameters $p=0.8$; “GLR”: GLR procedure; “ASR”: $T_{\mathrm{ASR}}\left(b\right)$ with $\mathsf{\Gamma }=\mathbb{R}$; “ACM”: $T_{\mathrm{ACM}}\left(b\right)$ with $\mathsf{\Gamma }=\mathbb{R}$. We run $10,000$ Monte Carlo trials to obtain each value. For each value, the standard deviation is less than one half of the value.

	$\mathit{n}=78$	$\mathit{n}=100$	$\mathit{n}=120$	$\mathit{n}=150$	$\mathit{n}=170$	$\mathit{n}=190$
CUSUM	473.11	2.06	2.00	2.00	2.00	2.00
GLR	2.00	1.96	1.27	1.00	1.00	1.00
ASR	8.64	6.39	5.08	3.92	3.36	2.94
ACM	8.67	6.37	5.07	3.88	3.32	2.94

. Our procedures are better than CUSUM procedure when n is small since the post-change parameters used in CUSUM procedure is far from the true parameter. Compared with the GLR procedure, our methods have a small performance loss, and the loss is almost negligible as n approaches to $d=190$.

Below are the specifications of Algorithm 1 in this case. For Bernoulli distribution with unknown parameter p, the natural parameter $\theta$ is equal to $log(p/(1-p\left)\right)$. Thus, we have $\mathsf{\Theta }=\mathbb{R}$, $\varphi \left(x\right)=x$, $dH\left(x\right)=1$, $\mathsf{\Phi }\left(\theta \right)=log(1+exp(\theta \left)\right)$, $\mu =exp\left(\theta \right)/(1+exp(\theta \left)\right)$ and ${\mathsf{\Phi }}^{\ast }\left(\mu \right)=\mu log\mu +(1-\mu )$$log(1-\mu )$.

4.4. Point Process Change-Point Detection: Poisson to Hawkes Processes

In this example, to illustrate the situation in Section 1.2, we consider a case where a homogeneous Poisson process switches to a Hawkes process (see, e.g., [12]); this can be viewed as a simplest case in Section 1.2 with one node. We construct ACM and ASR procedures. In this case, the MLE for the unknown post-change parameter cannot be found in close-form, yet ACM and ASR can be easily constructed and give reasonably good performance, although our theory no longer holds in this case due to the lack of i.i.d. samples.

The Hawkes process can be viewed as a non-homogeneous Poisson process where the intensity is influenced by historical events. The data consist of a sequence of events occurring at times $\{t_1,t_2,\dots ,t_n\}$ before a time horizon T: $t_i\le T$. Assume the intensity of the Poisson process is ${\lambda }_s,s\in (0,T)$ and there may exists a change-point $\kappa \in (0,T)$ such that the process changes. The null and alternative hypothesis tests are

$$\left\{\begin{array}{cc}{\mathsf{H}}_0:\hfill & {\lambda }_s=\mu ,0<s<T;\hfill \\ {\mathsf{H}}_1:\hfill & {\lambda }_s=\mu ,0<s<\kappa ,\hfill \\ & {\lambda }_s=\mu +\theta \sum _{\kappa <t_j<s}\phi (s-t_j),\kappa <s<T,\hfill \end{array}\right.$$

where $\mu$ is a known baseline intensity, $\theta >0$ is unknown magnitude of the change, $\phi \left(s\right)=\beta e^{-\beta s}$ is the normalized kernel function with pre-specified parameter $\beta >0$, which captures the influence from the past events. We treat the post-change influence parameter $\theta$ as unknown since it represents an anomaly.

We first use a sliding window to convert the event times into a sequence of vectors with overlapping events. Assume of size of the sliding window is L. For a given scanning time $T_i\le T$, we map all the events in $[T_i-L,T_i]$ to a vector $X_i={[t_{\left(1\right)},\dots ,t_{\left(m_i\right)}]}^{⊺}$, $t_{\left(i\right)}\in [T_i-L,T_i]$, where $m_i$ is the number of events falling into the window. Note that $X_i$ can have different length for different i. Consider a set of scanning times $T_1,T_2,\dots ,T_t$. This maps the event times into a sequence of vectors $X_1,X_2,\dots ,X_t$ of lengthes $m_1$, $m_2$, …, $m_t$. These scanning times can be arbitrary; here we set them to be event times so that there are at least one sample per sliding window.

For a hypothetical change-point location k, it can be shown that the log-likelihood ratio (between the Hawkes process and the Poisson process) as a function of $\theta$, is given by

$$\ell \left(\theta \right|X_i)=\sum _{t_q\in (T_i-L,T_i)}\log \left[\mu +\theta \sum _{t_j\in (T_i-L,t_q)}\beta e^{-\beta (t_q-t_j)}\right]-\mu L-\theta \sum _{t_q\in (T_i-L,T_i)}\left[1-e^{-\beta (T_i-t_q)}\right].$$

(19)

Now based on this sliding window approach, we can approximate the original change-point detection problem as the following. Without change, $X_1,\dots ,X_t$ are sampled from a Poisson process. Under the alternative, the change occurs at some time such that $X_1,\dots ,X_{\kappa }$ are sampled from a Poisson process, and $X_{\kappa +1},\dots ,X_t$ are sampled from a Hawkes process with parameter $\theta$, rather than a Poisson process. We define the estimator of $\theta$, for assumed change-point location $\kappa =k$ as follows

$${\widehat{\theta }}_{k,i}\triangleq {\widehat{\theta }}_{k,i}(X_k,\dots ,X_i)={\widehat{\theta }}_{k,i}(t_{\ell }\in [T_k,T_i])$$

(20)

Now, consider $k\in [i-w,i-1]$, and keep w estimators: ${\widehat{\theta }}_{i-w,i},\dots ,{\widehat{\theta }}_{i-1,i}$. The update for each estimator is based on stochastic gradient descent. By taking derivative with respect to $\theta$, we have

$$\frac{\partial \ell \left(\theta \right|X_i)}{\partial \alpha }=\sum _{t_q\in (T_i-L,T_i)}\frac{{\sum }_{t_j\in (T_i-L,t_q)}\beta e^{-\beta (t_q-t_j)}}{\mu +\theta {\sum }_{t_j\in (T_i-L,t_q)}\beta e^{-\beta (t_q-t_j)}}-\sum _{t_q\in (T_i-L,T_i)}\left[1-e^{-\beta (T_i-t_q)}\right],$$

Note that there is no close form expression for the MLE, which the solution to the above equation. We perform stochastic gradient descent instead

$${\widehat{\theta }}_{k,i+1}={\widehat{\theta }}_{k,i}-\gamma \frac{\partial \ell \left(\theta \right|X_{i+1})}{\partial \theta }{|}_{\theta ={\widehat{\theta }}_{k,i}},k=i-w+1,i-w,\dots ,i,$$

where $\gamma >0$ is the step-size. Now we can apply the ACM and ASR procedures, by using the fact that $f_{{\widehat{\theta }}_{k,t}}\left(X_{t+1}\right)/f_{{\theta }_0}\left(X_{t+1}\right)=\ell ({\widehat{\theta }}_{k,t}|X_{t+1})$ and calculating using (19).

Table 4. Point process change-point detection: EDD of ACM and ASR procedures for various values of true $\theta$; ARL of the procedure is controlled to be 5000 by selecting threshold via Monte Carlo simulation.

	$\mathit{\theta }=0.4$	$\mathit{\theta }=0.5$	$\mathit{\theta }=0.5$	$\mathit{\theta }=0.7$
ACM	33.03	27.75	20.39	16.16
ASR	38.59	24.96	20.17	13.91

shows the EDD for different $\alpha$. Here we choose the threshold such that ARL is 5000. We see that the scheme has a reasonably good performance, the detection delay decreases as the true signal strength $\theta$ increases.

5. Conclusions

In this paper, we consider sequential hypothesis testing and change-point detection with computationally efficient one-sample update schemes obtained from online mirror descent. We show that the loss of the statistical efficiency caused by the online mirror descent estimator (replacing the exact maximum likelihood estimator using the complete historical data) is related to the regret incurred by the online convex optimization procedure. The result can be generalized to any estimation method with logarithmic regret bound. This result sheds lights on the relationship between the statistical detection procedures and the online convex optimization.