#1 1

0 		 A Novel Sparse Bayesian Learning and Its Application to Fault Diagnosis for Multistation Assembly Systems
Jihoon Chunga, Bo Shenb, and Zhenyu (James) Konga

aGrado Department of Industrial and Systems Engineering, Virginia Tech, Blacksburg, US
bDepartment of Mechanical and Industrial Engineering,
New Jersey Institute of Technology, Newark, US

Sparse Bayesian Learning with Temporally Correlated  Solution Vectors and Prior Knowledge with Application to Fault Diagnosis in Multistation Assembly Systems
Author information is purposely removed for double-blind review

IISE Transactions  Template

This paper addresses the problem of fault diagnosis in multistation assembly systems. Fault diagnosis is to identify process faults that cause the excessive dimensional variation of the product using dimensional measurements. For such problems, the challenge is solving an underdetermined system caused by a common phenomenon in practice; namely, the number of measurements is less than that of the process errors. To address this challenge, this paper attempts to solve the following two problems: (1) how to utilize the temporal correlation in the time series data of each process error and (2) how to apply prior knowledge regarding which process errors are more likely to be process faults. A novel sparse Bayesian learning method is proposed to achieve the above objectives. The method consists of three hierarchical layers. The first layer has parameterized prior distribution that exploits the temporal correlation of each process error. Furthermore, the second and third layers achieve the prior distribution representing the prior knowledge of process faults. Then, these prior distributions are updated with the likelihood function of the measurement samples from the process, resulting in the accurate posterior distribution of process faults from an underdetermined system. Since posterior distributions of process faults are intractable, this paper derives approximate posterior distributions via Variational Bayes inference. Numerical and simulation case studies using an actual autobody assembly process are performed to demonstrate the effectiveness of the proposed method.

Keywords: Fault Diagnosis in Multistation Assembly Systems, Sparse Bayesian Learning, Temporal Correlation in Process Error, Prior Knowledge of Process Faults, Variational Bayes Inference.

$\lVert \textbf{x} \rVert_{p} $$\ell_{p}$ norm of the vector x.
$\lVert \text{A} \rVert_{F} $Frobenius norm of the matrix A.
$ \textbf{I}_{L}$Identity matrix with the size $L \times L$.

$ \text{A} \otimes \text{B}$Kronecker product of the matrices A and B.
$ \text{Vec(A)}$Vectorization of the matrix A.
$ \text{Tr(A)}$Trace of the matrix A.

$ \text{diag}\{ \alpha_{1},...,\alpha_{m} \}$Diagnoal matrix with principal diagonal elements being $\alpha_{1},...,\alpha_{m}$.

§ INTRODUCTION

Multistation assemblies represent the systems that perform assembly operations on multiple stations to assemble a final product. The quality of the final product relies on several factors, called key control characteristics (KCCs) or process errors [2]. The positioning accuracy of fixture locators is KCC in the multistation assembly [3]. Fixture locators clamp the parts during the assembly process, so their deviation from nominal may cause dimensional quality problems in the final product. Therefore, it is necessary to identify process errors that have mean shifts and/or large variance increases from their design specifications, namely, process faults [4]. Hence, fault diagnosis in multistation assemblies estimates the mean and variance of process errors, namely, the variations of fixture locators. This paper focuses on the process fault of the mean shift.

Because of the unduly cost and physical constraints, the dimensional variation of KCCs, namely, fixture locators, cannot be directly monitored using sensors in multistation assemblies [26]. Instead, key product characteristics (KPCs), namely, the key measurements from the final product, can be used to estimate the KCCs and, consequently, identify process faults among process errors. A fault-quality linear model of multistation process represents the relationship between KPCs and process errors as follows [21, 20]:

\begin{equation}\label{eq:1}
\textbf{y}=\Phi \textbf{x} +\textbf{v},
\end{equation}

where $\textbf{y}\in \mathbb{R}^{M \times 1} $represents $M$ dimensional measurements (i.e., KPCs), $ \textbf{x}\in \mathbb{R}^{N \times 1} $denotes $N$ process errors (i.e., KCCs), $\Phi \in \mathbb{R}^{M\times N}$ is a fault pattern matrix obtained from all the process information in the multistage process [13], and $\textbf{v} \in \mathbb{R}^{M \times 1}$ denotes the noise. Although the relationship between KPCs and KCCs is nonlinear in many manufacturing processes, a fault-quality linear model is used in general as the nonlinear relationship could be approximated to a linear model utilizing Taylor series expansion because the variations of KCCs are small, and the relationship between KCCs and KPCs is smooth and without sharp changes [26, 33, 12]. Since process errors indicate the mean deviations of the fixture locators, process faults refer to nonzero elements in $\textbf{x}$, namely, process errors with nonzero mean shifts [2].

Following the fault-quality model in Eq. (<ref>), some research on fault diagnosis for manufacturing systems has been investigated [12, 43]. All these methods assume that the number of measurements ($M$) is greater than the number of process errors ($N$) (i.e., $M>N$). However, this assumption may not always hold in actual manufacturing applications, as using an excessive number of sensors (measurements) will result in undue costs [2]. These approaches are unsuitable if this assumption is violated. This is because Eq. (<ref>) becomes an underdetermined system that results in the non-existence of a unique solution. To overcome the challenge, the sparse solution assumption [14] that x in Eq. (<ref>) has a minimal number of nonzero elements is required. In the context of the fault diagnosis problem, it denotes the sparsity of process faults in the fault-quality linear model. This is reasonable since it is likely to have a few process faults in practice [3].

Among the several sparse estimation methods, the Bayesian method called sparse Bayesian learning has received much attention recently because of its superior estimation performance [42].

Several studies have used sparse Bayesian learning for fault diagnosis in manufacturing systems [26, 2, 3, 28]. These studies successfully identified process faults by providing prior distribution of process errors (i.e., x in Eq. (<ref>)) to promote the sparsity of process faults. Especially, the work in [2] applied Bayesian learning to diagnose mean shift fault, which is the most relevant work to our study.  In [2], multiple KPCs samples have been used to estimate the mean deviation of fixture locators. Specifically, given the average of multiple KPCs samples, which is $\textbf{y}$ in Eq. (<ref>), the mean deviation of the fixture locators (namely, $\textbf{x}$ in Eq. (<ref>)) can be estimated. However, this work does not consider the characteristics of multiple KPCs samples where they are collected sequentially. In practice, time series data of each process error from the multiple KPCs samples may have a strong temporal correlation in the multistation assembly process. For example, the locator position of the fixture system could have a drifting due to wear [43], causing dimensional quality issues in the samples. If the process faults are not mitigated immediately, the effect of fixture deviation is auto-correlated in terms of time due to the degradation of wear of production tooling over time [33]. Similarly, the dimensional variability caused by machine-tool thermal distortions in the multistation processes is highly correlated between the samples assembled in a certain period [1]. Therefore, the temporal correlation of each process error causes product samples manufactured over a period of time to exhibit the same patterns of faults by a specific source of variation [29]. The relationship between the sequentially collected multiple KPCs samples and their process errors can be formulated as the following multiple measurements vectors (MMV) model [10] that extends from Eq. (<ref>):

\begin{equation}\label{eq:3}
\textbf{Y}=\Phi \textbf{X} +\textbf{V},
\end{equation}

where $\textbf{Y}=[\textbf{Y}_{\cdot 1},...,\textbf{Y}_{\cdot L}] \in \mathbb{R}^{M\times L}$ is a measurement matrix consisting of $L$  KPCs samples, and $\textbf{Y}_{\cdot i} \in \mathbb{R}^{M \times 1}$ is a vector that denotes the $i^{th}$ KPCs sample. $\textbf{X}=[\textbf{X}_{\cdot 1},...,\textbf{X}_{\cdot L}] \in \mathbb{R}^{N\times L}$ is a process error matrix, where $\textbf{X}_{\cdot i} \in \mathbb{R}^{N \times 1}$ is a vector that represents process errors (KCCs) of the $i^{th}$ KPCs sample. $\textbf{V}\in\mathbb{R}^{M\times L}$ is a noise matrix. Since $L$ KPCs samples have the same process faults, all columns of $\textbf{X}$ share the index of nonzero rows called support. It is called a common support assumption in the MMV model [10].  In addition,  $L$ elements in the $j^{th}$ row of $\textbf{X}$, representing the time series data of the $j^{th}$ process error,  are highly correlated as nonzero values if the $j^{th}$ process error is process fault.  However, the dynamic changes of process faults due to the complexity of the manufacturing process can easily violate the common support assumption in a large number of KPCs samples.

Therefore, this paper focuses on the small number of KPCs samples ($L$) to satisfy the common support assumption. Utilizing a small number of KPCs samples is also efficient in fault diagnosis of process faults in the multistation assembly process for time and cost reduction.

Beyond the common support assumption and temporal correlation in the time series data of each process error, utilizing prior knowledge of process faults is an additional way to improve the identification of sparse process faults [26]. Specifically, the prior knowledge regarding which process errors are more likely to be process faults than others. This knowledge can be obtained from domain-specific knowledge from practitioners or collected based on the fault diagnosis at the past time stamps. For example, the manufacturing engineers in an assembly line usually know that some fixture locators may malfunction more frequently than others based on their experiences. In practice, the prior knowledge provides only part of the actual process faults. In addition, the knowledge may contain some erroneous information as to process faults, which are actually not. Therefore, the prior knowledge provides partial and even erroneous information about process faults. However, utilizing the prior knowledge is expected to improve the identification of process faults if the correct and erroneous information can be properly distinguished.

In the sparse Bayesian learning  literature, there exist studies that consider the temporal correlation  of rows in matrix X in Eq. (<ref>) based on the common support assumption [30, 18], and the work that utilizes the partial with some erroneous prior knowledge of support to improve the performance of sparse estimation [16, 17], separately. However, these studies did not integrate these aspects to improve the sparse estimation.

To address this challenge, this paper aims to develop a novel sparse Bayesian hierarchical learning method that simultaneously utilizes the temporal correlation in the time series data of each process error, as well as prior knowledge, which may contain erroneous information. The new method, namely, support knowledge aided with temporally correlated process error sparse Bayesian learning (SA-TSBL), is proposed to achieve this objective. The contributions of this work are summarized as follows:

* From the methodological point of view, this paper proposes a novel sparse Bayesian learning that considers both temporal correlation in the time series data of each process error and prior knowledge of process faults to improve the sparse estimation. This method also derives an approximate posterior distribution of the sparse solution via Variational Bayes inference [32] to address the intractable computational challenge.

* From the application perspective, the proposed method is applied to fault diagnosis in the multistation assembly systems. The method mitigates the dimensional quality issues in the assembly operation by effectively identifying the fixture locators with excessive mean shifts. The effectiveness of the proposed method is validated in real-world simulation case studies that use an actual auto body assembly process.

The rest of this paper is organized as follows. A brief review of related research work is provided in Section <ref>. The proposed methodology is presented in Section <ref>, followed by numerical case studies to validate its effectiveness in Section <ref>. Section <ref> offers real-world case studies on fault diagnosis problems in the multistation assembly process. Finally, conclusions and future work are discussed in Section <ref>.

§ REVIEW OF RELATED WORK

The related existing studies of fault diagnosis in manufacturing systems are reviewed in Section <ref>. Then, the literature related to sparse Bayesian learning is provided in Section <ref>. Afterward, the research gaps in the current work are identified in Section <ref>.

§.§ Fault Diagnosis Methodologies in Manufacturing Systems

There has been a large body of research on fault diagnosis methodologies for manufacturing systems based on the fault-quality model in Eq. (<ref>). [23] developed a PCA-based orthogonal diagonalization strategy to transform the measurement data. It enabled the estimation of the variance of process errors in a multistation assembly system. [43] applied a mixed linear model to represent the relationship between measurements and process faults in a multistage machine process. Then, it used a maximum likelihood estimator for both mean shift and variance change to detect the process change. [7] presented a fault diagnosis method based on fault mapping procedures that combine principal component analysis (PCA) and pattern recognition.

All these approaches assume the number of measurements is greater than the number of the process errors (i.e., $M>N$ in Eq. (<ref>)), which may not be consistent with the industrial practice. However, if this assumption is violated, all of the approaches mentioned above are ineffective because the fault-quality linear model becomes an underdetermined system resulting in the non-existence of a unique solution.

To overcome an underdetermined system for fault diagnosis in the manufacturing system, sparse learning which has been actively used in fault diagnosis and detection in the manufacturing system can be considered. [19] developed a fault diagnosis method via dictionary learning and sparse representation-based classification for fault diagnosis in the manufacturing process.

[41] proposed a new non-convex penalty called the generalized logarithm penalty, which enables sparsity and reduces noise disturbance for bearing fault diagnosis.
[39] developed a dynamical multivariate functional data modeling via sparse subspace learning to detect process faults in the manufacturing processes. [11] proposed a group-sparsity learning approach for bearing fault diagnosis. In addition to fault diagnosis and detection in the general manufacturing process, sparse learning has been widely utilized to address the issue of an underdetermined system in the multistage assembly system. Specifically, sparse Bayesian learning has been actively utilized to incorporate the sparsity of process faults as the prior distribution. [3] proposed a fault diagnosis approach by combining the state-space model and the relevant vector machine to figure out process faults using the sparse estimate of the variance change of process errors. [28] developed a Bayesian variable selection-based method to identify both process faults and sensor faults in the assembly process. [2] proposed a spatially correlated sparse Bayesian learning to deal with the case when process errors have a spatial correlation. The work is based on the hypothesis that if one of the fixture locators deviates from its design specification, the neighboring locators are also expected to deviate. [26] presented a Bayesian approach for identifying variation sources in a multistage manufacturing process using the sparse variance component prior. The work focuses on the identification of process faults that have variance increases.

§.§ Sparse Bayesian Learning

Since [35] proposed Sparse Bayesian learning (SBL), many researchers have significantly extended it. [36] first introduced SBL to sparse estimation for the single measurement vector model in Eq. (<ref>). Then [37] extended it to the MMV model (Eq. (<ref>)), deriving the MSBL algorithm using the common support assumption. The advantage of SBL and MSBL is that global minimum of both methods are always the sparsest solution compared to the minimization based algorithms [34, 8], where global minimum is usually not the sparsest solution [42, 6]. In addition, sparse Bayesian learning have much fewer local minima than sparse learning with frequentist approaches such as the FOCUSS family [42, 10]. Based on the MMV model, several studies exploited the temporal correlation to improve the performance of sparse estimation. [42] presented a block SBL framework where a positive definite matrix captures the correlation structure of each row of X in Eq. (<ref>). [30] established a hierarchical Bayesian framework to model the temporally smooth signals using a multinomial distribution as the prior distribution of each row of X in Eq. (<ref>). [18] provided a Wishart distribution as the prior distribution to learn the temporal correlation from multiple measurement samples. Besides considering the temporal correlation, work utilizing prior knowledge of the support has been actively studied recently under the SBL framework. To utilize the prior knowledge in the Bayesian framework, [16, 27, 38] added one more hierarchical layer to conventional SBL to integrate the prior knowledge of support as the prior distribution.

Research Gap Analysis

The work summarized in Section <ref> investigates process faults in multistation assembly systems for quality assurance. Recently, the sparsity of process faults has been used via SBL to deal with the low dimensional measurements in fault diagnosis problems, which is common in industrial practice. However, there is a lack of efforts to identify process faults by considering temporal correlation in the time series data of each KCC and utilizing prior knowledge of process faults. Research efforts in Section <ref> introduce methods in SBL that consider the temporal correlation and prior knowledge of support. However, these two aspects were not integrated into the reported work. Therefore, this paper proposes a novel SBL method that simultaneously utilizes temporal correlation in the time series data of each KCC and prior knowledge of process faults to improve the identification of sparse process faults in multistation systems.

§ PROPOSED RESEARCH METHODOLOGY

This section proposes a novel sparse Bayesian hierarchical method: support knowledge aided temporally correlated process error SBL (SA-TSBL). The proposed SA-TSBL is described in Section <ref>, followed by Bayesian inference of the proposed method in Section <ref>.

§.§ Proposed Methodology

The proposed methodology is a sparse Bayesian hierarchical model using multiple KPCs samples that have the same process faults. The method considers the correlation in the time series data of each process error, and utilizes prior knowledge of process faults to improve the sparse estimation. To exploit the temporal correlation in the time series data of each process error, the proposed method transforms the MMV model in Eq. (<ref>) to the following block single measurement vector model [42].

\begin{equation}\label{eq:4}
\textbf{y}=\textbf{D} \textbf{x} +\textbf{v},
\end{equation}

where $\textbf{y}=$Vec$(\textbf{Y}^{\top})\in\mathbb{R}^{ML\times 1}, \textbf{D}=\Phi \otimes \textbf{I}_{L}, \textbf{x}=$Vec$(\textbf{X}^{\top})\in \mathbb{R}^{NL\times1}$. Assume noise vector v follows Gaussian distribution with zero mean and variance $\lambda$. Eq. (<ref>) can be rewritten as

\begin{equation}\label{eq:5}
\textbf{y}= [\phi_{1}\otimes\textbf{I}_{L},...,\phi_{N}\otimes\textbf{I}_{L}][\text{x}_{1}^{\top},...,\text{x}_{N}^{\top}]^{\top}+\textbf{v},
\end{equation}

where $\phi_{i}$ is the $i^{th}$ column in $\Phi$, and $\text{x}_{i}$ consists of the $i^{th}$ process errors of $L$ KPCs samples (i.e., $\text{x}_{i}^{\top}=(\tilde{x}_{i1}, \tilde{x}_{i2},...,\tilde{x}_{iL})$, where $\tilde{x}_{ij}$ denotes the element in the $i^{th}$ row and the $j^{th}$ column of matrix $\textbf{X}$ in Eq. (<ref>)). In other words, $\text{x}_{i}\in \mathbb{R}^{L \times 1}$ is the $i^{th}$ block of x in Eq. (<ref>), as illustrated in Figure <ref>.

Block structure of x.

Therefore, $K$ nonzero rows in X in Eq. (<ref>) are represented as $K$ nonzero blocks in x in the proposed method. It implies that there exist $K$ process errors that shifted from the design specification, namely, process faults. Since $\text{x}_{i}$ consists of the time series data of $i^{th}$ process error from $L$ KPCs samples, the proposed method exploits the correlation between the $L$ elements of $\text{x}_{i}$. Since noise vector v follows Gaussian distribution with zero mean and variance $\lambda$, the Gaussian likelihood is provided for the block model in Eq. (<ref>) as follows:

\begin{equation*}
p(\textbf{y}|\textbf{x};\lambda)\sim N(\textbf{Dx},\lambda\textbf{I}_{ML}).
\end{equation*}

The proposed method consists of the following three layers. The prior distribution in the first layer is provided to exploit the correlation of the time series data of each KCC. The second and third layers consist of prior distribution representing the prior knowledge of process faults among process errors. A graphical representation of the proposed method is shown in Figure <ref>.

Graphical representation of the proposed method. A circle indicates a random variable or hyperparameter that needs to be estimated. A dashed circle and square represent an observation and a constant, respectively.

In the first layer of the hierarchical model in Figure <ref>, the Gaussian distribution is provided as the prior distribution for process error (x) since it is the most common distribution considered for process errors in the literature and practice [26, 3]. The prior distribution is given by

\begin{equation*}
p(\textbf{x}|\boldsymbol{\alpha};\textbf{B}) \sim N(\textbf{0},\Sigma_{0}),
\end{equation*}

where $\Sigma_{0}$ is

\Sigma_{0}=\left[
\begin{array}{ccc}
\alpha_{1}^{-1}\text{B}_{1}^{-1}  & \cdots & 0 \\
\vdots & \ddots & \vdots \\
0 & \cdots & \alpha_{N}^{-1}\text{B}_{N}^{-1}
\end{array}
\right],

$\boldsymbol{\alpha}=\{ \alpha_{1},...,\alpha_{N}\}$, and $\textbf{B}=\{ \text{B}_1,...,\text{B}_N \}$. Since $\text{x}_{i}$ consists of the time series data of $i^{th}$ process error from $L$ KPCs samples, $\alpha_{i}^{-1}$ controls the sparsity of the $i^{th}$ process error ($\text{x}_{i}$). For example, when $\alpha_{i}^{-1}$ converges to zero, the associated block $\text{x}_{i}$ will be driven to zero [35]. $\text{B}_{i}$ is a positive definite matrix that captures the correlation of the time series data of $i^{th}$ process error. The proposed method assumes the independence between process errors (e.g., different fixture locators), which is a common assumption in the literature [26, 4, 28]. It results in the block diagonal matrix $\Sigma_{0}$.

The second layer in Figure <ref> specifies Gamma prior distributions over $\boldsymbol{\alpha}$. The distributions have an individual rate parameter $b_{i}$ for each parameter $\alpha_{i}$, namely,

\begin{equation}\label{eq:6}
\end{equation}

where $\boldsymbol{b}=\{b_{1},...,b_{N}\}$. Eq. (<ref>) effectively integrates the prior information of process faults into the SBL framework. Basically, $a$ and $b_{i}$ are set to be very small values (e.g., $10^{-4}$) to provide a large variance of a prior distribution over $\alpha_{i}$ in SBL [35]. It encourages the large values of $\alpha_{i}$, and promotes the $i^{th}$ process error as non-process faults (i.e., $\text{x}_{i}=\boldsymbol{0}$ in Eq. (<ref>)) [16]. Specifically, the marginal distribution of the time series data of $i^{th}$ process error, namely, $p(\text{x}_{i})=\int p(\text{x}_{i}|\alpha_{i}, \text{B})p(\alpha_{i}|a, b_{i})\text{d}\alpha_{i}$, follows multivariate student-t distribution which the probability is concentrated at zero [35]. Suppose set $P$ consists of the indexes of process errors that are likely to be process faults based on prior knowledge. Then, the corresponding rate parameters $\{b_{i}, i\in P \}$ are set as a relatively large value (e.g., 1) to provide a small variance of the prior distribution over $\alpha_{i}$ [16]. The prior distribution allows the small value of $\alpha_{i}$, and encourages the $i^{th}$ process error as process faults (i.e., $\text{x}_{i}\neq 0$ in Eq. (<ref>)). However, assigning a fixed value to $b_{i}$ in Eq. (<ref>) has limitations to deal with the situation when the set $P$ has some erroneous information of process faults.

To address this issue, the third layer in Figure <ref> assigns a prior distribution over $\{b_{i}, i\in P \}$ [27]. Gamma distribution is assigned because of the support of the variable. In addition, the distribution enables users to derive the closed form of the approximate posterior distribution of $\{b_{i}, i\in P \}$, which will be described in the following section. The prior distribution over $\{b_{i}, i\in P \}$ is provided as follows:

\begin{equation}\label{eq:7}
\text{Gamma}(b_{i}|p,q)=\Gamma(p)^{-1}q^{p}b_{i}^{p-1}e^{-qb_{i}}, i\in P.
\end{equation}

$p$ and $q$ in Eq. (<ref>) are specified to characterize prior distribution. These two values in the proposed method can be adjusted by the users to satisfy the following two conditions as follows.

* First, the values are set to provide a large mean value of the prior distribution of $\{b_{i}, i\in P \}$.

* Second, the values are determined to have a large variance in the prior distribution of $\{b_{i}, i\in P \}$.

The first condition encourages the $i^{th}$ process error as process faults based on prior knowledge. Compared to a small fixed value of $b_{i}=10^{-4}$ that promotes the large $\alpha_{i}$ in previous research [16, 35], the large mean value of the prior distribution of $\{b_{i}, i\in P \}$ lets the small $\alpha_{i}$. The small $\alpha_{i}$ provides a prior distribution that $i^{th}$ process error is likely to be process faults of the mean shift.

The second condition is designed to deal with the case when erroneous prior knowledge is provided. Assume incorrect prior knowledge that the $i^{th}$ process error may be a process fault is provided. If the variance of the prior distribution of $\{b_{i}, i\in P \}$ is small, the $i^{th}$ process error is likely to be misdiagnosed as process faults because of the large mean value of the prior distribution of $\{b_{i}, i\in P \}$ from the first condition. However, if the variance of the prior distribution of $\{b_{i}, i\in P \}$ is large, $\{b_{i}, i\in P \}$ is not significantly affected by the prior distribution. Instead, $\{b_{i}, i\in P \}$ is highly affected by the data itself, where the $i^{th}$ process error is not the process fault. This allows $\{b_{i}, i\in P \}$ to be learned as a small value from data and provides the sparsity prior to the $i^{th}$ process error [16].

For $\{b_{i}, i\in P^{c}\}$ in the proposed method, the parameters are set to a fixed small value of $10^{-4}$ as in the previous study [35].

§.§ Bayesian Inference of the Proposed Methodology

The proposed method in Section <ref> has several hidden variables that need to be estimated, that is, the process errors (x), the variable controlling the sparsity of process errors ($\boldsymbol{\alpha}$), and the variable used to provide prior knowledge of process faults ($b_{i}$, $\forall i\in P$). In addition, there are hyperparameters related to temporal correlations in the time series data of $i^{th}$ KCCs ($\text{B}_{i}$) and noise ($\lambda$), which also require estimation. To avoid too many parameters to being estimated, causing a challenging task in sparse estimation, the proposed method sets $\text{B}_{i}=\text{B}$ $(\forall i)$ [42].

However, the posterior distribution of hidden variables in the proposed method (Eq. (<ref>)) does not have a closed form because of the complexity of the proposed hierarchical model in Section <ref>. Specifically, the denominator in Eq. (<ref>) cannot be calculated as closed form. Let $b_{i}$ $(i\in P)$ as $\bar{b}$ for convenience.

\begin{align}\label{eq:80}
P(\textbf{x}, \boldsymbol{\alpha}, \bar{b} | \textbf{y}) & = \frac{P(\textbf{y}, \textbf{x}, \boldsymbol{\alpha},\bar{b})}{  \idotsint\limits_{\textbf{x}, \boldsymbol{\alpha}, \bar{b} } P(\textbf{y}, \textbf{x}, \boldsymbol{\alpha},\bar{b}) \, d\textbf{x} \, d\boldsymbol{\alpha} \, d\bar{b} }
\end{align}

To overcome this challenge, this paper derives approximate posterior distributions of hidden variables via Variational Bayes inference (VBI). Specifically, Variational Bayes Expectation Maximization (VBEM) [32] is utilized to estimate hidden variables and hyperparameters in the proposed method to identify mean shifts process faults. VBEM consists of (1) E-step: Variational Bayesian expectation step to estimate hidden variables $\textbf{x}, \boldsymbol{\alpha}, \bar{b}$ by approximating the posterior distribution of hidden variables; and (2) M-step: Variational Bayesian maximization step to update hyperparameters $\text{B}$ and $\lambda$ by maximizing the expected value of the logarithm of the complete likelihood.

Let $\boldsymbol{\theta}$ be a vector with all hidden variables in the proposed method (i.e., $\boldsymbol{\theta}= (\textbf{x}, \boldsymbol{\alpha} ,\bar{b}))$. VBI approximates the posterior distribution of $\boldsymbol{\theta}$, denoted as $q(\boldsymbol{\theta})$, by minimizing Kullback-Leibler (KL) divergence between $q(\boldsymbol{\theta})$ and the true posterior distribution, namely, $p(\boldsymbol{\theta}|\boldsymbol{y})$ (i.e., $D_{KL} (q(\boldsymbol{\theta})||p(\boldsymbol{\theta}|\boldsymbol{y})$)). $q(\boldsymbol{\theta})$ is factorized as

\begin{equation*}
\end{equation*}

by the mean-field approximation [9]. The approximate posterior distribution $q(\theta_{i})$, where $\theta_{i}$ is the $i^{th}$ element in the set $\boldsymbol{\theta}$ is derived as follows by minimizing the $D_{KL} (q(\boldsymbol{\theta})||p(\boldsymbol{\theta}|\boldsymbol{y})$ under the mean-field approximation.

\begin{equation}\label{eq:8}
\ln{q(\theta_{i})} = \mathbb{E}[\ln{p(\boldsymbol{y},\boldsymbol{\theta}})]_{\boldsymbol{\theta} \setminus \theta_{i}}+const,
\end{equation}

where $\mathbb{E}_{\boldsymbol{\theta} \setminus \theta_{i}}$ denotes the expectation taken with the set $\boldsymbol{\theta}$ without $\theta_{i}$. const can be obtained through normalization. Eq. (<ref>) is used in the following E-step of VBEM to approximate the posterior distributions of hidden variables.

E-step of VBEM: The posterior distributions of hidden variables that are related to process errors ($\textbf{x}$), sparsity of process errors ($\boldsymbol{\alpha}$), and prior knowledge of process faults ($\bar{b}$) are approximated by Eq. (<ref>), respectively, as follows.

\begin{align}\label{eq:9}
\ln{q(\textbf{x})} &= \langle \ln{p(\textbf{y},\textbf{x}, \boldsymbol{\alpha}, \bar{b})} \rangle_{q(\boldsymbol{\alpha})q(\bar{b})}+const  \nonumber \\
& =\langle \ln{p(\textbf{y}|\textbf{x};\lambda})p(\textbf{x}|\boldsymbol{\alpha};\textbf{B}) \rangle_{q(\boldsymbol{\alpha})}+const,
\end{align}

\begin{align}\label{eq:10}
\ln{q(\boldsymbol{\alpha})} &= \langle \ln{p(\textbf{y},\textbf{x}, \boldsymbol{\alpha}, \bar{b})} \rangle_{q(\textbf{x})q(\bar{b})}+const  \nonumber \\
& =\langle \ln{p(\textbf{x}|\boldsymbol{\alpha};\textbf{B}})p(\boldsymbol{\alpha}|\bar{b}) \rangle_{q(\boldsymbol{\alpha})}+const,
\end{align}

\begin{align}\label{eq:11}
\ln{q(\bar{b})} &= \langle \ln{p(\textbf{y},\textbf{x}, \boldsymbol{\alpha}, \bar{b})} \rangle_{q(\textbf{x})q(\boldsymbol{\alpha})}+const  \nonumber \\
& =\langle \ln{p(\boldsymbol{\alpha}|\bar{b}})p(\bar{b}) \rangle_{q(\boldsymbol{\alpha})}+const,
\end{align}

where $\langle \cdot \rangle$ indicates the expectation. Based on the statistical inference, the posterior distributions of hidden variables can be derived as

\begin{equation}\label{eq:12}
\end{equation}

\begin{equation}\label{eq:13}
\end{equation}

\begin{equation}\label{eq:14}
q(\bar{b})=\prod_{i\in P}\text{Gamma}(b_{i}|p+a,\langle \alpha_{i} \rangle+q).
\end{equation}

The expectations and moments of distributions in Eqs. (<ref>), (<ref>), and (<ref>) are

\begin{equation}\label{eq:15}
\mu_{\textbf{x}} = \frac{1}{\lambda}\Sigma_{\textbf{x}}\textbf{D}^{\top}\textbf{y},
\end{equation}

\begin{equation}\label{eq:16}
\langle \alpha_{i} \rangle = \frac{a+\frac{L}{2}}{\langle {b}_{i} \rangle}, \hspace{0.25em} (\forall i),
\end{equation}

\begin{equation}\label{eq:17}
\langle b_{i} \rangle = \frac{p+a}{q+\langle a_{i}\rangle },  \hspace{0.25em}(i\in P),
\end{equation}

\begin{equation}\label{eq:18}
\Sigma_{\textbf{x}} = (\frac{1}{\lambda}\textbf{D}^{\top}\textbf{D}+\langle \text{AB} \rangle)^{-1},
\end{equation}

where $\langle \text{AB} \rangle = \text{diag}[\langle \alpha_{1} \rangle \text{B},...,\langle \alpha_{N} \rangle\text{B}]$. The estimator of the time series data of $i^{th}$ process error ($\langle\text{x}_{i}\rangle$) can be obtained from $\mu_{\textbf{x}}$ in Eq. (<ref>) as follows.

\begin{equation}\label{eq:100}
\langle \text{x}_{i} \rangle = \mu_{\textbf{x}}((i-1)L+1:iL), i=1,...,N,
\end{equation}

Detailed derivations of Eqs. (<ref>), (<ref>), and (<ref>) are provided in the Appendices <ref>, <ref>, and <ref>, respectively.

M-step of VBEM: Temporal correlations in the time series data of each KCC (B) and noise ($\lambda$) are estimated in this step. Let $\tilde{\boldsymbol{\theta}}=\{ \text{B},\lambda \}$. Posterior distributions of $\textbf{x}, \boldsymbol{\alpha} ,\bar{b}$ obtained in Eqs. (<ref>), (<ref>), and (<ref>) are denoted as $q(\textbf{x};\tilde{\boldsymbol{\theta}}^{OLD}), q(\boldsymbol{\alpha};\tilde{\boldsymbol{\theta}}^{OLD})$, and $q(\bar{b};\tilde{\boldsymbol{\theta}}^{OLD})$ respectively. $\tilde{\boldsymbol{\theta}}$ can be updated by maximizing the complete likelihood as follows:

\begin{equation*}
\begin{split}
\tilde{\boldsymbol{\theta}}^{NEW} &= \argmaxB_{\tilde{\boldsymbol{\theta}}}\langle \ln{p(\textbf{y},\textbf{x},\boldsymbol{\alpha},\bar{b};\tilde{\boldsymbol{\theta}}}) \rangle_{q(\textbf{x};\tilde{\boldsymbol{\theta}}^{OLD})q(\boldsymbol{\alpha};\tilde{\boldsymbol{\theta}}^{OLD})q(\bar{b};\tilde{\boldsymbol{\theta}}^{OLD})} \\
&=\argmaxB_{\tilde{\boldsymbol{\theta}}}\langle \ln{p(\textbf{y}|\textbf{x};\tilde{\boldsymbol{\theta}}})p(\textbf{x}|\boldsymbol{\alpha};\tilde{\boldsymbol{\theta}}) \rangle_{q(\textbf{x};\tilde{\boldsymbol{\theta}}^{OLD})q(\boldsymbol{\alpha};\tilde{\boldsymbol{\theta}}^{OLD})}.
\end{split}
\end{equation*}

Let $Q(\tilde{\boldsymbol{\theta}}) = \langle \ln{p(\textbf{y}|\textbf{x};\tilde{\boldsymbol{\theta}}})p(\textbf{x}|\boldsymbol{\alpha};\tilde{\boldsymbol{\theta}}) \rangle_{q(\textbf{x};\tilde{\boldsymbol{\theta}}^{OLD})q(\boldsymbol{\alpha};\tilde{\boldsymbol{\theta}}^{OLD})}$, which results in

\begin{align}\label{eq:19}
Q(B,\lambda) &= \langle \ln{p(\textbf{y}|\textbf{x};\lambda)} \rangle_{q(\textbf{x};\tilde{\boldsymbol{\theta}}^{OLD})q(\boldsymbol{\alpha};\tilde{\boldsymbol{\theta}}^{OLD})}
+\langle \ln{p(\textbf{x}|\boldsymbol{\alpha};\text{B})} \rangle_{q(\textbf{x};\tilde{\boldsymbol{\theta}}^{OLD})q(\boldsymbol{\alpha};\tilde{\boldsymbol{\theta}}^{OLD})}  \nonumber \\
& =  \langle \ln{p(\textbf{y}|\textbf{x};\lambda)} \rangle_{q(\textbf{x};\tilde{\boldsymbol{\theta}}^{OLD})}
+\langle \ln{p(\textbf{x}|\boldsymbol{\alpha};\text{B})} \rangle_{q(\textbf{x};\tilde{\boldsymbol{\theta}}^{OLD})q(\boldsymbol{\alpha};\tilde{\boldsymbol{\theta}}^{OLD})}.
\end{align}

$\text{B}$ and $\lambda$ are estimated as Eqs. (<ref>) and (<ref>), respectively, by maximizing the Eq. (<ref>).

\begin{equation}\label{eq:22}
\text{B}= [\frac{1}{N}\sum_{i=1}^{N}\langle \alpha_{i} \rangle (\Sigma_{\text{x}_{i}} + \langle \text{x}_{i} \rangle \langle \text{x}_{i} \rangle^{\top})]^{-1}
\end{equation}

\begin{equation}\label{eq:25}
\lambda=\frac{[\lVert \textbf{y}-\textbf{D} \mu_{\textbf{x}} \rVert_{2}^{2} + \hat{\lambda} [NL - \text{Tr}(\Sigma_{\textbf{x}}\mathbb{E}_{q(\boldsymbol{\alpha};\tilde{\boldsymbol{\theta}}^{OLD})}(\Sigma_{0}^{-1}))]]}{ML}.
\end{equation}

Detailed derivations of Eqs. (<ref>) and (<ref>) are described in  Appendices <ref> and <ref>, respectively.

Algorithm <ref> shows the procedure of the proposed SA-TSBL method. Given the multiple KPCs samples ($\textbf{y}$) and fault pattern matrix ($\Phi$), the proposed method estimates the following variables and parameters in E and M steps, respectively.

* E-step: Process errors ($\mu_{\textbf{x}}$), the variable related to the sparsity of process errors ($ \boldsymbol{\alpha}$), and the variable that is used to provide prior knowledge of process faults ($\bar{b}$).

* M-step: Temporal correlations in the time series data of each KCC ($\text{B}$) and noise ($\lambda$).

These steps iterate until the estimator of process errors ($\mu_{\textbf{x}}$) is rarely updated, namely, $\lVert \mu_{\textbf{x}}^{t-1}-\mu_{\textbf{x}}^{t} \rVert_{\infty}<\gamma$, where $\lVert \cdot \rVert_{\infty}$ indicates infinity norm and $\gamma$ is a user-defined threshold (e.g., $\gamma=10^{-6}$). Then, the mean deviation of process errors, which is the output of Algorithm <ref>, are calculated by the following procedure.

* Step 1: The matrix $\tilde{\mu}_{\textbf{x}}$ is defined in which the $i^{th}$ row represents Eq. (<ref>) derived from $\mu_{\textbf{x}}$.

* Step 2: The average of each row of matrix $\tilde{\mu}_{\textbf{x}}$, that is, a vector  $\bar{\tilde{\mu}}_{\textbf{x}}$ is provided as the output of Algorithm <ref>, indicating the mean deviation of process errors.

Therefore, nonzero values in a vector $\bar{\tilde{\mu}}_{\textbf{x}}$ are process faults of the mean shifts.

Input: Multiple KPCs samples ($\textbf{y}$), Fault pattern matrix ($\Phi$).

Set $a=b_{i}\hspace{0.25em} (i \in P^{c})=10^{-4}$. $p$ and $q$ are set based on two conditions in Section <ref>.

Initialize  $\text{B}=\textbf{I}_{L}, \lambda=1,  b_{i} \hspace{0.25em} (i\in P) =1,,t=1.$

While  $\lVert \mu_{\textbf{x}}^{t-1}-\mu_{\textbf{x}}^{t} \rVert_{\infty}< \gamma $ do

E-step of VBEM:

Update $ \mu_{\textbf{x}}$ using Eq. (<ref>)

Update $ \boldsymbol{\alpha}$ using Eq. (<ref>)

Update $ \bar{b}$ using Eq. (<ref>)

M-step of VBEM:

Update $\text{B}$ using Eq. (<ref>)

Update $\lambda$ using Eq. (<ref>)

Output: Mean deviations of process errors $\bar{\tilde{\mu}}_{\textbf{x}}.$
Proposed SA-TSBL method

§ NUMERICAL CASE STUDIES

This section provides three scenarios to compare the performance between the proposed method and benchmark methods.

* Section <ref>  shows the performance evaluation by varying the temporal correlation. This study is to validate the effectiveness of strong temporal correlation in sparse estimation.

* Section <ref> provides the numerical study to investigate the impact of the number of KPCs samples (i.e., measurement samples) on the performance of process faults identification.

* Section <ref> illustrates the sparse estimation performance by varying the ratio between the number of measurements and process errors (i.e., the severity of the underdetermined systems).

All the numerical case studies consist of 100 independent trials. $p$ and $q$ in Eq. (<ref>) are determined as 1 and 0.1, respectively, in the case studies to satisfy the two conditions provided in Section <ref>. The code of the proposed SA-TSBL algorithm is implemented in Matlab 2017. The CPU of the computer used in this paper is an Intel® Core™ Processor i7-8750H.

The benchmark methods selected in this study are as follows, which are widely used in sparse Bayesian learning.

* MSBL proposed in [37] is a basic SBL method for the MMV model that assumes independence in the time series data of each KCC.

* T-MSBL proposed in [42] is a typical SBL method for the MMV model that considers temporal correlation in the time series data of each KCC.

* SA-MSBL proposed in [38] is the SBL method that considers prior knowledge of support in the MMV model. It assumes independence in the time series data of each KCC.

* SA-SBL proposed in [16] and [27] is the SBL method that considers prior knowledge of process faults in the single measurement vector (SMV) model. Like [2], the average of multiple KPCs samples is used to estimate the mean deviation of process errors.

Data Generations: The data generation process is summarized in Figure <ref>. The solution matrix $\textbf{X}_{true}\in \mathbb{R}^{N\times L}$ is randomly generated with $K$ nonzero rows. Indexes of the nonzero rows are randomly chosen in each trial. The nonzero rows in $\textbf{X}_{true}$ are generated as AR(1) process that initiates from the standard Gaussian distribution since AR(1) processes are sufficient to represent the temporal structure of the small number of measurement samples ($L$) [42]. The AR coefficient, defined as $\beta$, represents the temporal correlation. A dictionary matrix $\Phi \in \mathbb{R}^{M \times N}$ is constructed with columns drawn from the surface of a unit hyper-sphere uniformly [14]. Finally, the measurements matrix is built by $ \textbf{Y} =\Phi\textbf{X}_{true}+\textbf{V}$ in the final step of Figure <ref>, where $\textbf{V}$ is a Gaussian noise matrix with zero-mean [42]. The variance of the noise matrix is chosen to meet the determined value of the signal-to-noise ratio (SNR). SNR is defined by SNR(dB)$\vcentcolon= 20(\text{log}_{10} (\lVert \Phi \textbf{X}_{true} \rVert_{F} / \lVert \textbf{V} \rVert_{F}))$ [42].

Flowchart of the data generation process.

Performance Evaluation: Since the objective of this paper is to diagnose the mean deviation of process errors, the target of the proposed method is to accurately estimate $\bar{\textbf{X}}_{true}$, which is the row-wise mean of $\textbf{X}_{true}$.

Two performance measures are used in this paper. One is the failure rate defined in [27]. It measures the accuracy of detecting process fault, which are the nonzero rows in $\bar{\textbf{X}}_{true}$. Assume that the number of process faults $K$ is given. Then, the row indexes of the $K$ largest $\ell_2$-norms from $\bar{\tilde{\mu}}_{\textbf{x}}$ in Algorithm <ref> are identified. If the indexes are different from the indexes of nonzero rows in $\bar{\textbf{X}}_{true}$, it is considered a failed trial. The failure rate is the percentage of failed trials in the total trials.

The other performance measure is a normalized mean squared error (NMSE) that is defined by $\lVert \bar{\tilde{\mu}}_{\textbf{x}}-\bar{\textbf{X}}_{true} \rVert_{F}^{2} / \lVert \bar{\textbf{X}}_{true} \rVert_{F}^{2}$. Averages of failure rate and NMSE from 100 trials are used as performance measures.
(Performance Evaluations of SA-TSBL, SA-MSBL, and SA-SBL)

The SA-MSBL, SA-SBL, and the proposed SA-TSBL have several cases based on the set of prior knowledge (i.e., $P$ in Eq. (<ref>)) even in the same problem. Both correct and erroneous information of support exists in the set $P$ as prior knowledge. To differentiate between the correct and erroneous information of support, two subsets, namely, $P_{C}$ and $P_{E}$ are defined. $P_{C}$ consists of the indexes of nonzero rows obtained from prior knowledge among true nonzero rows in $\bar{\textbf{X}}_{true}$. In contrast, $P_{E}$ consists of the indexes of nonzero rows in prior knowledge but are actually zero in $\bar{\textbf{X}}_{true}$. The cardinalities of $P_{C}$ and $P_{E}$ in this paper are assumed to satisfy two conditions. The first condition is the cardinality of set $P_{C}$ and $P_{E}$ should be less than or equal to 75$\%$ and 50$\%$ of the number of nonzero rows in $\bar{\textbf{X}}_{true}$ (i.e., $K$), respectively. For example, if $K$ is 6, the cardinality of $P_{C}$ and $P_{E}$ are less than or equal to 4 and 3, respectively. The first condition illustrates the cardinality of prior knowledge is similar to $K$. The condition enables the performance evaluation of the proposed method in comprehensive situations. Since the cardinality of $P_{C}$ is less than or equal to 4, the prior knowledge is partial (i.e., $P_{C}$ misses some rows that are actually nonzero in $\bar{\textbf{X}}_{true}$) in all cases. In addition, the cases with a cardinality of $P_{E}$ greater than zero show the situations when the erroneous prior knowledge (i.e., $P_{E}$ contains the rows that are actually zero in $\bar{\textbf{X}}_{true}$) exist. The second condition is the cardinality of $P_{C}$ is greater than or equal to that of $P_{E}$. The second condition prevents erroneous prior knowledge from being dominant prior knowledge. Therefore, performances of these three methods are evaluated as the average of all cases with different cardinalities of $P_{C}$ and $P_{E}$, where each case consists of 100 trials.

§.§ Performance Evaluation in Various Temporal Correlations

This case study shows the performance of all methods in various temporal correlations under the noiseless case. The size of the dictionary matrix $\Phi$ is 8$\times$40. The number of nonzero rows of $\bar{\textbf{X}}_{true}$ ($K$), and measurement samples ($L$) are 6 and 3, respectively. Temporal correlation ($\beta$) varies among 0.1, 0.3, 0.5, 0.7, 0.9, and 0.99. The proposed method shows the best performance in all temporal correlations, as shown in Table <ref>.
Performance comparison with various temporal correlations ($\beta$).

5cFailure Rate     5cNMSE

$\beta$       0.1      0.3      0.6      0.9      0.99     0.1      0.3      0.6      0.9      0.99

T-MSBL       0.35     0.36     0.29     0.22     0.20     0.62     0.54     0.42     0.35     0.26

MSBL         0.54     0.57     0.73     0.88     0.98     0.58     0.59     0.84     1.02     1.25

SA-MSBL      0.98     0.98     0.98     0.98     0.99     0.75     0.78     0.82     0.89     0.97

SA-SBL       0.99     0.99     0.99     0.99     0.99     0.84     0.86     0.90     0.91     0.99

(Proposed)     0.27     0.27     0.23     0.18     0.16     0.45     0.40     0.34     0.28     0.22

The boxplot for (a) failure rate and (b) NMSE in all cases of SA-TSBL in various temporal correlations.

Specifically, the failure rate and NMSE of the proposed method and T-MSBL tend to decrease as $\beta$ increases since these two methods capture the temporal correlation. In contrast, both measures of MSBL and SA-MSBL show the opposite trend as the two methods assume independence among samples.

The performance evaluations of the proposed method with different cardinalities of $P_{C}$ and $P_{E}$ are represented as boxplots in Figure <ref>. Since the number of nonzero rows is 6, each correlation has 14 cases designed under the two conditions of Remark <ref>. In each boxplot, the horizontal line represents the minimum, the first quartile, median, the third quartile, and maximum value sequentially from the bottom, respectively. The minimum value denotes the smallest value within the 1.5 $\times$ interquartile range (IQR) below the first quartile. Similarly, the maximum value is defined as the largest value that is within the 1.5 $\times$ IQR above the third quartile [22]. A blue square in Figure <ref> shows the performance of the proposed method without any prior knowledge (i.e., $|P_{C}|=|P_{E}|=0$), and the black dot represents the average of all cases used as a performance evaluation measure of the proposed method. The trend of black dots shows the performance improvements of the proposed method as  temporal correlation ($\beta$) increases. In addition, boxplots in Figure <ref> represent all cases with partial and some erroneous prior knowledge achieve better performance in both measures than those without prior knowledge. Even in the case of $\beta=0.6$, a square is considered an outlier since it is located higher than
1.5 $\times$ IQR from the third quartile of the boxplot. The results show the effectiveness of prior knowledge in the proposed method, which successfully distinguishes the correct and incorrect prior knowledge in various temporal correlations.

Performance Evaluation in the Various Numbers of Measurement Samples

The objective of this case study is to compare the performance of all methods using the various numbers of measurement samples ($L$). The size of the dictionary matrix $\Phi$ is 7$\times$55, and the number of nonzero rows in $\bar{\textbf{X}}_{true}$ ($K$) is 4. SNR and temporal correlation ($\beta$) are set as 35dB, and 0.95, respectively.
The number of measurement samples varies from 2 to 4. Table <ref> shows that the proposed method achieves the best performance in all cases.

Performance comparison with various number of measurement samples ($L$).

3cFailure   Rate     3cNMSE

$L$            2              3             4             2          3          4

T-MSBL       0.59           0.17          0.13          0.75       0.19       0.13

MSBL         0.83           0.77          0.52          1.09       1.04       0.64

SA-MSBL      0.86           0.88          0.89          1.04       1.09       0.92

SA-SBL       0.88           0.89          0.89          1.06       1.05       0.92

(Proposed)      0.45           0.12          0.10          0.52       0.12       0.10

In addition, the results in Table <ref> illustrate that most of the MMV models achieve performance improvements as the number of measurement samples increases because of common support assumption [10].

The boxplot for (a) failure rate and (b) NMSE in all cases of SA-TSBL in various numbers of measurement samples.

Figure <ref> shows the performance evaluation of the proposed method with various sizes of $P_{C}$ and $P_{E}$. Especially, the performance of T-MSBL and the proposed method without prior knowledge (blue square in Figure <ref>) is similar. However, the proposed achieves better performance than T-MSBL by utilizing partial and some erroneous prior knowledge. Specifically, the performance of the proposed method improves from 4$\%$ to 31$\%$ of those without any prior knowledge (a relative improvement from square to dot in Figure <ref>) in the various numbers of measurement samples.

§.§ Performance Evaluation in Various Ratios between the Number of Measurements and Process Errors

Results in Table <ref> illustrate the performance of all methods by varying the ratio between the number of measurements and process errors (i.e., underdetermined ratio). In this study, the number of measurements ($M$) is fixed at 10, and the underdetermined ratio ($N/M$) is selected from 3, 5, 7, and 9 with 25dB for SNR, respectively. The number of nonzero rows in $\bar{\textbf{X}}_{true}$ ($K$) and measurement samples ($L$) are set as 4, 3, respectively. The temporal correlation ($\beta$) is 0.99. As the underdetermined ratio increases, it becomes more challenging to identify the sparse process faults. However, capturing temporal correlation and utilizing the partial and some erroneous prior knowledge enable the proposed method to achieve the best performance even in a high underdetermined ratio. Figure <ref> shows prior knowledge of support is still valuable in various underdetermined ratios. This study shows that the proposed algorithm can be applied to applications such as neuroimaging in that highly underdetermined systems exist.

Performance comparison with various underdetermined ratio ($N/M$).

4cFailure Rate                                  4cNMSE
$N/M$      3   5    7     9      3     5     7     9

T-MSBL      0.09      0.28      0.26      0.31      0.05      0.31      0.33      0.42

MSBL     0.29     0.56     0.64      0.69      0.30      0.61     0.77      0.85

SA-MSBL     0.42     0.66     0.73     0.82     0.35     0.59     0.65     0.73

SA-SBL     0.67     0.85     0.88     0.90     0.48     0.67     0.71     0.75

(Proposed)     0.08     0.21     0.22     0.24     0.04     0.21     0.28     0.29

The boxplot for (a) failure rate and (b) NMSE in all cases of SA-TSBL in various underdetermined ratio.

§ REAL-WORLD SIMULATION CASE STUDIES

An assembly model from a real auto body assembly process is utilized as a real-world case study. As shown in Figure <ref>, the assembled product is a floor pan of a car. It consists of four parts namely, the left floor pan, left bracket, right floor pan, and right bracket.

Floor-pan assembly model [3].

They are assembled in three stations as illustrated in Figure <ref>. During the assembly process, the parts are held by fixtures and part-mating features, which are the KCCs in this process [4]. KPCs are measured from four points, namely, M1, M2, M3, and M4, respectively, as shown in Figure <ref>. Each part has its own measurement place such as M1 on part 1, M2 on part 2, M3 on part 3, and M4 on part 4. These measurements can be measured in each station when the corresponding part has been assembled in the previous stations. For example, M4 on part 4 cannot be measured in station 1 because part 4 has not yet been assembled at station 1 [4]. In addition, KPCs are measured in three directions (X, Y, and Z) at each point. In this assembly process, there exists a total of 33 process errors, which are fixture locator dimensional errors [2]. The dimension of the fault pattern matrix is 12$\times$33, and is established based on the literature [4, 21, 24]. Since the variations of KCCs can be propagated to deviations of other KCCs in the subsequent stations, the transfer and accumulation of KCCs deviations between multiple stations are considered in the formulation of the fault pattern matrix $\Phi$ [13, 4]. The matrix is provided in Appendix <ref>. Since the number of measurements (12) is less than the number of process errors (33), it causes an underdetermined system in the fault quality linear model. It requires sparse estimation to identify process faults.

Floor-pan assembly procedure from three assembly stations [3].

To generate the multiple KPCs samples from the auto body assembly process, the time series data of each process fault is generated from AR (1) process initiated from the standard Gaussian distribution as in Section <ref>. The generated temporal correlated process faults ($\textbf{X}_{true}$) and fault pattern matrix provide the multiple KPCs samples ($\textbf{Y}$), as in Figure <ref>. In addition, performance evaluation measures,  benchmark methods, and the values of $p$ and $q$ in Eq. (<ref>) used in Section <ref> are still utilized in Section <ref>. Prior knowledge of process faults is also provided in the same way as Section <ref>. Sections <ref> and <ref> show the performance evaluation by varying the number of KPCs samples, and temporal correlation $\beta$, respectively.

Performance Evaluation in the Various Number of KPCs Samples

This case study aims to demonstrate the effectiveness of sparse estimation of process faults by varying the number of KPCs samples ($L$) when a strong correlation exists ($\beta$=0.99). The number of KPCs samples varies from 2 to 4, and SNR is 50dB. Three process faults ($K$) are determined among 33 process errors randomly. The proposed method shows the best performance in all various numbers of KPCs samples, as shown in Table <ref>.

Performance comparison by varying the number of KPCs samples.

3cFailure Rate     3cNMSE

KPCs samples                  2          3          4            2        3        4

T-MSBL         0.59       0.63       0.51      0.34     0.32     0.25

MSBL            0.59       0.62       0.50       0.36     0.33     0.24

SA-MSBL         0.57       0.62       0.52       0.40     0.40     0.34

SA-SBL         0.60       0.65       0.58       0.43     0.42     0.40
(Proposed)         0.48       0.54       0.38      0.30     0.30     0.21

All methods except SA-SBL generally tend to improve performance in both measures as the number of KPCs samples increases. However, the performances of all methods do not improve significantly as the number of KPCs samples increases, compared to Section <ref>. The reason could be the fault pattern matrix in this multistation assembly process, which is very structured compared to the random design matrix $\Phi$ in Section <ref>. The random matrix achieves an accurate sparse estimation based on some theoretical properties of the random matrix, such as low mutual coherence that measures the highest correlation between columns of $\Phi$ [5]. In contrast, if the mutual coherence is high, it causes inaccurate sparse estimation [15]. As shown in Appendix <ref>, the columns of fault pattern matrix $\Phi$ in the assembly process are highly correlated, and mutual coherence is 1, causing the challenging sparse estimation task than Section <ref>. Therefore, the results in this section are not significantly improved compared to Section <ref>, which uses a random matrix with low mutual coherence, even if the number of KPCs samples is increased.

Figure <ref> shows that utilizing prior knowledge of process faults is still effective in the proposed method for failure rate and NMSE, even in the structured design matrix $\Phi$. The property lets the proposed method obtain a more accurate sparse estimation than T-MSBL, MSBL that cannot incorporate the prior knowledge of process faults.

The boxplot for (a) failure rate and (b) NMSE in all cases of SA-TSBL in the various number of KPCs samples.

§.§ Performance Evaluation in Various Temporal Correlations

This case study presents the performance of all methods in various temporal correlations $\beta$, when there exist five KPCs samples ($L$). Temporal correlation varies among 0.0, 0.3, 0.6, 0.9, and 0.99. The study has three process faults ($K$) with noise level 80dB. Table <ref> shows the proposed method achieves the best performance, and the performance of other methods exhibit similar trends to the case in Section. <ref>. The performance of the proposed method and T-MSBL improves in general as $\beta$ increases, and MSBL and SA-MSBL show the opposite trends since they assume the independence in the time series data of each process error. Figure <ref> illustrates the prior knowledge of process faults is still effective in various temporal correlations even with the structured matrix $\Phi$.

Performance comparison by varying temporal correlations ($\beta$).

5cFailure Rate                                                  5cNMSE
$\beta$                  0.0               0.3               0.6               0.9               0.99              0.0               0.3               0.6               0.9               0.99

T-MSBL                0.62              0.59              0.56              0.52              0.54              0.19              0.21              0.21              0.19              0.20

MSBL                  0.59              0.55              0.58              0.61              0.63              0.18              0.20              0.21              0.24              0.30

SA-MSBL               0.57              0.58              0.55              0.57              0.55              0.29              0.30              0.29              0.31              0.35

SA-SBL                0.75              0.73              0.72              0.66              0.60              0.46              0.45              0.42              0.38              0.39

(Proposed)     0.52     0.52     0.48     0.47     0.40     0.17     0.20     0.18     0.16     0.16

The boxplot for (a) failure rate and (b) NMSE in all cases of SA-TSBL in various temporal correlations.

§ CONCLUSIONS

This paper proposes a novel sparse hierarchical Bayesian method, SA-TSBL, to effectively identify the sparse process faults in multistation assembly systems. The method identifies process faults by considering the temporal correlation of each KCC and utilizing partial with some erroneous prior knowledge of process faults. Since posterior distributions of process errors in the proposed method are computationally intractable, this paper derives approximate posterior distributions of process errors via Variational Bayes inference. The effectiveness of the proposed method is validated by both numerical cases and real-world simulation application that uses an actual auto body assembly system. Based on these studies, it is evident that the direct use of multiple KPCs samples in the proposed SA-TSBL is more effective to process mean shift identification than in a previous study using the average of multiple KPCs samples in the single measurement vector model (Eq. (<ref>)). This is because the proposed method can fully utilize multiple KPCs samples without any information loss. In addition, the results in case studies represent the proposed method achieves high performance in process faults estimation when the time series data of each KCC have a strong temporal correlation. Furthermore, utilizing the prior knowledge of process faults improves the performance of the proposed method in the case studies, even if the knowledge has some erroneous information. This is possible through the Bayesian framework in the proposed method that distinguishes between the correct and incorrect prior knowledge.

In this work, all fixture locators are assumed to be independent, which is a common assumption in the literature [26, 4, 28]. However, there exists a correlation between process errors in some manufacturing systems. Specifically, the fixture locators in the multistation assembly process have some spatial correlation with each other if the locators are in the composite tolerance mode [25]. Therefore, expanding the proposed SA-TSBL considering the spatial correlation among the process errors is a promising direction for future work. In addition, instead of utilizing the Gamma distribution in Eq. (<ref>) enforcing the sparsity of process error in the proposed method, other distributions such as  Laplacian distribution [31], and Generalized Pareto distribution [40] can also be considered as prior distribution since the probability of these distributions is concentrated at zero. Based on each of these distributions, deriving new hierarchical models considering both the temporal correlation of each process error and prior knowledge of process faults is also a valuable direction of future work. Finally, extending the proposed method to consider the different temporal correlations between the process errors is an additional direction for future work.

[1]
Abellan-Nebot, J. V., J. Liu, F. R. Subirón, and J. Shi (2012).
State space modeling of variation propagation in multistation
machining processes considering machining-induced variations.
Journal of Manufacturing Science and Engineering 134(2).

[2]
Bastani, K., B. Barazandeh, and Z. J. Kong (2018).
Fault diagnosis in multistation assembly systems using spatially
correlated bayesian learning algorithm.
Journal of Manufacturing Science and Engineering 140(3).

[3]
Bastani, K., Z. Kong, W. Huang, X. Huo, and Y. Zhou (2012).
Fault diagnosis using an enhanced relevance vector machine (rvm) for
partially diagnosable multistation assembly processes.
IEEE Transactions on Automation Science and Engineering 10(1), 124–136.

[4]
Bastani, K., Z. Kong, W. Huang, and Y. Zhou (2016).
Compressive sensing–based optimal sensor placement and fault
diagnosis for multi-station assembly processes.
IIE Transactions 48(5), 462–474.

[5]
Candès, E. J. and M. B. Wakin (2008).
An introduction to compressive sampling.
IEEE signal processing magazine 25(2), 21–30.

[6]
Candes, E. J., M. B. Wakin, and S. P. Boyd (2008).
Enhancing sparsity by reweighted $\ell_{1}$ 1 minimization.
Journal of Fourier analysis and applications 14(5),

[7]
Ceglarek, D. and J. Shi (1996).
Fixture failure diagnosis for autobody assembly using pattern
Journal of manufacturing science and engineering 118(1).

[8]
Chen, S. S., D. L. Donoho, and M. A. Saunders (2001).
Atomic decomposition by basis pursuit.
SIAM review 43(1), 129–159.

[9]
Cohn, I., T. El-Hay, N. Friedman, and R. Kupferman (2010).
Mean field variational approximation for continuous-time bayesian
The Journal of Machine Learning Research 11,

[10]
Cotter, S. F., B. D. Rao, K. Engan, and K. Kreutz-Delgado (2005).
Sparse solutions to linear inverse problems with multiple measurement
IEEE Transactions on Signal Processing 53(7),

[11]
Dai, J. and H. C. So (2021).
Group-sparsity learning approach for bearing fault diagnosis.
IEEE Transactions on Industrial Informatics 18(7),

[12]
Ding, Y., D. Ceglarek, and J. Shi (2002).
Fault diagnosis of multistage manufacturing processes by using state
space approach.
J. Manuf. Sci. Eng. 124(2), 313–322.

[13]
Ding, Y., D. Ceglarek, J. Shi, et al. (2000).
Modeling and diagnosis of multistage manufacturing processes: part i
state space model.
In Proceedings of the 2000 Japan/USA symposium on flexible
automation, pp.  23–26. Ann Arbor, MI.

[14]
Donoho, D. L. (2006).
For most large underdetermined systems of linear equations the
minimal $\ell_{1}$ norm solution is also the sparsest solution.
Communications on Pure and Applied Mathematics: A Journal Issued
by the Courant Institute of Mathematical Sciences 59(6), 797–829.

[15]
Elad, M. (2007).
Optimized projections for compressed sensing.
IEEE Transactions on Signal Processing 55(12),

[16]
Fang, J., Y. Shen, F. Li, H. Li, and Z. Chen (2015).
Support knowledge-aided sparse bayesian learning for compressed
In 2015 IEEE International Conference on Acoustics, Speech and
Signal Processing (ICASSP), pp.  3786–3790. IEEE.

[17]
Guo, Y., B. Sun, N. Li, and D. Fang (2017).
Variational bayesian inference-based counting and localization for
off-grid targets with faulty prior information in wireless sensor networks.
IEEE Transactions on Communications 66(3),

[18]
Han, N. and Z. Song (2018).
Bayesian multiple measurement vector problem with spatial structured
sparsity patterns.
Digital Signal Processing 75, 184–201.

[19]
Han, T., D. Jiang, Y. Sun, N. Wang, and Y. Yang (2018).
Intelligent fault diagnosis method for rotating machinery via
dictionary learning and sparse representation-based classification.
Measurement 118, 181–193.

[20]
Huang, W., J. Lin, M. Bezdecny, Z. Kong, and D. Ceglarek (2018).
Stream-of-variation modeling—part i: A generic three-dimensional
variation model for rigid-body assembly in single station assembly processes.
Journal of Manufacturing Science and Engineering 140(3).

[21]
Huang, W., J. Lin, Z. Kong, and D. Ceglarek (2007).
Stream-of-variation (sova) modeling ii: a generic 3d variation model
for rigid body assembly in multistation assembly processes.
Journal of Manufacturing Science and Engineering 129(4).

[22]
Hubert, M. and E. Vandervieren (2008).
An adjusted boxplot for skewed distributions.
Computational statistics & data analysis 52(12),

[23]
Kong, Z., D. Ceglarek, and W. Huang (2008).
Multiple fault diagnosis method in multistation assembly processes
using orthogonal diagonalization analysis.
Journal of manufacturing science and engineering 130(1).

[24]
Kong, Z., W. Huang, and A. Oztekin (2009).
Variation propagation analysis for multistation assembly process with
consideration of gd&t factors.
Journal of Manufacturing Science and Engineering 131(5).

[25]
Kong, Z., R. Kumar, S. Gogineni, Y. Zhou, and D. Ceglarek (2006).
Mode-based tolerance analysis in multi-station assembly using stream
of variation model.
Trans. NAMRI/SME 34, 469–476.

[26]
Lee, J., J. Son, S. Zhou, and Y. Chen (2020).
Variation source identification in manufacturing processes using
bayesian approach with sparse variance components prior.
IEEE Transactions on Automation Science and Engineering 17(3), 1469–1485.

[27]
Li, H., Y. Jiang, J. Fang, and M. Rangaswamy (2017).
Adaptive subspace signal detection with uncertain partial prior
IEEE Transactions on Signal Processing 65(16),

[28]
Li, S. and Y. Chen (2016).
A bayesian variable selection method for joint diagnosis of
manufacturing process and sensor faults.
IIE transactions 48(4), 313–323.

[29]
Liu, J. (2010).
Variation reduction for multistage manufacturing processes: a
comparison survey of statistical-process-control vs stream-of-variation
Quality and Reliability Engineering International 26(7), 645–661.

[30]
Luessi, M., S. D. Babacan, R. Molina, and A. K. Katsaggelos (2013).
Bayesian simultaneous sparse approximation with smooth signals.
IEEE Transactions on Signal Processing 61(22),

[31]
Park, T. and G. Casella (2008).
The bayesian lasso.
Journal of the American Statistical Association 103(482), 681–686.

[32]
Petersen, K. B., O. Winther, and L. K. Hansen (2005).
On the slow convergence of em and vbem in low-noise linear models.
Neural computation 17(9), 1921–1926.

[33]
Shi, J. (2022).
In-process quality improvement: Concepts, methodologies, and
IISE Transactions (just-accepted), 1–34.

[34]
Tibshirani, R. (1996).
Regression shrinkage and selection via the lasso.
Journal of the Royal Statistical Society: Series B
(Methodological) 58(1), 267–288.

[35]
Tipping, M. E. (2001).
Sparse bayesian learning and the relevance vector machine.
Journal of machine learning research 1(Jun),

[36]
Wipf, D. and S. Nagarajan (2009).
Solving sparse linear inverse problems: Analysis of reweighted
$\ell_1$ and $\ell_2$ methods.
In SPARS'09-Signal Processing with Adaptive Sparse Structured

[37]
Wipf, D. P. and B. D. Rao (2007).
An empirical bayesian strategy for solving the simultaneous sparse
approximation problem.
IEEE Transactions on Signal Processing 55(7),

[38]
Yu, D., Y. Guo, N. Li, and M. Wang (2019).
Sa-m-sbl: An algorithm for csi-based device-free localization with
faulty prior information.
IEEE Access 7, 61831–61839.

[39]
Zhang, C., H. Yan, S. Lee, and J. Shi (2021).
Dynamic multivariate functional data modeling via sparse subspace
Technometrics 63(3), 370–383.

[40]
Zhang, S., Y. Liu, X. Li, and G. Bi (2017).
Variational bayesian sparse signal recovery with lsm prior.
IEEE Access 5, 26690–26702.

[41]
Zhang, Z., W. Huang, Y. Liao, Z. Song, J. Shi, X. Jiang, C. Shen, and Z. Zhu
Bearing fault diagnosis via generalized logarithm sparse
Mechanical Systems and Signal Processing 167, 108576.

[42]
Zhang, Z. and B. D. Rao (2011).
Sparse signal recovery with temporally correlated source vectors
using sparse bayesian learning.
IEEE Journal of Selected Topics in Signal Processing 5(5), 912–926.

[43]
Zhou, S., Y. Chen, and J. Shi (2004).
Statistical estimation and testing for variation root-cause
identification of multistage manufacturing processes.
IEEE Transactions on Automation Science and Engineering 1(1), 73–83.

[1]
Abellan-Nebot, J. V., J. Liu, F. R. Subirón, and J. Shi (2012).
State space modeling of variation propagation in multistation
machining processes considering machining-induced variations.
Journal of Manufacturing Science and Engineering 134(2).

[2]
Bastani, K., B. Barazandeh, and Z. J. Kong (2018).
Fault diagnosis in multistation assembly systems using spatially
correlated bayesian learning algorithm.
Journal of Manufacturing Science and Engineering 140(3).

[3]
Bastani, K., Z. Kong, W. Huang, X. Huo, and Y. Zhou (2012).
Fault diagnosis using an enhanced relevance vector machine (rvm) for
partially diagnosable multistation assembly processes.
IEEE Transactions on Automation Science and Engineering 10(1), 124–136.

[4]
Bastani, K., Z. Kong, W. Huang, and Y. Zhou (2016).
Compressive sensing–based optimal sensor placement and fault
diagnosis for multi-station assembly processes.
IIE Transactions 48(5), 462–474.

[5]
Candès, E. J. and M. B. Wakin (2008).
An introduction to compressive sampling.
IEEE signal processing magazine 25(2), 21–30.

[6]
Candes, E. J., M. B. Wakin, and S. P. Boyd (2008).
Enhancing sparsity by reweighted $\ell_{1}$ 1 minimization.
Journal of Fourier analysis and applications 14(5),

[7]
Ceglarek, D. and J. Shi (1996).
Fixture failure diagnosis for autobody assembly using pattern
Journal of manufacturing science and engineering 118(1).

[8]
Chen, S. S., D. L. Donoho, and M. A. Saunders (2001).
Atomic decomposition by basis pursuit.
SIAM review 43(1), 129–159.

[9]
Cohn, I., T. El-Hay, N. Friedman, and R. Kupferman (2010).
Mean field variational approximation for continuous-time bayesian
The Journal of Machine Learning Research 11,

[10]
Cotter, S. F., B. D. Rao, K. Engan, and K. Kreutz-Delgado (2005).
Sparse solutions to linear inverse problems with multiple measurement
IEEE Transactions on Signal Processing 53(7),

[11]
Dai, J. and H. C. So (2021).
Group-sparsity learning approach for bearing fault diagnosis.
IEEE Transactions on Industrial Informatics 18(7),

[12]
Ding, Y., D. Ceglarek, and J. Shi (2002).
Fault diagnosis of multistage manufacturing processes by using state
space approach.
J. Manuf. Sci. Eng. 124(2), 313–322.

[13]
Ding, Y., D. Ceglarek, J. Shi, et al. (2000).
Modeling and diagnosis of multistage manufacturing processes: part i
state space model.
In Proceedings of the 2000 Japan/USA symposium on flexible
automation, pp.  23–26. Ann Arbor, MI.

[14]
Donoho, D. L. (2006).
For most large underdetermined systems of linear equations the
minimal $\ell_{1}$ norm solution is also the sparsest solution.
Communications on Pure and Applied Mathematics: A Journal Issued
by the Courant Institute of Mathematical Sciences 59(6), 797–829.

[15]
Elad, M. (2007).
Optimized projections for compressed sensing.
IEEE Transactions on Signal Processing 55(12),

[16]
Fang, J., Y. Shen, F. Li, H. Li, and Z. Chen (2015).
Support knowledge-aided sparse bayesian learning for compressed
In 2015 IEEE International Conference on Acoustics, Speech and
Signal Processing (ICASSP), pp.  3786–3790. IEEE.

[17]
Guo, Y., B. Sun, N. Li, and D. Fang (2017).
Variational bayesian inference-based counting and localization for
off-grid targets with faulty prior information in wireless sensor networks.
IEEE Transactions on Communications 66(3),

[18]
Han, N. and Z. Song (2018).
Bayesian multiple measurement vector problem with spatial structured
sparsity patterns.
Digital Signal Processing 75, 184–201.

[19]
Han, T., D. Jiang, Y. Sun, N. Wang, and Y. Yang (2018).
Intelligent fault diagnosis method for rotating machinery via
dictionary learning and sparse representation-based classification.
Measurement 118, 181–193.

[20]
Huang, W., J. Lin, M. Bezdecny, Z. Kong, and D. Ceglarek (2018).
Stream-of-variation modeling—part i: A generic three-dimensional
variation model for rigid-body assembly in single station assembly processes.
Journal of Manufacturing Science and Engineering 140(3).

[21]
Huang, W., J. Lin, Z. Kong, and D. Ceglarek (2007).
Stream-of-variation (sova) modeling ii: a generic 3d variation model
for rigid body assembly in multistation assembly processes.
Journal of Manufacturing Science and Engineering 129(4).

[22]
Hubert, M. and E. Vandervieren (2008).
An adjusted boxplot for skewed distributions.
Computational statistics & data analysis 52(12),

[23]
Kong, Z., D. Ceglarek, and W. Huang (2008).
Multiple fault diagnosis method in multistation assembly processes
using orthogonal diagonalization analysis.
Journal of manufacturing science and engineering 130(1).

[24]
Kong, Z., W. Huang, and A. Oztekin (2009).
Variation propagation analysis for multistation assembly process with
consideration of gd&t factors.
Journal of Manufacturing Science and Engineering 131(5).

[25]
Kong, Z., R. Kumar, S. Gogineni, Y. Zhou, and D. Ceglarek (2006).
Mode-based tolerance analysis in multi-station assembly using stream
of variation model.
Trans. NAMRI/SME 34, 469–476.

[26]
Lee, J., J. Son, S. Zhou, and Y. Chen (2020).
Variation source identification in manufacturing processes using
bayesian approach with sparse variance components prior.
IEEE Transactions on Automation Science and Engineering 17(3), 1469–1485.

[27]
Li, H., Y. Jiang, J. Fang, and M. Rangaswamy (2017).
Adaptive subspace signal detection with uncertain partial prior
IEEE Transactions on Signal Processing 65(16),

[28]
Li, S. and Y. Chen (2016).
A bayesian variable selection method for joint diagnosis of
manufacturing process and sensor faults.
IIE transactions 48(4), 313–323.

[29]
Liu, J. (2010).
Variation reduction for multistage manufacturing processes: a
comparison survey of statistical-process-control vs stream-of-variation
Quality and Reliability Engineering International 26(7), 645–661.

[30]
Luessi, M., S. D. Babacan, R. Molina, and A. K. Katsaggelos (2013).
Bayesian simultaneous sparse approximation with smooth signals.
IEEE Transactions on Signal Processing 61(22),

[31]
Park, T. and G. Casella (2008).
The bayesian lasso.
Journal of the American Statistical Association 103(482), 681–686.

[32]
Petersen, K. B., O. Winther, and L. K. Hansen (2005).
On the slow convergence of em and vbem in low-noise linear models.
Neural computation 17(9), 1921–1926.

[33]
Shi, J. (2022).
In-process quality improvement: Concepts, methodologies, and
IISE Transactions (just-accepted), 1–34.

[34]
Tibshirani, R. (1996).
Regression shrinkage and selection via the lasso.
Journal of the Royal Statistical Society: Series B
(Methodological) 58(1), 267–288.

[35]
Tipping, M. E. (2001).
Sparse bayesian learning and the relevance vector machine.
Journal of machine learning research 1(Jun),

[36]
Wipf, D. and S. Nagarajan (2009).
Solving sparse linear inverse problems: Analysis of reweighted
$\ell_1$ and $\ell_2$ methods.
In SPARS'09-Signal Processing with Adaptive Sparse Structured

[37]
Wipf, D. P. and B. D. Rao (2007).
An empirical bayesian strategy for solving the simultaneous sparse
approximation problem.
IEEE Transactions on Signal Processing 55(7),

[38]
Yu, D., Y. Guo, N. Li, and M. Wang (2019).
Sa-m-sbl: An algorithm for csi-based device-free localization with
faulty prior information.
IEEE Access 7, 61831–61839.

[39]
Zhang, C., H. Yan, S. Lee, and J. Shi (2021).
Dynamic multivariate functional data modeling via sparse subspace
Technometrics 63(3), 370–383.

[40]
Zhang, S., Y. Liu, X. Li, and G. Bi (2017).
Variational bayesian sparse signal recovery with lsm prior.
IEEE Access 5, 26690–26702.

[41]
Zhang, Z., W. Huang, Y. Liao, Z. Song, J. Shi, X. Jiang, C. Shen, and Z. Zhu
Bearing fault diagnosis via generalized logarithm sparse
Mechanical Systems and Signal Processing 167, 108576.

[42]
Zhang, Z. and B. D. Rao (2011).
Sparse signal recovery with temporally correlated source vectors
using sparse bayesian learning.
IEEE Journal of Selected Topics in Signal Processing 5(5), 912–926.

[43]
Zhou, S., Y. Chen, and J. Shi (2004).
Statistical estimation and testing for variation root-cause
identification of multistage manufacturing processes.
IEEE Transactions on Automation Science and Engineering 1(1), 73–83.

§ APPENDIX

In this section, the approximate posterior distributions of hidden variables which are $\textbf{x}, \boldsymbol{\alpha}$, and $\bar{b}$ are derived via Variational Bayes inference. In addition, the derivations of hyperparameters by maximizing the complete likelihood in Eq. (<ref>) are also described.

§.§ Inference for Eq. (<ref>)

Based on Eq. (<ref>), $\ln{q(\textbf{x}})=\langle \ln{p(\textbf{y}|\textbf{x};\lambda})p(\textbf{x}|\boldsymbol{\alpha};\text{B}) \rangle_{q(\boldsymbol{\alpha})}+const.$ Therefore,

\begin{align}
\ln{q(\textbf{x}}) &\propto \langle \ln{p(\textbf{y}|\textbf{x};\lambda})p(\textbf{x}|\boldsymbol{\alpha};\text{B}) \rangle_{q(\boldsymbol{\alpha})} \nonumber\\
& = \langle \ln{(N(\textbf{y}|\textbf{Dx},\lambda \textbf{I}_{ML})}\prod_{i=1}^{N}N(\text{x}_{i}|\textbf{0},(\alpha_{i}\text{B})^{-1}) \rangle_{q(\boldsymbol{\alpha})}  \nonumber\\
& \propto \langle -\frac{1}{2\lambda}(\textbf{y}-\textbf{Dx})(\textbf{y}-\textbf{Dx})^{\top}-\frac{1}{2}\textbf{x}[\text{AB}]\textbf{x} \rangle_{q(\boldsymbol{\alpha})} \nonumber
\end{align}

where $\text{AB}= \text{diag}[\alpha_{1}\text{B},...,\alpha_{N}\text{B}]$. Therefore, Eq. (<ref>) is derived, and $q(\textbf{x})$ follows Gaussian distribution as follows:

\begin{equation*}
q(\textbf{x})\sim N(\textbf{x}|\mu_{\textbf{x}},\Sigma_{\textbf{x}}),
\end{equation*}

where $\mu_{\textbf{x}}=\frac{1}{\lambda}\Sigma_{\textbf{x}}\textbf{D}^{\top}\textbf{y}$, $\Sigma_{\textbf{x}}=(\frac{1}{\lambda}\textbf{D}^{\top}\textbf{D}+\langle \text{AB} \rangle )^{-1}$.

§.§ Inference for Eq. (<ref>)

Based on Eq. (<ref>), $\ln{q(\boldsymbol{\alpha})}=\langle \ln{p(\textbf{x}|\boldsymbol{\alpha};\text{B})p(\boldsymbol{\alpha}|\bar{b})} \rangle_{q(\textbf{x})q(\bar{b})} +const$. Therefore,

\begin{align}
\ln{q(\boldsymbol{\alpha})} &\propto \langle \ln{p(\textbf{x}|\boldsymbol{\alpha};\text{B})p(\boldsymbol{\alpha}|\bar{b})}  \rangle_{q(\textbf{x})q(\bar{b})} \nonumber \\
&=\langle \ln{(\prod_{i=1}^{N}N(\text{x}_{i}|0,(\alpha_{i}\text{B})^{-1})  \prod_{i=1}^{N}\text{Gamma}(\alpha_{i}|a,b_{i})}\rangle_{q(\textbf{x})q(\bar{b})}  \nonumber\\
&\propto \langle \Sigma_{i=1}^{N}\ln{\alpha_{i}}^{\frac{L}{2}} - (\frac{\alpha_{i}}{2}\text{x}_{i}^{\top}\text{B}\text{x}_{i})+\ln{\alpha}_{i}^{a-1}-b_{i}\alpha_{i} \rangle_{q(\textbf{x})q(\bar{b})} \nonumber\\
& \propto \ln{\alpha}_{i}^{\frac{L}{2}+a-1}-\frac{\alpha_{i}}{2} \langle \text{x}_{i}^{\top}\text{B} \text{x}_{i} \rangle - \langle b_{i} \rangle \alpha_{i}, \nonumber
\end{align}

\begin{align}
\langle \text{x}_{i}^{\top}\text{B}\text{x}_{i} \rangle & = \langle \text{Tr}(\text{x}_{i}^{\top}\text{B}\text{x}_{i}) \rangle  \nonumber \\
& = \langle \text{Tr}(\text{x}_{i}^{\top}\text{x}_{i}\text{B}) \rangle \nonumber \\
& = \text{Tr}(\langle \text{x}_{i}^{\top}\text{x}_{i} \rangle \text{B}) \nonumber \\
&= \text{Tr}[( \Sigma_{\text{x}_{i}} + \langle \text{x}_{i} \rangle \langle \text{x}_{i}\rangle^{\top})\text{B}]. \nonumber
\end{align}

Therefore, Eq. (<ref>) is derived, and $q(\alpha_{i}), i\in P$ follows Gamma distribution as follows:

\begin{equation*}
q(\alpha_{i})=\text{Gamma}(\alpha_{i}|\frac{L}{2}+a,\frac{\text{Tr}[( \Sigma_{\text{x}_{i}} + \langle \text{x}_{i} \rangle \langle \text{x}_{i}\rangle^{\top})\text{B}]}{2} + \langle b_{i} \rangle).
\end{equation*}

For $q(\alpha_{i}), i\in P^{c}$ follows Gamma distribution as follows:

\begin{equation*}
q(\alpha_{i})=\text{Gamma}(\alpha_{i}|\frac{L}{2}+a,\frac{\text{Tr}[( \Sigma_{\text{x}_{i}} + \langle \text{x}_{i} \rangle \langle \text{x}_{i}\rangle^{\top})\text{B}]}{2} +  b_{i} ).
\end{equation*}

§.§ Inference for Eq. (<ref>)

Based on Eq. (<ref>), $\ln{q(\bar{b})}=\langle \ln{p(\boldsymbol{\alpha}|\bar{b})}p(\bar{b}) \rangle_{q(\textbf{x})q(\boldsymbol{\alpha})} +const.$ Therefore,

\begin{align}
\ln{q(\bar{b})}& \propto \langle \ln{p(\boldsymbol{\alpha}|\bar{b})}p(\bar{b}) \rangle_{q(\textbf{x})q(\boldsymbol{\alpha})} \nonumber\\
&= \langle \ln{\prod_{i\in P}\text{Gamma}(\alpha_{i}|a,b_{i})} \prod_{i\in P}\text{Gamma}(b_{i}|p,q) \rangle_{q(\textbf{x})q(\boldsymbol{\alpha})}   \nonumber \\
& \propto \Sigma_{i\in P }(-b_{i}\langle \alpha_{i} \rangle +a\ln{b_{i} + (p-1)\ln{b_{i}} -q b_{i}})  \nonumber \\
&=\Sigma_{i\in P}((p+a-1)\ln{b_{i}}-(q+\langle \alpha_{i} \rangle )b_{i}). \nonumber
\end{align}

Therefore Eq. (<ref>) is derived, and $q(b_{i}), \forall i \in P$ follows Gamma distribution as follows:

\begin{equation*}
q(b_{i})=\text{Gamma}(b_{i}|p+a,\langle \alpha_{i} \rangle+q), \forall i \in P.
\end{equation*}

§.§ Inference for Eq. (<ref>)

To estimate B, let the second term in Eq. (<ref>) as follows:

\begin{equation*}
Q(\text{B})=\langle \ln{p(\textbf{x}|\boldsymbol{\alpha};\text{B}}) \rangle_{{q(\textbf{x};\tilde{\boldsymbol{\theta}}^{OLD})q(\boldsymbol{\alpha};\tilde{\boldsymbol{\theta}}^{OLD})}}.
\end{equation*}

Let $\Gamma$=diag($\alpha_{1}^{-1},...,\alpha_{N}^{-1}$). It can be shown that

\begin{align}\label{eq:20}
\ln{p(\textbf{x}|\boldsymbol{\alpha};\text{B})} &= -\frac{1}{2}\ln{(|\Gamma^{-1}|^{L} |\text{B}^{-1} |^{N}   }) -\frac{1}{2}\textbf{x}^{\top}(\Gamma \otimes \text{B}) \textbf{x} \nonumber \\
&=-\frac{1}{2}\ln{|\Gamma^{-1}|^{L}}-\frac{1}{2}\ln{|\text{B}^{-1}|^{N}}-\frac{1}{2}\textbf{x}^{\top}(\Gamma \otimes \text{B}) \textbf{x}.
\end{align}

$Q(\text{B})$  can be calculated by taking the expectation to Eq. (<ref>). Then, taking derivative $Q(\text{B})$  with respect to $\text{B}$ leads to

\begin{equation}\label{eq:21}
\frac{\partial Q(\text{B})}{\partial \text{B}}=\frac{N}{2}\text{B}^{-1}-\frac{1}{2}\sum_{i=1}^{N}\langle \alpha_{i} \rangle (\Sigma_{\text{x}_{i}} + \langle \text{x}_{i} \rangle \langle \text{x}_{i} \rangle^{\top}).
\end{equation}

By letting Eq. (<ref>) equals to zero, B is estimated as follows:

\begin{equation*}
\text{B}= [\frac{1}{N}\sum_{i=1}^{N}\langle \alpha_{i} \rangle (\Sigma_{\text{x}_{i}} + \langle \text{x}_{i} \rangle \langle \text{x}_{i} \rangle^{\top})]^{-1}
\end{equation*}

§.§ Inference for Eq. (<ref>)

To estimate $\lambda$, let the first term in Eq. (<ref>) as follows:

\begin{equation*}
Q(\lambda)=\langle \ln{p(\textbf{y}|\textbf{x};\lambda)} \rangle_{q(\textbf{x};\tilde{\boldsymbol{\theta}}^{OLD})q(\boldsymbol{\alpha};\tilde{\boldsymbol{\theta}}^{OLD})}.
\end{equation*}

It can be shown that

\begin{align}
Q(\lambda)& \propto -\frac{ML}{2}\ln{\lambda}-\frac{1}{2\lambda}\mathbb{E}_{q(\textbf{x};\tilde{\boldsymbol{\theta}}^{OLD})q(\boldsymbol{\alpha};\tilde{\boldsymbol{\theta}}^{OLD})}\lVert \textbf{y}-\textbf{D} \textbf{x} \rVert_{2}^{2} \nonumber \\
& = -\frac{ML}{2}\ln{\lambda} -\frac{1}{2\lambda}\mathbb{E}_{q(\textbf{x};\tilde{\boldsymbol{\theta}}^{OLD})q(\boldsymbol{\alpha};\tilde{\boldsymbol{\theta}}^{OLD})}[\lVert \textbf{y}-\textbf{D}\mu_{\textbf{x}} \rVert_{2}^{2} + \text{Tr}(\Sigma_{\textbf{x}}\textbf{D}^{\top}\textbf{D})] \nonumber \\
& = -\frac{ML}{2}\ln{\lambda} -\frac{1}{2\lambda} [\lVert \textbf{y}-\textbf{D}\mu_{\textbf{x}} \rVert_{2}^{2} + \hat{\lambda} \mathbb{E}_{q(\boldsymbol{\alpha};\tilde{\boldsymbol{\theta}}^{OLD})} \text{Tr}(\Sigma_{\textbf{x}}(\Sigma_{\textbf{x}}^{-1} -\Sigma_{0}^{-1}))] \label{eq:23} \\
& = -\frac{ML}{2}\ln{\lambda} -\frac{1}{2\lambda} [\lVert \textbf{y}-\textbf{D}\mu_{\textbf{x}} \rVert_{2}^{2} + \hat{\lambda} [NL - \mathbb{E}_{q(\boldsymbol{\alpha};\tilde{\boldsymbol{\theta}}^{OLD})} \text{Tr}(\Sigma_{\textbf{x}}(\Sigma_{0}^{-1}))] \nonumber \\
& = -\frac{ML}{2}\ln{\lambda} -\frac{1}{2\lambda} [\lVert \textbf{y}-\textbf{D}\mu_{\textbf{x}} \rVert_{2}^{2} +  \hat{\lambda} [NL - \text{Tr}(\Sigma_{\textbf{x}}\mathbb{E}_{q(\boldsymbol{\alpha};\tilde{\boldsymbol{\theta}}^{OLD})}(\Sigma_{0}^{-1}))], \label{eq:24}
\end{align}

where Eq. (<ref>) follows Eq. (<ref>), and $\hat{\lambda}$ denotes the estimated $\lambda$ in the previous iteration. By setting the derivative of Eq. (<ref>) over $\lambda$ to zero, $\lambda$ can be estimated as follows:

\begin{equation*}
\lambda=\frac{[\lVert \textbf{y}-\textbf{D} \mu_{\textbf{x}} \rVert_{2}^{2} + \hat{\lambda} [NL - \text{Tr}(\Sigma_{\textbf{x}}\mathbb{E}_{q(\boldsymbol{\alpha};\tilde{\boldsymbol{\theta}}^{OLD})}(\Sigma_{0}^{-1}))]]}{ML}.
\end{equation*}

§.§ Fault pattern matrix TEXT in Section <ref>

Figure <ref> shows the fault pattern matrix $\Phi$ from previous researches [4, 20, 24] that used in Section <ref>.

Fault pattern matrix $\Phi$ in Section <ref> [2].


* Rombach et al. (2022a) Robin Rombach, Andreas Blattmann, Dominik Lorenz, Patrick Esser, and Björn Ommer.  High-resolution image synthesis with latent diffusion models.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 10684–10695, 2022a.
  * Rombach et al. (2022b) Robin Rombach, Andreas Blattmann, and Björn Ommer.  Text-guided synthesis of artistic images with retrieval-augmented diffusion models, 2022b.
  * Ryu and Ye (2022) Dohoon Ryu and Jong Chul Ye.  Pyramidal denoising diffusion probabilistic models, 2022.
  * Saharia et al. (2022) Chitwan Saharia, William Chan, Saurabh Saxena, Lala Li, Jay Whang, Emily Denton, Seyed Kamyar Seyed Ghasemipour, Burcu Karagol Ayan, S Sara Mahdavi, Rapha Gontijo Lopes, et al.  Photorealistic text-to-image diffusion models with deep language understanding.  _arXiv preprint arXiv:2205.11487_ , 2022.
  * Salimans and Ho (2022) Tim Salimans and Jonathan Ho.  Progressive distillation for fast sampling of diffusion models.  _arXiv preprint arXiv:2202.00512_ , 2022.
  * Salimans et al. (2016) Tim Salimans, Ian Goodfellow, Wojciech Zaremba, Vicki Cheung, Alec Radford, and Xi Chen.  Improved techniques for training gans.  _Advances in neural information processing systems_ , 29, 2016.
  * San-Roman et al. (2021) Robin San-Roman, Eliya Nachmani, and Lior Wolf.  Noise estimation for generative diffusion models, 2021.
  * Santos et al. (2023) Javier E Santos, Zachary R. Fox, Nicholas Lubbers, and Yen Ting Lin.  Blackout diffusion: Generative diffusion models in discrete-state spaces, 2023.
  * Sehwag et al. (2022) Vikash Sehwag, Caner Hazirbas, Albert Gordo, Firat Ozgenel, and Cristian Canton Ferrer.  Generating high fidelity data from low-density regions using diffusion models, 2022.
  * Shen et al. (2022) Max W. Shen, Ehsan Hajiramezanali, Gabriele Scalia, Alex Tseng, Nathaniel Diamant, Tommaso Biancalani, and Andreas Loukas.  Conditional diffusion with less explicit guidance via model predictive control, 2022.
  * Sheynin et al. (2022) Shelly Sheynin, Oron Ashual, Adam Polyak, Uriel Singer, Oran Gafni, Eliya Nachmani, and Yaniv Taigman.  Knn-diffusion: Image generation via large-scale retrieval, 2022.
  * Shih et al. (2023) Andy Shih, Suneel Belkhale, Stefano Ermon, Dorsa Sadigh, and Nima Anari.  Parallel sampling of diffusion models, 2023.
  * Skilling (1989) John Skilling.  The eigenvalues of mega-dimensional matrices.  In _Maximum Entropy and Bayesian Methods_ , pages 455–466. Springer, 1989.
  * Sohl-Dickstein et al. (2015) Jascha Sohl-Dickstein, Eric Weiss, Niru Maheswaranathan, and Surya Ganguli.  Deep unsupervised learning using nonequilibrium thermodynamics.  In _International Conference on Machine Learning_ , pages 2256–2265. PMLR, 2015.
  * Song et al. (2020a) Jiaming Song, Chenlin Meng, and Stefano Ermon.  Denoising diffusion implicit models.  _arXiv preprint arXiv:2010.02502_ , 2020a.
  * Song and Ermon (2019) Yang Song and Stefano Ermon.  Generative modeling by estimating gradients of the data distribution.  _Advances in Neural Information Processing Systems_ , 32, 2019.
  * Song and Ermon (2020) Yang Song and Stefano Ermon.  Improved techniques for training score-based generative models.  _Advances in neural information processing systems_ , 33:12438–12448, 2020.
  * Song and Kingma (2021) Yang Song and Diederik P Kingma.  How to train your energy-based models.  _arXiv preprint arXiv:2101.03288_ , 2021.
  * Song et al. (2020b) Yang Song, Sahaj Garg, Jiaxin Shi, and Stefano Ermon.  Sliced score matching: A scalable approach to density and score estimation.  In _Uncertainty in Artificial Intelligence_ , pages 574–584. PMLR, 2020b.
  * Song et al. (2020c) Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole.  Score-based generative modeling through stochastic differential equations.  _arXiv preprint arXiv:2011.13456_ , 2020c.
  * Song et al. (2021) Yang Song, Conor Durkan, Iain Murray, and Stefano Ermon.  Maximum likelihood training of score-based diffusion models.  _Advances in Neural Information Processing Systems_ , 34:1415–1428, 2021.
  * Song et al. (2023) Yang Song, Prafulla Dhariwal, Mark Chen, and Ilya Sutskever.  Consistency models.  2023\.
  * Stein et al. (2023) George Stein, Jesse C. Cresswell, Rasa Hosseinzadeh, Yi Sui, Brendan Leigh Ross, Valentin Villecroze, Zhaoyan Liu, Anthony L. Caterini, J. Eric T. Taylor, and Gabriel Loaiza-Ganem.  Exposing flaws of generative model evaluation metrics and their unfair treatment of diffusion models, 2023.
  * Su et al. (2022) Xuan Su, Jiaming Song, Chenlin Meng, and Stefano Ermon.  Dual diffusion implicit bridges for image-to-image translation.  _arXiv preprint arXiv:2203.08382_ , 2022.
  * Sun et al. (2023) Wujie Sun, Defang Chen, Can Wang, Deshi Ye, Yan Feng, and Chun Chen.  Accelerating diffusion sampling with classifier-based feature distillation.  In _2023 IEEE International Conference on Multimedia and Expo (ICME)_ , pages 810–815. IEEE, 2023.
  * Szegedy et al. (2016) Christian Szegedy, Vincent Vanhoucke, Sergey Ioffe, Jon Shlens, and Zbigniew Wojna.  Rethinking the inception architecture for computer vision.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 2818–2826, 2016.
  * Tabak and Turner (2013) Esteban G Tabak and Cristina V Turner.  A family of nonparametric density estimation algorithms.  _Communications on Pure and Applied Mathematics_ , 66(2):145–164, 2013.
  * Tabak and Vanden-Eijnden (2010) Esteban G Tabak and Eric Vanden-Eijnden.  Density estimation by dual ascent of the log-likelihood.  _Communications in Mathematical Sciences_ , 8(1):217–233, 2010.
  * Tachibana et al. (2021) Hideyuki Tachibana, Mocho Go, Muneyoshi Inahara, Yotaro Katayama, and Yotaro Watanabe.  It$\backslash$^$\\{$o$\\}$-taylor sampling scheme for denoising diffusion probabilistic models using ideal derivatives.  _arXiv preprint arXiv:2112.13339_ , 2021.
  * Ulhaq et al. (2022) Anwaar Ulhaq, Naveed Akhtar, and Ganna Pogrebna.  Efficient diffusion models for vision: A survey, 2022.
  * Vahdat et al. (2021) Arash Vahdat, Karsten Kreis, and Jan Kautz.  Score-based generative modeling in latent space.  _Advances in Neural Information Processing Systems_ , 34:11287–11302, 2021.
  * Van Den Oord et al. (2016) Aäron Van Den Oord, Nal Kalchbrenner, and Koray Kavukcuoglu.  Pixel recurrent neural networks.  In _International conference on machine learning_ , pages 1747–1756. PMLR, 2016.
  * Van Den Oord et al. (2017) Aaron Van Den Oord, Oriol Vinyals, et al.  Neural discrete representation learning.  _Advances in neural information processing systems_ , 30, 2017.
  * Vargas et al. (2023) Francisco Vargas, Will Grathwohl, and Arnaud Doucet.  Denoising diffusion samplers, 2023.
  * Verma et al. (2023) Prakhar Verma, Vincent Adam, and Arno Solin.  Variational gaussian process diffusion processes, 2023.
  * Vincent (2011) Pascal Vincent.  A connection between score matching and denoising autoencoders.  _Neural computation_ , 23(7):1661–1674, 2011.
  * Voleti et al. (2022) Vikram Voleti, Christopher Pal, and Adam Oberman.  Score-based denoising diffusion with non-isotropic gaussian noise models, 2022.
  * Wallace et al. (2023) Bram Wallace, Akash Gokul, Stefano Ermon, and Nikhil Naik.  End-to-end diffusion latent optimization improves classifier guidance, 2023.
  * Watson et al. (2021a) Daniel Watson, William Chan, Jonathan Ho, and Mohammad Norouzi.  Learning fast samplers for diffusion models by differentiating through sample quality.  In _International Conference on Learning Representations_ , 2021a.
  * Watson et al. (2021b) Daniel Watson, Jonathan Ho, Mohammad Norouzi, and William Chan.  Learning to efficiently sample from diffusion probabilistic models.  _arXiv preprint arXiv:2106.03802_ , 2021b.
  * Wu and la Torre (2022) Chen Henry Wu and Fernando De la Torre.  Unifying diffusion models’ latent space, with applications to cyclediffusion and guidance, 2022.
  * Wu et al. (2023) Zike Wu, Pan Zhou, Kenji Kawaguchi, and Hanwang Zhang.  Fast diffusion model, 2023.
  * Xiang et al. (2023) Chendong Xiang, Fan Bao, Chongxuan Li, Hang Su, and Jun Zhu.  A closer look at parameter-efficient tuning in diffusion models, 2023.
  * Xiao et al. (2021) Zhisheng Xiao, Karsten Kreis, and Arash Vahdat.  Tackling the generative learning trilemma with denoising diffusion gans.  _arXiv preprint arXiv:2112.07804_ , 2021.
  * Xu et al. (2022) Yilun Xu, Ziming Liu, Max Tegmark, and Tommi Jaakkola.  Poisson flow generative models.  _Advances in Neural Information Processing Systems_ , 35:16782–16795, 2022.
  * Xu et al. (2023a) Yilun Xu, Mingyang Deng, Xiang Cheng, Yonglong Tian, Ziming Liu, and Tommi Jaakkola.  Restart sampling for improving generative processes.  _arXiv preprint arXiv:2306.14878_ , 2023a.
  * Xu et al. (2023b) Yilun Xu, Ziming Liu, Yonglong Tian, Shangyuan Tong, Max Tegmark, and Tommi Jaakkola.  Pfgm++: Unlocking the potential of physics-inspired generative models.  _arXiv preprint arXiv:2302.04265_ , 2023b.
  * Xu et al. (2023c) Yilun Xu, Shangyuan Tong, and Tommi Jaakkola.  Stable target field for reduced variance score estimation in diffusion models, 2023c.
  * Yang et al. (2022a) Ling Yang, Zhilong Zhang, Yang Song, Shenda Hong, Runsheng Xu, Yue Zhao, Yingxia Shao, Wentao Zhang, Bin Cui, and Ming-Hsuan Yang.  Diffusion models: A comprehensive survey of methods and applications.  _arXiv preprint arXiv:2209.00796_ , 2022a.
  * Yang et al. (2022b) Xiulong Yang, Sheng-Min Shih, Yinlin Fu, Xiaoting Zhao, and Shihao Ji.  Your vit is secretly a hybrid discriminative-generative diffusion model, 2022b.
  * Yi et al. (2023) Mingyang Yi, Jiacheng Sun, and Zhenguo Li.  On the generalization of diffusion model, 2023.
  * Zhang et al. (2023a) Chenshuang Zhang, Chaoning Zhang, Mengchun Zhang, and In So Kweon.  Text-to-image diffusion model in generative ai: A survey.  _arXiv preprint arXiv:2303.07909_ , 2023a.
  * Zhang et al. (2023b) Chenshuang Zhang, Chaoning Zhang, Sheng Zheng, Mengchun Zhang, Maryam Qamar, Sung-Ho Bae, and In So Kweon.  A survey on audio diffusion models: Text to speech synthesis and enhancement in generative ai.  _arXiv preprint arXiv:2303.13336_ , 2, 2023b.
  * Zhang et al. (2023c) Han Zhang, Ruili Feng, Zhantao Yang, Lianghua Huang, Yu Liu, Yifei Zhang, Yujun Shen, Deli Zhao, Jingren Zhou, and Fan Cheng.  Dimensionality-varying diffusion process.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 14307–14316, 2023c.
  * Zhang and Chen (2021) Qinsheng Zhang and Yongxin Chen.  Diffusion normalizing flow.  _Advances in Neural Information Processing Systems_ , 34:16280–16291, 2021.
  * Zhang and Chen (2022) Qinsheng Zhang and Yongxin Chen.  Fast sampling of diffusion models with exponential integrator.  _arXiv preprint arXiv:2204.13902_ , 2022.
  * Zhang et al. (2022) Qinsheng Zhang, Molei Tao, and Yongxin Chen.  gddim: Generalized denoising diffusion implicit models.  _arXiv preprint arXiv:2206.05564_ , 2022.
  * Zheng et al. (2023) Hongkai Zheng, Weili Nie, Arash Vahdat, Kamyar Azizzadenesheli, and Anima Anandkumar.  Fast sampling of diffusion models via operator learning.  In _International Conference on Machine Learning_ , pages 42390–42402. PMLR, 2023.
  * Zheng and Zhou (2020) Huangjie Zheng and Mingyuan Zhou.  Act: Asymptotic conditional transport.  2020\.
  * Zheng et al. (2022) Huangjie Zheng, Pengcheng He, Weizhu Chen, and Mingyuan Zhou.  Truncated diffusion probabilistic models.  _stat_ , 1050:7, 2022.
  * Zou et al. (2023) Hao Zou, Zae Myung Kim, and Dongyeop Kang.  Diffusion models in nlp: A survey.  _arXiv preprint arXiv:2305.14671_ , 2023.