# Self-learning sparse PCA for multimode process monitoring

Jingxin Zhang, Donghua Zhou, , and Maoyin Chen This work was supported by
National Natural Science Foundation of China [grant numbers 62033008,
61873143]. (Corresponding authors: Donghua Zhou; Maoyin Chen)Jingxin Zhang is
with the Department of Automation, Tsinghua University, Beijing 100084, China
(e-mail: zjx18@mails.tsinghua.edu.cn). Donghua Zhou is with College of
Electrical Engineering and Automation, Shandong University of Science and
Technology, Qingdao 266000, China and also with the Department of Automation,
Tsinghua University, Beijing 100084, China (e-mail:
zdh@mail.tsinghua.edu.cn).Maoyin Chen is with the Department of Automation,
Tsinghua University, Beijing 100001, China and also with School of Automation
and Electrical Engineering, Linyi University, Linyi 276005, China (e-mail:
mychen@tsinghua.edu.cn). This paper has been submitted for IEEE Transactions
on Industrial Informatics for potential publication.

###### Abstract

This paper proposes a novel sparse principal component analysis algorithm with
self-learning ability for successive modes, where synaptic intelligence is
employed to measure the importance of variables and a regularization term is
added to preserve the learned knowledge of previous modes. Different from
traditional multimode monitoring methods, the monitoring model is updated
based on the current model and new data when a new mode arrives, thus
delivering prominent performance for sequential modes. Besides, the
computation and storage resources are saved in the long run, because it is not
necessary to retrain the model from scratch frequently and store data from
previous modes. More importantly, the model furnishes excellent
interpretability owing to the sparsity of parameters. Finally, a numerical
case and a practical pulverizing system are adopted to illustrate the
effectiveness of the proposed algorithm.

###### Index Terms:

Multimode process monitoring, sparse PCA, synaptic intelligence, self-learning

##  I Background

Multimode process monitoring is increasingly demanded and significant, as
industrial systems generally operate in varying modes due to raw materials,
changing load, etc [1, 2, 3]. Data from different modes have different
characteristics, such as mean value and variance [4, 5]. It is imperative to
investigate the effective manners for multimode processes [6, 7].

Marcos _et al_ [1] summarized the techniques for multimode processes and
divided the methods into two major categories, namely, single-model schemes
and multiple-model schemes. Single-model methods aim to find a transformation
to remove the multimodality features and then the fault is detected by a
decision function [3]. Multiple-model methods identify the mode and build the
monitoring model within each mode [8, 2, 9]. However, the information on all
modes should be complete, which is evidently impossible in real systems. When
a new mode arrives, the learned knowledge of previous modes may be overwritten
when training the same monitoring model for the current mode, which would lead
to a abrupt performance decrease. This phenomenon has been analyzed and
discussed in [10]. Overall, state-of-the-art methods aforementioned need to
retrain the monitoring model from scratch [10, 11].

Since the operating modes in practical systems appear successively, it is
consuming to progressively retain massive data and repeatedly retrain the
models. Therefore, it is important to establish a self-learning monitoring
model and update it continually when a new mode arrives. One strategy of
realizing self-learning ability is continual learning, where the model is
updated when new data are available. The technical core of continual learning
is to accommodate new information while preserving the acquired knowledge [12,
11]. However, there exists one longstanding challenge, namely, ‘catastrophic
forgetting’ issue, where the information of previous modes is overlapped by
new data and a new model based on new data may fail to monitor previous modes.
Various schemes have been developed to alleviate this issue and supply
excellent performance [12, 11, 13, 14]. Nevertheless, it is still scarce to
investigate the process monitoring techniques with self-learning ability [15,
16, 17].

Zhang _et al_ firstly focused on this research and illustrated the necessity
in [10], where a single monitoring model with continual learning ability was
investigated for successive modes. Elastic weight consolidation (EWC) [11] was
employed to solve the ‘catastrophic forgetting’ of traditional principal
component analysis (PCA) [10, 18], referred to as PCA-EWC, which retained the
significant information of previous modes by slowing down the changes of
certain influential parameters. Note that EWC estimates the importance measure
offline based on the point estimate of Fisher information matrix (FIM) [11].

Figure 1: Illustration of gradient decent optimization for different modes:
(a) The trajectory for one mode; (b) The trajectory when training the same
model on the second mode subsequently; (c) The trajectory when minimizing the
total loss from both modes (green) and gradients from each mode (red and blue)
[13].

Consider the poor interpretability of PCA, this paper investigates sparse PCA
(SPCA) with self-learning ability for multimode process monitoring, where the
model is updated based on the current model and new data. Instead of EWC, we
consider synaptic intelligence (SI), which calculates the importance matrix
along the entire learning trajectory [19, 20, 21], by computing the gradients
of loss and the parameter update. For convenience, the proposed SPCA with SI
is denoted as SPCA-SI. Compared with PCA-EWC [10], SPCA-SI utilizes $L_{1}$
regularization to enhance the sparsity of critical parameters, thus providing
better model interpretability. Besides, the importance measure by SI is easier
to estimate than EWC, because the gradients are usually available while FIM is
intractable. Moreover, sparse representation is also beneficial to reduce
catastrophic forgetting, as there are fewer model-sensitive parameters [13,
20, 14]. In this paper, we assume that the mode transition is accomplished in
a short time and the switching time is available.

The rest of this paper is organized below. Section II reviews SPCA and SI
briefly, and introduces the research problem. The proposed SPCA-SI algorithm
is elaborated in Section III and settled by accelerated proximal gradient
descent (APG) method. A novel $T^{2}$ statistic is proposed and the monitoring
procedure is summarized in Section IV. Besides, the influence of parameters
and computational complexity are discussed. The effectiveness of the proposed
approach is illustrated by a numerical case study and a practical coal
puzzling system in Section V. The concluding remark is given in Section VI.

##  II Preliminary

###  II-A Revisit of SPCA

Here we give another perspective of SPCA, where the projection vectors are
acquired one by one.

Given the dataset $\boldsymbol{X}\in R^{N\times m}$, $N$ is the number of
samples and $m$ is the number of variables. PCA aims to minimize the
reconstruction error, namely,

$\begin{array}[]{l}\min\;\left\|{\boldsymbol{X}-\boldsymbol{X}\boldsymbol{p}{\boldsymbol{p}^{T}}}\right\|_{F}^{2}\\\
s.t.\quad{\boldsymbol{p}^{T}}\boldsymbol{p}=1\end{array}$

where $\boldsymbol{p}\in R^{m}$ is the projection vector.

Consider the virtues of sparsity, $L_{1}$ regularization is adopted to enhance
interpretability and alleviate catastrophic forgetting simultaneously [13,
14]. Thus, the objective of SPCA is

$\begin{array}[]{l}\min\;\left\|{\boldsymbol{X}-\boldsymbol{X}\boldsymbol{p}{\boldsymbol{p}^{T}}}\right\|_{F}^{2}+\lambda\left\|\boldsymbol{p}\right\|_{1}\\\
s.t.\quad{\boldsymbol{p}^{T}}\boldsymbol{p}=1\end{array}$ (1)

where $\lambda$ is a regularization parameter. After $\boldsymbol{p}$ is
calculated, let
$\boldsymbol{X}=\boldsymbol{X}-\boldsymbol{X}\boldsymbol{p}{\boldsymbol{p}^{T}}$,
and repeat (1) until $l$ projection vectors are obtained. Here, the number of
principal components $l$ is determined by cumulative percent variance.

###  II-B Review of SI

The synaptic framework was proposed and detailed information has been
described in [19]. Here we overview the key points [19, 20].

For a learning process, we aim to seek for the optimal parameter
$\boldsymbol{\theta}$ given the objective function $J(\boldsymbol{\theta})$.
Gradient-based methods are effective manners to solve the optimization
problem. SI estimates the importance measure for each parameter along the
learning trajectory, which reflects the sensitivity of each parameter to the
loss.

The gradient is a conservative field, and the value of the integral along the
trajectory is equal to the difference between the end point and the start
point [19]. Consider an infinitesimal parameter update
$\boldsymbol{\delta}\left(k\right)$ at $k$th iteration, the change in loss is
approximated by

$J\left({\boldsymbol{\theta}\left(k\right)+\boldsymbol{\delta}\left(k\right)}\right)-J\left({\boldsymbol{\theta}\left(k\right)}\right)\approx\sum\nolimits_{i}{{g_{i}}\left(k\right){\delta_{i}}\left(k\right)}$

where $\boldsymbol{g}=\frac{{\partial J}}{{\partial{\boldsymbol{\theta}}}}$ is
the gradient,
${\delta_{i}}\left(k\right)=\theta_{i}\left(k\right)-\theta_{i}\left(k-1\right)$.
During the whole learning process, the change in loss over the entire
trajectory is calculated by

$\begin{aligned}
\sum\limits_{k}{\boldsymbol{g}\left(k\right)\boldsymbol{\delta}\left(k\right)}=&\sum\limits_{i}{\sum\limits_{k}{{g_{i}}\left(k\right){\delta_{i}}\left(k\right)}}\\\
=&-\sum\limits_{i}{{\varpi_{i}}}\end{aligned}$

More intuitively [21],

${\varpi_{i}}=\sum\limits_{k}{\left({{\theta_{i}}\left(k\right)-{\theta_{i}}\left({k-1}\right)}\right)\frac{{-\partial
J}}{{\partial{\theta_{i}\left(k\right)}}}}$ (2)

Then, the importance measure is normalized by [20]

${{\bar{\varpi}}_{i}}=\max\left({0,\frac{{{\varpi_{i}}}}{{{{\left({\Delta{\theta_{i}}}\right)}^{2}}+\zeta}}}\right)$
(3)

thus the regularization term and the loss shares the same unit.
$\Delta{\theta_{i}}{\rm{=}}\sum\nolimits_{k}{\left({{\theta_{i}}\left(k\right)-{\theta_{i}}\left({k-1}\right)}\right)}$
is the total change for each parameter and $\zeta$ is added to avoid ill-
conditioning issue.

###  II-C Problem reformulation

This paper proposes a self-learning process monitoring approach for successive
modes, which is built based on SPCA and the model is updated when a new mode
arrives. We take two modes to depict the research problem by Fig. 1 [13].

When training the monitoring model for mode $\mathcal{M}_{1}$, the
optimization issue is settled by gradient decent method and the trajectory of
SPCA is exhibited in Fig. 1. When a new mode $\mathcal{M}_{2}$ arrives,
traditional SPCA-based methods generally train the model subsequently, as
illustrated in Fig. 1. In this case, the learning of mode $\mathcal{M}_{2}$
leads to an overlap of the learned knowledge, which indicates that the
retrained model is not efficient for the previous mode $\mathcal{M}_{1}$. This
paper aims to accommodate new data by a continually updated model while
accumulating the learned knowledge, thus delivering brilliant performance for
two or more modes. As shown in Fig. 1, the total loss for both modes is
considered simultaneously and the optimal solution is an equilibrium between
the gradients of different modes [13].

##  III Methodology

In this section, we present the procedure of SPCA-SI for the first mode and
the objective is optimized by APG. Then, SPCA-SI is extended to more general
cases.

###  III-A SPCA-SI for the first mode

For the first mode $\mathcal{M}_{1}$, data $\boldsymbol{X}_{1}$ are collected.
We settle the issue (1) by augmented Lagrangian function:

$J=\left\|{\boldsymbol{X}_{1}-\boldsymbol{X}_{1}\boldsymbol{p}{\boldsymbol{p}^{T}}}\right\|_{F}^{2}+\lambda{\left\|\boldsymbol{p}\right\|_{1}}+\mu{\left({\boldsymbol{p}^{T}}\boldsymbol{p}-1\right)^{2}}$
(4)

where $\mu$ is the Lagrange parameter. (4) is nonconvex and nonsmooth, and can
not be directly settled by gradient-based methods.

APG is an effective technique to deal with this type of optimization issue and
employed in this paper [22]. We divide (4) into smooth part
$g\left(\boldsymbol{p}\right)$ and nonsmooth part
$h\left(\boldsymbol{p}\right)$, namely,

$g\left(\boldsymbol{p}\right)=\left\|{\boldsymbol{X}_{1}-\boldsymbol{X}_{1}\boldsymbol{p}\boldsymbol{p}^{T}}\right\|_{F}^{2}+\mu{\left({\boldsymbol{p}^{T}}\boldsymbol{p}-1\right)^{2}}$
(5) $h\left(\boldsymbol{p}\right)=\lambda{\left\|\boldsymbol{p}\right\|_{1}}$
(6)

###  III-B Solution with APG

The procedure of APG contains the gradient-based part of
$g\left(\boldsymbol{p}\right)$ and the proximal gradient.

For the smooth part $g\left(\boldsymbol{p}\right)$, (5) is equivalent to

$\begin{aligned}
g\left(\boldsymbol{p}\right)=&tr\left({\boldsymbol{X}_{1}^{T}}\boldsymbol{X}_{1}\right)+\mu+tr\left({\boldsymbol{p}{\boldsymbol{p}^{T}}\left({{\boldsymbol{X}_{1}^{T}}\boldsymbol{X}_{1}+\mu\boldsymbol{I}}\right)\boldsymbol{p}{\boldsymbol{p}^{T}}}\right)\\\
&-2tr\left({\boldsymbol{p}{\boldsymbol{p}^{T}}\left({{\boldsymbol{X}_{1}^{T}}\boldsymbol{X}_{1}+\mu\boldsymbol{I}}\right)}\right)\end{aligned}$

Thus, the gradients are

${\nabla_{\boldsymbol{p}}}g\left(\boldsymbol{p}\right)=\frac{{\partial
g}}{{\partial{\boldsymbol{p}}}}=\boldsymbol{p}{\boldsymbol{p}^{T}}\boldsymbol{G}_{1}\boldsymbol{p}+\boldsymbol{G}_{1}\boldsymbol{p}{\boldsymbol{p}^{T}}\boldsymbol{p}-2\boldsymbol{G}_{1}\boldsymbol{p}$
(7) ${\nabla_{\mu}}g\left(\boldsymbol{p}\right)=\frac{{\partial
g}}{{\partial\mu}}=\left({\boldsymbol{p}^{T}}\boldsymbol{p}-1\right)^{2}$ (8)

where
$\boldsymbol{G}_{1}=2\left({{\boldsymbol{X}_{1}^{T}}\boldsymbol{X}_{1}+\mu\boldsymbol{I}}\right)$.

For $h\left(\boldsymbol{p}\right)$, the proximal function is defined as [22,
23]

$\displaystyle{\boldsymbol{p}^{+}}=$
$\displaystyle\arg\mathop{\min}\limits_{\boldsymbol{z}}\frac{1}{{2t}}\left\|{\boldsymbol{z}-\left({\boldsymbol{p}-t\nabla_{\boldsymbol{p}}g\left(\boldsymbol{p}\right)}\right)}\right\|_{2}^{2}+h\left(\boldsymbol{z}\right)$
(9) $\displaystyle:=$ $\displaystyle
pro{\boldsymbol{x}_{h,t}}\left({\boldsymbol{p}-t\nabla_{\boldsymbol{p}}g\left(\boldsymbol{p}\right)}\right)$

Inspired by Adam [24], the learning rate $t$ is adaptively estimated to
accelerate convergence. At $k$th iteration, $t_{k}$ is calculated by

$\displaystyle{t_{k}}=$ $\displaystyle
f\left({\alpha,{t_{k-1}},\nabla{g_{k}}}\right)$ (10) $\displaystyle=$
$\displaystyle{\alpha\mathord{\left/{\vphantom{\alpha{\left({\sqrt{{{\left({{\tau_{2}}{t_{k-1}}+\left({1-{\tau_{1}}}\right){{\left\|{\nabla{g_{k}}}\right\|}^{2}}}\right)}\mathord{\left/{\vphantom{{\left({{\tau_{2}}{t_{k-1}}+\left({1-{\tau_{1}}}\right){{\left\|{\nabla{g_{k}}}\right\|}^{2}}}\right)}{\left({1-{\tau_{2}}}\right)}}}\right.\kern-1.2pt}{\left({1-{\tau_{2}}}\right)}}}+\varepsilon}\right)}}}\right.\kern-1.2pt}{\left({\sqrt{{{\left({{\tau_{2}}{t_{k-1}}+\left({1-{\tau_{1}}}\right){{\left\|{\nabla{g_{k}}}\right\|}^{2}}}\right)}\mathord{\left/{\vphantom{{\left({{\tau_{2}}{t_{k-1}}+\left({1-{\tau_{1}}}\right){{\left\|{\nabla{g_{k}}}\right\|}^{2}}}\right)}{\left({1-{\tau_{2}}}\right)}}}\right.\kern-1.2pt}{\left({1-{\tau_{2}}}\right)}}}+\varepsilon}\right)}}$

where $\nabla{g_{k}}$ is the corresponding gradient. $\alpha$, ${\tau_{1}}$
and ${\tau_{2}}$ are constants. $\varepsilon$ is added to avoid ill-
conditioning issue.

Algorithm 1 APG for optimization issue (4)

Input: Initialize $\boldsymbol{p}_{1}=\boldsymbol{p}_{0}$,
$\boldsymbol{z}_{1}=\boldsymbol{p}_{0}$, $t_{1}=t_{0}=0$,
${\boldsymbol{\bar{\omega}}}_{1}=\boldsymbol{0}$, $\tau_{1}=0.9$,
$\tau_{2}=0.999$, $\varepsilon=10^{-8}$, $t^{y}_{0}=10^{-4}$,
$t^{p}_{0}=10^{-4}$, $t^{\mu}_{0}=10^{-4}$, $\alpha^{p}=0.001$,
$\alpha^{\mu}=0.01$
Output: the optimal $\boldsymbol{p}$, and the importance measure
$\boldsymbol{\omega}$
for $k=1,2,3,\cdots$ do

  1. 1.

Update the projection vector:

    1. (a)

${\boldsymbol{y}_{k}}={\boldsymbol{p}_{k}}+\frac{{{t_{k-1}}}}{{{t_{k}}}}\left({{\boldsymbol{z}_{k}}-{\boldsymbol{p}_{k}}}\right)+\frac{{{t_{k-1}}-1}}{{{t_{k}}}}\left({{\boldsymbol{p}_{k}}-{\boldsymbol{p}_{k-1}}}\right)$

    2. (b)

${\boldsymbol{z}_{k+1}}=pro{x_{h,{t^{y}}}}\left({{\boldsymbol{y}_{k}}-{t^{y}_{k}}\nabla_{\boldsymbol{p}}g\left({{\boldsymbol{y}_{k}}}\right)}\right)$,
calculate
$t_{k}^{y}=f\left({{\alpha^{p}},t_{k-1}^{y},{\nabla_{\boldsymbol{p}}}g\left({{\boldsymbol{y}_{k}}}\right)}\right)$
by (10)

    3. (c)

${\boldsymbol{v}_{k+1}}=pro{x_{h,{t^{p}}}}\left({{\boldsymbol{p}_{k}}-{t^{p}_{k}}\nabla_{\boldsymbol{p}}g\left({{\boldsymbol{p}_{k}}}\right)}\right)$,
calculate
$t_{k}^{p}=f\left({{\alpha^{p}},t_{k-1}^{p},{\nabla_{\boldsymbol{p}}}g\left({{\boldsymbol{p}_{k}}}\right)}\right)$
by (10)

    4. (d)

${t_{k+1}}=\frac{{\sqrt{4{{\left({{t_{k}}}\right)}^{2}}+1}+1}}{2}$

    5. (e)

${\boldsymbol{p}_{k+1}}=\left\\{{\begin{array}[]{*{20}{c}}{{\boldsymbol{z}_{k+1}},}&{if\;J\left({{\boldsymbol{z}_{k+1}}}\right)\leq
J\left({{\boldsymbol{v}_{k+1}}}\right)}\\\
{{\boldsymbol{v}_{k+1}},}&{otherwise}\end{array}}\right.$

  2. 2.

Update $\mu$,
$\mu_{k+1}=\mu_{k}+t^{\mu}_{k}{\nabla_{\mu}}g\left(\boldsymbol{p}_{k+1}\right)$,
$t_{k}^{\mu}=f\left({{\alpha^{\mu}},t_{k-1}^{\mu},{\nabla_{\mu}}g\left({\boldsymbol{p_{k+1}}}\right)}\right)$

  3. 3.

Calculate the importance measure
${\boldsymbol{\bar{\omega}}}_{k+1}={\boldsymbol{\bar{\omega}}}_{k}-\left(\left(\nabla
g\left({{\boldsymbol{p}_{k+1}}}\right)\right)^{T}\odot\left({{\boldsymbol{p}_{k+1}}-{\boldsymbol{p}_{k}}}\right)^{T}\right)^{T}$

end for
Normalize $\bar{\boldsymbol{\omega}}$ by (13) and denote as
${{\boldsymbol{\omega}}}$

The proximal function $prox$ is defined and the proximal gradient is
calculated by the soft threshold [22]

$\begin{aligned}
pro{x_{h,t}}\left(\boldsymbol{p}\right)=&\arg\mathop{\min}\limits_{\boldsymbol{z}}\frac{1}{{2t}}\left\|{\boldsymbol{z}-\boldsymbol{p}}\right\|_{2}^{2}+\lambda{\left\|\boldsymbol{z}\right\|_{1}}\\\
=&{S_{\lambda t}}\left(\boldsymbol{p}\right)\end{aligned}$

The soft threshold ${S_{\lambda t}}\left(\boldsymbol{p}\right)$ has an
analytical solution [25]:

${\left[{{S_{\lambda
t}}}\right]_{i}}=\left\\{{\begin{array}[]{*{20}{l}}{{p_{i}}-\lambda
t,}&{{p_{i}}>\lambda t}\\\ {0,}&{\left|{{p_{i}}}\right|\leq\lambda t}\\\
{{p_{i}}+\lambda t,}&{{p_{i}}<\lambda t}\end{array}}\right.$ (11)

where $p_{i}$ is the $i$th element of $\boldsymbol{p}$.

According to (2), the importance measure is computed by

$\bar{\boldsymbol{\omega}}=\sum\limits_{k}\left(\left(\nabla
g\left({{\boldsymbol{p}_{k}}}\right)\right)^{T}\odot\left({{\boldsymbol{p}_{k}}-{\boldsymbol{p}_{k-1}}}\right)^{T}\right)^{T}$
(12)

where $\odot$ denotes the Khatri-Rao product and $\boldsymbol{p}_{k}$ is the
projection vector at $k$th iteration step. Accordingly, each element of
$\bar{\boldsymbol{\omega}}$ is normalized by

${{\omega}_{i}}=\max\left({0,\frac{{{\bar{\omega}_{i}}}}{{{{\left({\Delta{p_{i}}}\right)}^{2}}+\zeta}}}\right)$
(13)

where $\Delta{p_{i}}$ represents the total change, $1\leq i\leq m$. The
solution is summarized in Algorithm 1.

The procedure of SPCA-SI is summarized in Algorithm 2. For convenience, the
optimal projection matrix and the importance measure are denoted as
${\boldsymbol{P}}_{\mathcal{M}_{1}}$ and
${\boldsymbol{\Omega}}_{\mathcal{M}_{1}}$, respectively.

Algorithm 2 Procedure of SPCA-SI

0: data $\boldsymbol{X}$, $l$

0: The projection matrix $\boldsymbol{P}$, the importance measure
$\boldsymbol{\Omega}$

1: Initialize
$\boldsymbol{P}^{0}=\left[{\begin{array}[]{*{20}{c}}{{\boldsymbol{p}^{0}_{1}}}&\cdots&{{\boldsymbol{p}^{0}_{l}}}\end{array}}\right]={I_{m,l}}$,
$j=1$;

2: Scale $\boldsymbol{X}$ to zero mean and unit variance;

3: Let $\boldsymbol{p}_{0}={\boldsymbol{p}^{0}_{j}}$, solve (4) by APG as
summarized in Algorithm 1. The gradients are calculated by (7-8);

4: The optimal projection vector and importance measure are denoted as
${\boldsymbol{p}_{j}}$ and ${\boldsymbol{\omega}_{j}}$. Deflate
$\boldsymbol{X}$ as
$\boldsymbol{X}:=\boldsymbol{X}-\boldsymbol{X}\boldsymbol{p}_{j}{\boldsymbol{p}_{j}^{T}}$;

5: Let $j=j+1$, return to step 3 until $j>l$;

6:
${\boldsymbol{P}=\left[{\begin{array}[]{*{20}{c}}{{\boldsymbol{p}_{1}}}&\cdots&{{\boldsymbol{p}_{l}}}\end{array}}\right]}$,
${\boldsymbol{\Omega}=\left[{\begin{array}[]{*{20}{c}}{{\boldsymbol{\omega}_{1}}}&\cdots&{{\boldsymbol{\omega}_{l}}}\end{array}}\right]}$.

###  III-C SPCA-SI for multiple modes

When the mode $\mathcal{M}_{i}$ ($i\geq 2$) arrives, only data
$\boldsymbol{X}_{i}$ are available for training and data from previous modes
are not retained. SPCA-SI aims to learn the new mode while consolidating the
acquired knowledge of previous modes. As shown in Fig. 1, we need to minimize
the total loss of all modes, with the constraint that the loss functions of
previous trained modes are unavailable. To alleviate catastrophic forgetting,
drastic changes to influential parameters in the past should be avoided.
Therefore, a quadratic surrogate loss is introduced to approximate the total
loss of previous modes [19].

The model of SPCA-SI is updated based on the current model and new data. For
the $j$th projection vector ($1\leq j\leq l$), the objective is designed as:

$\displaystyle\min$
$\displaystyle\left\|{\boldsymbol{X}_{i}-\boldsymbol{X}_{i}\boldsymbol{p}{\boldsymbol{p}^{T}}}\right\|_{F}^{2}+{\left({\boldsymbol{p}-{\boldsymbol{p}_{\mathcal{M}_{i-1}}}}\right)^{T}}\bar{\boldsymbol{\Omega}}\left(\boldsymbol{p}-{\boldsymbol{p}_{\mathcal{M}_{i-1}}}\right)$
(14) $\displaystyle+\lambda\left\|\boldsymbol{p}\right\|_{1}$ $\displaystyle
s.t.$ $\displaystyle{\boldsymbol{p}^{T}}\boldsymbol{p}=1$

where ${\boldsymbol{p}_{\mathcal{M}_{i-1}}}$ is the $j$th column of
${\boldsymbol{P}_{\mathcal{M}_{i-1}}}$,
$\bar{\boldsymbol{\Omega}}=\gamma_{i}diag\left(\boldsymbol{\omega}\right)$,
$\boldsymbol{\omega}$ is the importance measure corresponding to
${\boldsymbol{p}_{\mathcal{M}_{i-1}}}$ and the $j$th column of
${{\mathord{\buildrel{\lower
3.0pt\hbox{$\scriptscriptstyle\smile$}}\over{\boldsymbol{\Omega}}}}_{\mathcal{M}_{i-1}}}$,
${{\mathord{\buildrel{\lower
3.0pt\hbox{$\scriptscriptstyle\smile$}}\over{\boldsymbol{\Omega}}}}_{\mathcal{M}_{i-1}}}=\sum\nolimits_{r=1}^{i-1}{{{\boldsymbol{\Omega}}_{\mathcal{M}_{r}}}}$,
and $\gamma_{i}$ is the weight which trades off previous versus current modes.
Similar to [11, 19], the additional regularization term is the quadratic
surrogate loss, which makes the optimal parameters of current mode close to
the previous one, with a small loss.

Similarly, the augmented Lagrangian function is depicted as

$\displaystyle J=$
$\displaystyle\left\|{\boldsymbol{X}_{i}-\boldsymbol{X}_{i}\boldsymbol{p}{\boldsymbol{p}^{T}}}\right\|_{F}^{2}+\lambda\left\|\boldsymbol{p}\right\|_{1}+\mu\left({\boldsymbol{p}^{T}}\boldsymbol{p}-1\right)^{2}$
(15)
$\displaystyle+{\left({\boldsymbol{p}-{\boldsymbol{p}_{\mathcal{M}_{i-1}}}}\right)^{T}}\bar{\boldsymbol{\Omega}}\left(\boldsymbol{p}-{\boldsymbol{p}_{\mathcal{M}_{i-1}}}\right)$

The smooth part $g\left(\boldsymbol{p}\right)$ and the corresponding gradient
are

$\displaystyle g\left(\boldsymbol{p}\right)=$
$\displaystyle\left\|{\boldsymbol{X}_{i}-\boldsymbol{X}_{i}\boldsymbol{p}\boldsymbol{p}^{T}}\right\|_{F}^{2}+{\left({\boldsymbol{p}-{\boldsymbol{p}_{\mathcal{M}_{i-1}}}}\right)^{T}}\bar{\boldsymbol{\Omega}}\left(\boldsymbol{p}-{\boldsymbol{p}_{\mathcal{M}_{i-1}}}\right)$
(16)
$\displaystyle+\mu{\left({\boldsymbol{p}^{T}}\boldsymbol{p}-1\right)^{2}}$
$\displaystyle{\nabla_{\boldsymbol{p}}}g\left(\boldsymbol{p}\right)=$
$\displaystyle\boldsymbol{p}{\boldsymbol{p}^{T}}\boldsymbol{G}_{i}\boldsymbol{p}+\boldsymbol{G}_{i}\boldsymbol{p}{\boldsymbol{p}^{T}}\boldsymbol{p}-2\boldsymbol{G}_{i}\boldsymbol{p}+2\bar{\boldsymbol{\Omega}}\left(\boldsymbol{p}-{\boldsymbol{p}_{\mathcal{M}_{i-1}}}\right)$
(17)

where
$\boldsymbol{G}_{i}=2\left({{\boldsymbol{X}_{i}^{T}}\boldsymbol{X}_{i}+\mu\boldsymbol{I}}\right)$.

The solution can refer to Algorithm 2. The optimization problem is (15), and
the gradients are calculated by (17) and (8). The optimal projection matrix
and importance measure are denoted as $\boldsymbol{P}_{\mathcal{M}_{i}}$ and
$\boldsymbol{\Omega}_{\mathcal{M}_{i}}$, respectively.

##  IV Monitoring model and Discussion

###  IV-A Monitoring statistics

Two statistics are designed to monitor the operating condition. To enhance the
monitoring performance for previous modes, the partial covariance information
of last mode is adopted to calculate $T^{2}$ statistic. For mode
$\mathcal{M}_{i}$ ($i\geq 1$),

${T^{2}}={\boldsymbol{x}^{T}}\boldsymbol{P}_{\mathcal{M}_{i}}{\boldsymbol{\Xi}_{\mathcal{M}_{i}}^{-1}}{\boldsymbol{P}_{\mathcal{M}_{i}}^{T}}\boldsymbol{x}$
(18)

where $\boldsymbol{x}\in\boldsymbol{X}_{i}$,
$\boldsymbol{P}_{\mathcal{M}_{i}}$ is the projection matrix,
$\boldsymbol{\Xi}_{\mathcal{M}_{i}}=\boldsymbol{P}_{\mathcal{M}_{i}}^{T}\left(\eta\frac{{{\boldsymbol{X}_{i}^{T}}\boldsymbol{X}_{i}}}{{N_{i}-1}}+(1-\eta)\boldsymbol{P}_{\mathcal{M}_{i-1}}{\boldsymbol{\Xi}}_{\mathcal{M}_{i-1}}{\boldsymbol{P}_{\mathcal{M}_{i-1}}^{T}}\right)\boldsymbol{P}_{\mathcal{M}_{i}}$,
and $N_{i}$ is the number of samples, $\eta$ trades off the previous versus
current modes with $0\leq\eta\leq 1$. $\boldsymbol{P}_{\mathcal{M}_{i-1}}$ and
${\boldsymbol{\Xi}}_{\mathcal{M}_{i-1}}$ are acquired from mode
$\mathcal{M}_{i-1}$, which represent the information of previous modes without
storing the original data. When $i=1$, let $\eta=1$,
$\boldsymbol{\Xi}_{\mathcal{M}_{1}}=\frac{{{\boldsymbol{X}_{1}^{T}}\boldsymbol{X}_{1}}}{{N_{1}-1}}$.
If $i>1$, $\eta$ is estimated by the importance of previous modes.
Correspondingly, the squared prediction error (SPE) is calculated by

$SPE={\boldsymbol{x}^{T}}\left({I-\boldsymbol{P}_{\mathcal{M}_{i}}{\boldsymbol{P}_{\mathcal{M}_{i}}^{T}}}\right)\boldsymbol{x}$
(19)

The thresholds are calculated by kernel density estimation (KDE) [26] and the
confidence level is $99\%$. Once one statistic is beyond its threshold, a
fault is detected. The training procedure is summarized in Algorithm 3.

Algorithm 3 Off-line training phase of SPCA-SI

1: For the mode $\mathcal{M}_{1}$, collect data $\boldsymbol{X}_{1}$;

2: Normalize $\boldsymbol{X}_{1}$ to zero mean and unit variance;

3: Perform traditional PCA on $\boldsymbol{X}_{1}$ and calculate the number of
principal components $l$;

4: Solve the issue (4) by Algorithm 2, acquire
$\boldsymbol{P}_{\mathcal{M}_{1}}$ and
${\boldsymbol{\Omega}}_{\mathcal{M}_{1}}$. The gradients are calculated by
(7-8);

5: Calculate statistics by (18-19) and thresholds by KDE;

6: For the mode $\mathcal{M}_{i}$ ($i\geq 2$), collect data
$\boldsymbol{X}_{i}$;

7: Scale $\boldsymbol{X}_{i}$ to zero mean and unit variance;

8: Solve the issue (15) by Algorithm 2, acquire
$\boldsymbol{P}_{\mathcal{M}_{i}}$ and
${\boldsymbol{\Omega}}_{\mathcal{M}_{i}}$. The gradients are calculated by (8)
and (17) ;

9: Calculate statistics by (18-19) and thresholds by KDE.

###  IV-B Discussion

####  IV-B1 Parameter setting

SPCA-SI has three regularization parameters. $\lambda$ affects the sparsity
and SPCA is transformed to PCA when $\lambda=0$. Thus, SPCA-SI is equivalent
to PCA-SI. Similarly, PCA-SI offers continual learning ability but the
parameters are not sparse.

Then, we explain how the continual learning ability is influenced by
$\gamma_{i}$ and $\eta$. Two extreme cases are given as an example. When
$\gamma_{i}=0$ and $\eta=1$, SPCA-SI is equivalent to traditional SPCA and
information of previous modes is forgotten catastrophically (similar to Fig.
1). When $\gamma_{i}\to\infty$ and $\eta=0$, the information of mode
$\mathcal{M}_{1}$ is completely preserved while the knowledge of subsequent
modes is not learned. In other cases, information of different modes is
memorized and beneficial to monitor multiple modes simultaneously.

####  IV-B2 Computational complexity

The computational complexity focuses on Algorithms 1 and 2. $k_{total}$ is the
total number of iterations. Here we use $flam$ to reflect the complexity. The
calculation of $\boldsymbol{G}_{i}$ needs $\frac{1}{2}N_{i}m^{2}+2m$ flam. For
mode $\mathcal{M}_{1}$, $10m^{2}+38m+14$ flam is required for each iteration
in Algorithm 1. Thus, the computational complexity is
$(10m^{2}+38m+14)k_{total}+l(\frac{1}{2}N_{1}m^{2}+3N_{1}m+2m)$ flam for
training. For mode $\mathcal{M}_{i}$ ($i\geq 2$), Algorithm 1 needs
$10m^{2}+50m+14$ flam in total for each iteration. The training phase needs
$(10m^{2}+50m+14)k_{total}+l(\frac{1}{2}N_{i}m^{2}+3N_{i}m+2m)$ flam.

##  V Case study

This section adopts two case studies to illustrate the effectiveness of SPCA-
SI. Recursive PCA (RPCA) [27] and improved mixture of probabilistic PCA
(IMPPCA) [26] are employed for comparison. For IMPPCA, the mode is
automatically identified by membership degree. The fault detection rate (FDR)
and false alarm rate (FAR) are employed to evaluate the monitoring
performance.

TABLE I: Comparative scheme for numerical case | Methods |  Training sources |  Model label |  Testing sources
---|---|---|---|---
Situation 1 | SPCA | Mode 1 | A | Mode 1
Situation 2 | SPCA-SI | Model A+Mode 2 | B | Mode 2
Situation 3 | SPCA-SI | - | B | Mode 1
Situation 4 | SPCA | Mode 2 | C | Mode 2
Situation 5 | SPCA | - | C | Mode 1
Situation 6 | RPCA | Modes 1,2 | D | Mode 1
Situation 7 | RPCA | - | D | Mode 2
Situation 8 | IMPPCA | Modes 1,2 | E | Mode 1
Situation 9 | IMPPCA | - | E | Mode 2

###  V-A Numerical case

The following numerical case is adopted:

$\left[{\begin{array}[]{*{20}{c}}{{x_{1}}}\\\ {{x_{2}}}\\\ {{x_{3}}}\\\
{{x_{4}}}\\\ {{x_{5}}}\\\ {{x_{6}}}\\\ {{x_{7}}}\\\
{{x_{8}}}\end{array}}\right]=\left[{\begin{array}[]{*{20}{c}}{0.55}&{0.82}&{0.94}\\\
{0.23}&{0.45}&{0.62}\\\ {-0.61}&{0.62}&{0.41}\\\ {0.49}&{0.79}&{0.89}\\\
{0.89}&{-0.92}&{0.06}\\\ {0.76}&{0.74}&{0.35}\\\ {0.46}&{0.28}&{0.81}\\\
{-0.02}&{0.41}&{0.01}\end{array}}\right]\left[{\begin{array}[]{*{20}{c}}{{s_{1}}}\\\
{{s_{2}}}\\\ {{s_{3}}}\end{array}}\right]+\boldsymbol{e}$

where the noise $\boldsymbol{e}$ follows Gaussian distribution with
$e_{i}\sim\mathbb{N}(0,0.001),i=1,\cdots,8$. Sequential data are generated
successively from two modes:
$\bullet$ Mode 1: $s_{1}\sim\mathbb{U}([-10,-9.7])$,
$s_{2}\sim\mathbb{N}(-5,1)$, and $s_{3}\sim\mathbb{U}([2,3])$;
$\bullet$ Mode 2: $s_{1}\sim\mathbb{U}([-6,-5.7])$,
$s_{2}\sim\mathbb{N}(-1,1)$, and $s_{3}\sim\mathbb{U}([3,4.2])$;
where $\mathbb{U}([-10,-9.7])$ represents the uniform distribution between
$-10$ and $-9.7$, and so on.

We generate 1000 normal samples to train the monitoring model and 1000 samples
for fault detection, including the first 500 normal samples and the last 500
faulty samples through the following scenarios:
$\bullet$ Fault 1: step fault of $x_{3}$, $x_{3}=x_{3}^{*}+0.08$;
$\bullet$ Fault 2: step fault of $x_{6}$, $x_{6}=x_{6}^{*}+0.08$;
$\bullet$ Fault 3: slope drift of $x_{1}$, $x_{1}=x_{1}^{*}+0.001(k-500)$;

where $500\leq k\leq 1000$, $x_{1}^{*}$, $x_{3}^{*}$ and $x_{6}^{*}$ are
normal.

The simulation scheme is designed to illustrate the effectiveness and
superiorities of SPCA-SI, as summarized in Table I. Note that ‘-’ indicates
that there is no need to retrain the model and the current monitoring model is
adopted for fault detection. The first five situations are utilized to verify
that SPCA-SI alleviates the catastrophic forgetting of traditional SPCA and
furnishes self-learning ability. Specifically, the results of Situations 1 and
4 testify the effectiveness of SPCA for a mode. The monitoring results of
Situations 2 and 3 are used to prove that SPCA-SI is able to monitor two modes
simultaneously by a continually updated model, which assimilates new data when
a new mode appears. Situation 5 is designed to show that SPCA fails to monitor
the previous mode and the learned knowledge is overwritten by new information.
For situations 6-9, RPCA and IMPPCA are compared to illustrate the
superiorities of SPCA-SI further.

() Situation 2

() Situation 3

() Situation 4

() Situation 5

() Situation 6

() Situation 7

() Situation 8

() Situation 9

Figure 2: Monitoring charts of Fault 1 TABLE II: FDR ($\%$) and FAR ($\%$) for numerical case Fault type | Fault 1 | Fault 2 | Fault 3
---|---|---|---
Indexes | FDR | FAR | FDR | FAR | FDR | FAR
Situation 1 | 100 | 7.4 | 100 | 2.6 | 96.8 | 0
Situation 2 | 100 | 6.6 | 100 | 2.2 | 91.0 | 0
Situation 3 | 98.6 | 2.4 | 99.6 | 1.0 | 90.6 | 4.6
Situation 4 | 100 | 8.4 | 100 | 2.6 | 96.4 | 0
Situation 5 | 100 | 93.4 | 100 | 90 | 98.8 | 65.2
Situation 6 | 100 | 98.4 | 100 | 98.4 | 100 | 98.4
Situation 7 | 100 | 33.2 | 100 | 100 | 100 | 100
Situation 8 | 100 | 2.2 | 100 | 2.6 | 95.2 | 1.4
Situation 9 | 100 | 2.2 | 100 | 1.8 | 96.0 | 2.0

Take Fault 1 as an example to interpret the monitoring consequences in detail,
as described in Fig. 2. As the monitoring charts for Mode 1 and Mode 2 are
similar for SPCA, the simulation chart of Situation 1 is not listed. In Figs.
2-2, SPCA-SI enables to detect the fault in Modes 1 and 2 accurately by model
B, which is updated based on the existing model A and data from Mode 2. The
model C fails to detect the fault in Mode 1, as depicted in Fig. 2. The FDR is
$100\%$ and FAR is $93.4\%$, which indicates that the learned knowledge of
Mode 1 is forgotten when training the model C (similar to Fig. 1). According
to Situations 1-5, SPCA-SI alleviates the catastrophic forgetting of SPCA and
provides self-learning ability for successive modes. For Situation 6-7, RPCA
is unable to track the system changes and distinguish the novelty from
normality. In Figs. 2-2, IMPPCA detects the fault accurately and the FDRs are
$100\%$.

The monitoring results of three faults are summarized in Table II. The
analysis aforementioned is equally applied to Fault 2 and Fault 3. Traditional
SPCA forgets the significant features of previous modes when training the same
model sequentially and fails to deliver prominent performance. SPCA-SI is
capable of monitoring two modes simultaneously based on a continually updated
model and settles the catastrophic forgetting of SPCA. RPCA fails to
distinguish between the normal modes and faults. IMPPCA detects the fault
accurately and needs to be retrained from scratch when new modes appear. Thus,
it consumes much more storage space and computational resource than SPCA-SI.
Overall, SPCA-SI with self-learning ability is superior to IMPPCA and RPCA, as
it can monitor multiple modes accurately, and the computation and storage
resources are saved in the long run.

Figure 3: Schematic diagram of the coal pulverizing system

###  V-B Pulverizing system process monitoring

The 1000-MW ultra-supercritical thermal power plant is increasingly popular
due to potential economic benefits and low pollution [10]. This paper
investigates the coal pulverizing system in Zhoushan Power Plant in China,
which provides high quality pulverized coal for boiler. As depicted in Fig. 3,
it contains coal feeder, coal mill, rotary separator, raw coal hopper and
stone coal scuttle. In practical systems, the types of coal and unit load
would change frequently, thus generating successive operating modes.

In this paper, we focus on two typical faults, namely, abnormality from outlet
temperature (Fault 4) and rotary separator (Fault 5). The information is
summarized in Table III. Note that the number of training samples and testing
samples are shorted for NoTrS and NoTeS, respectively. Nine key variables are
selected by prior knowledge. Assume that when new modes appear, the system
operates under the normal condition at the preliminary stage.

To illustrate the effectiveness and self-learning ability of SPCA-SI, 17
situations are designed in Table IV. Three different modes are considered in
this case. Similar to the numerical case, Situations 1, 4 and 9 are utilized
to illustrate the effectiveness of SPCA for a mode. Situations 2, 3, 6, 7 and
8 are designed to verify that SPCA-SI can monitor several modes
simultaneously, where new data are assimilated while preserving the learned
knowledge. Situations 5, 10 and 11 are employed to show the catastrophic
forgetting issue of SPCA. For Situations 12-14, RPCA is adopted to track the
successive modes. Take Situation 12 as an example, Modes 1-3 appear
sequentially and Mode 1 occurs again, and the fault occurs in Mode 1. For
Situations 15-17, IMPPCA is adopted to monitor three modes and the training
data are required to be complete.

TABLE III: Data information of the coal pulverizing system Fault type |  Mode number | NoTrS | NoTeS |  Fault location
---|---|---|---|---
Fault 4 |  $1$ | 2160 | 2880 | 909
$2$ | 1080 | 1080 | 533
$3$ | 1440 | 1440 | 626
Fault 5 |  $1$ | 2880 | 1080 | 806
$2$ | 720 | 720 | 352
$3$ | 2880 | 2160 | 134
TABLE IV: Simulation scheme for the pulverizing system | Methods |  Training sources |  Model label |  Testing sources
---|---|---|---|---
Situation 1 | SPCA | Mode 1 | A | Mode 1
Situation 2 | SPCA-SI | Model A+Mode 2 | B | Mode 2
Situation 3 | SPCA-SI | - | B | Mode 1
Situation 4 | SPCA | Mode 2 | C | Mode 2
Situation 5 | SPCA | - | C | Mode 1
Situation 6 | SPCA-SI | Model B+Mode 3 | D | Mode 3
Situation 7 | SPCA-SI | - | D | Mode 1
Situation 8 | SPCA-SI | - | D | Mode 2
Situation 9 | SPCA | Mode 3 | E | Mode 3
Situation 10 | SPCA | - | E | Mode 1
Situation 11 | SPCA | - | E | Mode 2
Situation 12 | RPCA | Modes 1, 2, 3 | F | Mode 1
Situation 13 | RPCA | - | F | Mode 2
Situation 14 | RPCA | - | F | Mode 3
Situation 15 | IMPPCA | Modes 1, 2, 3 | H | Mode 1
Situation 16 | IMPPCA | - | H | Mode 2
Situation 17 | IMPPCA | - | H | Mode 3

The monitoring results of Fault 4 and Fault 5 are summarized in Table V. Take
Fault 4 as an instance to explain the results detailedly. SPCA can detect the
fault in Mode 1 accurately, but the FAR is $5.62\%$. When a new mode 2
arrives, the model is updated based on the model A and the newly collected
data. Thus, the monitoring model B is able to monitor the two successive modes
simultaneously, and the FDRs are higher than $99\%$. Besides, the FAR of
Situation 3 is lower than that of Situation 1, which indicates that
information of Mode 2 is beneficial to reduce the false alarms of Mode 1. The
monitoring model C enables to monitor two modes. But the FAR of Mode 1 is
$6.5\%$ and higher than that of SPCA-SI. When the new mode 3 appears, the
proposed SPCA-SI trains the model D based on the model B and the current data
collected. It enables to monitor the three modes simultaneously and the FDRs
approach to $100\%$. Besides, the FARs are the lowest among all situations. It
is revealed that SPCA-SI can preserve partial significant information of
trained modes, which is advantageous to monitor other similar modes. RPCA
fails to detect the fault in three modes and the FARs approximate to $100\%$.
Although IMPPCA is capable of monitoring modes 1 and 2, the FAR for Mode 3 is
$41.12\%$. IMPPCA fails to deliver desired expert level monitoring
performance.

Owing to the paper length limitation, we just select 8 representative
monitoring charts of Fault 4, as exhibited in Fig. 4. The simulation results
of SPCA-SI are mainly listed. Two charts of RPCA and IMPPCA are selected as
comparison. The analysis of Fault 4 also applies to Fault 5. Note that the FAR
of IMPPCA for Mode 2 is $51.28\%$.

According to Table IV, SPCA is vulnerable to catastrophic forgetting issue and
fails to monitor multiple modes based on the same model. SPCA-SI provides
self-learning ability and the model is updated when a new mode arrives, which
enables it to monitor multiple modes accurately. Moreover, the learned
knowledge of previous modes is preserved continually and the model-sensitive
parameters are fewer than PCA. Similar to numerical case, RPCA is incapable of
separating normal modes and faults. Besides, IMPPCA is unable to monitor three
modes accurately as the FARs are more than $20\%$. When a new mode appears, we
need to store data and retrain the model from scratch, which costs
considerable resources and energy. In conclusion, SPCA-SI outperforms other
comparative methods in consideration of detection accuracy and demanding
resources in the long term.

TABLE V: FDR ($\%$) and FAR ($\%$) for the practical case Fault type | Fault 4 | Fault 5
---|---|---
Indexes | FDR | FAR | FDR | FAR
Situation 1 | 99.95 | 5.62 | 100 | 0
Situation 2 | 99.45 | 0 | 100 | 5.98
Situation 3 | 99.95 | 4.07 | 100 | 0
Situation 4 | 99.45 | 0 | 100 | 15.38
Situation 5 | 99.95 | 6.5 | 100 | 0
Situation 6 | 100 | 0.32 | 93.49 | 0
Situation 7 | 99.95 | 1.54 | 100 | 0
Situation 8 | 99.45 | 0 | 100 | 13.96
Situation 9 | 100 | 0.48 | 92.75 | 0
Situation 10 | 99.95 | 75.77 | 100 | 0
Situation 11 | 100 | 100 | 100 | 94.87
Situation 12 | 100 | 99.45 | 100 | 100
Situation 13 | 100 | 100 | 100 | 100
Situation 14 | 100 | 100 | 100 | 100
Situation 15 | 100 | 6.61 | 100 | 3.23
Situation 16 | 100 | 6.95 | 100 | 51.28
Situation 17 | 100 | 41.12 | 98.96 | 0

() Situation 2

() Situation 3

() Situation 6

() Situation 7

() Situation 8

() Situation 11

() Situation 14

() Situation 17

Figure 4: Monitoring charts of Fault 4

##  VI Conclusion

This paper presented a novel SPCA-SI method with self-learning ability for
monitoring successive modes. The importance measure of variables is evaluated
by SI along the learning trajectory. The acquired knowledge of previous modes
is accumulated and the model is updated when new data are available, thus
delivering excellent performance for successive modes. The optimization issue
is settled by APG and the learning rate is adaptively determined to accelerate
convergence. Besides, the influence of parameters is discussed and different
methods can be converted by specific parameter setting. SPCA-SI furnishes
excellent model interpretability, as the critical parameters are sparse.
Moreover, a novel $T^{2}$ statistic is presented, where the significant
information of previous modes is consolidated further. The effectiveness of
the proposed method has been illustrated by a numerical case and a practical
industrial system.

This proposed method requires the similarity among different modes and prior
information about mode switching time. In future, we’ll investigate the
numerous and diverse modes, with the mode switching time identified
automatically.

## References

  * [1] M. Qui$\tilde{n}$ones-Grueiro, A. Prieto-Moreno, C. Verde, and O. Llanes-Santiago, “Data-driven monitoring of multimode continuous processes: A review,” _Chemometrics and Intelligent Laboratory Systems_ , vol. 189, pp. 56–71, 2019.
  * [2] B. Wang, Z. Li, Z. Dai, N. Lawrence, and X. Yan, “Data-driven mode identification and unsupervised fault detection for nonlinear multimode processes,” _IEEE Transactions on Industrial Informatics_ , vol. 16, no. 6, pp. 3651–3661, 2020.
  * [3] K. Zhang, K. Peng, and J. Dong, “A common and individual feature extraction-based multimode process monitoring method with application to the finishing mill process,” _IEEE Transactions on Industrial Informatics_ , vol. 14, no. 11, pp. 4841–4850, 2018.
  * [4] L. Zhou, J. Zheng, Z. Ge, Z. Song, and S. Shan, “Multimode process monitoring based on switching autoregressive dynamic latent variable model,” _IEEE Transactions on Industrial Electronics_ , vol. 65, no. 10, pp. 8184–8194, 2018.
  * [5] Y. Liu, J. Zeng, J. Bao, and L. Xie, “A unified probabilistic monitoring framework for multimode processes based on probabilistic linear discriminant analysis,” _IEEE Transactions on Industrial Informatics_ , vol. 16, no. 10, pp. 6291–6300, 2020.
  * [6] Y. Jiang and S. Yin, “Recent advances in key-performance-indicator oriented prognosis and diagnosis with a MATLAB toolbox: DB-KIT,” _IEEE Transactions on Industrial Informatics_ , vol. 15, no. 5, pp. 2849–2858, 2019.
  * [7] J. Zhang, M. Chen, H. Chen, X. Hong, and D. Zhou, “Process monitoring based on orthogonal locality preserving projection with maximum likelihood estimation,” _Industrial & Engineering Chemistry Research_, vol. 58, no. 14, pp. 5579–5587, 2019.
  * [8] X. Peng, Y. Tang, W. Du, and F. Qian, “Multimode process monitoring and fault detection: A sparse modeling and dictionary learning method,” _IEEE Transactions on Industrial Electronics_ , vol. 64, no. 6, pp. 4866–4875, 2017.
  * [9] J. Shang, D. Zhou, M. Chen, H. Ji, and H. Zhang, “Incipient sensor fault diagnosis in multimode processes using conditionally independent bayesian learning based recursive transformed component statistical analysis,” _Journal of Process Control_ , vol. 77, pp. 7–19, 2019.
  * [10] J. Zhang, D. Zhou, and M. Chen, “Monitoring multimode processes: a modified PCA algorithm with continual learning ability,” _arXiv:2012.07044_ , 2020\.
  * [11] J. Kirkpatrick, R. Pascanu, N. Rabinowitz, J. Veness, G. Desjardins, A. A. Rusu, K. Milan, J. Quan, T. Ramalho, and A. Grabska-Barwinska, “Overcoming catastrophic forgetting in neural networks,” _Proceedings of the National Academy of Sciences of the United States of America_ , vol. 114, no. 13, pp. 3521–3526, 2017.
  * [12] R. Aljundi, K. Kelchtermans, and T. Tuytelaars, “Task-free continual learning,” in _2019 IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2019, pp. 11 254–11 263.
  * [13] R. Hadsell, D. Rao, A. A. Rusu, and R. Pascanu, “Embracing change: Continual learning in deep neural networks,” _Trends in Cognitive Sciences_ , vol. 24, no. 12, pp. 1028–1040, 2020.
  * [14] J. Yoon1, E. Yang, J. Lee, and S. J. Hwang, “Lifelong learning with dynamically expandable networks,” in _International Conference on Learning Representations_ , 2018.
  * [15] X. Xu, H. Yang, C. Lian, and J. Liu, “Self-learning control using dual heuristic programming with global Laplacian eigenmaps,” _IEEE Transactions on Industrial Electronics_ , vol. 64, no. 12, pp. 9517–9526, 2017\.
  * [16] F. Ye, K. Chakrabarty, Z. Zhang, and X. Gu, “Self-learning and adaptive board-level functional fault diagnosis,” in _The 20th Asia and South Pacific Design Automation Conference_ , 2015, pp. 294–301.
  * [17] J. Feldmann, N. Youngblood, C. D. Wright, H. Bhaskaran, and W. H. P. Pernice, “All-optical spiking neurosynaptic networks with self-learning capabilities.” _Nature_ , vol. 569, no. 7755, pp. 208–214, 2019.
  * [18] J. Zhang, D. Zhou, and M. Chen, “Monitoring nonstationary processes based on recursive cointegration analysis and elastic weight consolidation,” _arXiv:2101.08579_ , 2021.
  * [19] F. Zenke, B. Poole, and S. Ganguli, “Continual learning through synaptic intelligence,” in _Proceedings of the 34th International Conference on Machine Learning_ , vol. 70, 2017, pp. 3987–3995.
  * [20] N. Y. Masse, G. D. Grant, and D. J. Freedman, “Alleviating catastrophic forgetting using context-dependent gating and synaptic stabilization,” _Proceedings of the National Academy of Sciences_ , vol. 115, no. 44, pp. E10 467–E10 475, 2018.
  * [21] G. M. van de Ven, H. T. Siegelmann, and A. S. Tolias, “Brain-inspired replay for continual learning with artificial neural networks.” _Nature Communications_ , vol. 11, no. 1, p. 4069, 2020.
  * [22] H. Li and Z. Lin, “Accelerated proximal gradient methods for nonconvex programming,” in _Proceedings of the 28th International Conference on Neural Information Processing Systems - Volume 1_ , ser. NIPS’15. MIT Press, 2015, p. 379–387.
  * [23] A. Beck and M. Teboulle, “Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems,” _IEEE Transactions on Image Processing_ , vol. 18, no. 11, pp. 2419–2434, 2009.
  * [24] D. P. Kingma and J. L. Ba, “Adam: A method for stochastic optimization,” in _ICLR 2015 : International Conference on Learning Representations 2015_ , 2015.
  * [25] N. K. Dhingra, S. Z. Khong, and M. R. Jovanovic, “The proximal augmented Lagrangian method for nonsmooth composite optimization,” _IEEE Transactions on Automatic Control_ , vol. 64, no. 7, pp. 2861–2868, 2019.
  * [26] J. Zhang, H. Chen, S. Chen, and X. Hong, “An improved mixture of probabilistic PCA for nonlinear data-driven process monitoring,” _IEEE Transactions on Cybernetics_ , vol. 49, no. 1, pp. 198–210, 2019.
  * [27] W. Li, H. Yue, S. Valle-Cervantes, and S. Qin, “Recursive PCA for adaptive process monitoring,” _Journal of Process Control_ , vol. 10, no. 5, pp. 471–486, 2000.


[a,b]Phuong Nguyen

# Performance Optimization of Baryon-block Construction in the Stochastic LapH
Method

Ben Hörz

###### Abstract

Implementations of measurement kernels in high-level Lattice QCD frameworks
enable rapid prototyping, but can leave hardware capabilities significantly
underutilized. This is an acceptable tradeoff if the time spent in unoptimized
routines is generally small. The computational cost of modern spectroscopy
projects however can be comparable to or even exceed the cost of generating
gauge configurations and computing solutions of the Dirac equation. One such
key kernel in the stochastic LapH method is the computation of baryon blocks;
we discuss several implementation strategies and achieve a 7.2x speedup over
the current implementation on a system with Intel® Xeon® Platinum 8358
processors, formerly Ice Lake.

##  1 Introduction

With the advent of modern spectroscopy methods for multi-hadron systems [1, 2,
3], ever more complicated physical systems are coming into reach. While a
variety of two-meson systems including with several coupled channels have been
investigated (see [4] for a plenary review at this conference), and more
recently even studies of three-meson systems have started appearing [5, 6, 7,
8, 9, 10, 11, 12] (also subject of a topical review this year [13]), the
situation for systems with baryons is comparably less advanced. Hadron
interactions involving baryons are however among the fundamental building
blocks for an understanding of nuclear physics rooted in QCD (for a review,
see for instance [14]). Examples include two-nucleon interactions, which –
when determined at sufficiently light pion mass – can serve as a validation
system to establish the reliability of these kinds of lattice QCD calculations
[15, 16, 17, 18]; three-nucleon interactions, which are difficult to access
phenomenologically and hence present a great opportunity for lattice QCD to
provide useful data; a variety of meson-baryon systems, for instance nucleon-
pion scattering in the isospin $I=3/2$ channel featuring the $\Delta(1232)$
resonance [19, 20, 21].

Even though modern spectroscopy methods present no conceptual difficulties
generalizing to systems involving baryons, there are a few challenges in
practice. Baryonic systems typically suffer from a worse signal-to-noise ratio
than purely mesonic systems, requiring larger amounts of statistics to obtain
meaningful results. In addition to necessitating more correlation-function
samples, every individual sample tends to be more computationally expensive
compared to the mesonic sector due to the increased number of quark fields. At
the correlator-construction level, algorithms eliminating redundant
computations have been devised to alleviate the proliferation of Wick
contractions [22, 23, 24, 5].

This work is concerned with improving the efficiency of the computation of
baryon functions in the stochastic LapH method [2]. In the stochastic LapH
framework, baryon blocks are rank-three tensors in dilution indices, carrying
additional labels identifying the baryon operator (flavor, spin, hadron
momentum) as well as the three noises used for the stochastic estimate of the
quark propagators in the LapH subspace. Correlators are then computed through
tensor contractions over dilution indices of those baryon functions as
governed by Wick’s theorem. This beneficial property of the stochastic LapH
method – the evaluation of complicated multi-hadron correlation functions is
reduced to tensor contractions involving blocks representing the constituent
hadrons – enables the re-use of baryon blocks for a wide variety of physical
systems. The optimizations presented in this work are hence immediately
applicable to a breadth of calculations involving baryons.

This contribution is organized as follows: section 2 defines the baryon-block
kernel and discusses its computational characteristics, section 3 contrasts
several implementation strategies, and benchmark results are presented in
section 4.

##  2 Baryon blocks in the stochastic LapH method

The stochastic LapH method [2] is a stochastic variant of distillation [1]
which avoids the $V^{2}$ scaling with the spatial simulation volume $V$ by
stochastically estimating the quark propagator projected into the LapH (or
distillation) subspace. A useful quantity for the computation of multi-hadron
correlation functions involving baryons is the (single-site) baryon function
defined per time slice of the simulation volume,

$\displaystyle
B_{d_{1}d_{2}d_{3}}^{(\vec{p},\Lambda,\mu,\eta)}=c^{(\vec{p},\Lambda,\mu)}_{\alpha\beta\gamma}\sum_{\vec{x}}\mathrm{e}^{-\mathrm{i}\,\vec{p}\vec{x}}\epsilon_{abc}q^{(\eta_{1},d_{1})}_{\alpha
a\vec{x}}q^{\prime(\eta_{2},d_{2})}_{\beta
b\vec{x}}q^{\prime\prime(\eta_{3},d_{3})}_{\gamma c\vec{x}},$ (1)

with color indices $a$, $b$, $c$, and the summation runs over the sites of a
three-dimensional time slice of the lattice. The relevant combinations of spin
indices $\alpha$, $\beta$, $\gamma$ are selected according to the group-
theoretical projection coefficients
$c^{(\vec{p},\Lambda,\mu)}_{\alpha\beta\gamma}$ for a given hadron momentum
$\vec{p}$, irrep $\Lambda$ and irrep row $\mu$ [25]. The quark fields $q$,
$q^{\prime}$, $q^{\prime\prime}$, which have been projected into the LapH
subspace, are obtained by repeatedly solving the Dirac equation with diluted
stochastic sources identified by a noise label $\eta$ and dilution index
$d_{1}=1,\dots,N_{\mathrm{dil}}$. Crucially, those solutions of the Dirac
equation for a given noise need only be computed once, and can then be cheaply
stored on disk and used for various multi-hadron projects by reconstructing
the smeared quark fields from their coefficients $Q$ in the basis of
eigenvectors of the three-dimensional gauge-covariant Laplacian $\phi^{(l)}$,
$l=1,\dots,N_{\mathrm{ev}}$, which is used to define the LapH subspace,

$\displaystyle q^{(\eta,d)}_{\alpha
a\vec{x}}=\sum_{l=1}^{N_{\mathrm{ev}}}Q^{(\eta,d)}_{\alpha
l}\phi^{(l)}_{a\vec{x}},$ (2)

and similarly for $q^{\prime}$ and $q^{\prime\prime}$, which differ only in
their coefficients $Q$, but with the same eigenvectors $\phi$.

Both the group-theoretical projection involving the spin indices and the
bookkeeping of noise indices in (1) are handled by the calling application,
leaving

$\displaystyle
B_{d_{1}d_{2}d_{3}}^{(\vec{p})}=\sum_{\vec{x}}\mathrm{e}^{-\mathrm{i}\,\vec{p}\vec{x}}\epsilon_{abc}Q^{(1)}_{d_{1}l_{1}}Q^{(2)}_{d_{2}l_{2}}Q^{(3)}_{d_{3}l_{3}}\phi^{(l_{1})}_{a\vec{x}}\phi^{(l_{2})}_{b\vec{x}}\phi^{(l_{3})}_{c\vec{x}},$
(3)

where summation over the eigenvector indices $l_{1}$, $l_{2}$ and $l_{3}$ is
implied, as the computational kernel, which is called many times for different
noise and spin combinations, i.e. different quark field coefficients, but with
the same momentum set and Laplacian eigenvectors.

Depending on the sizes of $N_{\mathrm{ev}}$ and $N_{\mathrm{dil}}$, as well as
how many times the kernel (3) is called per time slice, one of the following
two different approaches is preferable.

For moderate values of $N_{\mathrm{ev}}$, (3) can be efficiently computed
using a two-step procedure: During the setup stage, the mode-triplets

$\displaystyle
T_{l_{1}l_{2}l_{3}}^{\vec{p}}=\sum_{\vec{x}}\mathrm{e}^{-\mathrm{i}\,\vec{p}\vec{x}}\epsilon_{abc}\phi^{(l_{1})}_{a\vec{x}}\phi^{(l_{2})}_{b\vec{x}}\phi^{(l_{3})}_{c\vec{x}},$
(4)

which are spin-, noise- and flavor-blind, are computed and kept in memory. All
lattice-sized objects in (3) have then been consumed, and the baryon function
for a given set of quark-field coefficients can be computed by tensor-
contracting them onto the precomputed mode-triplet. Those tensor contractions
can be performed with high performance, so the majority of the runtime tends
to be associated with the initial setup phase, which needs to be amortized
over many kernel invocations. The major drawback of this mode-triplet approach
is the need to keep one $N_{\mathrm{ev}}^{3}$-sized object per momentum in
memory111Based on the symmetries of (4), only ${N_{\mathrm{ev}}\choose 3}$
elements of a mode-triplet are independent. Exploting that symmetry with a
sparse storage scheme however complicates the subsequent tensor contractions
of quark-field coefficients onto the mode-triplet..

For large number of eigenvectors $N_{\mathrm{ev}}$, a more economical approach
is to first reconstruct the quark fields from the coefficients as per (2), and
subsequently perform the reduction (3) over sets of lattice-sized objects. The
quark-field reconstruction can be efficiently implemented using matrix-matrix
multiplication and reduces the complexity of subsequent lattice-sized
reductions to $N_{\mathrm{dil}}^{3}$ (rather than $N_{\mathrm{ev}}^{3}$ for
the mode-triplet approach), which however must be performed for every kernel
invocation.

In view of the requirements of baryon calculations in large volumes – such as
the E250 ensemble [26] generated by the CLS effort [27, 28], where employing
the mode-triplet approach is not feasible222First results in the mesonic
sector presented in [29] used the analogous mode-doublet approach for meson
construction, which is still affordable due to its slightly weaker
$N_{\mathrm{ev}}^{2}$ scaling. For the $N_{\mathrm{ev}}=1536$ employed in that
work, the mode triplet on the other hand occupies $1.8\,\mathrm{TB}$ of memory
already for the moderate number of momenta $N_{\mathrm{mom}}=33$, clearly
making this approach impractical. – the goal of this work is to provide an
efficient implementation of (3), which utilizes the great compute capabilities
of modern hardware by exploiting the $N_{\mathrm{dil}}^{3}$ compute complexity
with only linear-in-$N_{\mathrm{dil}}$ memory traffic.

##  3 Implementation details

Algorithm 1 Cache blocking algorithm with pseudo code

1:$q_{1},q_{2},q_{3},phase$

2:$baryon$

3:

4:! $Block{D_{i}}\leftarrow N_{D_{i}}\,/\,Bsize{D_{i}}$ ($i=1,2,3$)

5:! $BlockX\leftarrow N_{X}\,/\,Bsize{X}$

6:function BaryonConstruct($q_{1},q_{2},q_{3},phase$)

7: for $Block{D_{1}}$ do in parallel

8: for $Block{D_{2}}$ do in parallel

9: tmpBuf $\leftarrow 0.$

10: for each $Block{X}$ do

11: for $d_{1}\leftarrow 1$ to $Bsize{D_{1}}$ do

12: for $d_{2}\leftarrow 1$ to $Bsize{D_{2}}$ do

13: for $x\leftarrow 1$ to $Bsize{X}$ do

14: $diq(d_{1},d_{2},:,x)\leftarrow q_{1}(\tilde{d_{1}},:,x)\times
q_{2}(\tilde{d_{2}},:,x)$

15: end for

16: end for

17: end for

18: for each $Block{D_{3}}$ do

19: for $d_{3}\leftarrow 1$ to $Bsize{D_{3}}$ do

20: for each $diq_{i}$ in $diq$ do

21: for $x\leftarrow 1$ to $Bsize{X}$ do

22: $singlet(d_{1},d_{2},d_{3},x)\leftarrow diq_{i}\times
q_{3}(\tilde{d_{3}},:,x)$space

23: end for

24: end for

25: end for

26: tmpBuf(:) $\leftarrow$ tmpBuf(:) \+ $singlet\times phase(:)$
$\triangleright$ Intel® MKL JIT GEMM

27: end for

28: end for

29: $baryon\leftarrow$ tmpBuf

30: end for parallel

31: end for parallel

32: return $baryon$

33:end function

Cache Blocking in $N_{D_{1}}$, $N_{D_{2}}$, $N_{D_{3}}$, $N_{X}$

The quark-field reconstruction can be performed efficiently using matrix-
matrix multiplication, for which highly optimized implementations are
available for all hardware architectures. Hence, in the following section, we
focus on optimizing the baryon-block calculation given the reconstructed quark
fields $q_{1}$, $q_{2}$, $q_{3}$.

The optimized algorithm is shown in Algorithm 1. Typically, the number of
requested hadron momenta $N_{\mathrm{mom}}$ is much smaller than the number of
allowed momenta (e.g. $33\ll 64^{3}$); therefore using a fast Fourier
transform is not beneficial. Hence, the phase factor
$e^{-\mathrm{i}\vec{p}\vec{x}}$ for the momentum projection in (3) can be
precomputed and re-used for several kernel invocations.

Cache blocking techniques are employed in conjunction with an appropriate data
layout to optimize data locality. Blocking is implemented both in the spatial
indices $x$ and the three dilution indices $d_{1},d_{2},d_{3}$. The blocking
in $x$ allows the kernel to exploit the available inherent input reuse. For
example, each block of $q_{1}$ input can be kept in the cache and re-used for
different $diq$ calculations with different $q_{2}$ since the input size is
small enough to stay in the cache (ll. 8-14). The blocking in
$d_{1},d_{2},d_{3}$ enables the kernel to keep the intermediate data (diq,
singlet, tmpBuf) in cache and use it for subsequent calculations (ll. 16-22).
A suitable data memory layout is (in row-major convention)
$N_{\mathrm{dil}}\times N_{\mathrm{BlockX}}\times N_{\mathrm{color}}\times
N_{\mathrm{BsizeX}}$ for the input $q_{1},q_{2},q_{3}$, ensuring that the data
is accessed contiguously in $x$ for each each color component. Furthermore,
this data layout stores the BlockX-sized chunks for the three colors of
$q_{1},q_{2},q_{3}$ adjacently, enhancing spatial locality in the calculations
of diq and singlet.

The small matrix-matrix multiplication in l. 23 to compute
$\textit{tmpBuf(:)}+=singlet\times phase(:)$, utilizes the Intel® Math Kernel
Library (Intel® MKL) with just-in-time (JIT) code generation for small
matrices with
$m=\mathrm{BsizeD}_{1}\times\mathrm{BsizeD}_{2}\times\mathrm{BsizeD}_{3}$,
$n=N_{\mathrm{mom}}$ and $k=\mathrm{BsizeX}$. These matrices – $singlet$ and
tmpBuf – should be sufficiently small to remain in the cache. General-purpose
GEMM implementations are typically optimized targeting larger matrix sizes.
Thus, for the small matrix-matrix multiplication required here, the Intel®
Math Kernel Library (Intel® MKL) with JIT compilation is used to generate
target microarchitecture code for the kernel which is optimized for small-
sized complex-valued matrix multiplication problems. As the matrix sizes are
fixed, the kernel can be produced once and then called many times, amortizing
the JIT compilation overhead.

Parallelization for multiple threads is achieved by distributing the work in
the loops over dilution-index blocks, (ll. 4-5) . In practice, the loops over
$\mathrm{BlockD}_{1}$ and $\mathrm{BlockD}_{2}$ are collapsed into a joint
iteration space and parallelized for multiple threads with OpenMP. In this
approach, there are no data dependencies since computations of each thread are
fully independent. Furthermore, as each block works on an identical amount of
data, the load is expected to be well-balanced between threads. Lastly,
parallelization is implemented at an outer level, encompassing plenty of work
per loop trip to amortize the OpenMP runtime overhead for instance for thread
scheduling.

##  4 Performance results

Figure 1: Performance of the previous and optimized kernel normalized to the
single-thread performance of the previous implementation. The optimized kernel
outperforms the baseline by up to 6.8x for the test with fixed frequency
(left) and by up to 7.2x for the test when turbo boost is enabled (right).

The implementation of Algorithm 1 is evaluated on a test system with two
Intel® Xeon® Platinum 8358 processors @ 2.60 GHz for a moderately large
problem size of $L=64$ with $N_{\mathrm{dil}}=64$ dilution indices per quark
field and number of requested momenta $n_{\mathrm{mom}}=33$. Figure 1 shows
the performance of the optimized kernel compared to the previous
implementation using 1 to 64 cores.

The block sizes are tunable parameters which generally depend on the target
architecture and should be tuned experimentally. For our test node, the block
sizes yielding the best performance are $\mathrm{BsizeX}=32$,
$\mathrm{BsizeD}_{1}=4$, $\mathrm{BsizeD}_{2}=8$, $\mathrm{BsizeD}_{3}=16$.
This can be understood as aligning the memory footprint of each single
intermediate object with the available cache hierarchy. For instance, the
$diq$ array stores $\mathrm{BsizeD}_{1}\times\mathrm{BsizeD}_{2}\times
N_{\mathrm{colors}}\times\mathrm{BsizeX}=3072$ complex double-precision
values, occupying 48 kB of memory, thus fitting perfectly into the L1 cache.
Similarly, the $singlet$ has a size of 256 kB which fits into the L2 cache.

When run on a single core, the optimized kernel is 1.8x and 2.5x faster than
the baseline in tests with and without turbo boost, respectively. With in-
depth profiling, the superior performance of the optimized implementation can
be traced back to better use of the memory system, as expected. The optimized
kernel operates at 1.6x and 3x higher arithmetic intensity in Read and Write,
respectively. In addition, the memory access patterns are also significantly
improved such that fewer cycles are spent on load and store operations in the
optimized kernel in comparison to the baseline (42% reduction in loads and 70%
reduction in stores). The optimized kernel is indicated to be L1-bound instead
of DRAM-bound. As a result, with a performance of 37.7 DP GFlops/s, the kernel
reaches 54% of the theoretical peak performance.

Figure 2: Strong-scaling behavior of the optimized kernel. Left: With fixed
frequency, the kernel achieves almost perfect parallel efficiency for a 62.5x
speedup with 64 threads. Right: With turbo enabled, the kernel scales well up
to eight cores. Beyond eight cores, the variable frequency gets throttled and
the decrease in parallel efficiency matches the decrease in frequency,
implying that the optimized kernel is compute-bound.

The improvement of the optimized kernel becomes particularly apparent in
multithreaded runs, where it outperforms the baseline by 1.8x to 6.8x at fixed
frequency and 2.5x to 7.2x for tests with turbo boost on (Figure 1). This
boost in performance can be understood as being due to improved temporal
locality. As long as threads progress at a comparable rate, adjacent threads
working on a different BlockD2 but sharing the same BlockD1 access the same
input data, which may be served from the cache hierarchy.

Figure 2 shows the strong-scaling behavior of the optimized kernel at fixed
clock frequency, achieving almost perfect parallel efficiency at a speedup of
62.5x when using 64 cores. In addition, the parallel efficiency is above 0.98
for all tests, implying that the kernel scales almost perfectly within a node.

For production runs with turbo boost enabled the clock frequency can increase
up to the boost frequency (3.6 GHz on our test system). While performance in
absolute terms is slightly better than at fixed frequency, the strong-scaling
parallel efficiency deteriorates, showing a 39.6x speedup with 64 cores
(Figure 2). The loss in parallel efficiency is due to frequency throttling.
While for one thread the average frequency is $3.285\,\mathrm{GHz}$, it drops
to $2.219\,\mathrm{GHz}$ when running with 32 threads. The frequency ratio
$0.67$ matches the parallel efficiency, indicating that the optimized kernel
is compute-bound and indeed limited by the frequency throttling.

Figure 3: Scaling performance of the optimized kernel with respect to the
problem size. The kernel scales as expected to within a few percent both with
the spatial volume and the number of dilution indices.

The scalability of the optimized kernel with respect to the problem size, of
importance in view of ever-increasing simulation volumes, is shown in Figure
3. As is evident from (3), the computational cost scales as
$\mathcal{O}(L^{3})$ and $\mathcal{O}(N_{\mathrm{dil}}^{3})$ with the spatial
volume and number of dilution indices, respectively. The optimized kernel
scales as expected to within a few percent as a function of the problem size
both with the spatial volume as well as the number of dilution indices.

##  5 Summary

We have presented an optimized implementation of the kernel computing baryon
blocks in the stochastic LapH method, achieving an up to 7.2x speedup over the
previos implementation. Exploiting the high arithmetic intensity of (3) by
blocking in dilution and spatial indices, and performing the momentum
projection with a JIT-compiled microkernel provided by the Intel® Math Kernel
Library (Intel® MKL), we achieve good single-core performance on a test system
with Intel® Xeon® Platinum 8358 processors. Parallelizing over blocks of
dilution indices using multithreading, we also observe good scalability all
the way to the maximum number of cores per socket333Depending on memory
requirements, this parallelization can straightforwardly be supplemented by a
standard domain decomposition over MPI ranks in order to enable scaling beyond
one socket or node at the expense of an $N_{\mathrm{dil}}^{3}$-sized MPI
reduction at the end of kernel execution.

This optimized implementation has been upstreamed into the chroma_laph
measurement suite and is ready for use in production runs on large lattice
volumes.

## Acknowledgments

We acknowledge useful conversations with Colin Morningstar as well as our
Intel colleagues Christoph Bauinger and Carsten Uphoff. Intel and the Intel
logo are trademarks of Intel Corporation or its subsidiaries.

## References

  * [1] Hadron Spectrum collaboration, _A Novel quark-field creation operator construction for hadronic physics in lattice QCD_ , _Phys. Rev. D_ 80 (2009) 054506 [0905.2160].
  * [2] C. Morningstar, J. Bulava, J. Foley, K.J. Juge, D. Lenkner, M. Peardon et al., _Improved stochastic estimation of quark propagation with Laplacian Heaviside smearing in lattice QCD_ , _Phys. Rev. D_ 83 (2011) 114505 [1104.3870].
  * [3] W. Detmold, D.J. Murphy, A.V. Pochinsky, M.J. Savage, P.E. Shanahan and M.L. Wagman, _Sparsening algorithm for multihadron lattice QCD correlation functions_ , _Phys. Rev. D_ 104 (2021) 034502 [1908.07050].
  * [4] L. Liu, _Hadron Spectroscopy and Interactions from Lattice QCD_ , _PoS_ LATTICE2022 (2022) 448.
  * [5] B. Hörz and A. Hanlon, _Two- and three-pion finite-volume spectra at maximal isospin from lattice QCD_ , _Phys. Rev. Lett._ 123 (2019) 142002 [1905.04277].
  * [6] T.D. Blanton, F. Romero-López and S.R. Sharpe, _$I=3$ Three-Pion Scattering Amplitude from Lattice QCD_, _Phys. Rev. Lett._ 124 (2020) 032001 [1909.02973].
  * [7] M. Mai, M. Döring, C. Culver and A. Alexandru, _Three-body unitarity versus finite-volume $\pi^{+}\pi^{+}\pi^{+}$ spectrum from lattice QCD_, _Phys. Rev. D_ 101 (2020) 054510 [1909.05749].
  * [8] C. Culver, M. Mai, R. Brett, A. Alexandru and M. Döring, _Three pion spectrum in the $I=3$ channel from lattice QCD_, _Phys. Rev. D_ 101 (2020) 114507 [1911.09047].
  * [9] M. Fischer, B. Kostrzewa, L. Liu, F. Romero-López, M. Ueding and C. Urbach, _Scattering of two and three physical pions at maximal isospin from lattice QCD_ , _Eur. Phys. J. C_ 81 (2021) 436 [2008.03035].
  * [10] Hadron Spectrum collaboration, _Energy-Dependent $\pi^{+}\pi^{+}\pi^{+}$ Scattering Amplitude from QCD_, _Phys. Rev. Lett._ 126 (2021) 012001 [2009.04931].
  * [11] A. Alexandru, R. Brett, C. Culver, M. Döring, D. Guo, F.X. Lee et al., _Finite-volume energy spectrum of the $K^{-}K^{-}K^{-}$ system_, _Phys. Rev. D_ 102 (2020) 114523 [2009.12358].
  * [12] T.D. Blanton, A.D. Hanlon, B. Hörz, C. Morningstar, F. Romero-López and S.R. Sharpe, _Interactions of two and three mesons including higher partial waves from lattice QCD_ , _JHEP_ 10 (2021) 023 [2106.05590].
  * [13] F. Romero-López, _Topical plenary on multi-particle interactions from lattice QCD_ , _PoS_ LATTICE2022 (2022) 440.
  * [14] C. Drischler, W. Haxton, K. McElvain, E. Mereghetti, A. Nicholson, P. Vranas et al., _Towards grounding nuclear physics in QCD_ , _Prog. Part. Nucl. Phys._ 121 (2021) 103888 [1910.07961].
  * [15] A. Francis, J.R. Green, P.M. Junnarkar, C. Miao, T.D. Rae and H. Wittig, _Lattice QCD study of the $H$ dibaryon using hexaquark and two-baryon interpolators_, _Phys. Rev. D_ 99 (2019) 074505 [1805.03966].
  * [16] J.R. Green, A.D. Hanlon, P.M. Junnarkar and H. Wittig, _Weakly bound $H$ dibaryon from SU(3)-flavor-symmetric QCD_, 2103.01054.
  * [17] B. Hörz et al., _Two-nucleon S-wave interactions at the $SU(3)$ flavor-symmetric point with $m_{ud}\simeq m_{s}^{\rm phys}$: A first lattice QCD calculation with the stochastic Laplacian Heaviside method_, _Phys. Rev. C_ 103 (2021) 014003 [2009.11825].
  * [18] S. Amarasinghe, R. Baghdadi, Z. Davoudi, W. Detmold, M. Illa, A. Parreno et al., _A variational study of two-nucleon systems with lattice QCD_ , 2108.10835.
  * [19] C.W. Andersen, J. Bulava, B. Hörz and C. Morningstar, _Elastic $I=3/2p$-wave nucleon-pion scattering amplitude and the $\Delta$(1232) resonance from Nf=2+1 lattice QCD_, _Phys. Rev. D_ 97 (2018) 014506 [1710.01557].
  * [20] G. Silvi et al., _$P$ -wave nucleon-pion scattering amplitude in the $\Delta$(1232) channel from lattice QCD_, _Phys. Rev. D_ 103 (2021) 094508 [2101.00689].
  * [21] J. Bulava, A.D. Hanlon, B. Hörz, C. Morningstar, A. Nicholson, F. Romero-López et al., _Elastic nucleon-pion scattering at $m_{\pi}=200~{}{\rm MeV}$ from lattice QCD_, 2208.03867.
  * [22] T. Doi and M.G. Endres, _Unified contraction algorithm for multi-baryon correlators on the lattice_ , _Comput. Phys. Commun._ 184 (2013) 117 [1205.0585].
  * [23] W. Detmold and K. Orginos, _Nuclear correlation functions in lattice QCD_ , _Phys. Rev. D_ 87 (2013) 114512 [1207.1452].
  * [24] J. Günther, B.C. Toth and L. Varnhorst, _Recursive approach to determine correlation functions in multibaryon systems_ , _Phys. Rev. D_ 87 (2013) 094513 [1301.4895].
  * [25] C. Morningstar, J. Bulava, B. Fahy, J. Foley, Y.C. Jhang, K.J. Juge et al., _Extended hadron and two-hadron operators of definite momentum for spectrum calculations in lattice QCD_ , _Phys. Rev. D_ 88 (2013) 014511 [1303.6816].
  * [26] D. Mohler, S. Schaefer and J. Simeth, _CLS 2+1 flavor simulations at physical light- and strange-quark masses_ , _EPJ Web Conf._ 175 (2018) 02010 [1712.04884].
  * [27] M. Bruno et al., _Simulation of QCD with N ${}_{f}=$ 2 $+$ 1 flavors of non-perturbatively improved Wilson fermions_, _JHEP_ 02 (2015) 043 [1411.3982].
  * [28] M. Bruno, T. Korzec and S. Schaefer, _Setting the scale for the CLS $2+1$ flavor ensembles_, _Phys. Rev. D_ 95 (2017) 074504 [1608.08900].
  * [29] S. Paul, A.D. Hanlon, B. Hörz, D. Mohler, C. Morningstar and H. Wittig, _I=1 $\pi$-$\pi$ scattering at the physical point_, _PoS_ LATTICE2021 (2022) 551 [2112.07385].