# Adaptive Partially-Observed Sequential Change Detection and Isolation

Xinyu Zhao1, Jiuyun Hu1, Yajun Mei2, Hao Yan1
1School of Computing and Augmented Intelligence
Arizona State University,
2School of Industrial and Systems Engineering,
Georgia Institute of Technology
The authors gratefully acknowledge the support from NSF DMS 1830363, and CMMI
1922739

###### Abstract

High-dimensional data has become popular due to the easy accessibility of
sensors in modern industrial applications. However, one specific challenge is
that it is often not easy to obtain complete measurements due to limited
sensing powers and resource constraints. Furthermore, distinct failure
patterns may exist in the systems, and it is necessary to identify the true
failure pattern. This work focuses on the online adaptive monitoring of high-
dimensional data in resource-constrained environments with multiple potential
failure modes. To achieve this, we propose to apply the Shiryaev–Roberts
procedure on the failure mode level and utilize the multi-arm bandit to
balance the exploration and exploitation. We further discuss the theoretical
property of the proposed algorithm to show that the proposed method can
correctly isolate the failure mode. Finally, extensive simulations and two
case studies demonstrate that the change point detection performance and the
failure mode isolation accuracy can be greatly improved.

Keywords: Shiryaev–Roberts procedure, multi-arm bandit, sequential change-
point detection, adaptive sampling, multiple failure modes

##  1 Introduction

Nowadays, most industrial applications are instrumented with hundreds or
thousands of sensors due to the advancement in sensing technology. Real-time
process monitoring and fault diagnosis are among the benefits that can be
gained from effective modeling and analysis of the produced high-dimensional
streaming data. Classical researches for process monitoring of high-
dimensional streaming data focus on a fully observable process, which means at
each sampling time point, all the variables can be observed for analysis (Yan
et al., 2018). However, it is often infeasible to acquire measurements of all
these sensing variables in real time due to limited sensing resources, sensing
capacity, sensor battery, or other constraints such as system transmission
bandwidth, memory, storage space, and processing speed in modern industrial
applications (Liu et al., 2015). Furthermore, under change detection and
isolation setting, we assume that the engineered systems that are being
studied have several distinct failure modes and patterns but do not know which
failure mode may occur beforehand. Overall, this paper focuses on change point
detection under resource-constrained environments with multiple potential
failure modes.

The first motivating example is in the hot forming process (Li and Jin, 2010)
as shown in Fig. 1(a). There are five sensing variables in the system: the
final dimension of workpiece ${\bf X}_{1}$, the tension in workpiece ${\bf
X}_{2}$, material flow stress ${\bf X}_{3}$, temperature ${\bf X}_{4}$, and
blank holding force ${\bf X}_{5}$. These five variables can be represented as
a Bayesian network, as shown in Fig. 1(a). For example, if we know that the
change of ${\bf X}_{4}$ and ${\bf X}_{5}$ are the two major failure sources in
the system. If ${\bf X}_{4}$ changes, $({\bf X}_{1},{\bf X}_{2},{\bf
X}_{3},{\bf X}_{4})$ will also change. Furthermore, If ${\bf X}_{5}$ changes,
only $({\bf X}_{1},{\bf X}_{2},{\bf X}_{5})$ will change. Therefore, different
failure modes may affect a different subset of sensors differently.

(a) Hot foaming process

(b) 3D printing example

Figure 1: Examples of Complex data in Various Industrial Applications Left
figure shows an example of a hot foaming process. Right figure shows an
example of monitoring the thermal images in additive manufacturing.

Another example comes from in-situ hot-spots detection in the laser powder bed
fusion (LPBF) process in the metal additive manufacturing process. A thermal
camera is often used to monitor the stability of the process while the product
is being produced on a layer-by-layer basis. Here, detecting the hot-spots
early is crucial for further product quality control. Fig. 1(b) show an
example of such hot-spots from the thermal camera. Given that the anomaly or
hot-spots can only occur on the edge/corner of the scanning path, multiple
failure modes can be defined.

There are a few challenges of sequential change-point detection under the
sampling constraint: 1) From the previous examples, the failure mode
distribution can be quite complicated. For example, in the hot foaming
process, as shown in Fig. 1(a), we aim to detect the failure mode with the
weakly conditional dependency on the graph; In the laser powder bed fusion
process, as shown in Fig. 1(b), we aim to detect the spatially clustered hot-
spots. 2) Even though we assume that we have prior knowledge of different
potential failure modes, we do not know which failure mode may occur in the
system. The main challenge is to balance the exploration of all potential
failure modes and the exploitation to focus on the most probable failure mode.
A conceptual illustration of the proposed algorithm is provided in Fig. 2. The
illustration example has shown an example that the sampled points are
performed on the 2D spatial domain. The sampling patterns at time $t_{1}$ and
$t_{2}$ focus on exploration for all failure modes and the sampling patterns
at $t_{3},\cdots,t_{n}$ focus on exploitation for failure mode 3. In general,
it is hard to decide when the algorithms should change to exploitation or
which failure mode they should focus on. Finally, given that the multiple
failure modes have quite complex shapes and distributions, the exploration and
exploitation among these modes are often quite challenging.

Figure 2: Conceptual Illustration of the Balance of Exploration and
Exploitation; The sampling patterns at $t_{1},t_{2}$ focus on the exploration
of all failure modes. The sampling pattern at $t_{3},\cdots,t_{n}$ focuses on
the exploitation of the failure mode 3.

There are also many works focusing on change-point detection under resources
constraint. Most of the existing works are proposed based on the ”local
monitoring and global decision” framework, which focuses on monitoring each
data stream independently using local monitoring statistics and then fusing
these local monitoring statistics together via a global decision framework.
For example, Liu et al. (2015) proposed scalable and efficient algorithms for
adaptive sampling for online monitoring. The method introduced a compensation
parameter for the unobserved variables to increase the chance of exploring
them. Recently, Zhang and Mei (2020) proposed to combine the powerful tools of
the multi-arm bandit problem for efficient real-time monitoring of HD
streaming data. However, these works either assume the data stream is
independent or cannot take advantage of the failure mode information in some
systems, which fails to monitor and identify the correct failure pattern. For
a complete literature review of monitoring of high-dimensional streaming data,
please see Section 2.

To generalize the sequential change-point detection framework to both detect
and identify the correct failure modes, change detection and isolation
literature has been proposed in the literature, which also inspires this
research. The change-point detection and isolation often assume that there are
a set of pre-defined post-change distributions. The goal is not only to detect
the change with the shortest detection delay but also to identify which change
mode occurs in the system. For example, Chen et al. (2020) proposed a Bayesian
method to decide on a procedure to identify both the change point as well as
the correct change mode. For a complete literature review of change detection
and isolation, please see Section 3.3. However, these works typically assume
that the data is fully observed, which cannot be applied to partially observed
data.

To address the challenge of multiple failure modes and partially observed
data, we propose a novel Multiple Thompson Sampling Shiryaev-Roberts-Pollak
(MTSSRP) Method by a modified ”local monitoring and global decision
framework”. As far as the authors know, this is the first work that discusses
the adaptive sampling framework for failure mode detection and isolation.
Unlike the literature on the monitoring of HD streaming data, where the local
monitoring statistics are defined at each individual sensor, we propose to
define the local statistics for each individual failure mode. This enables the
proposed MTSSRP to take advantage of the failure mode information, which is
very important in the high-dimensional space, given that failure mode
information can significantly reduce the search space since there are
unlimited ways that change may occur in the high-dimensional space. To
quantify the uncertainty of unobserved sensing variables for different failure
modes, we propose to apply the Shiryaev-Robert (SR) procedure for sequential
change point detection on the failure mode level.

Furthermore, to balance the exploration and exploitation, we will borrow the
idea from Multi-arm Bandit (MAB). MAB aims to sequentially allocate a limited
set of resources between competing ”arms” to maximize their expected gain,
where the reward function for each arm is not known. MAB provides a principled
way to balance exploration and exploitation. To apply MAB for change point
detection under the sampling constraint, we propose to use the SR statistics
of the selected failure modes as the reward function in the Multi-arm Bandit
(MAB) problem (Zhang and Mei, 2020). However, different from (Zhang and Mei,
2020), the selection of arm is on the sensor level, where the SR statistics is
defined on the failure-mode level. For high-dimensional data, specifying the
joint distribution of high-dimensional data can be very challenging.
Therefore, this paper will explore spatial structures for defining the failure
mode distributions as shown in 4.7. This paper also discussed that with the
independence assumption of the distribution variable, the computational
efficiency can be greatly improved.

The paper is organized as follows. In Section 2, we will review the existing
literature on change detection and isolation. We will also discuss works on
process monitoring with resource constraints. We further introduce our
proposed method and then discuss its property in Section 3. Then, we apply the
proposed approach to both the simulated data and evaluate its performance and
compare the existing methods in Section 5. Furthermore, we apply the proposed
method to two real cases in Section 6, respectively. Concluding remarks are
given in Section 7.

##  2 Literature Review

In this section, we will provide a more detailed review of statistical process
control or sequential change point detection methods. We will briefly classify
the methods for the following four categories: monitoring of independent HD
streaming data, monitoring of functional data or profile monitoring, process
monitoring with the resource constraint, change-point detection, and
isolation.

In the first category, monitoring the HD streaming data has often been treated
as monitoring the multiple independent univariate data streams. There are two
distinct frameworks for monitoring independent data streams in recent years.
First, the ”global monitoring” framework focused on directly designing the
global monitoring statistics for the process monitoring or change point
detection for high-dimensional data (Xie and Siegmund, 2013; Wang and Mei,
2015; Cho and Fryzlewicz, 2015; Chan et al., 2017). However, the global
monitoring framework is typically computationally inefficient for high-
dimensional data. Second, the ”local monitoring and global decision” framework
focus on monitoring each data stream independently using local monitoring
statistics and then fusing these monitoring statistics together via the global
statistics (Mei, 2010, 2011). The benefit is that these methods are typically
computationally efficient and can be scalable to high-dimensional data.
However, these methods are often limited to the independent data stream.
Finally, this framework is targeted for the case of fully observed data, which
may not be applicable under the resource constraint.

In the second category, profile monitoring techniques have been proposed to
tackle the complex spatial correlation structures. Dimensionality reduction
techniques, such as principal component analysis (PCA), are widely used.
Various types of alternatives such as multivariate functional PCA (Paynabar et
al., 2016), multi-linear PCA (Grasso et al., 2014), and tensor-based PCA (Yan
et al., 2015) are proposed. On the other hand, non-parametric methods based on
local kernel regression (Qiu et al., 2010) and splines (Chang and Yadama,
2010) are developed. To monitor the non-smooth waveform signals, a wavelet-
based mixed effect model is proposed in (Paynabar and Jin, 2011). However, for
both PCA-based methods and non-parametric methods, they typically assume that
the change alternative is not known. To utilize the anomaly structures, smooth
sparse decomposition methods have been proposed and utilize two sets of basis
functions, the background basis and anomaly basis, to represent the spatial
structures of the background and anomaly, which have been applied to smooth
profiles (Yan et al., 2017, 2018) and waveform profiles (Yue et al., 2017).
However, all the profile monitoring techniques assume that the complete
measurements are given and cannot be applied for HD data with partial
observations.

In the third category, many existing works focus on the change point detection
with the sampling constraint. Here, we will briefly classify the existing
monitoring methods with the sampling constraint into two categories,
monitoring the i.i.d data stream and monitoring the correlated data stream.
For monitoring the i.i.d data stream with the sampling constraint, Liu et al.
(2015) proposed a top-R-based adaptive sampling strategy as a combination of
random sampling in the in-control state and fixed sampling in the out-of-
control state. Another work by (Zhang and Mei, 2020) converts the problem into
a MAB framework and adaptively selects the sensors with Thompson Sampling.
Recent work by (Gopalan et al., 2021) provides an information-theoretic lower
bound for the detection delay. However, this method can be applied to multi-
dimensional problems, but cannot be applied to the case with multiple complex
failure-mode (i.e., after-change) distributions.

However, due to the i.i.d assumption, these methods might not be suitable for
data with complex distributions in reality. To deal with this problem, Xian et
al. (2018) proposed an adaptive sampling strategy that can handle the
correlated data generated from a multinomial distribution. For monitoring
correlated data streams with the sampling constraint, these methods can be
classified into monitoring data generated from the Bayesian Network and
spatial profile. For example, Liu et al. (2013) and Liu and Shi (2013)
proposed a sensor allocation strategy according to a Bayesian Network to
detect changes with multivariate $T^{2}$ control chart. Another work discussed
the problem when there is a spatial correlation among sensors and proposed a
spatial-adaptive sampling strategy to focus on suspicious spatial
clusters(Wang et al., 2018). However, these methods either consider the data
stream with spatial correlation (Wang et al., 2018; Ren et al., 2020; Gómez et
al., 2022) or modeled by the Bayesian Network structures (Liu et al., 2013;
Liu and Shi, 2013), which fails to apply to the problem with general failure
mode distributions as discussed in this paper.

Finally, there is a large amount of work focused on the case when there are
multiple failure modes, and it is necessary to identify the true failure while
detecting the changes. The problem of sequential change detection with
multiple failure modes is usually called change detection and isolation. The
goal is to find the best decision procedure that can control the false alarm
rate as well as the false isolation probability. The problem is of importance
since it is common in different applications like fault diagnosis, process
monitoring, and object identification (Nikiforov et al., 1993; Willsky, 1976;
Malladi and Speyer, 1999). The major works in change detection and isolation
can be categorized into Bayesian and non-Bayesian directions. Nikiforov (1995)
proposed a change detection/isolation framework as an extension of Lorden’s
results (Lorden, 1971) which follows non-Bayesian schema. Another two works
formulate the problem into a Bayesian version, which considers the change
point as a random variable (Chen et al., 2020; Malladi and Speyer, 1999).

##  3 Proposed Methodology

In this section, we will first describe the problem formulation of partially
observed multi-mode change detection based on high-dimensional (HD) streaming
data with sampling control in Section 3.1. We will review some relevant
methodology in Section 4. We will describe the proposed MTSSRP methodology in
Section 4.4. We will prove important properties of the proposed algorithms
about the average run length and failure mode isolation guarantee in Section
4.5. We will give the guidelines to select the tuning parameters of the
proposed MTSSRP method in Section 4.6. Finally, we will give a discussion and
several guidelines on selecting the failure mode distributions in Section 4.7.

###  3.1 Problem Formulation and Background

Suppose we are monitoring data stream $X_{j,t}$ for $j=1,\cdots,p$ and
$t=1,2,\cdots,T$, where $p$ is the number of dimensionality in the system and
$T$ is the monitored time length. We assume that the data streams follow joint
distribution $f_{0}$ before change as ${\bf X}_{t}\sim f_{0}$ for $t<\nu$. At
some unknown change time $\nu\in\\{1,2,\cdots,\\}$, an undesirable event
occurs and causes an abrupt change of the data stream into one or few failure
modes. For example, the after-change distributions $f_{k}\in\mathcal{F}$ can
be anyone from a family of distributions
$\mathcal{F}=\\{f_{1},\cdots,f_{K}\\}$. In another word, after the change,
${\bf X}_{t}\sim f_{k},k=1,\cdots,K$ for $t>\nu$. In other words, we do not
assume that we know which failure mode occurs in the system. Following the
change detection and isolation framework, we do assume that a single failure
mode $f_{k}$ may occur after the change. We will discuss the case with
multiple failure mode. Here, $f_{0},f_{1},\cdots f_{K}$ are the joint
distribution for all sensing variables. The sensing variables can also also be
correlated. Finally, in practice, we can set the magnitude of the joint
distribution as the interested magnitude of the change to be detected.

Furthermore, we assume that given the resource constraint, it is not possible
to observe all the data streams. For the partially observed data with sampling
control, the set of the observed data is denoted by
$\mathbf{y}_{t}=\\{X_{j,t},j\in C_{t}\\}$. Here, $C_{t}$ is the set of
observed sensor indices at time $t$, which can be selected online. In other
words, we can define $a_{j,t}$ as the binary variable denoting whether the
variable $j$ is observed at time $t$, $C_{t}=\\{j:a_{j,t}=1\\}$. Finally, the
sampling constraint is represented by $\sum_{j=1}^{p}a_{j,t}=q$ at each time
$t=1,2\cdots$, which means at each time $t$, only $q$ sensing variables can be
observed from all $p$ variables.

The objective of this paper is to design an efficient adaptive sampling
algorithm and the change point detection algorithm to automatically distribute
sensing resources according to the knowledge of the system failure modes such
that the change can be detected quickly as soon as it occurs and the
corresponding failure mode can be identified accurately while maintaining the
false alarm constraint.

##  4 Review of Relevant Methodology

In Section 4, we will review the formulations of the relevant methodology in
detail. Section 4.1 reviews the Bayesian decision framework, namely the
Shiyaev-Roberts (SR) procedure for sequential change-point detection methods.
Section 4.2 reviews the extension of the SR procedure for HD data monitoring.

###  4.1 Review of Shiyaev-Roberts statistics for uni-variate Sequential
Change Point Detection

We will first review the Bayesian decision approach for the sequential change-
point detection by Shiryaev (1963). Consider a univariate streaming data ${\bf
X}_{1},{\bf X}_{2}\cdots,$ where the distribution of the data changed from
$f_{0}$ to $f_{1}$ at some unknown time $\nu$. In another word, We assume that
${\bf X}_{1},\ldots,{\bf X}_{\nu}\overset{i.i.d}{\sim}f_{0}$ and ${\bf
X}_{\nu+1},\ldots\overset{i.i.d}{\sim}f_{1}$. The goal is to propose a
statistical decision policy to determine the stopping time $T$, where $T=t$
implies that a change has happened at time $t$.

Here, the goal is to find a decision policy such that the detection delay can
be minimized under the false alarm constraint. A Bayesian formulation on the
statistical decision policy has been proposed, where the true change point
$\nu$ is assumed to has geometric distribution as $P(\nu=t)=p(1-p)^{t-1}$.
When $p\rightarrow 0$, (Shiryaev, 1963) proposed a statistics
$R_{t}=\sum_{j=1}^{t}\prod_{i=j}^{t}\frac{f_{1}({\bf X}_{i})}{f_{0}({\bf
X}_{i})}$, which can be computed recursively as

$R_{t}=(R_{t-1}+1)\frac{f_{1}({\bf X}_{t})}{f_{0}({\bf X}_{t})}.$ (1)

For some pre-specified constant $A$, the decision to detect a change point in
distribution can be made as $T_{A}=\inf_{t}\\{R_{t}\geq A\\}.$ Pollak (1985,
1987) proved that this change point detection rules has an asymptotic minimax
property, which is to minimize
$\sup_{1\leq\nu<\infty}\mathbb{E}(T_{A}-\nu|T_{A}\geq\nu)$ under the
constraint that $\mathbb{E}(T_{A}|\nu=\infty)\geq B$ and $B\to\infty$.

Despite the statistical efficiency, SR procedure is defined for the change
point detection for univariate data, which cannot be used to monitor the high-
dimensional data.

###  4.2 Review of Thompson Sampling with the SR statistics

Zhang and Mei (2020) combined the Thompson Sampling with the SR statistics to
monitor high-dimensional with adaptive sampling control. Suppose for any time
$t$, the data is $p$ dimensional. Before change time $\nu$, the $p$th
dimension is $\overset{i.i.d}{\sim}f_{0}$, and after change time $\nu$, the
$p$th dimension is $\overset{i.i.d}{\sim}f_{1}$. Zhang and Mei (2020) proposed
a statistic for each dimension with

$R^{\prime}_{j,t}=\begin{cases}(R^{\prime}_{j,t-1}+1)\frac{f_{1}(X_{j,t})}{f_{0}(X_{j,t})}&j~{}{\rm
observed}\\\ R^{\prime}_{j,t-1}+1&j~{}{\rm unobserved}\end{cases}$ (2)

The decision of a change point of distribution is a threshold of top $r$
statistics, which is

$T=\inf\\{t\geq 1:\sum_{j=1}^{r}R^{\prime}_{j,t}\geq A\\}$

To determine which dimensions to observe, Zhang and Mei (2020) implemented
Thompson Sampling, where a random variable $R^{\prime}_{j,0}$ is sampled from
a prior distribution $G$ and added to $R^{\prime}_{j,t}$ with some
coefficients. The $q$ dimensions with the largest result are observed at time
$t+1$. Even though this algorithm provided a theoretical understanding of the
compensation terms in (Liu et al., 2015) from the Bayesian perspective, it is
still focused on the case where the failures occurred on individual sensor
level, which can not be applied to HD data with complex failure mode
distributions.

###  4.3 Review of Change Point Detection and Isolation

We will briefly review the change-point detection and isolation framework by
Chen et al. (2020). Consider a streaming data ${\bf X}_{1},{\bf X}_{2}\cdots,$
where the distribution of the data changed from $f_{0}$ to $f_{k}$ at some
unknown time $\nu$. In another word, We assume that ${\bf X}_{1},\ldots,{\bf
X}_{\nu}\overset{i.i.d}{\sim}f_{0}$ and ${\bf
X}_{\nu+1},\ldots\overset{i.i.d}{\sim}f_{k}$, where
$f_{k}\in\mathcal{F}=\\{f_{1},\cdots,f_{K}\\}$ is a known pre-defined set of
after-change distributions.

The goal of a change-point detection and isolation algorithm is to compute a
terminal pair $\delta=(T,\hat{k})$ online based on the observations ${\bf
X}_{1},{\bf X}_{2},\cdots$, where $T$ is the alarm time at which
$\hat{k}$-type change is identified. The goal is to propose a statistical
decision policy to determine the stopping time $T_{A}$, where $T_{A}=t$
implies that a change has happened at time $t$. Most of the change-point
detection framework will use pre-specified constant $A$, the decision to
detect a change point in distribution can be made as
$T_{A}=\inf_{t}\\{R_{t}\geq A\\}.$

Similar to the change point detection problem, we would like to minimize the
detection delay as $\sup_{1\leq\nu<\infty}\mathbb{E}(T_{A}-\nu|T_{A}\geq\nu)$.

Besides, there are typically two types of constraints: 1) under the case that
there is no change, the average detection time should be larger than a
threshold $B$ as $\mathbb{E}(T_{A}|\nu=\infty)\geq B$. 2) If a true change
mode $k$ happens, the false isolation rate, which is defined as:
$\max\limits_{1\leq k\leq K}\mathbb{P}^{k}\\{\hat{k}\neq k|T\geq\nu\\}$ should
also be small. However, these methods assume that the data is fully
observable, which cannot be used to partially observed data.

###  4.4 Proposed Algorithm

In Section 4.4, we will introduce the proposed methodology with the following
major steps, monitoring statistics update, change point detection decision,
failure mode isolation, and planning for adaptive sampling. The overall
framework is shown in Fig. 3, and the detailed steps are as follows:

  1. 1.

Monitoring statistics update: We first construct the monitoring statistics of
partially observed HD streaming data for each failure mode based on the SR
procedure. The detailed step is discussed in Section 4.4.1.

  2. 2.

Change point detection decision: According to the updated monitoring
statistics for each failure mode, a top-R statistic is used to conduct the
global decision. We will then raise a global alarm if the process has gone out
of control and decide which failure mode has occurred. The detailed step is
discussed in Section 4.4.2.

  3. 3.

Planning for adaptive sampling: If the change is not detected, we will update
the sampling layout dynamically according to the historical observations. To
achieve this, we propose to borrow the Thompson sampling idea to decide the
next sampling layout, where the data is randomly sampled from the identified
failure mode distribution. The detailed step is discussed in Section 4.4.3.
Furthermore, the optimization algorithm to solve this planning and optimal
sampling decision is discussed in Section 4.4.4. The selected sampling
patterns will be used to update the monitoring statistics recursively.

  4. 4.

Failure mode isolation: Finally, if the change is detected, we will isolate
and identify the true failure mode in the system.

Figure 3: Procedure of the proposed method

####  4.4.1 Recursive Monitoring Statistics Update

In this subsection, we will discuss the proposed method of constructing the
Shiyaev-Roberts (SR) statistics for each failure mode $k\in\\{1,\cdots,K\\}$
with missing observations. Here, we denote the local SR statistics at time $t$
as $R_{k,t}$. Here, $\mathbf{y}_{t}$ is the set of observed data streams. We
will follow the same rule of updating the local statistics $r_{k,t}$.

$R_{k,t}=(R_{k,t-1}+1)\frac{\tilde{f}_{C_{t},k}(\mathbf{y}_{t})}{\tilde{f}_{C_{t},0}(\mathbf{y}_{t})}.$

Here, $\tilde{f}_{C_{t},k}$ is the joint distribution of the observed data
$\mathbf{y}_{t}$ at set $C_{t}$. For computational efficiency and stability,
it is recommended to use the $r_{k,t}=\log R_{k,t}$, updated as

$r_{k,t}=\log(\exp(r_{k,t-1})+1)+\log\frac{\tilde{f}_{C_{t},k}(\mathbf{y}_{t})}{\tilde{f}_{C_{t},0}(\mathbf{y}_{t})}.$

We set $R_{k,0}=0$ initially and update the statistics accordingly.

####  4.4.2 Detection Decision and Failure Mode Isolation

We will combine the local statistics for each failure mode to construct a
global stopping time. Here, if we know that the system only contains one
failure mode, we propose to use the largest of the monitoring statistics
$r_{k,t}$ to trigger the alarm. For example, if $r_{k,t}$ is larger than a
control limit $A$, to raise the alarm.

$T=\inf\\{t\geq 1:\max\limits_{k}r_{k,t}\geq A\\}$ (3)

However, if we know that there are multiple failure modes in the system, the
summation of the top $r_{k,t}$ statistics can be used to trigger the alarm.

$T=\inf\\{t\geq 1:\sum_{k=1}^{K_{s}}r_{(k),t}\geq A\\}$ (4)

Finally, to isolate the most probable failure mode when the change is
detected, we propose to use the monitoring statistics with the largest index
$\hat{k}$, computed as

$\hat{k}=\arg\max\limits_{k}r_{k,T}.$ (5)

####  4.4.3 Planning for Adaptive Sampling

In this section, we present an efficient method to plan and select the best
sampling pattern to observe at the next time point. Suppose that we have
observed $\mathbf{y}_{1},\cdots,\mathbf{y}_{t-1}$ now and the goal is to
determine $\\{a_{j,t}\\}$ at next time $t$, which is a binary variable
denoting whether variable $j$ is observed or not at time $t$. Inspired by the
MAB, we propose to maximize the reward function, defined by the monitoring
statistics of the top few selected failure modes. More specifically, we
propose to use the summation of the SR statistics of the top-$K_{s}$ failure
modes as the reward function, where $K_{s}$ is a pre-defined parameter to
balance the exploration and exploitation. In other words, we can define the
reward function as $S_{t}=\sum_{k=1}^{K_{s}}r_{(k),t}$, where $r_{(k),t}$ is
the rank of the statistics such as $r_{(1),t}\geq\cdots\geq r_{(K_{s}),t}$.
One specific challenge is that to compute the reward function $S_{t}$ for
planning, we need to compute $r_{(k),t}$, which requires $\mathbf{x}_{t}$ to
be fully observed. However, given that we are still at time $t-1$ yet and data
$x_{t}$ has not been observed yet, it is impossible to compute and optimize
$S_{t}$ for the planning problem.

To solve this planning problem, we propose to optimize a sampled version of
the monitoring statistics $S_{t}$, defined as
$\tilde{S}_{t}=\sum_{k=1}^{K_{s}}a_{k,t}r_{(k),t}$. Borrowing from the
Thompson sampling algorithm, we would like to use the sampled version of
$\mathbf{x}_{t}$, denoted as $\mathbf{\tilde{x}}^{k}_{t}$ as

$\tilde{r}_{k,t}=\log(\exp{(r_{k,t-1})}+1)+\log\frac{f_{k}(\tilde{\mathbf{x}}^{k}_{t})}{f_{0}(\tilde{\mathbf{x}}^{k}_{t})},$
(6)

where $\tilde{\mathbf{x}}^{k}_{t}\sim f_{k}$ is sampled from the $k^{th}$
failure mode. Finally, one can optimize $a_{j,t}$ by maximizing the sampled
version of $S_{t}$, denoted as $\tilde{S}_{t}$ as

$\displaystyle\max_{a_{j,t}}\tilde{S}_{t}\quad\text{subject to}$
$\displaystyle\sum_{j}a_{j,t}=q,a_{j,t}=\\{0,1\\}$ (7)

Finally, we will discuss how to solve the optimization (7) in Section 4.4.4.

####  4.4.4 Optimization for Planning

In general, optimizing (7) is often challenging. If the problem dimension $p$
and the number of selected sensing variables $q$ are small, we can enumerate
all the $\tbinom{p}{q}$ possible combinations of the sampling layouts.
However, given the time complexity is $O\left(\tbinom{p}{q}\right)$,
enumeration of all possible combinations is not feasible for large $p,q$.

Here, we will first present the closed-form solution for a special case of the
proposed algorithm, where the joint failure mode distribution
$f_{k}(\mathbf{x}_{t})=\prod_{j}f_{j,k}(x_{t,j})$ can be approximated by
independent but not necessarily identical distributions in each dimension $j$.
Notice that if the possible failure modes in each dimension $j$ is finite, the
total number of possible failure modes for all dimensions is finite. Here, we
would like to derive the analytical solution to optimize (7) under this
setting in Proposition 1.

######  Proposition 1.

If $f_{k}(\mathbf{x}_{t})=\prod_{j}f_{j,k}(x_{t,j})$ for $k=0,\cdots,K$, the
set $C_{t}$ in (7) can be solved by selecting the indices of the largest $q$
of $s_{t,j}$, denoted as $s_{t,(1)},s_{t,(2)},\cdots,s_{t,(q)}$. Here
$s_{t,j}$ is defined as

$s_{t,j}=\sum_{k=1}^{K_{s}}\log\frac{f_{j,(k)}(\tilde{x}^{k}_{j,t})}{f_{j,0}(\tilde{x}^{k}_{j,t})}$
(8)

$s_{t,(j)}$ is the order statistics, defined as $s_{t,(1)}\geq
s_{t,(2)}\geq\cdots\geq s_{t,(q)}\geq s_{t,(q+1)}\geq\cdots\geq s_{t,(p)}$.

We would like to mention that the computation of (8) in Proposition 1 is
actually very efficient. To compute each $s_{t,j}$, it requires the summation
of $K_{s}$ terms, which is of $O(K_{s})$ complexity. To compute all $p$
sensing variables at each time $t$, it requires only $O(pK_{s})$ complexity at
each time to decide the best sampling layout. Here, the limitation is that
$f_{k,j}$ is assumed to be independent over different data dimension $j$.
However, we find that even the distribution of each failure mode is not
independent, this approximation can still achieve a pretty reasonable
solution.

In this paper, we will only focus on the monitoring of continuous variables
and assume that the data follows a normal distribution $f_{k}\sim
N(\boldsymbol{\mu}_{k},\Sigma_{k})$. However, as derived in the proposed
framework, this method can be generalized to other distributions quite easily.
Finally, as mentioned in Proposition 1, if we will further assume that
$\Sigma_{k}$ is diagonal as
$\Sigma_{k}=diag(\sigma_{k,1}^{2},\cdots\sigma_{k,p}^{2})$, we can derive a
simpler formula for $s_{t,j}$ in Proposition 2.

######  Proposition 2.

Given that $f_{k,j}=N(\mu_{k,j},\Sigma_{k})$, where
$\Sigma_{k}=diag(\sigma_{k,1}^{2},\cdots\sigma_{k,p}^{2})$. We can derive
$s_{t,j}=\sum_{k=1}^{K_{s}}(\frac{1}{\sigma_{k,j}^{2}}(\tilde{x}^{k}_{j,t}-\mu_{k,j})^{2}-\frac{1}{\sigma_{k,0}^{2}}(\tilde{x}^{k}_{j,t}-\mu_{k,j})^{2})$,
$\tilde{x}_{t}\sim f_{k}$,

We would like to emphasize that the independence assumption of each failure
mode distribution is actually not required for the proposed algorithm. It is
only useful to derive the closed-form solution in solving (7). If the spatial
dimension is not independent for different failure modes, the proposed
planning procedure can still be optimized without the assumptions by
approximating the optimal solution. Here, we propose a greedy algorithm to
detect the $a_{j,t}$ sequentially. The detailed step is given as follows.
First, we can select the first sensing variable $j_{1}$ to optimize the
sampled version $\tilde{S}_{t}$ by
$j_{1}=\arg\max\tilde{S}_{t},\sum_{j}a_{j,t}=1,a_{j,t}=\\{0,1\\}.$ After
$j_{1}$ is decided, we would like to choose the second sensing variable
$j_{2}$ by $j_{2}=arg\max\tilde{S}_{t},\sum
a_{j,t}=2,a_{j_{1},t}=1,a_{j,t}=\\{0,1\\}$. We will continue the procedure
until $j_{q}$ is selected. In conclusion, the set of the observed sensing
index is given as $C_{t}=\\{j_{1},\cdots,j_{q}\\}$. Given that we only need to
enumerate all $p$ dimensions in each of the $q$ iterations, the time
complexity can be reduced to $O(pq)$. Despite the efficiency, the greedy
forward selection strategy usually does not produce a global optimal solution.

Here, we would like to highlight the major difference of the proposed method
compared to the existing literature of monitoring of the i.i.d data stream
such as (Zhang and Mei, 2020): 1) the number of failure modes $K$ does not
need to be the same as the dimensionality $p$ of the data stream; 2) The
normal data distribution $f_{0}$ does not need to be i.i.d according to each
dimension as $x_{j}\overset{i.i.d}{\sim}f(x),\text{for all }j$; 3) For each
failure mode, it can include overlapping sets of sensing variables.

Finally, we would like to point out a special version of the proposed method
and how it links to (Zhang and Mei, 2020), if we are interested in monitoring
the i.i.d data stream with the focus of detecting the change of each
individual sensing variable.

######  Proposition 3.

For before change $H_{0}:x_{j}\overset{i.i.d}{\sim}f(x),\text{for all }j$.
After change, for $j^{th}$ failure mode, where $j\in\\{1,\cdots,p\\}$, only 1
distribution changed the distribution to $g(x)$ as $x_{j}{\sim}g(x)$, where
the rest $x_{j^{\prime}}\overset{i.i.d}{\sim}f(x),j^{\prime}\neq j$ still
follows the pre-change distribution. The proposed algorithm will result in the
sampled updating rule as in (Zhang and Mei, 2020):

$R_{j,t}=\begin{cases}\frac{g(\tilde{x}_{j,t})}{f(\tilde{x}_{j,t})}(R_{j,t-1}+1)&a_{j,t}=1\\\
R_{j,t-1}+1&a_{j,t}=0\end{cases}$ (9)

Proposition 3 shows a special case for the proposed algorithm, which assumes
the change only affected a few data streams and the algorithms try to identify
the change with the resources constraint. Under the current setting, the
proposed algorithm will become another sampled version of the TSSRP algorithm.
Many previous works, including (Liu et al., 2015; Zhang and Mei, 2020) have
studied this setting. However, the proposed algorithm can be generalized into
any other joint distributions of different failure modes.

###  4.5 Properties of the Proposed Algorithm

Here, we will prove two important properties of the algorithms about the bound
of the average run length and the failure mode isolation in Theorem 4 and
Theorem 5, respectively.

######  Theorem 4 (Average Run Length).

Let $T=\inf\\{t\geq 1:r_{(1),t}\geq A_{2}\\}$. Then we have that under the
null hypothesis where no changes occur, $\mathbb{E}T\geq A_{1}/K$,
$\mathbb{E}T=O(A_{1})$, where $A_{1}=e^{A_{2}}$.

Theorem 4 provides a lower and upper bound for the Average Run Length if no
changes occur. Theorem 4 provides us the guidance to select conservative upper
and lower bounds of the control limit $A$. Specifically, $K*ARL$ can serve as
the upper bound in the bisection search to speed up the threshold choosing
procedure.

######  Theorem 5 (Failure Mode Isolation).

Assume ${\bf X}_{1}\ldots{\bf X}_{\nu}\sim f_{0}$, ${\bf X}_{\nu+1}\ldots\sim
f_{k}$. All distributions of failure modes are continuous, and the KL
divergence of the distributions of failure mode $l$ and the true failure mode
$k$ follows $0<KL(f_{k}\|f_{l})<\infty$, and $\mathrm{Var}_{x\sim
f_{k}}[\log\frac{f_{k}}{f_{l}}]<\infty$, $\text{for all }l\neq k$. The
probability that $P(r_{k,t}>r_{l,t})\rightarrow 1$ as $t\to\infty$.

Theorem 5 provides the behavior of the largest SRP statistics when time goes
to infinity under the alternative hypothesis (where the failure mode $k$
occurs). It shows that the adaptive sampling algorithm will always be able to
isolate the true failure mode $k$ if $t\rightarrow\infty$.

In some cases, there might be multiple failure modes happening at different
time points after the change point $\nu$, i.e. for some $t_{1}>\nu$, ${\bf
X}_{t_{1}}\sim f_{k}$; for some $t_{2}>\nu$, ${\bf X}_{t_{2}}\sim f_{l}$.
Corollary 6 is proved.

######  Corollary 6.

Assume ${\bf X}_{1}\ldots{\bf X}_{\nu}\sim f_{0}$. Let $\mathcal{K}=\\{k:{\bf
X}_{t}\sim f_{k}\textit{ for some }t>\nu\\}$ be the set of true failure modes.
If we further assume that the support of different failure modes are non-
overlapping, then for any $k\in\mathcal{K}$ and $l\notin\mathcal{K}$, we have
$\lim_{t\to\infty}P(r_{k,t}>r_{l,t})=1$.

Corollary 6 assumes that different failure modes are not overlapping with each
other, we can prove in Corollary 6 that SRP statistics of the true failure
modes will be larger compared to those of the other potential failure modes.
We would like to point out that Corollary 6 is not always true for failure
mode distributions that are potentially overlapped with each other. For
example, if half of the data after the change follows $f_{1}$ and the other
half after the change follows $f_{2}$, it might be possible that the true
failure mode identified would be $f_{3}=\frac{1}{2}(f_{1}+f_{2})$.

###  4.6 Choice of Parameters

Here, we will present practical guidelines for tuning parameter selection.
Given that the number of sensing variables $q$ typically depends on the
available sensing resources in the particular applications, we only need to
select the following parameters: the number of top-R selected failure modes
$K_{s}$, the control limit threshold $A$, and the failure mode distribution
$f_{k}$ and $f_{0}$.

Choice of the number of observed failure modes $K_{s}$: First, the number of
selected failure modes for the monitoring statistics should be smaller than
the total number of potential failure modes. Ideally, $K_{s}$ should be chosen
as large as the total number of true failure modes in the system. In practice,
we found that increasing $K_{s}$ to be more than the true number of failure
modes in the system would lead the algorithm to explore more potential failure
modes or increase the exploration power. However, if $K_{s}$ is too large, the
algorithm is not able to focus on the actual failure modes, which decreases
the exploitation power.

Choice of threshold $A$: The choice of control limit $A$ can be determined by
the In-control ARL (or $ARL_{0}$). If $A$ is large, $ARL_{0}$ will also
increase. In practice, we can set an upper bound of the $A$ by utilizing the
Theorem 4 and then use the binary search algorithm to find the best $A$ for a
fixed $ARL_{0}$.

###  4.7 Selection of failure mode distribution $f_{k}$

Finally, the complex joint distributions of different failure modes also bring
significant computational complexity, which will be addressed in this paper.

Selecting the failure mode distribution is very important to achieve better
change detection and isolation performance. However, it is very challenging to
provide accurate failure modes definition for high-dimension data without any
domain knowledge. In this work, we mainly focus on detecting mean-shift in
high dimension data. We further discuss how to define the failure modes based
on our knowledge of the high-dimension data. Overall, there are two strategies
for choosing the most appropriate failure mode distribution $f_{k}$. 1) If
there is prior knowledge about the failure mode distributions, we can set the
distribution according to the prior knowledge. For example, if we know that
the hot spots are clustered, each failure mode distribution can be assumed as
the mean-shift of the IC distribution $f_{0}$ with an individual B-spline or
Gaussian kernel basis. If we know that the post-change distribution is sparse,
a simple way is to set the failure mode distribution as the mean shift of each
individual sensor. 2) If we do not know the failure mode distributions, we can
collect some samples for each failure mode and use these samples to estimate
the failure mode distribution.

##  5 Simulation Study

Here, we will evaluate the proposed method in a simulation study. We will
start with the simulation setup for single failure mode in Section 5.1 with
two different scenarios: the non-overlapping case and the overlapping case.
Then, we will evaluate the proposed algorithm in these two scenarios. To test
the robustness of the proposed algorithm in the case when multiple failure
modes coexist, we also perform the sensitivity analysis to evaluate the
performance in 5.2

###  5.1 Simulation Study for a Single Failure Mode

Here, we will discuss the two scenarios for the simulation setup. We are
trying to distinguish multiple failure modes by whether these failure modes
have overlapping support. For example, for the first ”non-overlapping” case,
we assume that different failure modes $f_{i}$ and $f_{j}$ have non-
overlapping support.

####  5.1.1 Scenario 1: the non-overlapping case

We will first discuss the non-overlapping case of the proposed method. In the
simulation, we let the data dimension $p=1000$, the number of failure mode
$K=50$, each time the algorithm will select $q=10$ sensors at each time. Here,
we assume the normal data or in control (IC) data follows $f_{0}=N(0,I)$.
After the change, the out-of-control (OC) data have $K$ failure modes, where
$f_{k}=N(\boldsymbol{\mu}_{k},I)$ and
$\boldsymbol{\mu}_{k}=\sum_{j=sk}^{(s+1)k}e_{j}$.
$e_{j}=(0,\cdots,1,\cdots,0)$, and only $j^{t}h$ element is 1. Here, three
failure modes have been selected after the change. In other words, if we
organize the 1000 data streams into a $50\times 20$ image, each failure
pattern would be each row of pixels which can be visualized in Fig. 7(a).
Therefore, different failure modes are not overlapped, given they contain
different sensing variables.

(a) All 50 Non-overlapping potential failure modes

(b) All 49 overlapping potential failure modes

Figure 4: Failure Modes for Overlapping and Non-overlapping Cases

Finally, we assume that we are only accessible to $10$ out of $p=1000$ data
streams to observe. At each time step, we adaptive select data streams to
monitor the whole process. We want to detect the correct failure mode as soon
as possible.

####  5.1.2 Scenario 2: the overlapping case

We will discuss the second scenario, where the failure patterns are generated
as small spatial clusters. In this case, we set the data as 2-D images with
size $30\times 30$ with total dimension $p=900$. Here, the failure modes are
generated using B-spline basis with $7$ knots in both $x$ and $y$ directions.
As shown in Fig. 7(b), we end up with $7^{2}=49$ potential failure modes.
After the change happens, we randomly select a failure mode as the true
failure mode. In this situation, some failure modes might overlap with each
other, which will be more challenging for the algorithm to isolate the real
changes. Finally, during the monitoring process, we can select 10 out of
$p=900$ data streams adaptively to observe online.

####  5.1.3 Simulation Result

Here, we will compare the proposed MTSSRP with the following benchmark
methods: 1) TSSRP method (Zhang and Mei, 2020), which is introduced in detail
in Appendix. 2) TRAS method Liu et al. (2015), where the local CUSUM
statistics are used for each individual data stream and later fused together
via the Top-R rule. To show the upper-bound and lower-bound performance, we
will also add three simple alternatives: 1) Random, where we randomly select
$q$ sensors at each time step with the top-r statistics by monitoring each
sensor individually. 2) Oracle, where we have not only access to all the data
streams but also the failure mode distribution information using the
monitoring statistics as MTSSRP. 3) MRandom, where we apply the same
monitoring statistics as MTSSRP, which considers the failure mode distribution
information in the monitoring statistics, but we randomly select the sensors
at each time step. We evaluate the proposed method with two metrics which are
detection delay and failure isolation accuracy.

First, we would like to compare the detection delay of the proposed method and
all the benchmark methods. Here, we set the in-control average run length
(i.e., denoted as $ARL_{0}$) for all methods as $200$ and compare their out-
of-control ARL or average detection delay (i.e., denoted as $ARL_{1}$) with
1000 replications. Here, we will compute the $ARL_{1}$ for different change
magnitudes ($\delta$ = $0.5$, $0.8$) in Table 1. We also compared the
$ARL_{1}$ and Isolation Accuracy from different magnitude $\delta$ ($0.1$ to
$0.8$) in Fig. 5. From the results, we can see that the proposed MTSSRP has
better $ARL_{1}$ compared to other benchmark methods. MTSSRP performs much
better than MRandom, which validates the efficiency of the proposed sampling
strategy. The advantage of MTSSRP over TSSRP shows that considering the
failure mode information can greatly improve the performance. We further
compare the isolation accuracy as shown in Fig. 8(c) and Fig. 8(d). It can be
seen that the proposed MTSSRP achieves better performance than others. It can
also reach pretty high accuracy when $\delta$ is greater than 0.6.

Table 1: Average Run Length and Failure Mode Isolation Accuracy for single
failure

_Case_ nonoverlap overlap _Change Magnitude_ $\delta=0.5$ $\delta=0.8$
$\delta=0.5$ $\delta=0.8$ _Metrics_ $ARL_{1}$ Accuracy $ARL_{1}$ Accuracy
$ARL_{1}$ Accuracy $ARL_{1}$ Accuracy _Oracle_ 13.71(11.56) 0.99(0.04)
2.52(1.08) 1.0(0.0) 13.68(12.05) 1.0(0.05) 2.57(1.43) 1.0(0.03) _Competing
Methods_ MTSSRP 112.05(66.96) 0.75(0.43) 26.49(17.09) 0.99(0.03) 111.16(66.01)
0.75(0.43) 26.66(17.71) 0.98(0.14) TSSRP 124.66(58.2) 0.55(0.50) 61.19(30.67)
0.81(0.39) 134.51(58.38) 0.48(0.5) 60.5(35.92) 0.81(0.39) TRAS 164.84(49.07)
0.46(0.50) 75.53(42.51) 0.98(0.13) 161.75(51.48) 0.47(0.5) 84.08(50.73)
0.93(0.25) MRandom 178.1(41.49) 0.28(0.45) 104.55(48.71) 0.92(0.28)
181.67(39.91) 0.23(0.42) 121.16(54.09) 0.8(0.4) Random 199.97(1.07) -
199.66(5.08) - 200.0(0.0) - 199.95(1.61) -

To understand how the proposed algorithm balances the exploration and
exploitation automatically, we would like to plot both the SR statistics for
each failure mode (i.e., in red) and when this particular failure mode has
observed sensors (i.e., in black dot) for both the potential failure mode
(i.e., failure mode doesn’t happen in this run) and the true failure mode in
Fig. 6 together. Here, the failure mode with the observed sensors can be
defined as that there are observed sensing variables located in the non-zero
location of the mean-shift of that particular failure mode. From Fig. 6, it is
clear that the statistics for potential is quite small compared to the
statistics for the true failure mode. From Fig. 6(a), we can also observe that
the SR statistics will grow naturally if this particular failure mode is not
observed. This will encourage the sensing variables to be allocated to this
particular failure mode eventually. Furthermore, if the particular failure
mode is observed where no change occurs, the monitoring statistics will drop
significantly, indicating that this failure mode has dropped significantly. On
the other hand, from Fig. 6(b), we can clearly see that after time $100$ when
the change occurs, the true failure mode statistics will grow significantly as
long as that particular failure mode is observed.

(a) $ARL_{1}$ in Non-overlapping Cases

(b) $ARL_{1}$ in Overlapping Cases

(c) Accuracy in Non-overlapping Cases

(d) Accuracy in Overlapping Cases

Figure 5: Out-of-control Average Run Length ($ARL_{1}$) and Failure isolation
accuracy for Different Change Magnitude $\delta$

(a) Potential Failure Mode

(b) True Potential Failure Mode

Figure 6: SR statistics for the True Failure Mode. The change happens at time
$t=100$. Left figure shows the monitoring statistics for the potential failure
mode, where the monitoring statistics is small. Right figure shows the
monitoring statistics for the true failure mode and the monitoring statistics
increase dramatically at time $t=100$.

###  5.2 Simulation Study for Multiple Failure Modes

We further design another simulation study to evaluate the performance of the
proposed methods when multiple failure modes happen together. Fig. 7
illustrates the potential failure modes and generated multiple failure modes.

(a) All 50 Non-overlapping potential failure modes

(b) All 49 overlapping potential failure modes

(c) Non-overlapping true failure modes

(d) Overlapping true failure modes

Figure 7: Generated Potential and True Failure Mode Patterns for Overlapping
and Non-overlapping Cases Table 2: Average Run Length and Failure Mode
Isolation Accuracy for multiple failure

_Case_ nonoverlap overlap _Change Magnitude_ $\delta=0.5$ $\delta=0.8$
$\delta=0.5$ $\delta=0.8$ _Metrics_ $ARL_{1}$ Accuracy $ARL_{1}$ Accuracy
$ARL_{1}$ Accuracy $ARL_{1}$ Accuracy _Oracle_ 6.05(0.12) 0.91(0.29)
1.78(0.02) 0.93(0.26) 5.28(0.12) 0.99(0.09) 1.56(0.02) 1.0(0.0) _Competing
Methods_ MTSSRP 57.26(1.32) 0.91(0.29) 14.27(0.22) 0.93(0.26) 53.73(1.36)
0.93(0.25) 14.38(0.23) 0.97(0.16) TSSRP 77.78(1.26) 0.85(0.36) 36.32(0.53)
0.88(0.26) 79.15(1.35) 0.89(0.31) 32.12(0.54) 0.92(0.26) TRAS 104.18(1.67)
0.43(0.50) 38.1(0.56) 0.91(0.29) 102.81(1.71) 0.48(0.50) 39.34(0.68)
0.93(0.25) MRandom 146.09(1.64) 0.59(0.49) 64.18(0.82) 0.92(0.26) 149.29(1.69)
0.57(0.50) 69.63(1.06) 0.98(0.15) Random 199.96(0.03) - 198.6(0.29) -
199.93(0.05) - 198.44(0.32) -

(a) $ARL_{1}$ for the Non-overlapping Cases

(b) $ARL_{1}$ for Overlapping Cases

(c) Accuracy for the Non-overlapping Cases

(d) Accuracy for Overlapping Cases

Figure 8: Out-of-control Average Run Length ($ARL_{1}$) and Failure isolation
accuracy for Different Change Magnitude $\delta$ for multiple failure modes

Similar to the experiment in the single failure mode. we will compare the
proposed MTSSRP with the following benchmark methods: 1) TSSRP method (Zhang
and Mei, 2020). 2) TRAS method Liu et al. (2015), where the local CUSUM
statistics are used for each individual data stream and later fused together
via the Top-R rule. To show the upper-bound and lower-bound performance, we
will also add three simple alternatives: 1) Random, where we randomly select
$q$ sensors at each time step with the top-r statistics by monitoring each
sensor individually. 2) Oracle, where we have not only access to all the data
streams but also the failure mode distribution information using the
monitoring statistics as MTSSRP. 3) MRandom, where we apply the same
monitoring statistics as MTSSRP, which considers the failure mode distribution
information in the monitoring statistics, but we randomly select the sensors
at each time step. We evaluate the proposed method with two metrics which are
detection delay and failure isolation accuracy. Here, we will compute the
$ARL_{1}$ for different change magnitudes ($\delta$ = $0.5$, $0.8$) in Table
2. We also compared the $ARL_{1}$ and Isolation Accuracy from different
magnitude $\delta$ ($0.1$ to $0.8$) in Fig. 8

We then show the sampling point distributions before the change and after the
change for both scenarios. To visualize the sampling distribution better, we
have aggregated 1000 generate sampled IC data and OC data in Fig. 9(a) and
Fig. 9(b), respectively for the non-overlapping failure modes and Fig. 9(c)
and Fig. 9(d) for the overlapping failure modes. It is clear that before the
change, the sampling distribution is pretty random. After the change, most of
the points will gather around the three true failure modes.

Finally, we showed the heatmap of the identified top failure modes from MTSSRP
method at $\delta=4$ in Fig. 10. It shows the top most likely failure modes
(columns) at a different time (rows) from our proposed algorithm. The actual
change point of the data is at $t=100$. Column $0$ shows that we identify a
failure pattern at around time $100$, which is consistent with the change
time. Prior to time $100$, the failure mode patterns are actually quite
random. Interestingly, we find that after the algorithm finds the first
failure pattern, the algorithm will continue to search for other underlying
failure patterns. At around time 120, we further detect the other two failure
patterns in our experiments in both overlapping and non-overlapping cases.
Another interesting behavior is that the other failure modes (except the top
3) are quite random in the non-overlapping case but not as random in the
overlapping case. The reason is that in the non-overlapping failure modes,
given that potential failure modes have no overlapping with the true failure
modes, there are no particular orders on the potential failure modes. However,
in the overlapping case, besides the 3 true failure modes, there are also some
potential failure modes overlapped with the true failure modes, which will be
selected. In conclusion, this behavior shows that the proposed algorithm is
able to detect the true failure modes, as indicated by the Theorem 5.

(a) IC Sampling Point Distribution for Non-overlapping Failure Modes

(b) OC Sampling Point Distribution for Non-overlapping Failure Modes

(c) IC Sampling Point Distribution for Overlapping Failure Modes

(d) OC Sampling Point Distribution for Overlapping Failure Modes

Figure 9: Examples of Simulated Data and Sampled Point Distribution for Both
the Overlapping and Non-overlapping Failure Modes; It is clear that the OC
sampling point distributions are the same as the true failure mode patterns.

(a) Non-overlapping case

(b) Overlapping case

Figure 10: Failure Mode Selection history of MTSSRP. Here, X-axis refers to
the ranking of the failure mode, and Y-axis refers to the time. All the
changes happen at time 100.

##  6 Case Study

In this section, we will evaluate the proposed MTSSRP algorithm in the laser
powder bed fusion process monitoring. We will evaluate the performance of the
MTSSRP and compare it with the state-of-the-art benchmark methods.

Figure 11: Hot-spots Detection in LPBF

###  6.1 Hot-spots Detection in Laser Powder Bed Fusion Process

Here, we will implement the proposed algorithm into the hot-spots detection in
the process monitoring in the Laser Powder Bed Fusion (LPBF) process. A $300$
fps video sequence was acquired during the realization of one layer of the
part by using the setup shown in Fig. 11, which consists of a thermal camera
mounting outside the LPBF chamber monitoring the hot-spots events. The
observed image is of size $121\times 71$ pixels. Previous studies showed that
the occurrence of local over-heating conditions might yield geometrical
distortions (Yan et al., 2020; Colosimo and Grasso, 2018). The hot-spots
caused by the formation of solidified balls will cause the local heat
accumulation and inflate from one layer to another. Therefore, the overall
goal of this study is to detect such hot-spots quickly. For more details about
the setup of this experiment and some preliminary works related to this
dataset, please refer to (Grasso et al., 2017; Colosimo and Grasso, 2018; Yan
et al., 2020). The dataset is also publicly available at
http://doi.org/10.6084/m9.figshare.7092863.

In this example, it is not easy to obtain the failure mode data beforehand,
and therefore, we rely on domain knowledge to define the failure mode
distribution. First, we know that the hot-spots must be in the scanning path.
Second, we know that the hot-spots must be locally clustered. Therefore, we
define the failure modes as each individual B-spline basis overlapped with the
printing regions. In this dataset, there are four different events starting
from 77, 94, 150, and 162. From Table 3, the proposed algorithm can detect the
change at time 79, 95, 152, and 162 with only 200 sensing variables out of
$8591$ sensing variables. In comparison, TSSRP can detect all four events but
with a much larger delay. However, TRAS can only detect Event 4, which fails
to detect the first three events. The original image frame, the sampling
patterns, and the selected failure modes at the detected time for these four
events are shown in Fig. 12. From Fig. 12, we can observe that the algorithm
can quickly converge the sampled points to the true hot-spots location at the
upper left corner.

Figure 12: Original Thermal Images, Sampling Patterns and Detected Failure Modes Table 3: Detection Time in the LPBF Process |  | Time of first signal
---|---|---
_Event Times_ | Event 1 | Event 2 | Event 3 | Event 4
_Actual Change Time_ | 77 | 94 | 150 | 162
_Competing Methods_ | MTSSRP | 79 | 95 | 152 | 162
| TSSRP | 80 | 99 | 156 | 164
| TRAS | - | - | - | 165

###  6.2 Tonnage Signal Monitoring

We will evaluate the proposed methodology to monitor the tonnage signals
collected in a multi-operation forging process, where four strain gauge
sensors on four columns of the forging machine measure the exerted tonnage
force of the press uprights as shown in Fig. 13. This results in the tonnage
profiles in each cycle of operation.

Figure 13: Tonnage Signal Monitoring

The data contains 305 in-control profiles collected under the normal
production condition and each abnormal production condition have 68 out-of-
control profiles for three different failure modes. 10 samples from each
failure mode is shown in Fig. 14. The four channel tonnage profiles results in
4804 dimensions in total.

In our experiment, since we do not have the prior knowledge about the failure
mode,. Therefore, we utilize 20% of the samples for training of the failure
modes by assuming that the data under each failure mode follows the Gaussian
distribution with diagonal covariance matrix approximation. We only select 500
out of 4804 data streams for the adaptive sampling. Here, we have conducted a
simulation study where the normal samples are randomly drawn from the $305$
in-control profiles. For out-of-control profiles, we have conducted three
scenarios, where each scenarios the out-of-control profiles are randomly
sampled from 68 out-of-control profiles with replacement. For all the methods,
we selected the threshold as the $95\%$ for the IC samples and compare the
out-of-control ARL (denoted as $ARL_{1}$) with $500$ replications for all this
three scenarios with the proposed MTSSRP and benchmark methods TSSRP and TRAS.

Figure 14: Tonnage Data

The result is shown in Table 4. From Table 4, we can conclude that failure
mode 3 is fairly easy to detect and all methods can achieve $ARL_{1}$ around
$1$. For failure mode 1, it is more similar to the normal dataset. The
proposed MTSSRP can achieve $ARL_{1}$ of $1.96$, which is at least half of the
$ARL_{1}$ of other methods such as TSSRP and TRAS. Finally, the most
challenging case is the failure mode 2, where the difference of failure mode 2
and the normal data is almost neglectable, as seen from Fig. 14. The proposed
MTSSRP can achieve $ARL_{1}$ around $7.46$, which is much smaller than the
$ARL_{1}$ of TSSRP and TRAS. In Fig. 15, the monitoring statistics for three
failure production conditions have been shown and our method is able to
identified the true failure mode correctly.

(a) Change of Failure mode 1

(b) Change of failure mode 2

(c) Change of failure mode 3

Figure 15: Monitoring Statistics for Each Failure Mode. The monitoring statistics of the true failure mode increase significantly after the changed time $t=305$.  Table 4: Out-of-control Average Run Length for Tonnage Signals | Out-of-control Average Run Length
---|---
_Failure Modes_ | Mode 1 | Mode 2 | Mode 3
_Competing Methods_ | MTSSRP | 1.96(0.05) | 7.46(0.26) | 1.02(0.01)
| TSSRP | 4.2(0.15) | 28.35(0.34) | 1.03(0.02)
| TRAS | 6.0(0.6) | 45.85(0.44) | 1.05(0.01)

###  6.3 COVID-19 Cases Detection for Hotspot Detetion

To better understand the COVID-19 status, different types of testing resource
is typically distributed into different regions. Centers for Disease Control
and Prevention (CDC) has classified the testing for COVID-19 into the
following two categories: 1) Diagnostic testing is intended to identify
current infection in individuals and is performed when a person has signs or
symptoms consistent with COVID-19, or is asymptomatic, but has recently known
or suspected exposure to COVID-19. 2) Screening tests are recommended for
unvaccinated people to identify those who are asymptomatic and do not have
known, suspected, or reported exposure to COVID-19. Screening helps to
identify unknown cases so that measures can be taken to prevent further
transmission. CDC (2022).

Overall, the screening test is very useful in randomly distributed test in
some underdeveloped areas to identify unknown cases so that measures can be
taken to prevent further transmission. However, in the Screening test, the
decision maker may have limited sampling resources. Therefore, adaptive
sampling techniques can be used to decide which region to sample at each time
epoch based on the collected testing results from the previous iterations.

In the case study, we will use the real COVID-19 test report data from Johns
Hopkins University Center for Systems Science and Engineering (JHU CSSE) (Dong
et al. (2020)). The dataset is available on
https://github.com/CSSEGISandData/COVID-19. More specifically, we will use the
confirmed COVID-19 cases from all 39 counties in Washington state. The source
of daily positive cases in Washington State is the Department of Health
(https://www.doh.wa.gov/Emergencies/COVID19). The time-series data is updated
daily around 23:59 (UTC). We will use the data from Jan 23, 2020, to Sep 13,
2020, a total of 235 days, as an illustration. On each day, the confirmed
cases are recorded in all counties in the United States. Overall, in the case
study, we assume at each time epoch, the state government will focus on doing
screening tests on Yakima County and Okanogan County over the $39$ counties.
In this example, it is important to detect the outbreak in each individual
county. Therefore, we set the failure mode as the mean shift of the infection
rate of each individual county. The outbreak time for Yakima and Okanogan are
118 and 169, as shown in Table. 5. The proposed method is able to detect the
outbreak rapidly. TSSRP method can also detect the outbreak with some delay,
while TRAS method cannot detect the outbreak. We further present the sampling
pattern shown in Fig. 17 and Fig. 18. As shown in Fig. 17, the algorithm
started to focus on certain counties starting around 140. Fig. 18 presents the
aggregated sampling frequency for each county. It can be seen that the
algorithm samples all the counties evenly during the in-control phase, as
shown in Fig. 18(a). During the out-of-control phase, as shown in Fig. 18(b),
it started to focus on the true hotspot counties such as Yakima and Okanogan,
which helps us identify the outbreak quickly in those counties.

(a) Monitoring Statistics for Hotspot Detection

(b) OMonitoring Statistics for different counties

Figure 16: Monitoring Statistics for Covid-19 Case Figure 17: Sampling Pattern

(a) Sampling frequency during in control phase

(b) Sampling frequency during out control phase

Figure 18: Sampling Pattern for Covid-19 Case Table 5: Detection Time for Covid hotspot detection | Time of first signal
---|---
Location |  | Yakima County
---
| Okanogan County
---
Infection Rate $>$ 0.01 | 118 | 169
Competing Methods | MTSSRP | 138 | 170
TSSRP | 173 | 183
TRAS | - | -

##  7 Conclusion

Online change detection of high-dimensional data under multiple failure modes
is an important problem in reality. In this paper, we propose to borrow the
concept from Bayesian change point detection and MAB to adaptively sample
useful local components, given the distributions of the multiple failure
modes. Our proposed algorithm can balance between exploration of all possible
system failure modes or exploitation of the most probable system failure mode.
Furthermore, we also studied the properties of the proposed methods and showed
the proposed algorithm could isolate the correct failure mode. Our simulation
and case study show that the proposed algorithm, by considering the failure
mode information, can significantly reduce the detection delay.

##  Appendix A Proof of Proposition 1

###### Proof.

It is worth noting that under the assumption that distribution is independent,
$\tilde{r}_{k,t}$ can be derived as

$\displaystyle\tilde{S}_{t}$
$\displaystyle=\sum_{k=1}^{K_{s}}\tilde{r}_{(k),t}$
$\displaystyle=\sum_{k=1}^{K_{s}}(\log(\exp(r_{(k),t-1})+1)+\sum_{j\in
C_{t}}\log\frac{f_{j,(k)}(\tilde{x}^{k}_{j,t})}{f_{j,0}(\tilde{x}^{k}_{j,t})})$
$\displaystyle=\sum_{k=1}^{K_{s}}(\log(\exp(r_{(k),t-1})+1))+\sum_{j\in
C_{t}}\sum_{k=1}^{K_{s}}\log\frac{f_{j,(k)}(\tilde{x}^{k}_{j,t})}{f_{j,0}(\tilde{x}^{k}_{j,t})})$
$\displaystyle=C_{0}+\sum_{j\in C_{t}}s_{j}$

Therefore, $\tilde{S}_{t}$ can be derived as
$\tilde{S}_{t}=\sum_{k=1}^{K_{s}}\tilde{r}_{(k),t}=C_{0}+\sum_{j\in
C_{t}}s_{j}$. Here,
$s_{j}=\sum_{k=1}^{K_{s}}\log\frac{f_{j,(k)}(\tilde{x}^{k}_{j,t})}{f_{j,0}(\tilde{x}^{k}_{j,t})}$
and $C_{0}=\sum_{k=1}^{K_{s}}(\log(\exp(r_{k,t-1})+1))$ is a con-stance. To
optimize $\tilde{S}_{t+1}$, we can always select the largest $q$ from $s_{j}$
after ranking $s_{(1)}\geq s_{(2)}\geq\cdots\geq s_{(q)}\geq
s_{(q+1)}\geq\cdots\geq s_{(p)}$.

∎

##  Appendix B Proof of Proposition 2

Therefore,
$\log\frac{f_{k,j}(X_{j,t})}{f_{0,j}(X_{j,t})}=(\frac{1}{\sigma_{k,j}^{2}}(\tilde{x}^{k}_{j,t}-\mu_{k,j})^{2}-\frac{1}{\sigma_{k,0}^{2}}(\tilde{x}^{k}_{j,t}-\mu_{k,j})^{2})$.
Therefore, we know that

$s_{j}=\sum_{k=1}^{K_{s}}\log\frac{f_{j,(k)}(\tilde{x}^{k}_{j,t})}{f_{j,0}(\tilde{x}^{k}_{j,t})}=\sum_{k=1}^{K_{s}}(\frac{1}{\sigma_{k,j}^{2}}(\tilde{x}^{k}_{j,t}-\mu_{k,j})^{2}-\frac{1}{\sigma_{k,0}^{2}}(\tilde{x}^{k}_{j,t}-\mu_{k,j})^{2}).$

##  Appendix C Proof of Proposition 3

###### Proof.

Here, If sensor $k$ is not selected,
$\frac{{f}_{k}(y_{t})}{{f}_{0}(y_{t})}=1$. Therefore, $R_{k,t}=R_{k,t-1}+1$.
This implies that if sensor $k$ is not observed, the corresponding $R_{k,t}$
will increase by $1$. Furthermore, if sensor $k$ is selected, we will update
$R_{k,t}$ as
$R_{k,t}=\frac{g(\tilde{x}_{k,t})}{f(\tilde{x}_{k,t})}(R_{k,t-1}+1)$. ∎

##  Appendix D Proof of Theorem 4

###### Proof.

Let $T_{1}=\inf\\{t\geq 1:R_{(1),t}\geq A_{1}\\}$, $T_{2}=\inf\\{t\geq
1:r_{(1),t}\geq A_{2}\\}$. Since the event $\\{R_{(1),t}\geq A_{1}\\}$ is
equivalent to the event $\\{r_{(1),t}\geq\log A_{1}\\}$, we only need to
consider the behavior for $T_{1}$.

For an upper bound of $T_{1}$, we define $T^{\prime}=\inf\\{t:R_{k,t}\geq
A_{1}\\}$ for some $k$. The theorem from Pollak (1987) is valid in
multivariate cases and thus
$\mathbb{E}T_{1}\leq\mathbb{E}T^{\prime}=O(A_{1})$. We show that
$\sum_{k=1}^{K}R_{k,t}-Kt$ is a martingale under the null hypothesis.

$\displaystyle\mathbb{E}[\sum_{k=1}^{K}R_{k,t+1}-K(t+1)|\mathcal{F}_{t}]$
$\displaystyle=\mathbb{E}\left[\sum_{k=1}^{K}(R_{k,t}+1)\frac{\tilde{f}_{C_{t+1},k}({\bf
X}_{C_{t+1},t+1})}{\tilde{f}_{C_{t+1},0}({\bf
X}_{C_{t+1},t+1})}-K(t+1)|\mathcal{F}_{t}\right]$
$\displaystyle=\sum_{k=1}^{K}(R_{k,t}+1)\mathbb{E}\left[\frac{\tilde{f}_{C_{t+1},k}({\bf
X}_{C_{t+1},t+1})}{\tilde{f}_{C_{t+1},0}({\bf
X}_{C_{t+1},t+1})}|\mathcal{F}_{t}\right]-K(t+1)$

For any $k$, we have that

$\displaystyle E\left[\frac{\tilde{f}_{C_{t+1},k}({\bf
X}_{C_{t+1},t+1})}{\tilde{f}_{C_{t+1},0}({\bf
X}_{C_{t+1},t+1})}|\mathcal{F}_{t}\right]$
$\displaystyle=\int_{x_{1},\ldots,x_{n}}\frac{\tilde{f}_{C_{t+1},k}({\bf
X}_{C_{t+1},t+1})}{\tilde{f}_{C_{t+1},0}({\bf
X}_{C_{t+1},t+1})}f_{0}(x_{1},\ldots,x_{n})$
$\displaystyle=\int_{x_{1},\ldots,x_{n}}\tilde{f}_{C_{t+1},k}({\bf
X}_{C_{t+1},t+1})f_{C_{t+1},0}(x_{unobserved}|{\bf X}_{C_{t+1},t+1})$
$\displaystyle=\int_{x_{observed}}\tilde{f}_{C_{t+1},k}({\bf
X}_{C_{t+1},t+1})\int_{x_{unobserved}}f_{C_{t+1},0}(x_{unobserved}|{\bf
X}_{C_{t+1},t+1})$ $\displaystyle=1$

Thus,
$\mathbb{E}[\sum_{k=1}^{K}R_{k,t+1}-K(t+1)|\mathcal{F}_{t}]=\sum_{k=1}^{K}R_{k,t}-Kt$
and $\sum_{k=1}^{K}R_{k,t}-Kt$ is a martingale. Since we have that
$\liminf_{t\to\infty}\int_{T_{1}>t}|\sum_{k=1}^{K}R_{k,t}-Kt|{\rm
d}\mathbb{P}_{\infty}=0$, $\sum_{k=1}^{K}R_{k,t}-Kt$ is uniformly integrable.
Thus
$\mathbb{E}_{\infty}(\sum_{k=1}^{K}R_{k,T_{1}}-KT_{1})=\mathbb{E}_{\infty}(\sum_{k=1}^{K}R_{k,0})=0$.
Thus, $\mathbb{E}KT_{1}=\mathbb{E}_{\infty}\sum_{k=1}^{K}R_{k,T_{1}}\geq
A_{1}$.
This completes the part for $T_{1}$. When $A_{2}=\log A_{1}$, the events
$\\{R_{(1),t}\geq A_{1}\\}$ and $\\{r_{(1),t}\geq\log A_{1}\\}$ are
equivalent. Therefore, the behavior of $T_{2}$ follows the same rule.

∎

##  Appendix E Proof of Theorem 5

###### Proof.

From (Pollak, 1987), we have

$\log R_{k,\nu+r}=\sum_{t=\nu+1}^{\nu+r}Z_{k,t}+\log
R_{k,\nu}+\sum_{i=0}^{r-1}\log[1+\frac{1}{R_{k,\nu+i}}]$

The increment of $R_{k,t}$ will be almost $Z_{k,t}$ where
$Z_{k,t}=\log(\frac{\tilde{f}_{C_{t},k}({\bf
X}_{C_{t},t})}{\tilde{f}_{C_{t},0}({\bf X}_{C_{t},t})})$. For any $l\neq k$,

$\lim_{t\to\infty}\mathbb{P}(R_{l,t}\leq
R_{k,t})=\lim_{T\to\infty}\mathbb{P}(\sum_{t=\nu+1}^{T}Z_{l,t}-Z_{k,t}\leq
b),$

where $b$ is some constant.

$\displaystyle\mathbb{P}(\sum_{t=\nu+1}^{T}Z_{l,t}-Z_{k,t}\leq b)$
$\displaystyle=\mathbb{P}\left(\sum_{t=\nu+1}^{T}\log(\frac{\tilde{f}_{C_{t},l}({\bf
X}_{C_{t},t})}{\tilde{f}_{C_{t},k}({\bf X}_{C_{t},t})})\leq b\right)$

${\bf X}_{C_{t}}$ is a sample of $q$ variables from $p$ dimensional data
$\mathbf{{\bf X}}_{t}$. There are $\binom{p}{q}$ kinds of different samples.
We divide the variables $\log(\frac{\tilde{f}_{C_{t},l}({\bf
X}_{C_{t},t})}{\tilde{f}_{C_{t},k}({\bf X}_{C_{t},t})})$ by the categories of
different sample results, and denote them by
$Y_{j,t}=\log(\frac{\tilde{f}_{C_{t},l}({\bf
X}_{C_{t},t})}{\tilde{f}_{C_{t},k}({\bf X}_{C_{t},t})})$,
$j=1,2,\ldots,\binom{p}{q}$ and different $j$ represents different selection
of $C_{t}$. We notice that when $j$ or $t$ is different, it has to be sampled
from different times. Thus, they are independent. Since we consider the
situation when $T$ goes to infinity, there exists at least one selection that
is observed infinitely many times. We first consider the case that there is
only one selection observed infinitely many times. Without loss of generality,
let this selection be the first one. Let the times that other selections of
data is observed be $c_{2},\ldots,c_{\binom{p}{q}}$ respectively. Denote
$Y_{j,t}=\log(\frac{f_{l,j}(X_{t,j})}{f_{k,j}(X_{t,j})})$. When $j=1$,
$t=1,2,\ldots$; when $j=2,\ldots,\binom{p}{q}$,
$t=1,2,\ldots,c_{\binom{p}{q}}$. It suffices to show
$\mathbb{P}\left(\sum_{t=1}^{\infty}Y_{1,t}+\sum_{t=1}^{c_{2}}Y_{2,t}+\cdots+\sum_{t=1}^{c_{\binom{p}{q}}}Y_{\binom{p}{q},t}\leq
b\right)=1$ for any constant $b$. Since $Y_{1,t}$ is i.i.d with negative
expectation, we have that for any constant $b$,
$\mathbb{P}\left(\sum_{t=1}^{\infty}Y_{1,t}\leq b\right)=1$
For any $\varepsilon>0$, let $d_{j}$ be the constant s.t.
$\mathbb{P}\left(\sum_{t=1}^{c_{j}}Y_{j,t}<d_{j}\right)>1-\frac{\varepsilon}{\binom{p}{q}-1}$.

$\displaystyle\mathbb{P}\left(\sum_{t=1}^{\infty}Y_{1,t}+\sum_{t=1}^{c_{2}}Y_{2,t}+\cdots+\sum_{t=1}^{c_{\binom{p}{q}}}Y_{\binom{p}{q},t}\leq
b\right)$
$\displaystyle>\mathbb{P}(\sum_{t=1}^{c_{2}}Y_{2,t}<d_{2})\mathbb{P}(\sum_{t=1}^{c_{3}}Y_{3,t}<d_{3})\cdots\mathbb{P}(\sum_{t=1}^{c_{\binom{p}{q}}}Y_{{\binom{p}{q}},t}<d_{\binom{p}{q}})\mathbb{P}\left(\sum_{t=1}^{\infty}Y_{1,t}\leq
b-d_{2}-\ldots-d_{\binom{p}{q}}\right)$
$\displaystyle>(1-\frac{\varepsilon}{\binom{p}{q}-1})^{\binom{p}{q}-1}\mathbb{P}\left(\sum_{t=1}^{\infty}Y_{1,t}\leq
b-d_{2}-\ldots-d_{\binom{p}{q}}\right)>1-\varepsilon$

Notice that this is true for any $\varepsilon>0$. We have

$\mathbb{P}\left(\sum_{t=1}^{\infty}Y_{1,t}+\sum_{t=1}^{c_{2}}Y_{2,t}+\cdots+\sum_{t=1}^{c_{\binom{p}{q}}}Y_{{\binom{p}{q}},t}\leq
b\right)=1.$

Using similar method, we can also prove the case when there are more than one
dimension observed infinitely many times. Thus,
$\lim_{t\to\infty}\mathbb{P}(R_{l,t}\leq R_{k,t})=1$. Since $l$ is chosen
arbitrarily, we complete the proof. ∎

##  Appendix F Proof of Corollary 6

###### Proof.

We extend the notation from Appendix E. For any $k\in\mathcal{K}$ and
$l\notin\mathcal{K}$, consider

$\displaystyle\mathbb{P}(\sum_{t=\nu+1}^{T}Z_{l,t}-Z_{k,t}\leq b)$
$\displaystyle=\mathbb{P}\left(\sum_{t=\nu+1}^{T}\log(\frac{\tilde{f}_{C_{t},l}({\bf
X}_{C_{t},t})}{\tilde{f}_{C_{t},k}({\bf X}_{C_{t},t})})\leq b\right)$

∎

Since $k\in\mathcal{K}$, there exists some $t_{0}$ s.t. ${\bf X}_{t}\sim
f_{k}$. Since the support of $f_{k}$ and $f_{l}$ are different,
$\log(\frac{\tilde{f}_{C_{t_{0}},l}({\bf
X}_{C_{t_{0}},t_{0}})}{\tilde{f}_{C_{t_{0}},k}({\bf
X}_{C_{t_{0}},t_{0}})})=-\infty$. For other times that ${\bf X}_{t}\nsim
f_{k}$, we have an undefined term of $\log\frac{0}{0}$. Since this term can be
interpreted as the log-likelihood difference at time $t$ and at this time,
neither $l$ and $k$ are selected, we may assume $\log\frac{0}{0}=0$ under this
assumption. Therefore,

$\mathbb{P}(\sum_{t=\nu+1}^{T}Z_{l,t}-Z_{k,t}\leq b)=1$

## References

  * CDC (2022) CDC.  Overview of testing for sars-cov-2 (covid-19), 2022.  URL https://www.cdc.gov/coronavirus/2019-ncov/hcp/testing-overview.html.
  * Chan et al. (2017) Hock Peng Chan et al.  Optimal sequential detection in multi-stream data.  _The Annals of Statistics_ , 45(6):2736–2763, 2017.
  * Chang and Yadama (2010) Shing I. Chang and Srikanth Yadama.  Statistical process control for monitoring non-linear profiles using wavelet filtering and b-spline approximation.  _International Journal of Production Research_ , 48(4):1049–1068, 2010.
  * Chen et al. (2020) Jie Chen, Wenyi Zhang, and H Vincent Poor.  A bayesian approach to sequential change detection and isolation problems.  _IEEE Transactions on Information Theory_ , 67(3):1796–1803, 2020.
  * Cho and Fryzlewicz (2015) Haeran Cho and Piotr Fryzlewicz.  Multiple-change-point detection for high dimensional time series via sparsified binary segmentation.  _Journal of the Royal Statistical Society: Series B: Statistical Methodology_ , pages 475–507, 2015.
  * Colosimo and Grasso (2018) Bianca M Colosimo and Marco Grasso.  Spatially weighted pca for monitoring video image data with application to additive manufacturing.  _Journal of Quality Technology_ , 50(4):391–417, 2018.
  * Dong et al. (2020) Ensheng Dong, Hongru Du, and Lauren Gardner.  An interactive web-based dashboard to track covid-19 in real time.  _The Lancet infectious diseases_ , 20(5):533–534, 2020.
  * Gómez et al. (2022) Ana María Estrada Gómez, Dan Li, and Kamran Paynabar.  An adaptive sampling strategy for online monitoring and diagnosis of high-dimensional streaming data.  _Technometrics_ , 64(2):253–269, 2022.
  * Gopalan et al. (2021) Aditya Gopalan, Braghadeesh Lakshminarayanan, and Venkatesh Saligrama.  Bandit quickest changepoint detection.  _Advances in Neural Information Processing Systems_ , 34:29064–29073, 2021.
  * Grasso et al. (2014) M Grasso, BM Colosimo, and M Pacella.  Profile monitoring via sensor fusion: the use of pca methods for multi-channel data.  _International Journal of Production Research_ , 52(20):6110–6135, 2014.
  * Grasso et al. (2017) Marco Grasso, Vittorio Laguzza, Quirico Semeraro, and Bianca Maria Colosimo.  In-process monitoring of selective laser melting: spatial detection of defects via image data analysis.  _Journal of Manufacturing Science and Engineering_ , 139(5), 2017.
  * Li and Jin (2010) Jing Li and Jionghua Jin.  Optimal sensor allocation by integrating causal models and set-covering algorithms.  _IIE Transactions_ , 42(8):564–576, May 2010\.
  * Liu and Shi (2013) Kaibo Liu and Jianjun Shi.  Objective-oriented optimal sensor allocation strategy for process monitoring and diagnosis by multivariate analysis in a Bayesian network.  _IIE Transactions_ , 45(6):630–643, June 2013\.
  * Liu et al. (2013) Kaibo Liu, Xi Zhang, and Jianjun Shi.  Adaptive sensor allocation strategy for process monitoring and diagnosis in a bayesian network.  _IEEE Transactions on Automation Science and Engineering_ , 11(2):452–462, 2013.
  * Liu et al. (2015) Kaibo Liu, Yajun Mei, and Jianjun Shi.  An Adaptive Sampling Strategy for Online High-Dimensional Process Monitoring.  _Technometrics_ , 57(3):305–319, July 2015.
  * Lorden (1971) Gary Lorden.  Procedures for reacting to a change in distribution.  _The Annals of Mathematical Statistics_ , pages 1897–1908, 1971.
  * Malladi and Speyer (1999) Durga P Malladi and Jason L Speyer.  A generalized shiryayev sequential probability ratio test for change detection and isolation.  _IEEE Transactions on Automatic Control_ , 44(8):1522–1534, 1999.
  * Mei (2010) Yajun Mei.  Efficient scalable schemes for monitoring a large number of data streams.  _Biometrika_ , 97(2):419–433, 2010.
  * Mei (2011) Yajun Mei.  Quickest detection in censoring sensor networks.  In _2011 IEEE International Symposium on Information Theory Proceedings_ , pages 2148–2152. IEEE, 2011.
  * Nikiforov et al. (1993) I Nikiforov, V Varavva, and V Kireichikov.  Application of statistical fault detection algorithms to navigation systems monitoring.  _Automatica_ , 29(5):1275–1290, 1993.
  * Nikiforov (1995) Igor V Nikiforov.  A generalized change detection problem.  _IEEE Transactions on Information theory_ , 41(1):171–187, 1995.
  * Paynabar and Jin (2011) Kamran Paynabar and Jionghua Jin.  Characterization of non-linear profiles variations using mixed-effect models and wavelets.  _IIE transactions_ , 43(4):275–290, 2011.
  * Paynabar et al. (2016) Kamran Paynabar, Changliang Zou, and Peihua Qiu.  A change-point approach for phase-i analysis in multivariate profile monitoring and diagnosis.  _Technometrics_ , 58(2):191–204, 2016.
  * Pollak (1985) Moshe Pollak.  Optimal detection of a change in distribution.  _The Annals of Statistics_ , pages 206–227, 1985.
  * Pollak (1987) Moshe Pollak.  Average run lengths of an optimal method of detecting a change in distribution.  _The Annals of Statistics_ , pages 749–779, 1987.
  * Qiu et al. (2010) Peihua Qiu, Changliang Zou, and Zhaojun Wang.  Nonparametric profile monitoring by mixed effects modeling.  _Technometrics_ , 52(3), 2010.
  * Ren et al. (2020) Haojie Ren, Changliang Zou, Nan Chen, and Runze Li.  Large-scale datastreams surveillance via pattern-oriented-sampling.  _Journal of the American Statistical Association_ , pages 1–15, 2020\.
  * Shiryaev (1963) Albert N Shiryaev.  On optimum methods in quickest detection problems.  _Theory of Probability & Its Applications_, 8(1):22–46, 1963.
  * Wang et al. (2018) Andi Wang, Xiaochen Xian, Fugee Tsung, and Kaibo Liu.  A spatial-adaptive sampling procedure for online monitoring of big data streams.  _Journal of Quality Technology_ , 50(4):329–343, October 2018.
  * Wang and Mei (2015) Yuan Wang and Yajun Mei.  Large-scale multi-stream quickest change detection via shrinkage post-change estimation.  _IEEE Transactions on Information Theory_ , 61(12):6926–6938, 2015.
  * Willsky (1976) Alan S Willsky.  A survey of design methods for failure detection in dynamic systems.  _Automatica_ , 12(6):601–611, 1976.
  * Xian et al. (2018) Xiaochen Xian, Andi Wang, and Kaibo Liu.  A Nonparametric Adaptive Sampling Strategy for Online Monitoring of Big Data Streams.  _Technometrics_ , 60(1):14–25, January 2018.
  * Xie and Siegmund (2013) Yao Xie and David Siegmund.  Sequential multi-sensor change-point detection.  In _2013 Information Theory and Applications Workshop (ITA)_ , pages 1–20. IEEE, 2013.
  * Yan et al. (2015) H. Yan, K. Paynabar, and J. Shi.  Image-based process monitoring using low-rank tensor decomposition.  _IEEE Transactions on Automation Science and Engineering_ , 12(1):216–227, 2015.
  * Yan et al. (2017) Hao Yan, Kamran Paynabar, and Jianjun Shi.  Anomaly detection in images with smooth background via smooth-sparse decomposition.  _Technometrics_ , 59(1):102–114, 2017.
  * Yan et al. (2018) Hao Yan, Kamran Paynabar, and Jianjun Shi.  Real-Time Monitoring of High-Dimensional Functional Data Streams via Spatio-Temporal Smooth Sparse Decomposition.  _Technometrics_ , 60(2):181–197, April 2018.
  * Yan et al. (2020) Hao Yan, Marco Grasso, Kamran Paynabar, and Bianca Maria Colosimo.  Real-time detection of clustered events in video-imaging data with applications to additive manufacturing.  _arXiv preprint arXiv:2004.10977_ , 2020.
  * Yue et al. (2017) Xiaowei Yue, Hao Yan, Jin Gyu Park, Zhiyong Liang, and Jianjun Shi.  A wavelet-based penalized mixed-effects decomposition for multichannel profile detection of in-line raman spectroscopy.  _IEEE Transactions on Automation Science and Engineering_ , 15(3):1258–1271, 2017.
  * Zhang and Mei (2020) Wanrong Zhang and Yajun Mei.  Bandit change-point detection for real-time monitoring high-dimensional data under sampling control.  _arXiv preprint arXiv:2009.11891_ , 2020.