Year: 2020

Abstract:
We consider the quickest change detection problem where both the parameters of preand postchange distributions are unknown, which prevent the use of classical simple hypothesis testing. Without additional assumptions, optimal solutions are not tractable as they rely on some minimax and robust variant of the objective. As a consequence, change points might be detected too late for practical applications (in economics, health care or maintenance for instance). Other approaches solve a relaxed version of the problem through the use of particular probability distributions or the use of domain knowledge. We tackle this problem in the more complex Markovian case and we provide a new scalable approximate algorithm with near optimal performance that runs in O(1).

heading: 1 INTRODUCTION
text: Quickest Change Detection (QCD) problems arise naturally in settings where a latent state controls observable signals (Basseville et al., 1993). In biology, it is applied in genomic sequencing (Caron et al., 2012) and in reliable healthcare monitoring (Salem et al., 2014). In industry, it finds application in faulty machinery detection (Lu et al., 2017; Martí et al., 2015) and in leak surveillance (Wang et al., 2014). It also has environmental applications such as traffic-related pollutants detection (Carslaw et al., 2006).
Any autonomous agent designed to interact with the world and achieve multiple goals must be able to detect relevant changes in the signals it is sensing in order to adapt its behaviour accordingly. This is particularly true for reinforcement learning based agents in multi-task settings as their policy is conditioned to some task parameter (Gupta et al., 2018; Teh et al., 2017). In order for the agent to be truly autonomous, it needs to identify the task at hand according to the environment requirement. For example, a robot built to assist cooks in the kitchen should be able to recognise the task being executed (chopping vegetables, cutting meat, ..) without external help to assist them efficiently. Otherwise, the agent requires a higher intelligence (one of the cooks for instance) to control it (by stating the task to be executed). In the general case, the current task is unknown and has to be identified sequentially from external sensory signals. The agent must track the changes as quickly as possible to adapt to its environment. However, current solutions for the QCD problem when task parameters are unknown, either do not scale or impose restrictive conditions on the setting. (i.i.d. observations, exponential family distributions, partial knowledge of the parameters, etc.).
In this paper, we construct a scalable algorithm with similar performances to optimal solutions. For this purpose, we use the change detection delay under known parameters as a lower bound for the delay in the unknown case. This improves our estimations of the parameters and thus improves our change point detection. We consider the case where the data is generated by some Markovian processes as in reinforcement learning. We assess our algorithm performances on synthetic data generated using distributions parameterised with neural networks in order to match the complexity level of real life applications. We also evaluate our algorithm on standard reinforcement learning environment.

heading: 2 QUICKEST CHANGE DETECTION PROBLEMS
text: Formally, consider a sequence of random observations (Xt) where each Xt belongs to some observation space X (say, an Euclidean space for simplicity) and is drawn from fθt(.|Xt−1), where
the parameter θt belongs to some task parameter space Θ and {fθ, θ ∈ Θ} is a parametric probability distribution family (non trivial, in the sense that all fθ are different). The main idea is that, at almost all stages, θt+1 = θt but there are some “change points” where those two parameters differ. Let us denote by tk the different change points and, with a slight abuse of notations, by θk the different values of the parameters. Formally, the data generating process is therefore: Xt ∼ ∑K k=0 fθk(.|Xt−1)1tk≤.<tk+1(t).
The overarching objective is to identify as quickly as possible the change points tk and the associated parameters θk, based on the observations (Xt). Typical procedures propose to tackle iteratively the simple change point detection problem, and for this reason we will focus mainly on the simpler setting of a single change point, where K = 2, t0 = 0, t1 = λ and t2 = ∞, where λ is unknown and must be estimated. For a formal description of the model and the different metrics, we will also assume that the parameters (θ0, θ1) are drawn from some distribution F over Θ. As a consequence, the data generating process we consider boils down to the following system 1:
θ0, θ1 ∼ F and { Xt+1 ∼ fθ0(.|Xt) if t ≤ λ Xt+1 ∼ fθ1(.|Xt) if t > λ . (1)

heading: 2.1 CRITERIA DEFINITIONS
text: As mentioned, the objective is to detect change points as quickly as possible while controlling the errors. There exist different metrics to evaluate algorithm performances; they basically all minimise some delay measurements while keeping the rate of type I errors (false positive change detection) under a certain level. We will describe later on different existing definitions of these criteria (type I error and delay). In order to evaluate a probability of error and an expected delay, we obviously need to define relevant probability measures first. Traditionally, there are two antagonistic ways to construct them: the MIN-MAX and the BAYESIAN settings (Veeravalli & Banerjee, 2014).
First, we denote by Pn (resp. En) the data probability distribution (resp. the expectation) conditioned to the change point happening at λ = n. This last event happens with some probability µ(λ = n) – with the notation that bold characters designate random variables–, and we denote by Pµ (resp. Eµ) the data probability distribution (resp. the expectation), integrated over µ, i.e., for any event Ω, it holds that Pµ(Ω) = ∑ n µ(λ = n)Pn(Ω). In the following, we describe the major existing formulations of the QCD problem where the goal is to identify an optimal stopping time τ :
BAYESIAN formulation: In this formulation, the error is the Probability of False Alarms (PFA) with respect to Pµ. The delay is evaluated as the Average Detection Delay (ADD) with respect to Eµ.
PFA(τ ) = Pµ(τ < λ) and ADD(τ ) = Eµ[(τ − λ)|τ > λ] In this setting, the goal is to minimise ADD while keeping PFA below a certain level α (as in Shiryaev formulation (Shiryaev, 1963)). Formally, this rewrites into:
(SHIRYAEV) { να = arg minτ∈∆α ADD(τ ) ∆α = {τ : PFA(τ ) < α} (2)
Min-Max formulation: The MIN-MAX formulation disregards prior distribution over the change point. As a consequence, the error is measured as the False Alarm Rate (FAR) with respect to the worst case scenario where no change occurs (P∞). As for the delay, two possibilities are studied: the Worst Average Detection Delay (WADD) and the Conditional Average Detection Delay (CADD). WADD evaluates this delay with respect to the worst scenario in terms of both the change point and the observations. CADD is a less pessimistic evaluation as it only considers the worst scenario in terms of change point. Mathematically they are defined as:
FAR(τ ) = 1 E∞[τ ]
and {
WADD(τ ) = supn esssupXn En[(τ − n)+|Xn] CADD(τ ) = supnEn[(τ − n)|τ > λ]
where Xn designates all observation up to the nth one. In this setting, the goal is to minimise either WADD (Lorden formulation (Lorden et al., 1971)); or CADD (Pollak formulation (Pollak, 1985)) while keeping FAR below a certain level α. Formally, these problems are written as follows:
(LORDEN) { να = arg minτ∈∆α WADD(τ ) ∆α = {τ : FAR(τ ) < α} and (POLLAK) { να = arg minτ∈∆α CADD(τ ) ∆α = {τ : FAR(τ ) < α}
(3)

heading: 3 TEMPORAL WEIGHT REDISTRIBUTION
text: Under known parameters, optimal solutions for the QCD consist in computing some statistics, denoted by Sθ0,θ1n along with a threshold Bα. The stopping time να is the first time S θ0,θ1 n exceeds Bα. We provide a more detailed description in Appendix A.1.

heading: 3.1 ASYMPTOTIC BEHAVIOUR OF THE SOLUTION TO THE BAYESIAN FORMULATION
text: The SHIRYAEV algorithm (Shiryaev, 1963) is asymptotically optimal in the i.i.d. case (Tartakovsky & Veeravalli, 2005). This result is extended to the non i.i.d. case under the following assumption:
Hypothesis 1 Given fθ0 and fθ1 , there exists q ∈ R and r ∈ N such that for any k ∈ N:
1 n ∑k+n t=k fθ1 (Xt+1|Xt) fθ0 (Xt+1|Xt) r-quickly−→ n→+∞ q (4)
Under Hypothesis 1, and with exponential or heavy tail prior, the SHIRYAEV algorithm is asymptotically optimal (Tartakovsky & Veeravalli, 2005). In addition, the moments of ADD satisfy the following property for all m ≤ r (Tartakovsky & Veeravalli, 2005):
Eµ[(να − λ)m|τ > λ] α→0∼ Eµ[(νSα − λ)m] α→0∼ ( | log(α)| q+d )m (5)
where νSα is the SHIRYAEV algorithm stopping rule and d = − limn→∞ logP(λ > n+ 1)/n. Our approach to the QCD problem under unknown parameters relies on the asymptotic behaviour from Equation (5) up to the second order (r = 2). We use it as a lower bound for the detection delay, providing a natural segmentation to our data when approximating the parameters. We will discuss this more in details in what follows.

heading: 3.2 RATIONAL FOR TEMPORAL WEIGHT REDISTRIBUTION
text: In the following, we denote the true parameters by θ∗0 and θ ∗ 1 and the corresponding stationary distributions by Π∗0 and Π ∗ 1 (that are well defined for irreducible Markov Chain). The purpose of this section is to devise an algorithm learning the parameters θ∗0 and θ ∗ 1 using the asymptotic behaviour of the detection delay under known parameters. The delay predicted in Equation (5) is a performance lower bound in our setting (otherwise it would contradict the optimality of the SHIRYAEV algorithm). Intuitively, the idea is to learn θ∗1 using the last | log(α)| q+d observation and to learn θ ∗ 0 using previous observation. This simple technique happens to be efficient and computationally simple.
Parameters optimisation: The optimal Sθ0,θ1t – called the SHIRYAEV statistic – has a specific form (see Equation (12) in Appendix) and the QCD problem under unknown parameters is equivalent to:
να = arg mint{t|S θ0,θ1 t > Bα} where
{ θ0 = arg minθ f0(θ, θ
1) θ1 = arg maxθ f1(θ 0, θ) , (6)
and (θ0, θ1) are random initialisation parameters, f0(θ0, θ1) = E0 [
log fθ1 (Xt+1|Xt) fθ0 (Xt+1|Xt)
] , f1(θ0, θ1) =
E1 [
log fθ1 (Xt+1|Xt) fθ0 (Xt+1|Xt)
] , and Ek∈{0,1} is the expectation with respect to the probability distribution
Pk(Xt, Xt+1) = Π ∗ k(Xt)fθ∗k(Xt+1|Xt). In fact, the solution να to Equation (6) is the optimal stopping time under the true parameters. This is a consequence of the following Lemma 1:
Lemma 1 θ∗0 and θ∗1 verifies the following for any θ0, θ1:
θ∗0 = arg minθ f0(θ, θ1) and θ ∗ 1 = arg maxθ f1(θ0, θ) (7)
In order to solve the equivalent QCD problem of Equation (6), a good approximation of the functions f0 and f1 given the observations (Xt)n0 is required. This implies the ability to distinguish pre-change samples from post-change samples and this is precisely why optimal solutions are intractable. They compute the statistics for any possible change point and consider the worst case scenario. Previous approximate solutions simplify the setting by using domain knowledge to infer the pre-change
distribution and by using a sliding window to evaluate the post-change distribution. The window size w is an irreducible cost in the expected delay.
With known parameters, optimal algorithms have an average detection delay proportional to the error threshold α and to the inverse of the KL divergence between the pre and post change distributions. Given the optimal stopping time να, it’s possible to evaluate the posterior distribution of the change point P(λ = t|να = n), which in turn is a good classifier of the pre and post change observation:
P(Xt ∼ fθ0 |να = n) ∝ P(λ > t|να = n) ; P(Xt ∼ fθ1 |να = n) ∝ P(λ < t|να = n) . (8) The objective of this section is to exploit this classifier to construct a family of functions that approximate well f0 and f1. Consider the following family of functions:
fn0 (θ0, θ1) = E n 0
[ log(
fθ1 fθ0
)(Xτ+1|Xτ ) ] and fn1 (θ0, θ1) = E n 1 [ log(
fθ1 fθ0
)(Xτ+1|Xτ ) ] ,
where both the observations Xt and the indicies τ are random variables. The observations are generated using Equation (1) with a finite horizon t2. The indicies are sampled with respect to τ ∼ P(λ > τ |να = n)/ ∑t2 i=0P(λ > i|να = n) in En0 while they are sampled with respect to
τ ∼ P(λ < τ |να = n)/ ∑t2 i=0P(λ < i|να = n) in En1 .
The family of functions fnk are a re-weighted expectation of the log-likelihood ratio of the observations using P(λ|να = n). Basically, under the assumption that the optimal detection happens at t = n, we associate to each observation Xτ a weight proportional to P(λ > τ |να = n) when estimating f0 (and respectively proportional to P(λ < τ |να = n) when estimating f1). In addition, this family of functions is practical as we can approximate asymptotically P(λ|να = n) using the theoretical delay behaviour of the SHIRYAEV algorithms presented in Equation (5).
Lemma 2 For any given parameters (θ0, θ1), the following convergences hold:
If λ =∞, then: limn,t2→∞ |fn0 (θ0, θ1)− f0(θ0, θ1)| = limn,t2→∞ |fn1 (θ0, θ1)− f0(θ0, θ1)| = 0 If λ <∞, then: {
limn,t2→∞ |fn0 (θ0, θ1)− f1(θ0, θ1)| = 0 limn,t2→∞ |fn1 (θ0, θ1)− f1(θ0, θ1)| = 0 limt2→∞ |f να 1 (θ0, θ1)− f1(θ0, θ1)| = 0
For any integer n ∈ N, if t2 = λ+ n, then: limλ,t2→∞ |f να 0 − f0| = 0.
A major implication of Lemma 2 is that with enough observations, fνα0 and f να 1 will eventually converge to the functions f0 and f1. In fact, f t0 approaches f0 up to να and then start degrading (as it ends up converging to f1) whereas f t1 becomes a better and better approximation of f1 around να and improves asymptotically. As such, we are going to consider the following problem as proxy to the original QCD problem (Equation (6)).
ν̂α = arg mint{t|S θt0,θ t 1 t > Bα} where { θt0 = arg minθ f t 0(θ, θ t−1 1 )
θt1 = arg maxθ f t 1(θ t−1 0 , θ)
. (9)
The rational behind this choice is that around the change point, Sθ t 0,θ t 1 t converges to S θ∗0 ,θ ∗ 1
t . Having access to a perfect evaluation of the functions f0 and f1 guarantees the fastest possible detection. In fact, να – the optimal stopping time under known parameters – is a lower bound to the detection delay in this setting. In practice, only the first t observation are accessible to evaluate f t0 and f t 1. This introduces approximation errors to the estimates θt0 and θ t 1, thus delaying the detection. However, this is not detrimental to the evaluation of the stopping time. Indeed:
Before να: θt0 is a good estimator of θ∗0 and θt1 is a bad estimator of θ∗1 . In fact, the fraction of prechange observations used to learn θ∗1 is of the order of min(1, λ/t). This helps maintaining a low log-likelihood ratio (equivalently maintain Sθ t 0,θ t 1
t below the cutting threshold), as the estimation of θ∗1 will converge to a parameter close to θ ∗ 0 .
After να: θt1 is a good approximation of θ∗1 , but θt0 is a noisy estimation of θ∗0 (as post-change observations are used to learn it). The fraction of the data generated using θ∗1 but used to learn θ∗0 is proportional to t−να t . This favours – up to a certain horizon – higher log-
likelihood ratio estimates (equivalently an incremental behaviour of the sequence Sθ t 0,θ t 1
t ). We will discuss further the performance issues of solving Equation (9) in the following.
Distribution Approximation: In order to exploit the previously introduced results, we need a good approximation of P(λ|να = n). If the observation satisfy Hypothesis 1, then the moments of P(λ|να = n) up to the rth order satisfy Equation (5). When r =∞, this is equivalent to the Hausdorff moment problem, which admits a unique solution. The generalised method of moments can be used to approximate this distribution. In the remaining of this paper, we will restrict ourselves to the case where r = 2. There is an infinite number of distributions that satisfy Equation (5) in this case; however, given a parametric family, we can compute analytically a solution of the problem so we consider the Logistic distribution. For a given θ0, θ1, we approximate P(λ|να = n) with fθ0,θ1n = Logistic(µ, s) where µ = n−
| log(α)| DKL(fθ0 |fθ1 )+d and
s = √
3 | log(α)|π(DKL(fθ0 |fθ1 )+d) . Verifying that fθ0,θ1n satisfies Equation (5) for m ≤ 2 is straightforward.
As a consequence, P(τ > t|τ ∼ fθ0,θ1n ) is a fair approximation of P(λ > t|να = n). Limitations: Solving the problem given by Equation (9) unfortunately yields poor performances for extreme values of the error threshold α. This is due to a degraded log likelihood estimation both before λ and as t→∞. In fact, the log-likelihood ratios Lt and L∗t defined as:
Lt(Xt+1|Xt) = log fθt1 fθt0 (Xt+1|Xt) and L∗t (Xt+1|Xt) = log fθ∗1 fθ∗0 (Xt+1|Xt),
satisfy the following lemma:
Lemma 3 The log-likelihood ratio Lt and L∗t admit the following asymptotic behaviours:
lim λ→∞
Eλ[L ∗ λ] < 0 but Lλ a.s−→ λ→+∞ 0 ; lim t→∞ Eλ[L ∗ t |λ <∞] > 0 but Lt a.s−→ t→+∞ 0 (10)
Practically, this means that before the change point, the log likelihood ratio is over-estimated (as it should be smaller than 0), while after the change point, the log likelihood ratio is eventually under-estimated (as it should exceed 0). This is a consequence of Lemma 2. In fact, for t λ (respectively for t ≤ λ), both functions f t0 and f t1 are approaching f1 (respectively f0). Thus in both cases θt0 and θ t 1 converge to the same value. However, combined with the behaviour of f t 0 and f t 1 around να from Lemma 2, we can expect the Kullback–Leibler (KL) divergence DKL(fθt0‖fθt1) to converge to 0 before and after the change point while peaking around the optimal stopping time. We exploit this observation to mitigate the problems highlighted with Lemma 3.
Annealing: The under-estimation of L∗t after the change point is due to the use of post change observations (drawn from fθ∗1 ) when estimating θ ∗ 0 . In practice, this is particularly problematic when the error rate must be controlled under a low threshold α (equivalently a high cutting threshold Bα) or when the pre and post change distributions are very similar. When t approaches the optimal stopping time (t ≈ να), the minimising argument of f t0 is converging to θ∗0 . As t grows, the approximation f t0 is degraded, while f t 1 is becoming a better approximation, thus θ t 1 starts converging to θ ∗ 1 . This leads to an increase in the KL divergence. Our proposition is to keep optimising the version of f t0 that coincides with this increase (which in turn is fνα0 : the best candidate to identify θ ∗ 0). As a fix, we propose: 1- to shift the probability P(λ > τ |νSα = n) by replacing it with P(λ > τ −∆|νSα = n) (where ∆ is the ‘delay’ of the detection with respect to the Shiryaev stopping time: (n−νSα )+), and 2- anneal the optimisation steps of θt0 as the KL divergence increases. These tweaks correct the objective function used to learn θ∗0 and stop its optimisation when the observations start to deviate from the learned pre-change distribution. The delay can be formally introduced by replacing f t0 with f t−∆ 0 in Equation (9). In practice, ∆ is a delay heuristic that increases as the KL divergence increases. This reduces the noise when learning θ∗0 . The second idea is to use a probability p0 = 1−∆ of executing a gradient step when learning the pre-change parameter. This anneals the optimisation of θ∗0 .
Penalisation: The over-estimation of L∗t before the change point, is due to the exclusive use of pre change observation when estimating θ∗1 . This is particularly problematic for application where a high error threshold is tolerable (equivalently a low cutting threshold for the used statistic). This is also an issue when the pre and post change distribution are very different. As a solution, we penalise the log-likelihood using the inverse KL divergence. Before the change point, both parameters converge to θ∗0 : this means that the KL divergence between the estimated distributions is around 0. Inversely, after the change point, if the optimisation of θt0 is annealed, the KL divergence between the estimated parameters must hover around DKL(fθ∗0‖fθ∗1 ). Formally, penalising the log likelihood can be seen as using L̂t when computing the statistics Sn where L̂t = Lt − c/DKL(fθt0‖fθt1) for some real c > 0.
We provide in the appendix an ablation analysis as well as experimental proof of the effectiveness of both annealing and penalisation.
Algorithm: In practice, only the first t observations are available in order to evaluate the parameters (θ∗0 , θ ∗ 1). Hence, we design a stochastic gradient descent (SGD) based algorithm to solve Equation (9). Consider the following loss functions, where I is a set of time indices:
Ln0 (I, θ0, θ1) = ∑ t∈I P(τ < t|τ ∼ fθ0,θ1n ) log( fθ1 fθ0
)(Xt+1|Xt) Ln1 (I, θ0, θ1) = ∑ t∈I P(τ > t|τ ∼ fθ0,θ1n ) log( fθ1 fθ0 )(Xt+1|Xt) (11)
The quantity Ltk is an estimator of f tk using the observations (Xt)t∈I . The rational previously discussed implies that by re-weighting random samples from the observations we can approximate the functions fk using the familly of functions f tk.
A good approximation of θt0 and θ t 1 can be obtained with few SGD steps. These parameters are used to update the SHIRYAEV or the CUSUM statistic. The change point is declared once the threshold Bα is exceeded. The proposed procedure is described in Algorithm 1. In addition to the classical SGD hyper-parameters (number of epochs, gradient step size, momentum, . . . ), we propose to control the penalisation and the annealing procedure using the coefficients c, ≥ 0 respectively.
Algorithm 1 Temporal Weight Redistribution
1: Procedure: TWR(θ0, θ1, (Xt)t2t=0, Ne, Bα, , c) 2: Initialise S ← 0; θ0, θ1 ← θ0, θ1; ∆← 0; and p0 ← 1 3: for t ∈ [1, t2] do 4: for e ∈ [1, Ne] do 5: Randomly sample I: a set of indices from [0, t] 6: θ0 ← gradient update rule(θ0,∇θ0Lt−∆0 (I, θ0, θ1)) with probability p0 7: θ1 ← gradient update rule(θ1,−∇θ1Lt1(I, θ0, θ1)) 8: end for 9: S ← statistic update rule(S, fθ1fθ0 (Xt+1|Xt)− c DKL(fθ0 |fθ1 )
10: If S > Bα then declare change point and break 11: If DKL(fθ0 |fθ1) > D̄ then ∆← ∆ + 1 and p0 ← p0 − 12: D̄ ← t−1t D̄ + 1 tDKL(fθ0 |fθ1) 13: end for

heading: 4 RELATED WORKS
text: Existing solutions for the QCD problem under known parameters can be extended to the case of unknown parameters using generalised likelihood ratio (GLR) statistics. For example, an optimal solution for the Lorden formulation is a generalisation of the CUSUM algorithm using the GLR statistic for testing the null hypothesis of no change-point, based on (Xi)n0 , versus the alternative hypothesis of a single change-point prior to n. However, given n observation, this solution runs in O(n2). Others approached the problem from a Bayesian perspective by keeping track of the run length posterior distribution (the time since the last change point). Since in practice this distribution is highly peaked, the Bayesian Online Change Point Detection (BOCPD) algorithm (Adams & MacKay, 2007) is namely a good alternative in the i.i.d. setting. In fact by pruning out run lengths with probabilities below or by considering only the K most probable run lengths, the algorithm runs respectively in O(n/ ) and O(Kn) without critical loss of performance. Extensions of BOCPD algorithm to the case of non i.i.d observations, namely the Gaussian Process Change Point Models run inO(n4) (Saatchi et al., 2010). If pruning is applied, the complexity can be reduced toO(nR2/ ) or O(nKR2) where R is the typical unpruned max run length (Saatchi et al., 2010). To avoid this complexity issue, multiple approximate solutions that run in O(1) have been proposed over the years. However, a common practice (Nikiforov, 2000; Li et al., 2009; Unnikrishnan et al., 2011; Singamasetty et al., 2017; Molloy & Ford, 2018) is to assume that pre-change parameters are known as this is not a major limitation in many applications (fault detection, quality control...) where the pre-change corresponds to nominal behaviour of a process and can be characterised offline. A sub-optimal solution (Nikiforov, 2000) is to track in parallel Sθ0,θ1n over a subset of possible post-change distributions and consider the worst case scenario. Somewhat Similarly to our approach,
adaptive algorithms (Li et al., 2009; Singamasetty et al., 2017) attempt to improve this procedure by tracking Sθ0,θ̄1n where θ̄1 = arg maxθ f1(θ0, θ). Different from proceeding works, our contribution in this paper is an approximate procedure that can handle the case of unknown pre- and post- change parameters, runs in O(1), and have near optimal performances.

heading: 5 EXPERIMENTAL RESULTS
text: Detection delay on synthetic data: The distribution family fθ used experimentally to simulate data given the pre/post change parameters and the previous observations relies on two deep neural network (φµ, φσ) used to evaluate the parameters of a Gaussian distribution such as: fθ(.|x) = N (φµ(θ, x), φσ(θ, x)). We compare our algorithm to three other possibilities: the optimal loglikelihood ratio (an oracle having access to the true parameters), the generalised log-likelihood ratio (GLR) and the adaptive log-likelihood ratio where the pre change parameters are learned offline (a fraction of pre-change observations are used to learn θ0).
We also introduce an alternative metric to the average detection delay: the regretted detection delay. An oracle with access to the true parameters will predict change as fast as the optimal algorithm under known parameters. In order to put performances into perspective, for a given cutting threshold B, the regretted average detection delayRB(a) of an algorithm, is the additional delay with respect to the oracle. If we denote the oracle’s optimal stopping time νB and the algorithm stopping time νaB , then: RB(a) = Eλ[(νaB − νB)+|νaB ≥ νB ≥ λ]. This formulation of the regret is conditioned to νaB ≥ νB for the same reason that the measurement of the average detection delay is conditioned to νB ≥ λ. An algorithm that predicts change faster than the optimal one under known parameters, is over-fitting the noise in the data to some extent. For this reason, when comparing performances, there is no true added value in detecting changes faster than the optimal algorithm in the known parameter setting.
We evaluate the performances using the SHIRYAEV update rule for the statistic (solution to the BAYESIAN setting). We provide in the appendix the analysis using the CUSUM update rule (solution to the MIN-MAX setting) with the same conclusions. We evaluate the regretted detection delay in the case where X = R10, Θ = R10, and (φµ, φσ) are 5-layer deep, 32-neurons wide neural network. We use a synthetic example where the KL divergence between preand post- change distribution is DKL(fθ∗0 |fθ∗1 ) = .3. This allows us to view performances in difficult settings. For higher KL divergence values, change detection becomes easier. We provide in the appendix the same analysis for the case whereDKL(fθ∗0 |fθ∗1 ) = 3. In Figure 1, we provide the average regret overs 500 simulations. Each simulated trajectory is 1000 long with a change point at the 500th observation. We use the full pipeline of our algorithm with
a penalisation coefficient (c = 0.01), and we allow the adaptive algorithm to exploit 10% of the pre change observations. Our algorithm (the blue curve) achieves comparable performance to the GLR statistic (without requiring intensive computational resources) and achieves better and more stable results compared to the adaptive algorithm (without requiring domain knowledge).
Multi-task reinforcement learning: Consider a task space Θ and for each task θ ∈ Θ the associated Markov decision processMθ = {S,A,P,Rθ, γ} where S is the state space, A is the action space, P is the environment dynamics, Rθ is the task specific reward function and γ is the discount factor. We evaluate algorithms using the copycat agent scenario: a main agent is solving a set of unknown tasks (θi)Ki=0 and switches from task to task at unknown change points (ti) K i=0. We
denote θ∗t = ∑ i θi1ti≤t<ti+1(t) the task being solved over time. Our goal is to construct a copycat agent predicting the main agent’s tasks online through the observation of the generated state action (st, at) ∈ S ×A. This scenario has many different real life applications: Users of any website can be viewed as optimal agents solving iteratively a set of tasks (listening to sets of music genre, looking for different types of clothes, ..). The state space corresponds to the different web-pages and the action space is the set of buttons available to the user. We argue that predicting the task being solved online is an important feature for behaviour forecasting and content recommendation.
In the following let xt = (st, at) denote the observations and let π∗θ denote the optimal policy of Mθ. The copycat problem is a QCD problem where the observations are drawn from fθ(xt+1|xt) = P(st+1|st, at)π∗θ(at+1|st+1). However this implies that we have the optimal policy π∗θ . In practice it is sufficient to either have access to the reward functions Rθ or a set of task labelled observations. In both cases the problem of learning π∗θ is well studied. WhenRθ is known, then coupling Hindsight Experience Replay (Andrychowicz et al., 2017) with state of the art RL algorithms such as Soft Actor Critic (Haarnoja et al., 2018) can solve a wide range of challenging multi-task problems (Plappert et al., 2018). When a history of observations are available, we can use generative adversarial inverse reinforcement learning (Ho & Ermon, 2016) or Maximum Entropy Inverse Reinforcement Learning (Gleave & Habryka, 2018; Yu et al., 2019) to learn π∗θ .
We evaluate our algorithm performances on the FetchReach environment from the gym library. Each time the main agent achieves the desired goal, a new task is sampled. We run this experiment for 500 time step. Temporal Weight Redistribution (TWR) is used in order to evaluate the tasks online. Algorithm 1 is easily adapted to the multiple change point setting: each time a change point is detected, we reinitialise the parameters ∆ and p0 and we set the task parameters to the previous post change parameters. As such, we use the running estimate θt0 overtime to approximate the main agent task θ∗t . Experimental results are reported in Figure 2.
We used the TWR algorithm with an annealing parameter = 0.1 and a penalisation coefficient c = 0.01, and the CUSUM statistic with a cutting threshold Bα = 50. We designate the task change points of the main agent with dashed black lines and the detected change points with a dashed red lines. On Figure 2a, we observe that the log-likelihood ratio (LLR) log(fθt1/fθt0)(xt) starts off negative, stabilises around 0 and becomes strictly positive after the actual change point. Once the change detected, we reset the TWR and the LLR falls down to a negative value. This trend is consistently reproduced with each new task. The main agent have a small probability of making a random action, this causes the LLR to peak sometimes before the change point. However this doesn’t cause the statistic to exceed the threshold Bα. We provide on Figure 2b the estimation error of θ∗t . The blue curve is the error of the TWR estimate ‖θt0 − θ∗t ‖. As a baseline we estimate the task using the last 10 observation: θ̂t = arg maxθ ∑10 i=1 log(fθ(xt−i). The error of this estimator ‖θ̂t − θ∗t ‖ is plotted in green. The prediction of the TWR copycat agent are clearly more reliable. Additional experimental results for the Fetchpush and the PickAndPlace environments as well as an analysis of the impact of using learned policies when solving the copycat problem are provided in the appendix.

heading: 6 CONCLUSION
text: We studied the QCD problem under unknown parameters for Markovian process. Extending our techniques to handle higher order Markov chains, (i.e., for some window parameter w > 0, fθ depends of Xt−w:t rather than just the last observation) is trivial, thus covering a wide range of natural processes. We provide a tractable algorithm with near optimal performances that relies on the asymptotic behaviour of optimal algorithms in the known parameter case. This mitigates the performance-cost dilemma providing the best of both worlds: a scalable procedure with low detection delay. Empirically, we were able to outperform previous approximate procedure (adaptive algorithms) and approach the performances of prohibitively expensive ones (GLR). Interesting future direction include taking into account the delay caused by the generalisation error.

heading: A SOLUTION FOR THE QCD PROBLEM UNDER KNOWN PARAMETERS
text: In this section, we assume that both the parameters θ0 and θ1 in Equation (1) are known. We present in the following the main optimality existing results.
Min-Max formulation: In the i.i.d. case (i.e. when fθ(.|x) = fθ(.)), the CUSUM algorithm (Lorden et al., 1971) is an optimal solution (Moustakides et al., 1986; Ritov, 1990) for the Lorden formulation (3). However, even though strong optimality results hold, optimising WADD is pessimistic. With the delay criterion base on CADD (Pollak formulation), algorithms based on the SHIRYAEVROBERTS statistic are asymptotically within a constant of the best possible performance (Pollak, 1985; Tartakovsky et al., 2012). In the non i.i.d. case, state of the art methods (Lai, 1998) are modified versions of the CUSUM algorithm that converge asymptotically to a lower bound of CADD (and thus a lower bound to WADD).
Bayesian formulation: The SHIRYAEV algorithm (Shiryaev, 1963) is asymptotically optimal in the i.i.d. case. This result is extended to the non i.i.d. case under Hypothesis 1 (Tartakovsky & Veeravalli, 2005).

heading: A.1 ALGORITHMS
text: In this section, we formally introduce the algorithms used to solve the change point problems under known parameters. A common trait of these algorithms is the computation of a statistic Sθ0,θ1n and the definition of a cutting threshold Bα. The stopping time να is the first time the statistic exceeds the threshold. The value Bα is chosen such that PFA (respectively FAR) of the associated stopping time does not exceed the level α. In the case of the BAYESIAN formulation, the threshold is simplified into Bα = 1−αα .
The SHIRYAEV Algorithm: The general formulation of the statistic Sθ0,θ1n , is the likelihood-ratio of the test of H0 : λ ≤ n versus H1 : λ > n given the observations Xt≤n. In the general case, with a prior P(λ = k) = πk, this statistic writes as:
Sθ0,θ1n :=
∑ k≤n πk ∏k−1 t=1 fθ0(Xt+1|Xt) ∏n+1 t=k fθ1(Xt+1|Xt)∑
k>n πk ∏n−1 t=1 fθ0(Xt+1|Xt)
(12)
In the case of geometric prior with parameter ρ, the statistic is simplified into:
Sθ0,θ1n = 1 (1−ρ)n ∑n k=1(1− ρ)k−1 ∏n t=k Lt(θ0, θ1) where Lt(θ0, θ1) = fθ1 (Xt+1|X t)
fθ0 (Xt+1|Xt)
The statistic Sθ0,θ1n can be computed recursively under this simplification:
Sθ0,θ10 = 0 and S θ0,θ1 n+1 = 1+Sθ0,θ1n 1−ρ Ln+1(θ0, θ1)
The SHIRYAEV-ROBERT Algorithm: The statistic to be computed in this case can be seen as an extension of the one in the SHIRYAEV algorithm when ρ = 0. The recursive formulation is indeed:
Sθ0,θ10 = 0 and S θ0,θ1 n+1 = (1 + S θ0,θ1 n )Ln+1(θ0, θ1)
The CUSUM Algorithm: The relevant statistic is defined as: Sθ0,θ1n = maxk≤n ∑n t=k log(Lt), and can be computed recursively by:
Sθ0,θ10 = 0 and S θ0,θ1 n+1 = (S θ0,θ1 n + log(Ln+1)) +

heading: B R-QUICK CONVERGENCE HYPOTHESIS
text: Hypothesis 1 is not too restrictive, and r-quick convergence conditions were previously used to establish the asymptotic optimality of sequential hypothesis tests for general statistical models (Draglia et al., 1999; Lai, 1976; 1981).
Using the law of large numbers, this hypothesis is satisfied for q = EXt∼fθ1 [ fθ1 (Xt) fθ0 (Xt) ] in the i.i.d. case. In the (ergodic) Markovian case, the convergence is achieved for irreducible Markov chains due to their ergodicity.
However, the speed of convergence of Equation (4), depends on the speed of convergence to the stationary distribution of the Markov chain. Even though for any r > 0 construction of Markov chains that satisfy assumption 4 is possible (Draglia et al., 1999), providing a guarantee of this speed is not a trivial question (Hairer, 2010).

heading: C ALGORITHMIC ANALYSIS
text: C.1 IMPACT ANALYSIS (ANNEALING AND PENALISATION)
In order to justify experimentally the use of annealing and penalisation to solve the issues highlighted with Lemma 3, we simulate simple examples where the phenomena are amplified to their extreme.
Annealing: We consider a one dimensional signal. The KL divergence between the pre and post change distribution is 0.06. As discussed in the paper, the issue is that after the change point we start using observations from the post change distribution to learn the pre change one. This leads to θt0 and θt1 converging to the post change parameter θ ∗ 1 . In fact, without the annealing procedure, the KL divergence DKL(fθt0‖fθt1) (the blue curve, in the top right subplot in Figure 3) ends up collapsing to 0. This has the problematic consequence of slowing down the statistics. In fact, both the SHIRYAEV and the CUSUM statistics (the blue curves in the lower subplots of Figure 3) start adjacent to the optimal statistic (in red) and slowly degrade in quality. On the other hand, implementing the annealing procedure -as described in the paper with a step size = 0.01- solves this issue. The KL divergence (green curve) hovers around the true value and the computed statistics are almost within a constant of the optimal one.
Penalising: We consider in this case a one dimensional signal with a KL divergence of 2. The problem we are analysing in this setting is the exclusive use of pre change observations when learning the parameters before the change point. The optimal log-likelihood ratio (LLR) before the change point is −2 (the red curve in the top right subplot of Figure 4) while the learned one (in blue) is
around 0. This is due to θt0 and θ t 1 converging to the pre change parameter θ ∗ 0 . This over-estimation of the LLR leads to a higher estimate of the SHIRYAEV statistic. In the lower left subplot of Figure 4, we have an increase of 100 fold compared to the optimal value. Penalising the LLR, using a coefficient c = 0.01, solves this issue without being detrimental to the post change performance. In fact the penalised LLR (the green curve) is strictly negative before the change point and has a similar values to the optimal ones afterwards. The associated SHIRYAEV statistic (the green curve) is less prone to detect false positives as it has the same behaviour as the optimal one. This safety break (penalisation) comes with a slight drawback for higher cutting thresholds as it induces an additional delay due to the penalising component. Interestingly, the CUSUM statistic does not have the same behaviour in this example. This is explained with the nature of the LORDEN formulation. As it minimises a pessimistic evaluation of the delay (WADD), the impact of a LLR of 0 is less visible. However, on average, the same phenomena occurs. This will be observed in the ablation analysis, where on average, there was no notable differences between the different formulations.

heading: C.2 ABLATION ANALYSIS (ANNEALING AND PENALISATION)
text: Since the algorithm integrates different ideas to mitigate the issues discussed in the paper, an ablation study is conducted to understand their contribution. We simulate the case where X = R10, Θ = R10, and (φµ, φσ) are 5-layer deep, 32-neurons wide neural network. The penalisation coefficient is fixed to 0.01 and the annealing step is fixed to 0.02. In Figure 5, the regretted detection delay is evaluated using 500 simulations. We consider a mid-ground complexity case, with a KL divergence of 1.5, which is a more realistic scenario in real life situations.
We provide the regret with respect to both the optimal performance under known parameters for both the BAYESIAN setting (SHIRYAEV statistic on the left) and the MIN-MAX setting (CUSUM statistic on the right). We split the analysis of the curves to three scenarios: The low, the mid and the high cutting threshold cases. They correspond respectively to 3 categories of risk inclinations levels: the high, the mid and the low.
The high risk inclination: In the presented case, this correspond to a cutting threshold B smaller than 1020 in the shiryaev algorithm and than 50 in the CUSUM algorithm. This can be associated to household applications and relatively simple settings with high reaction time. In this case, the
simplest version of the algorithm is sufficient to achieve relatively small regret values. In fact, the green curve (associated to the configuration without annealing) presents a relatively low regret with respect to the optimal delay.
The low risk inclination: In the presented case, this correspond toB higher than 1050 in the shiryaev algorithm and than 100 in the cusum algorithm. This can be associated to critical applications where precision outweighs the benefice of speed such as energy production. In this case, configurations without annealing (presented in green here) have poor performance. The best configuration is to perform annealing but without penalisation of the log likelihood.
The intermediate risk inclination: This coincides with the remaining spectrum of cases. For example, when the application is critical but requires high reaction time, or when the application has no reaction time constraints but doesn’t require extreme safety measures. In this setting, we observe that the best performance require some amount of penalisation of the log likelihood. The choice of the coefficient should be guided by the KL divergence between the pre and post change distributions, with c = 0 when the change is subtle (small KL).

heading: C.3 DIVERGENCE MEASURES BASED ALGORITHMIC VARIANT
text: In order to solve the QCD problem under unknown parameters, we propose in this paper to approximate the log-likelihood ratio L∗t efficiently. We achieved this by optimising a surrogate of the KL divergence between the estimated parameters (θ0t , θ 1 t ) and the true parameters (θ ∗ 0 , θ ∗ 1). Algorithm 1 devises a theoretically grounded weighting technique that approximates expectations under pre/post change distributions. For this reason, other divergence measure based loss functions can be used instead of the proposed ones in Equation 11. In particular, consider the following loss functions:
Ln0 (I, θ0, θ1, g) = ∑ t∈I P(τ < t|τ ∼ fθ0,θ1n )g( fθ1 fθ0
)(Xt+1|Xt) Ln1 (I, θ0, θ1, g) = ∑ t∈I P(τ > t|τ ∼ fθ0,θ1n )g( fθ1 fθ0 )(Xt+1|Xt)
For k ∈ {0, 1}, the loss function Lnk (I, θ0, θ1, g) approximates Ek [ g( fθ1 fθ0 ) ] . Under the assumptions that g is an increasing function with g(1) = 0 and h(x) = x.g(x) is convex, Ln0 becomes a proxy of the h-divergence measure Dh(fθ∗1‖fθ∗0 ) (respectively Dh(fθ∗0‖fθ∗1 ) for L n 1 ). This formulation is coherent with the QCD objective as it minimises Lt under the pre-change distribution and maximises it under the post change distribution and we can use it as a valid change point detection algorithm. The obtained parameters (θ0t , θ 1 t ) are the furthest apart with respect to the associated h-divergence. The KL divergence (associated with g(x) = log(x)) is a natural fit for our setting as it simultaneously ensures that the evaluated parameters converge to (θ∗0 , θ ∗ 1).
We keep using the same hyper-parameters from the ablation analysis in order to asses the ADD of the TWR algorithm using different divergence measures. The experimental results are reported in Figure 6. In addition to the KL-based version (dark blue curve), we consider the case where g(x) = √ x− 1 (purple curve) and the case where g(x) = (x− 1). log(x) (cyan curve). We also provide the adaptive (yellow) and GLR (green) average detection delay as a baseline.
The KL-divergence based approach provides comparable ADD to the GLR, outperforming both the adaptive algorithm and the other variant of the TWR.

heading: D PROOFS OF TECHNICAL RESULTS
text: Proof of Lemma 1: We start by observing the following: f0(θ, θ1) = E0 [ log( fθ1 fθ∗0 ) ] + E0 [ log( fθ∗0 fθ ) ] f1(θ0, θ) = E1 [
log( fθ∗1 fθ0
) ] − E1 [ log( fθ∗1 fθ ) ]
Notice that both quantities E0 [
log( fθ1 fθ∗0
) ] and E1 [
log( fθ∗1 fθ0
) ] are constants with respect to the param-
eter θ, and that:  E0 [ log( fθ∗0 fθ ) ] = DKL(fθ∗0 |fθ) E1 [ log( fθ∗1 fθ ) ] = DKL(fθ∗1 |fθ)
where DKL is the Kullback–Leibler divergence in the i.i.d case and the Kullback–Leibler divergence rate in the Markov chains case, i.e.
DKL(fθ0 |fθ1) =
{ ∫ x fθ0(x) log( fθ0 fθ1
)(x)∫ x Π0(x) ∫ y fθ0(y|x) log( fθ0 fθ1 )(y|x)
When well defined, DKL(fθ0 |fθ1) is strictly positive except when fθ0 = fθ1 . This allows us to conclude that: {
arg minθ f0(θ, θ1) = arg minθDKL(fθ∗0 |fθ) = θ ∗ 0 arg maxθ f1(θ0, θ) = arg maxθ −DKL(fθ∗1 |fθ) = θ ∗ 1
This concludes the proof of the Lemma.
Proof of Lemma 2: We denote µ0 the initial distribution according to which X0 is sampled and µt the distribution of Xt. We re-write the functions fnk as:
fn0 (θ0, θ1) = ∑t2 t=0
P(t<λ|να=n)∑t2 i=0 P(i<λ|να=n) E0,t[log( fθ1 fθ0 )(X|Xp)]
fn1 (θ0, θ1) = ∑t2 t=0
P(t>λ|να=n)∑t2 i=0 P(i>λ|να=n) E1,t[log( fθ1 fθ0 )(X|Xp)]
where Ek,t is taken with respect to (Xp, X) ∼ µt(Xp)× (fθ∗0 (X|Xp)1t<λ + fθ∗1 (X|Xp)1t>λ). It’s important to notice for what follows that P(t < λ|να = n) (respectively P(t > λ|να = n)) is
decreasing (increasing) with respect to t. We also highlight that for any finite set I , we have the following: ∑
t∈I P(t<λ|να=n)∑t2 i=0 P(i<λ|να=n) n,t2→∞−−−−−→ 0∑ t∈I
P(t>λ|να=n)∑t2 i=0 P(i>λ|να=n) n,t2→∞−−−−−→ 0
Consider the case where λ =∞. As Xt is a irreducible Markov chain, we have that
||Π∗0 − µt||T.V t→∞−−−→ 0
where ||.||T.V is the total variation norm. This means that for any small value > 0, there exist T > 0 such that:
∀t ≥ T ||Π∗0 − µt||T.V ≤ Thus, as the observation Xt are bounded, we have for any continuous function f , there exist T > 0
such that: ∀t ≥ T |Ek,t[f(X|Xp)]− E0[f(X|Xp)]| ≤ This implies that when λ =∞ then f∞0 and f∞1 converges to f0. In fact for any , there exist N and T2 such that for all n > N :
∀t2 ≥ T2 ∑ t<T
P(t>λ|να=n)∑t2 i=0 P(i>λ|να=n) ≤
∀t2 ≥ T2 ∑ t<T
P(t<λ|να=n)∑t2 i=0 P(i<λ|να=n)
≤ ∀t2 ≥ t ≥ T |Ek,t[f(X|Xp)]− E0[f(X|Xp)]| ≤
The case where λ is finite can be deduced by considering the sequence Xt>λ. All the observation are sampled according to fθ∗1 , and thus f ∞ 0 and f ∞ 1 converges to f1.
This result is also valid for fνα1 . In fact, as P(t > λ|να = n) is increasing with respect to t, then for any fixed horizon H we have:∑
t<H
P(t > λ|να = n)∑t2 i=0P(i > λ|να = n) t2→∞−−−−→ 0.
As this remains true for H = T , we obtain that fνα1 converges to f1.
For a given integer n, if t2 = λ+ n, the convergence of fνα0 to f0 is achieved as λ→∞. In order to establish this result, we need to prove that ∑t2 t=λ
P(t<λ|να)∑t2 i=0 P(i<λ|να) is decreasing with respect to λ.
This is a consequence of P(t < λ|να) being decreasing with respect to t. In fact:∑t2 t=λP(t < λ|να) ≤ ∑t2 t=λP(λ < λ|να) = nP(λ < λ|να)∑t2
t=0P(t < λ|να) ≥ ∑λ t=0P(λ < λ|να) = λP(λ < λ|να)
as a consequence the following inequality holds: t2∑ t=λ P(t < λ|να)∑t2 i=0P(i < λ|να) ≤ n λ λ→∞−−−−→ 0.
As such, for any > 0, there exists An > 0 such that: ∀λ ≥ An ∑t2 t=λ
P(t<λ|να)∑t2 i=0 P(i<λ|να) ≤
∀λ ≥ T ∑T t=0
P(t<λ|να)∑t2 i=0 P(i<λ|να)
≤ ∀λ ≥ t ≥ T |E0,t[f(X|Xp)]− E0[f(X|Xp)]| ≤
Proof of Lemma 3: We have fλk converges to f0 as λ → ∞ (respectively f tk converges to f1 as t → ∞). Given that θ∗0 minimises f0 (respectively θ∗1 minimises f1), then θλk and θtk satisfy the following:
θλ0 , θ λ 1 a.s−→ λ→+∞ θ∗0 and θ t 0, θ t 1 a.s−→ t→+∞ θ∗1 .
As such, we obtain the following result:
Lλ a.s−→
λ→+∞ 0 and Lt a.s−→ t→+∞ 0
On the other hand we have the following:
limλ→∞Eλ[L ∗ λ] = EXλ−1∼Π∗0 ,Xλ∼fθ∗0 [L∗λ] = −DK.L(fθ∗0‖fθ∗1 ) < 0 limt→∞Eλ[L
∗ t ] = EXt−1∼Π∗1 ,Xt∼fθ∗1 [L∗λ] = DK.L(fθ∗0‖fθ∗1 ) > 0
This concludes the proof of the lemma.

heading: E ADDITIONAL EXPERIMENTAL RESULTS

heading: E.1 HYPER-PARAMETER SELECTION
text: The algorithm we designed requires the selection of a few hyper-parameters in order to run properly. In this section we address the issue of tuning them.
The first set of parameters are for optimisation purpose, and thus we advise selecting them according to the complexity of the probability distribution family {fθ, θ ∈ Θ}. In all the experimental settings we used a batch size of 32, a number of epoch equal to 25 and a gradient step size of 0.001. The initial task parameters θ0 and θ1 are chosen randomly unless stated otherwise.
As for the cutting threshold, the penalisation and the annealing coefficients, they depend on the KL divergence between the pre- and post- change distribution (thus the task sampling distribution F ) and on the mixing times of these distribution. Fine tuning using grid search is the most efficient way to identify suitable candidates.
From our experimental inquiries, we advise a low penalisation coefficient (e.g. c = 0.001) unless the threshold Bα is smaller than 100 for the SHIRYAEV statistic and smaller than 10 for the CUSUM statistic. The annealing parameter depends on the mixing time of the Markov chain. In fact reflects to what extent we keep learning the pre-change parameter once it’s similar to the post-change one. In practice we found out that a small coefficient (e.g. = 0.01) is advised when the observations converge slowly to the stationary distribution. On the other hand, when convergence occurs in few observations, a higher coefficient is a safer choice.

heading: E.2 COPYCAT AGENT PROBLEM: COMPLEMENTARY RESULTS
text: In this section we further the analysis of the copycat agent problem. The Fetch set of environments provides a good mix of complexity and interpretability. The ability of the TWR algorithm to detect changes and to solve the copycat agent problem in these environments is proof of it’s robustness.

heading: SHARED BEHAVIOURS
text: When different tasks are associated to a shared behavior, it becomes difficult to evaluate the parameter θ∗t . This is true for the push and the PickandPlace environments. In both cases, the goal is to control a robotic hand (through its joints’ movements) in order to interact with an object until it reaches a final position (either push it there or pick it and then place it there).
The task space in these problems is the set of possible final positions. No matter what the task is, the main agent is going to reach to the object. The observations associated to this intermediate process of reaching are probably going to yield a poor estimation of the final position. However as the agent start to interact with the object, it starts to become clear what task is being executed. For this reason, an artificial change is detected between the reaching behaviour and the interaction one when using the TWR algorithm. Even though this detection is unwanted (in the sense that there was no actual change of behavior) it allows to reduce the estimation error ‖θt0 − θ∗t ‖. In this section, the experiment we run consists of performing three randomly sampled tasks, each over a window of 20 observations. We estimate the running task θ∗t using the TWR algorithm and a maximum log-likelihood estimator (MLE) over both the last 5 and 10 observations. We average the errors over 500 trajectories and report the results in Figure 7. In both environments estimating the task with the running pre-change parameter of the TWR algorithm outperforms MLE. The gap grows narrower however in the PickAndPlace setting as the shared behaviour last longer. This makes estimations more difficult.
We also notice that the first task estimation of TWR is better than the MLE estimator despite the fact that all the observations are generated with a single parameter. This strange result is due to the detection of an artificial change after the "shared" behaviour. optimising only observations that are clearly correlated to a particular task, leads to better estimations. As for the gap in the first estimation, this is explained by a built-in agnostic behaviour in the TWR algorithm. When the number of observation is less than the ADD of the SHIRYAEV algorithm, the pre-change weight term P(τ < t|τ ∼ fθ t 0,θ t 1
n ) of Ln0 is extremely small. and thus the estimation remains close to the initialisation. In addition, the tow target functions Lt0 and Lt1 are adversarial as they kind of optimise the opposite of each other. This leads to smaller gradient updates. On the other hand, MLE tend to commit to a particular estimate using observations from the ’shared’ behaviour phase. This leads generally to a heavily over-fitting estimate.
The agnostic behaviour is more pronounced in Figure 7b (PickAndPlace) as the robotic hand picks the object exactly the same way no mater what task is being solved. As long as there is few evidence of a distinctive behaviour occurring, the TWR estimate remains close to the initialisation. In fact the tasks here are positions and the error is the distance between the true and predicted coordinates. In the PickAndPlace scenario, the target is randomly sampled in a cube of edge length 0.3, in the Push scenario, it is sampled in a square with the same edge length. with random initialisation we get an average error of 0.1× √ 3 and 0.1× √ 2 respectively with a random initialisation. This coincide with the TWR error (blue curve) in both environments given few observations proving that the learned task remained close to the initialisation.
In Figure 7a, the improvement due to the artificial detection of change between the shared behaviour and the task specific one is marked around the 5th observation. In fact once the robotic hand starts pushing the object in a particular way, the estimation error starts to decay faster.

heading: PERFORMANCE COST OF LEARNING THE POLICY
text: Estimating the running task θ∗t in the copycat problem with either the TWR algorithm or the MLE requires access to the probability distribution π∗θ . However, in real life situations, we can only learn an approximation of this distribution through historical observations.
In this section we evaluate experimentally the impact of using inverse reinforcement learning to construct an approximate policy π̂θ and using it as a substitute for π∗θ . We consider the Reach environment where the task is to move the robotic hand to a particular position. The main agent execute 5 different tasks, sampled randomly, each over a window of 20 time steps. We evaluate the running task with TWR and MLE using the last 5 observations using the actual policy π∗θ and the learned one π̂θ. We average the estimation error over 100 simulations and we provide the experimental results in Figure 8. The approximate policy is constructed using the GAIL algorithm and a data-set of task labelled observations generated using π∗θ .
As established before, the TWR estimates using the true policy (in blue) outperform MLE based approaches. This trends remains true on average when using π̂θ. However, the true policy based MLE (in green) converges to a better estimation in the last observations before the change point, while the GAIL based TWR estimate (in red) as well as the GAIL based MLE (in orange) seem to be unable
to improve. This is the cost of using a policy approximation. The additional error is built in due to estimation error of the learned policy.

heading: E.3 COMPLEMENTARY PERFORMANCE ANALYSIS
text: In this section we provide complementary experimental results to the ones introduced in the paper. We keep using the same hyper-parameters (i.e. X = R10, Θ = R10, and (φµ, φσ) are 5-layer deep, 32-neurons wide neural network).

heading: SINGLE CHANGE POINT PERFORMANCE
text: In Figure 9, we provide the average regret overs 500 simulations in both the BAYESIAN (on the left) and the MIN-MAX (on the right) settings. Each simulated trajectory is 1000 observation long with a change point at the 500th observation. We use the full pipeline of our algorithm with a penalisation coefficient (c = 0.01), and we allow the adaptive algorithm to exploit 10% of the pre change observations. Our algorithm (the blue curve) achieves comparable performance to the GLR statistic (without requiring intensive computational resources) and achieves better and more stable results compared to the adaptive algorithm (without requiring domain knowledge).
In all experiments and for all considered change detection algorithms, we observe a spike and a high variance of the average detection delay for low cutting threshold. This is a direct consequence of the definition of Bα. In fact a low cutting threshold implies a high tolerance of the type I error. This implies that the optimal detection time is no longer a reliable estimate.

heading: MULTIPLE CHANGE POINTS PERFORMANCE
text: We reconsider in this section the multiple change points problem introduced early on in the paper. If the change points are sufficiently separated (i.e. tk − tk−1 is big enough), they can be seen as independent single change point problems. However there is no good reason to believe that this is the case of real life applications. For this reason, the detection delay of the first change point t1 will probably reduce the accuracy at which we estimate the parameter θ1, this in turn will affect the accuracy of estimating the log-likelihood for the next change point t2. This behaviour will keep on snowballing until we reach a breaking point at which we will miss a change point. Some attempts to deal with this issue have been made in the past to improve the convergence rate of the GLR algorithm in order to reduce the impact of this problem (Ross, 2014). These approaches remain computationally extensive and thus are not considered in our comparison setting.
We provide however an analysis of the impact of the number of pre change observations on the average detection delay in the single change point setting. We keep using the same hyper-parameters where X = R10, Θ = R10, and (φµ, φσ) are 5-layer deep, 32-neurons wide neural network. We use parameters achieving a KL divergence of 1.5 between the pre and post change parameters, and
evaluate the average detection delay over 500 simulations where t2 = λ + 250. By varying λ we simulate different multiple change point scenarios where the accumulated delay exhausts most of the available observations. The first thing to notice is that our algorithms (in the blue curve) has similar average detection delay to the GLR (green curve) in all presented cases. When we have sufficient observations (λ = 250), all algorithms will have reasonable performances compared to the optimal delay. However, for smaller horizon, the adaptive algorithm (the orange curve) is predicting change way before the optimal algorithm. This is explained with the bad estimation of the pre change distribution due to a reduced sample size. In fact for cutting threshold up to 1025, the adaptive algorithm is predicting the change point before it happens. In addition, we observe that our approach varies less than the adaptive one as it’s less prone to fit statistical fluctuations.
We provide an analysis of the impact of the number of pre change observations on the average detection delay in both the BAYESIAN (on the left) and the MIN-MAX (on the right) settings. We use parameters achieving a KL divergence of 1.5 between the pre and post change parameters, and evaluate the average detection delay over 500 simulations where t2 = λ + 250. By varying λ we simulate different multiple change point scenarios where the accumulated delay exhausts most of the available observations. There is no notable variation in the algorithms behaviour when using either the SHIRYAEV or the CUSUM statistic. As for the impact of the number of available pre-change observation before the change point (λ), the results presented in Figure 10 illustrate how our algorithm is less prone than adaptive procedures to over fit the observation when their number is limited. In fact the blue curve (our algorithm) is closer than the adaptive one (orange curve) to the GLR (green curve) performances. All approaches end up having a constant regret with respect to the optimal detection given enough observations (λ = 250).

heading: PERFORMANCE GIVEN A NOMINAL BEHAVIOUR
text: In this section, we keep using the same hyper-parameters from the previous experiments. However we consider the setting where we have access to the pre-change parameter (which correspond to a nominal behaviour). Beside the ADD of the GLR and TWR statistics (computed with no prior knowledge of the pre-change parameter), we computed the ADD of the adaptive and TWR statistics using the prior knowledge of the true pre-change distribution.
The first notable thing is that the adaptive algorithm performances (yellow curve) did not improve by much compared to the case where the parameter were estimated using 10% of the available observations. In fact TWR statistics outperforms it even with no prior knowledge of the nominal behaviour (dark blue curve). The second notable observation, is that given the pre-change parameter, TWR statistic performance (light blue) is within a constant of the GLR statistic. These observations confirm again that using the asymptotic behaviour of the Shriyaev detection delay is a very powerful approximation

Review summaries:
This paper studies the quickest change detection for Markovian data, when both the parameters of pre- and post-change distributions are unknown. The main contribution is a scalable algorithm that sequentially estimates the unknown parameters and plug-in to classical detection schemes to get the stopping rule. A notable feature is that this is a joint estimation and detection framework. And the authors incorporate several tools, like SGD, annealing, penalization, into the detection task, which turns out to have good performance compared with existing benchmarks. 

The author(s) propose a quickest change detection technique under known parameter scenario. They use a Markovian dynamics to generate the pre and post change-point distributions and use a Shirayev test-statistic based on the asymptotic behavior of the optimal delay under know parameters. The proposed methodology is validated on synthetic data and a multitask reinforcement learning example. There are several issues which restricts the paper to reach an optimal level. These are highlighted below