Large Multistream Data Analytics for Monitoring and Diagnostics in Manufacturing Systems
Samaneh Ebrahimi Chitta Ranjan Kamran Paynabar

 




 
 


H. Milton Stewart School of Industrial and Systems Engineering,

          Georgia Institute of Technology, Atlanta, GA, USA




Large Multistream Data Analytics for Monitoring and Diagnostics in Manufacturing Systems
    Kamran Paynabar
    December 30, 2023
========================================================================================






The high-dimensionality and volume of large scale multistream data has inhibited significant research progress in developing an integrated monitoring and diagnostics (M&D) approach. This data, also categorized as big data, is becoming common in manufacturing plants.
In this paper, we propose an integrated M&D approach for large scale streaming data. We developed a novel monitoring method named Adaptive
Principal Component monitoring (APC) which adaptively chooses PCs
that are most likely to vary due to the change for early detection.
Importantly, we integrate a novel diagnostic approach, Principal Component Signal Recovery (PCSR), to enable a streamlined SPC. This diagnostics approach draws inspiration from
Compressed Sensing and uses Adaptive Lasso for identifying the sparse
change in the process. We theoretically motivate our approaches and
do a performance evaluation of our integrated M&D method through simulations
and case studies.





§ INTRODUCTION 


Recently, the problem of process monitoring and diagnosis using a large scale multi-stream data has become an active research area in statistical process control (SPC). The reason is two-fold: first, sensing technologies have enabled fast measurement of a large number of process variables, resulting in large data streams, and, second, conventional multivariate methods such
as Hotelling's T^2, MEWMA, and MCUSUM ( <cit.>)
are not scalable in terms of the computational time and the
detection power. 

An example of large streams can be found in gas turbine
systems used for power generation. In these systems, the performance
of the confined combustion process is being monitored using hundreds of
sensors measuring temperature, vibration, pressure, etc., in different
chambers and segments of the turbine. Early detection of any changes
in the system, followed by the diagnosis of the faulty variables is
necessary to avoid imminent blowout that leads to relighting the combustor and costly shutdowns. 

Another application of large streaming data is in image-based process monitoring in which each pixel of an image
can be considered as a single data stream. For example, in a rolling
process where a set of rollers are used to reduce the cross-section
of a long steel bar by applying compressive forces, the quality of
produced bars can be inspected by a vision system that is set up to
take images of the bar surface at short time intervals. A sample of
such an image is shown in Figure <ref>. In this image, each row contains 512 pixels, and each pixel can be considered as a variable resulting in an high-dimensional correlated data stream.




Despite its importance, existing SPC literature lacks a scalable integrated
M&D approach using large data streams. In the following, we will discuss the
existing work on monitoring and diagnostics, their shortcomings—especially 
their lack of integrability—, and motivate our approach.



  §.§.§ Monitoring


Conventional multivariate monitoring charts are effective on small or moderate data streams. However, their performance deteriorates as the number of data streams increases. 
To address the high-dimensionality issue, more recent works have focused
on employing variable selection techniques to reduce the dimensionality
by removing the variables that are less susceptible to the process
change. Examples of the variable-selection-based method include <cit.>.
However, these methods are not scalable and generally require intensive
 computation as the dimension grows. Moreover, most of
these methods are difficult to interpret for a process
engineer. 

There is a section of scalable multivariate monitoring methods that are
developed based on the assumption that data streams are independent
and the change is sparse (only a small subset of variables is affected by a process change). 
For example, <cit.>
assumed that exactly one variable changes at a time and
proposed an approach based on the maximum of CUSUM statistics from
each individual data stream. <cit.>
developed robust monitoring scheme based on the sum of (the top-r)
local CUSUM statistics assuming all variables are independent and
measurable. 

For the case that variables are not easily or efficiently
measurable, <cit.> presented TRAS (top-r based adaptive
sampling), which is an adaptive sampling strategy that uses the sum
of top r local CUSUM statistics for monitoring. The sparsity
assumption is generally valid in practice as a change or fault often
affects only a small subset of variables. However, although theoretically
and computationally appealing, the independence assumption is typically
unrealistic.

To address the dependency and high-dimensionality issues, Principal
Component Analysis (PCA) has been widely used for monitoring multivariate
data streams. PCA is a well-known projection technique that transforms
dependent data to uncorrelated features known as Principal Component
(PC) scores. Often, only a few PC scores that explain the most variation
of original data are used for monitoring in a dimension
reduction ( <cit.>).
However, monitoring top PCs with the highest variance may not always be a right approach.

As an example, consider a bivariate normal distribution given in Figure <ref>, 
in which PC_1 represents the direction of the eigenvector
with larger eigenvalue, and the red arrow indicates the direction
of change in the mean of the distribution. As can be seen from the
figure, the effect of the change on both PC-scores is the same. However,
the fact that PC1 constitutes the most of the process variance makes
it less sensitive to the change compared with PC2. In other words,
small changes may be masked by the variation present in top PC scores,
hence becoming undetectable. 



Other PC selection criteria for process
monitoring include the variance of reconstruction error (VRE) approach
by <cit.>, and the fault signal to noise ratio
(SNR) by <cit.> and <cit.>. The
VRE method selects a subset of PCs which minimizes fault reconstruction
error, while the fault SNR method is based on fault detection sensitivity.
The main disadvantage of these methods is that they require prior knowledge of the fault direction. 

In this paper, we present an adaptive
PC Selection (APC) approach based on hard-thresholding for selecting
the set of PC scores that are most susceptible to an unknown change.
Unlike the top-r-PCs, in our approach the number of features
may vary at each sampling time. Also the PCs are adaptively selected based on the observed sample and its standardized distance from the in-control mean. 
Additionally, the proposed APC approach does not require any prior knowledge
about a fault or change direction, which makes it more universally
applicable.



  §.§.§ Diagnostics


Another long-standing issue with PCA-based monitoring methods is the
lack of diagnosability. This is because that the PC scores used as
monitoring features are linear combinations of original measurements. Therefore, if a PC score initiates an out-of-control alarm, it is difficult to attribute it to any specific process variable(s). Interpretation and decomposition of these additive statistics are often theoretically difficult and/or computationally expensive
in high-dimensional data streams. 

For diagnostics on a PC-based monitoring, one
common approach is the use of contribution plots that specifies the
contribution of each variable to the out of control statistic  <cit.>.
Contribution plots are popular because of their ease of implementation
and their ability to work without any a priori knowledge. However,
correct isolation with contribution plots is not guaranteed for multiple
sensor faults  <cit.>. 

To overcome this problem, hierarchical contribution plots was proposed  <cit.>.
However, it will perform poorly if the initial partitioning is not correct.
Moreover, in the context of high-dimensional data, these methods become difficult
to interpret and are also computationally expensive.


For the purpose of diagnosis in multivariate control charts with original
measurements, <cit.> and  <cit.>
proposed variable selection techniques. Both methods optimize a penalized
likelihood function for multivariate normal observations to identify
the subset of altered variables. The ℒ_1-penalized regression
method of <cit.> provides more computational
advantages in implementation. <cit.> combined Bayesian
Information Criterion (BIC) with penalization techniques to assist
the fault localization process and suggested an Adaptive Lasso-based
diagnostic procedure. However, these methods assume that the change
point is already known and focus only on diagnosis. 

Additionally,
they cannot easily be integrated with a PCA-based monitoring approach.
To address these shortcomings, we propose a new diagnostics approach
that seamlessly integrates with our proposed PCA-based monitoring
method. The developed approach draws inspiration from Compressed Sensing and uses Adaptive Lasso to identify the shifted variables.
In this paper, we focus on detecting mean shifts and we assume
the shift is sparse. As mentioned earlier, this is a reasonable assumption
because in real-world usually only a small number of variables change.

The major contributions in this paper are, a. countering the traditional
view of top-PC scores as the best for process monitoring, and
proposing an adaptive PC selection approach as an alternative; and
b) proposing a new diagnostics approach that integrates with the proposed
PCA-based monitoring framework. 
An overview of the proposed Monitoring and Diagnostics (M&D) approach
is shown in Figure <ref>.















§ INTEGRATED PCA-BASED MONITORING AND DIAGNOSTICS








 §.§ Background


PCA is a linear transformation widely used for dimension reduction
and generating uncorrelated features. Suppose 
a p-dimensional data stream denoted as X={𝐱^(t):𝐱^(t)∈ℝ^p; t=1,2,…}
is collected at sampling time t. Without loss of generality, assume
the data streams are centered (zero mean) with a covariance matrix Σ.
By applying PCA, this set of correlated observations can be converted
into a set of linearly uncorrelated variables known as principal
component scores. The PC scores can be computed by 𝐲=A^T𝐱,
where A∈ℝ^p× p is the matrix of eigenvectors
of Σ and 𝐲∈ℝ^p are the PCs.
Also, it can be shown that var(𝐲_j)=λ_j and
cov(𝐲_j,𝐲_k)=0, ∀(j,k)∈{1,…,p}, j≠ k.

For ease of interpretation, the eigenvectors in A are arranged
such that their corresponding eigenvalues are in decreasing order,
i.e. λ_1≥λ_2≥…≥λ_p. This
ordering will be further referred throughout the paper for developing
our methodology. In this paper, we call the principal scores corresponding
to higher and lower eigenvalues as top-PCs and bottom-PCs, respectively.

In most conventional PCA-based methods, top-k PCs are selected
for monitoring because they contain more process information. This
approach, however, may not always result in an appropriate set of monitoring variables.
To illustrate this, we synthesized in-control samples of correlated
data from a multivariate normal distribution with a dimension of 500 (=p) and μ=0, σ=0.1 followed by out-of-control samples. In the out-of-control data, a random 10% set of variables are shifted by 0.05σ. We perform PCA on the data and
monitor all PCs separately.

Figure <ref> shows
the control charts for top 5 and bottom 5 of PCs. As shown in the
figure (left), the top PCs fail to detect the change. As mentioned
earlier, this is due to the fact that top PCs have large variances,
which make them insensitive to small shifts in the mean. On the other
hand, the bottom PCs with smaller variances are more sensitive and
can detect the small shift at the time it occurs, i.e., t=50. This
experiment was repeated several times, and each time similar results
were found. 

This shows that depending on the direction and the size of a
change, the traditional approach of selecting top PCs may severely
underperform. Therefore, it is imperative to develop a PC selection
approach that can adaptively select the set of most sensitive PCs and does not depend on the a priori knowledge about the direction
of the change.






 §.§ Adaptive PC selection (APC) for Process Monitoring 


In this section, we propose an adaptive PC Selection approach for
process monitoring. This approach selects and monitors a set
of PCs that show a higher deviation from a known in-control state.
Suppose, the in-control observations follow 𝐱_𝑡∼ N(0,Σ);t<τ,
and at an unknown time τ a mean shift occurs such that 𝐱_𝑡∼ N(μ,Σ);t⩾τ,
where, μ is a non-zero sparse vector, and the process covariance
is assumed to remain constant. Therefore, given the eigenvector matrix
A∈ℝ^p× p, the PC scores after the process change
will become 𝐲_𝑡=A^T𝐱_𝑡∼ N(A^Tμ,Λ),
where Λ=diag(λ_1,λ_2,…λ_p). 

The
standardized expected shift magnitude along the j^th PC can be
obtained by δ_j=||μ||cosθ_j/√(λ_j);j=1,2,...,p,
where θ_j is the angle between the shift direction and the
j^th PC. As can also be seen from Figure <ref>, this can imply that
a PC closer to the shift direction (i.e. smaller θ) will capture
a larger shift magnitude. Moreover, if θ is similar for two
PCs, then the one with the smaller variance will be more sensitive
to the change. Therefore, to take both measures into account, we work
with standardized PC score, denoted by ỹ_tj=y_tj√(λ_j),
that contains both magnitude and sensitivity information.







We choose the EWMA statistic for monitoring as it is more sensitive
to small changes and includes the information of previous samples.
The EWMA statistic, denoted by z_tj, is defined as z_tj=γỹ_tj+(1-γ)z_(t-1)j; t=1,2,...; j=1,2,...,p,
where z_0=0, and γ∈[0,1] is a weight. Under the in-control
process, z_tj∼ N(0,σ^2=γ/1-γ).
Consequently, its squared standardized value follows a Chi-squared
distribution with one degree of freedom, i.e., d_tj=(z_tj/√(γ/1-γ))^2∼χ_(1)^2.

When the process is out-of-control, depending on the direction of
the mean shift, a few d_tj values will inflate, while the rest
are slightly affected (or unaffected) by the mean shift. To increase
the detection power of the monitoring procedure, these PCs should
be filtered out. For this purpose, following <cit.>, we use a soft-thresholding operator to define the following aggregated monitoring statistics,


    R_t=∑_j=1^p(d_tj-ν)_+,


where the operator (·)_+=max{ 0,·}, and
ν is the threshold value selected based on a desired significance
level of χ^2 test. We monitor the R_t statistic and raise
an alarm if, R_t>R_0, where R_0 is the threshold level
found for a desired in-control ARL using simulations.

Selection of Control Limit (R_0)

To determine an appropriate value of R_0, we require the
distribution of monitoring statistic R. To find the distribution, we first specify the moments
of thresholded values in Proposition <ref>.

 
 
If d_tj∼χ_1^2, then their soft thresholded values d̃_tj=(d_tj-ν)_+ follows a bimodal truncated χ_1^2 distribution, with the following moments

E(d̃_tj)=E((d_tj-ν)_+)=1/Γ(0.5)[Γ(0.5,ν/2+e^-ν/2√(2ν)]-ν P(χ_1^2>ν)

E(d̃^2_tj)=E((d_tj-ν)_+^2)=1/Γ(0.5)[3Γ(0.5,ν/2+e^-ν/2√(2ν)(3+ν)]-2ν E(d̃)-ν^2P(χ_1^2>ν)


Proof is provided in Appendix <ref>. 

Using the Central Limit Theorem, R   N(pμ_d̃, √(p)σ_d̃).
Hence, R_0 for a desired type I error, α, is

    R_0=pμ_d̃+√(p)σ_d̃Φ^-1(1-α),


where Φ is the inverse normal cdf. To validate this approach and the normal
approximation, we perform simulations in Section <ref>.
The results show that the empirical α obtained by this approach
is very close to the true α.



 §.§ PC-based Signal Recovery (PCSR) Diagnosis Methodology


In monitoring high-dimensional data streams, apart from quick detection of changes,
precise fault diagnosis to identify accountable variables is extremely
crucial. Diagnosis aims at isolating
the shifted variables, which will help identify and eliminate the root
causes of a problem. However, despite its importance, very
few diagnostic methods exist for high-dimensional data streams that is 
integrable with a PCA-based monitoring.

To that end, we propose a diagnostics approach that seamlessly
integrates with the proposed PCA-based monitoring for large data streams.
Inspired by Compressed Sensing (CS), we develop an adaptive lasso
formulation to recover the variables responsible for an out-of-control
alarm. We assume that only the process mean has shifted and the
shift is sparse. 

In CS, a high-dimensional sparse original signal can be reconstructed
from noisy transformed observations by finding solutions to an englishunderdetermined
linear system. In other words, given a set of observations 𝐲,
and a transformation (sensing) matrix Υ, a sparse unknown
original signal μ can be recovered from 𝐲=Υμ+ϵ,
where ϵ denotes the random errors. 

The outcome of a PC monitoring method can be formulated similarly
to identify the shifted process variables. Without loss of generality,
we suppose the process has mean 0 during in-control that
changes to a sparse mean μ during the out-of-control of state.
Therefore, the out-of-control observations follow 𝐱=μ+ϵ,
where ϵ∼(0,Σ). Consequently, the out-of-control
PC scores are,


    𝐲=A𝐱=Aμ+ϵ̃,


where, ϵ̃=Aϵ is the noise in the
PC domain, with zero mean and covariance of Λ= diag(λ_1,λ_2,…λ_p).

Looking at Eq. <ref>, we can notice its similarity
with a compressed sensing problem. In Eq. <ref>,
the eigenvector matrix, A, and the principal scores, y, are
known, and we wish to estimate the shifted mean μ^ when an
out-of-control situation is detected after monitoring. <cit.>
and <cit.> showed that a least squares objective
function with L_1 penalty, also known as lasso, can be used to
estimate the sparse vector μ. Since lasso estimates in general
are not consistent, we use adaptive lasso <cit.>
to build our diagnosis model. Specifically, 


    μ̂=argmin_μ{||𝐲-Aμ||_l_2^2+r∑_j=1^pw_j|μ_j|,


where r is a nonnegative regularization parameter and w=1/μ̂_OLS is the data dependent weight vector.


One problem in solving Eq.<ref> is that the covariance matrix
of ϵ' is not homogeneous. The variance heterogeneity may affect
the estimation performance. To address this issue, we apply the following
transformation to get constant variances for all error terms. 

    y^*=Λ^-1/2y,     A^*=Λ^-1/2A,    ϵ^*=Λ^-1/2ϵ^'.


Consequently Eq. <ref> is transformed to 𝐲^*=A^*μ+ϵ^*,
where ϵ^*∼(0,I) with I
as a p dimensional identity matrix. The updated adaptive lasso
formulation is given by.


    μ̂=argmin_μ{||𝐲^*-A^*μ||_l_2^2+r∑_j=1^pw_j|μ_j|}


Where w_j=1/|μ̂_j|, and μ̂ is a
root p-consistent estimator to μ, e.g., μ̂=μ̂_ols.
The optimization problem in Eq. <ref> can be solved
using various optimization algorithms, such as gradient descent, proximal
descent, and LARS. In our implementation, we used the gradient descent
method. 

After finding the solution, the set of variables whose corresponding
estimated μ_j is non-zero is considered as the altered variables.
It should be noted that according to Theorem 2 in <cit.>,
the estimated mean is consistent, loosely meaning that when A^*
has large dimension (i.e. large p), the non-zero components of
μ are correctly identified. This implies that the larger the
number of data streams the higher is the likelihood of correct
diagnosis. See Appendix <ref> for more details. 

To determine the value of parameter r, one can use Bayesian
Information Criterion (BIC) (<cit.>) and choose
the r value that results in the smallest BIC value. The reason
behind choosing BIC is that it can determine the true sparse model
if the true model is included in the candidate set (<cit.>).
Since, in the diagnosis problem the objective is to detect the nonzero
elements (shifted variables) rather than estimation of the out-of-control
mean, BIC is a proper criterion for diagnosis <cit.>. 



§ EXPERIMENTAL ANALYSIS 


In this section, first we validate Proposition <ref> using
simulations. Afterwards, we study the performance of the proposed
monitoring-diagnostic method in change detection and in terms of quick
detection of mean shifts and identification of altered variables.
For all the experiments, we simulate data streams such that in-control data follows a multivariate normal distribution N(0,Σ) and the out-of-control is N(μ_1,Σ), μ_1 is sparse.
We carry out the simulations for different levels and types of shifts and the covariance structure and compare the results with existing methods as benchmarks.



 §.§ Validation of Proposition <ref> for Choosing Control Limits

 

To validate Proposition <ref>,
we perform two sets of experiments. In the first experiment, we generate
d_j∼χ_1^2 for j=1,⋯,p, and R_t=∑_j=1^p(d_tj-ν)_+
for t=1,⋯,1000, similar to Eq. <ref>. Then
we calculate R_0 using Eq. <ref> for desired α=0.05.
We calculate the empirical Type I error, denoted by α̃,
as the fraction of times R_t's pass the control limit R_0.
We perform this experiment for different values of p and ν
and replicate each scenario 1000 times. Finally, we report the average
empirical Type I errors in table <ref>.

In the second experiment, first we simulate x_t for
t=1,⋯,1000 as a p dimensional normal distribution random
variables with random covariance matrix and zero mean. Given the eigenvector
matrix A, we calculate its PC scores, and its
corresponding EWMA statistic using γ=0.4. 
Consequently, its squared standardized value are calculated as d_tj=(z_tj/√(γ/1-γ))^2.

Here, for each observation we define R_t=∑_j=1^p(d_tj-ν)_+,
we repeat this procedure 1000 times. Similar to previous experiment,
we calculate R_0 using Eq. <ref> for desired α=0.05.
We calculate the empirical Type I error, α̃, as the
fraction of times R_t's pass the control limit R_0. We replicate
each (p, ν) scenario 1000 times and we report the average empirical
type I errors in table <ref>.

As can be seen from Table <ref>-<ref>, 
as p increases, the empirical
Type I error approaches to its true value α=0.05. Moreover,
for large p, the result is less sensitive to the choice of the
threshold value, ν. Hence, it shows the validity of the proposed
approach for finding control limits. 

Note that the main difference
between these studies is the independence of R_t's. In the first
study, R_t's are independently generated, whereas in the second
study, R_t's are calculated using EWMA statistics, which are not
independent. The larger bias in the results of the second study is
mainly because the monitoring statistic values are autocorrelated.
However, for very large p (e.g, p>5000), this difference is
negligible. For smaller p, we would suggest using a Monte Carlo
simulation to determine the control limits.








 §.§ Monitoring Methods Analysis


In this section, we conduct various simulations to validate the performance
of the proposed monitoring method based on the Average Run Length
(ARL) and its standard error for different magnitudes of shifts. Specifically,
the following scenarios are considered: 


  I. Random covariance structure and random shift: To generate
the random covariance matrix, we use the Wishart distribution with
diagonal entries equal to 1. To generate a sparse mean shift, we randomly
select 20% of the process variables and shift them by δ. 

  II. Block diagonal covariance: This scenario mimics the situations
where each data stream is correlated with only a subset of the data
streams. The covariance matrix used in this scenario has K=12 blocks,
denoted as B_k, k=1,...,K. Each block B_k is a
random semi-positive definite matrix generated from a Wishart distribution.
To generate out-of-control data, we shift the mean of some of the variables
that belong to only one of the blocks (B_k), by δ , i.e.,
μ_j=[ δ j∈ B_k; 0 j∉ B_k ]. 










In each scenario, p data streams are generated. We run the simulations for, p= 100, p= 1000, and p= 10,000 to evaluate the performances in different dimensions.
We apply the proposed monitoring method, APC, and compare it with
three existing methods: 

a) T_new by <cit.>.
This monitoring method is based on a goodness-of-fit test of the local
CUSUM statistics from each data stream. 

b) TRAS by <cit.>.
Top-r based adaptive sampling (TRAS) is an adaptive sampling strategy
that uses the sum of top-r local statistics for monitoring. Since
this method works only for independent variables, we will implement
it on PCs rather than original data, 


c) Traditional PCA-based monitoring. In this approach, the selected
number of components to retain in the model is based on the cumulative
percentage of variance (CPV) equal to 90. Control charts are constructed by using
the Hotelling's T^2 statistic and the Q statistic (<cit.>).

To detect an out-of-control condition, the control limits are set
such that the ARL for in-control observations is equal to 200 (this
corresponds to a significance level of 0.005). Each control limit
is calculated through 1,000 replications. The results are shown in
Figures <ref>-<ref>. 

For p=100, as shown in Figure <ref>(a),
the proposed APC markedly outperforms the other benchmark methods. Even for shifts as small as δ=0.1σ, the ARL for APC is 4.06. This is about fourteen times smaller than the second best method, which is TRAS with ARL equal to 56. Moreover, for shifts δ≥0.1σ, APC detects the shift almost instantly (i.e., ARL_1=1). 

The results for Scenario II in Figure <ref>(b) also show that APC is superior to others, especially for moderate and large shifts. As an example, for a shift with the magnitude of
δ=0.25σ, APC's ARL is 17.67, while this values for the
best benchmark (TRAS) is 61.37.  As expected, the out-of-control ARL values for all methods in Scenario II is larger than those in Scenario I. 

For higher dimensions, the APC's ability to detect shifts becomes even better while the other method's performances stay the same or deteriorate. This shows that as dimension grows, the shifts are easier to be captured in PC scores that are in the direction of the shift.


To summarize, this study indicates that the APC method outperforms other methods for detecting
small values of shifts. Also, as dimension grows APC works better in detecting a change promptly. This can be attributed to the adaptive nature of the proposed monitoring statistic. 





















 §.§ Diagnosis Analysis


In this section, in addition to Scenarios I and II presented in the previous section, we add another scenario (Scenario III) with an autoregressive covariance
matrix, i.e., ρ_ij=|0.5|^(i-j) for variables i,j. This covariance matrix impose the variables close to each other to have higher covariance and as the variables go farther in the matrix, their covariance becomes smaller.
We validate the performance of the proposed diagnosis method for different
percentages of shifted variables (PS) as well as their shift magnitudes
(δ), using the following performance measures: (a) false negative
percentage (%FN), defined as the percentage of the number of variables
that are not detected over the number of all faulty variables; (b)
false positive percentage (%FP), defined as the percentage of the ratio of number of variables that are mistakenly detected as faulty
over the number of all not-faulty variables; (c) parameter selection
score (PSS), defined as the total number of variables that are labeled
incorrectly (either as faulty or not-faulty); and (d) F1-Score, defined
as the harmonic average of the precision and recall, and indicates
our overall performance combining the FP and FN. 

For FN, FP and PSS measures, the smaller the value, the better
the performance, whereas for F1-score the higher the better. We compare
the performance of our proposed method with the Lasso-based  diagnosis approach proposed by <cit.> called LEB. LEB
is an LASSO-based diagnostic approach approach for
diagnosis of sparse changes using BIC and the adaptive LASSO variable selection. The comparison results for p=100  are shown in Tables <ref>-<ref>, 
and Figures <ref>-<ref> for scenarios I, II, and III, respectively.  

In these tables PS denotes the percentage
of shifted process variables. As shown in Table <ref>
and Figure <ref>, under scenario I, when shift
occurs in a random set of variables with a random covariance matrix,
our proposed PCSR performs better than the LEB (<cit.>) for most of the cases
except for the case with PS=10% and small shifts (i.e., δ=0.7σ
). Even in this case, PCST is close to LEB. However,
for larger shifts, PCSR outperforms LEB. For instance, for a shift equal to 1σ and percentage of shifted variables equal to 10% the F1 accuracy using PCSR is equal to 0.9881 while it is equal to 0.9363 for LEB. 

In Scenario II, PCSR
clearly outperforms the LEB method. For example for a shift equal to 0.7σ
on 10% of variables, PCSR's F1-score is 0.6802 while, the LEB F1 score is 0.3725 (see Table <ref> and Figure <ref>). 

Also, the results in scenario III indicate the superior performance of that our method (see Table <ref> and Figure <ref>). For instance, for a  0.5σ shift on 25% of variables, PCSR's F1-score is 0.7173 and 0.4648 for LEB. 

These results show that for random non-sparse covariance matrices, the LEB method and PCSR method performs almost similarly. However, for sparse covariance matrices such as a block covariance or an autoregressive covariance, PCSR clearly outperforms LEB method. These sparse occurance of covariance matrices are very common in real world. This is because of the fact that in many situations, each data stream is correlated only with a small group of other data streams, but is not correlated with all other data streams collected in the system. Hence, a method that can detect the changes in such systems is necessary and more appropriate for real-world applications.

In sum, the results of the simulation study show the effectiveness
of our method in identifying the set of altered variables and its superiority over the current state-of-the-art.











































§ CASE STUDY 


In this section, we apply the proposed monitoring and diagnosis methods
on two case studies, a) defect detection in a steel rolling process,
and b) quality monitoring of wine. Additionally, we compare our results
with the existing methods.



 §.§ Defect detection in Steel Rolling Process


Early detection of process shifts in a rolling process is necessary
to avoid damage to products and reduce manufacturing costs. Rolling
is a high-speed process that makes its monitoring particularly challenging.
In this study, we show that the PCA-based method can effectively
detect anomalies and damages imprinted on a steel bar after rolling.
The dataset we consider here, includes images of size 128×512
pixels of the surface of rolled bars collected by a high-speed camera
<cit.>. Of the 100 images, the first 50 images are
in-control. One example of the image of rolling data for in-control vs
out-of-control process is shown in Figure <ref>.


We use this data to simulate an image with in-control observations
in the first 126 rows and out-of-control observations in the remaining
72 rows. The generated image is presented in Figure <ref>.
Also, we crop the image at the right end to avoid the non-informative dark segment of the image. Hence, our generated picture is of the size
of 198×300. In this study, each row of an image (a vector
of 300×1 ) is treated as an observation, creating a multi-stream data with the size of 300. As the picture shows, for out-of-control
observations, some small black lines, indicating anomalies, emerge
at the left part of the frame. We are interested to see whether our
monitoring approach can detect this change, and whether the diagnosis
approach can determine the changed pixels.






in each image there are streams of 128 strips, and with 100 images,
we will have 12,800 (128×100) data stream in total. Of the
12,800 data sample, we use the first 5,120 in-control data (the first
40 images) to obtain the control limits. The control limits are found
as explained in Sec. <ref> to
obtain in-control ARL of 200. 


Monitoring. We apply our proposed APC method for monitoring the process. We use
the first 70 in-control data (the first 70 rows of the image) to obtain
the control limits. The control limits are determined according to
the procedure explained in Sec. <ref>
to achieve the in-control ARL of 200. The resulting control chart
is shown in Figure <ref>. As can be seen
from the figure, after the change point, our monitoring statistic
instantly inflates and raises an out-of-control alarm by the first
observation after the change. Furthermore, to compare its performance
with existing state-of-the-art methods, we report the run
lengths (the number of observations before the change is detected)
for each method in Table <ref>. In the results, APC
has the smallest run-length (RL). This implies APC is the fastest in detecting
the change in comparison to other benchmarks.


The result for our APC monitoring method is shown in Figure <ref>.
Also comparing the RL after the change for APC, LPC, Conventional
PCA, Zou's method and Hao's method is provided in <ref>.
As this result shows LPC is fastest in detecting the change. 






Diagnosis. To check the performance of our PCSR method,
we performed diagnosis using our method vs LEB method on the out of
control data. The phase-1 data is used as the ground truth (sample size 70), and
25 out-of-control observations are used to detect the changed pixels
in the generated image. The area selected as out-of-control for each
method as well as the in-control and out-of-control images are shown
in Figure <ref>. The identified pixels are shown in black and the remaining
unchanged pixels are shown in white in Figure <ref>(c)
and (d), respectively. 

As the results show, PCSR method clearly detects
the changed pixels in the image with no false detection. Note that
although LEB can identify the changed pixels, it generates a few false
detection areas. 





 §.§ Wine Quality Monitoring


In this section, we demonstrate the efficacy of our proposed methodology
by applying it to a real dataset from a white wine production process.
The data is taken from the UCI data repository [http://archive.ics.uci.edu/ml/datasets/Wine+Quality].
The data has 4898 observations obtained between May 2004 to February
2007 for the purpose of improving the quality of Portuguese Vinho
Verde wine. The collected data has eleven variables named as fixed
acidity, volatile acidity, citric acid, residual sugar, chlorides,
free sulfur dioxide, total sulfur dioxide, density, PH, sulphates
and alcohol. An additional (manually annotated) quality
variable is available that will be used as ground truth for the
wine quality. This variable ranges between 0 (very bad) and 10 (very
excellent), and is provided based on sensory analysis <cit.>.

Our objective is to monitor the wine quality using
the variables and diagnose the shifted variables, if there is a shift.
We perform the APC study on this dataset. Similar to <cit.>'s
study on this data,
we focus on a subset of the data in which the quality variable is
either 6 or 7. The observations with the quality of 7 are considered as
acceptable while the rest are unacceptable, hence out-of-control. The quality variable is the ground truth that is used to gauge the performance of our monitoring—the monitoring should raise an out-of-control alarm as soon as the quality variable is going down from 7 to 6. When the alarm is raised, our diagnosis approach should be able to pinpoint the actual shifted process variables.

Monitoring. Overall there are 880 observations with the quality equal to 7 of
which 830 observations are used for phase I monitoring. Also, we set
the control limits (R_0 in APC method) such that ARL for in-control
observation is 1000. To do the comparison, we implement our method
along with the existing methods shown in Sec. <ref>.
All parameters in the methods are set to achieve in-control ARL of
1000 so that the methods are comparable.

For phase II monitoring, we use the remaining 40 points with the quality
of 7 followed by observations with the quality of 6. The goal is to investigate
how fast and accurately our monitoring algorithm detects the change
point in comparison to the existing methods.

The results are shown in Table <ref> and Figure <ref>. As shown in Table  <ref>,
the APC monitoring method detects the change after 11 observations.
On the other hand, other methods took more than twice as many observations to detect the change.





Diagnosis. 
Among the eleven variables, four were determined
as shifted variables by PCSR, viz. residual sugar, chlorides, density, and alcohol. 
The LEB method selected
chlorides, density and alcohol as the shifted variables. 



§ CONCLUSION 


In this paper, we proposed an SPC framework for high-dimensional data
streams that seamlessly integrates monitoring and diagnostics. We
proposed a new PCA-based monitoring approaches, viz. Adaptive PC Selection
(APC) monitoring. We first negated the common belief that the high-PCs
(principal components with highest variances) should be used for
monitoring, and then, showed that monitoring adaptively selected
PCs will be more effective. 
Using simulations, we showed that adaptively selected PCs outperforms
other benchmark methods for different types of covariance matrix structures
and types of shifts. Moreover, in all the stated scenarios, the conventional
approach of monitoring high-PC was shown to have poorer performance.


In the diagnosis module, we first discussed the challenge in finding
the shifted variables after a PCA-based monitoring procedure. The
challenge lies in isolating the process variable from the signaling PC. 
To address this, we used the CS principle to formulate an adaptive Lasso estimation
to detect the shifted variables. This formulation takes the eigenvectors
and principal components (after a shift) as inputs and yields the
process variables that caused the shift. Our experimental validations
showed that the proposed PCSR performs significanlty better than the current state-of-the-art.


Furthermore, we showed the practical applicability and validity of
our methods via real-world case studies. The first case study was
on defect detection in a steel rolling process, in which we found
that the proposed APC detects the shift faster than all the other
methods. Moreover, the PCSR diagnosis approach detects the change
pixels better than the existing method with fewer false positives. In another case study, we
monitored wine quality and diagnosed the shift. 
Our monitoring approach was again faster and our diagnosis approach
could find an additional shifted process variable, sugar, that
was undetected by the existing diagnostics approach. 

In this paper, we
have focused on monitoring and diagnosing the mean shifts. While the
developed APC can potentially be used to detect shifts in covariances,
further research is required to extend the PCSR diagnostics approach
to the covariance matrix monitoring. 


 





§ FIRST AND SECOND MOMENTS OF THE THRESHOLDED STATISTIC


The first moment of d̃ is calculated as,

    E(d̃)   =E((d_tj-ν)_+)=E((d_tj-ν)|d_tj>ν)P(d_tj>ν)
       =E((d_tj|d_tj>ν)p(d_tj>ν)-ν P(d_tj>ν)
       =∫_ν^∞d P(d)-ν P(χ_1^2>ν)
       =∫_ν^∞{x 1/√(2)Γ(0.5)x^-1/2e^-x/2} dx-ν P(χ_1^2>ν)
       =1/Γ(0.5){Γ(0.5,ν/2)+e^-ν/2√(2ν)}-ν P(χ_1^2>ν)


To calculate the second moment of d̃,

    E(d̃^2)   =E((d_tj-ν)_+^2)=E((d_tj-ν)_+^2)|d_tj>ν) P(d_tj>ν)
    
    =   E(d_tj^2|d_tj>ν) P(d_tj>ν)-2ν E(d_tj|d_tj>ν) P(d_tj>ν)+ν^2 P(d_tj>ν) Eq. <ref>→
    
    =   ∫_ν^∞{x^2 1/√(2)Γ(0.5)x^-1/2e^-x/2} dx-2ν{E(d̃)+ν P(χ_1^2>ν)}+ν^2P(χ_1^2>ν)
    
    =   ∫_ν^∞{x^2 1/√(2)Γ(0.5)x^-1/2e^-x/2} dx-2ν E(d̃)-ν^2P(χ_1^2>ν)
    
    =   1/Γ(0.5)[3Γ(0.5,ν/2)+e^-ν/2√(2ν)(3+ν)]-2ν E(d̃)-ν^2P(χ_1^2>ν)




§ CONSISTENCY OF THE DIAGNOSIS METHOD


To prove the consistency of our diagnosis model, we use the derivations in <cit.>. To show the consistency in adaptive lasso  used the following conditions,  

noise have independent identical distribution with mean 0 and variance σ^2


for observation matrix X,and number of observations n, 1/nX^T X→C. Where C is a positive definite matrix.


Condition <ref> is valid for Eq. <ref>. As we showed, after the transformations, the model noise has iid distribution with variance equal to 1.
To show the validity of condition <ref> we need to show that this consition holds for A^* instead of X. We use lemma <ref> and its proof to show it.


For A^* given in Eq. <ref> , C=1/pA^*^T A^* is a positive definite matrix



 The proof is a follows: 

    A^*   =Λ^-1/2A 
    C=1/pA^*^T A^*= 1/p A^T    Λ^-1/2Λ^-1/2 A ⇒C= A^T Λ^-1A


To show that C is positive definite matrix it suffice to show that x^T C x>0, ∀ x≠ 0

    x^T C x   =x^T A^T Λ^-1A x 
       =z^T Λ^-1  z
       =∑_i=1^p 1/λ_i z_i^2

Where λ_i is the pc score i. Since pc scores are positive, hence:

    x^T C x>0, ∀ x≠ 0



Since our model holds the above conditions, we can now show the consistency of our model, as follows:


Suppose that λ in Eq. <ref> varies with p. If λ_p√(p)→ 0, and λ_p →inf, then the adaptive lasso estimate must satisfy the following:


  * Consistency in variable selection: lim_p P(S^*=S)=1 

  * Asymptotic Normality: √(p)(μ̂^*_S-μ̂_S→_d N(0,σ^')

Where σ^'=C_11^-1





for proof of Theorem, please refer to Theorem 2 in  <cit.> 

