混合概率PCA：一种用于非线性过程监控的故障诊断方法

混合概率PCA （mixture of probabilistic principal component analysis, MPPCA）

1. 引言
2. 潜变量模型和PCA
3. 概率PCA
4. 混合概率PCA
5. 局部模型数量确定

引言

作为一个坚守主成分分析(principal component analysis, PCA) 的学渣，虽然了解大量关于PCA的拓展方法，但依然觉得PCA还是有许多可以研究的内容。本文浅谈一种改进的混合概率PCA 方法，将概率PCA拓展到非线性异常监测中。此外，该方法也适合于多工况异常监测，工况可自动辨识。
本文参考：Zhang, J., Chen, H., Chen, S., & Hong, X. (2019). An Improved Mixture of Probabilistic PCA for Nonlinear Data-Driven Process Monitoring. IEEE Transactions on Systems, Man, and Cybernetics, 49(1), 198–210.
正文链接
预印版链接

潜变量模型和PCA

考虑下列模型
$\mathbf{t}=\mathbf{y}(\mathbf{x};\mathbf{w})+\mathbf{\xi }$,
其中，观测数据为$\mathbf{t}\in {\Re }^d$, 潜变量为$\mathbf{x}\in {\Re }^q$，$\mathbf{w}$为模型参数，$\mathbf{\epsilon }$为独立噪声。线性化模型为：
$\mathbf{t}=\mathbf{W}\mathbf{x}+\mathbf{\mu }+\mathbf{\xi }$。
假设，$\mathbf{x}\sim \mathrm{N}(0,\mathbf{I})$, $\mathbf{\xi }\sim \mathrm{N}(0,\mathbf{\Psi })$， $\mathbf{\Psi }\in {\Re }^{d\times d}$为对角矩阵，$\mathbf{\mu }\in {\Re }^d$为输出均值向量，$\mathbf{W}\in {\Re }^{d\times q}$为载荷矩阵。 在高斯假设下，上述线性模型产生的观测数据服从高斯分布，即 $\mathbf{t}\sim \mathrm{N}(\mathbf{\mu },\mathbf{C})$ ， $\mathbf{C}=\mathbf{\Psi }+\mathbf{W}{\mathbf{W}}^{\mathrm{T}}\in {\Re }^{d\times d}$.

观测数据集 { t n } n = 1 N \{{\boldsymbol t}_n \}_{n=1}^N {tn​}n=1N​，线性模型的矩阵表达为

$\mathbf{T}=\mathbf{W}\mathbf{X}+\mathbf{u}{1}^{\mathrm{T}}+\mathbf{\Xi }$,

其中， $\mathbf{T}=[{\mathbf{t}}_1,...,{\mathbf{t}}_N]\in {\Re }^{d\times N}$, $\mathbf{X}=[{\mathbf{x}}_1,...{\mathbf{x}}_N]\in {\Re }^{q\times N}$, $\mathbf{\Xi }=[{\mathbf{\xi }}_1,...{\mathbf{\xi }}_N]\in {\Re }^{d\times N}$。均值为$\mathbf{u}=\frac{1}{N}{\sum }_{n=1}^N{\mathbf{t}}_n$，协方差矩阵$\mathbf{S}=\frac{1}{N}(\mathbf{T}-\mathbf{u}{1}^{\mathrm{T}})(\mathbf{T}-\mathbf{u}{1}^{\mathrm{T}}{)}^{\mathrm{T}}$。

概率PCA

PCA的主要目标是寻找到低维空间的投影向量，同时保证重构误差最小。
给定噪声$\mathbf{\epsilon }\sim \mathrm{N}(0,{\sigma }^2\mathbf{I})$, 则通过以下公式揭示特定$\mathbf{x}$在$\mathbf{t}$-空间上的概率分布
$p(\mathbf{t}|\mathbf{x})=(2\pi {\sigma }^2{)}^{-d/2}\exp \left\{-\frac{1}{2{\sigma }^2}{\Vert \mathbf{t}-\mathbf{W}\mathbf{x}-\mathbf{\mu }\Vert }^2\right\}$
$\mathbf{x}$上的高斯先验概率可以定义为

$p(\mathbf{x})=(2\pi {)}^{-q/2}\exp \left\{-\frac{1}{2}{\mathbf{x}}^{\mathrm{T}}\mathbf{x}\right\}$

$\mathbf{t}$ 的边缘分布为

$\begin{array}{l}p(\mathbf{t})=\int p(\mathbf{t}|\mathbf{x})p(\mathbf{x})d\mathbf{x}\\ =(2\pi {)}^{-d/2}{\left|\mathbf{C}\right|}^{-1/2}\exp \left\{-\frac{1}{2}{(\mathbf{t}-\mathbf{\mu })}^{\mathrm{T}}{\mathbf{C}}^{-1}(\mathbf{t}-\mathbf{\mu })\right\}\end{array}$

其中. 模型协方差为 $\mathbf{C}={\sigma }^2\mathbf{I}+\mathbf{W}{\mathbf{W}}^{\mathrm{T}}$
在 Bayesian理论下， 给定$\mathbf{t}$, $\mathbf{x}$ 的后验分布为

$p(\mathbf{x}|\mathbf{t})=\mathrm{e}\mathrm{x}\mathrm{p}\{-\frac{1}{2}{\left\{\mathbf{x}-{\mathbf{M}}^{-1}{\mathbf{W}}^{\mathrm{T}}(\mathbf{t}-\mathbf{\mu })\right\}}^{\mathrm{T}}({\sigma }^{-2}\mathbf{M})\left\{\mathbf{x}-{\mathbf{M}}^{-1}{\mathbf{W}}^{\mathrm{T}}(\mathbf{t}-\mathbf{\mu })\right\}\}\times (2\pi {)}^{-\mathrm{q}/2}{\left|{\sigma }^{-2}\mathbf{M}\right|}^{1/2}$

其中，后验协方差为 ${\sigma }^2{\mathbf{M}}^{-1}={\sigma }^2({\sigma }^2\mathbf{I}+{\mathbf{W}}^{\mathrm{T}}\mathbf{W}{)}^{-1}$. $\mathbf{M}\in {\Re }^{q\times q}$, $\mathbf{C}\in {\Re }^{d\times d}$ .
对数似然函数为
L = ∑ n = 1 N ln ⁡ { p ( t n ) }      = − N 2 { d ln ⁡ ( 2 π ) + ln ⁡ ∣ C ∣ + t r ( C − 1 S ) } \begin{array}{l} L = \sum\limits_{n = 1}^N {\ln \left\{ {p\left( {{{\boldsymbol t_n}}} \right)} \right\}} \\ \;\;{\kern 1pt} = - \frac{N}{2}\left\{ {d\ln \left( {2\pi } \right) + \ln \left|{\boldsymbol C } \right| + tr\left( {{{\boldsymbol C}^{ - 1}}{\boldsymbol S}} \right)} \right\} \end{array} L=n=1∑N​ln{p(tn​)}=−2N​{dln(2π)+ln∣C∣+tr(C−1S)}​

当$\mathbf{W}$的列跨越数据的主子空间时，对数似然最大。解析解也可以通过$\mathbf{S}$的特征分解和噪声方差${\sigma }^2$的估计（基于$\mathbf{S}$的最小特征值）获得。或者，可以使用EM算法的迭代方法来生成以下完整数据对数似然

L c = ∑ n = 1 N ln ⁡ { p ( t n , x n ) } = ∑ n = 1 N ln ⁡ { ( 2 π σ 2 ) − d / 2 exp ⁡ { − 1 2 σ 2 ∥ t n − W x n − μ ∥ 2 } ( 2 π ) − q / 2 exp ⁡ { − 1 2 x n T x n } } L_c=\sum_{n=1}^{N} \ln\left\{p({\boldsymbol t}_n, {\boldsymbol x}_n ) \right\}\\ = \sum_{n=1}^{N} \ln\left\{{(2\pi {\sigma ^2})^{ - d/2}}\exp \left\{ { - \frac{1}{{2{\sigma ^2}}}{{\left\| {\boldsymbol t}_n - {\boldsymbol W} {\boldsymbol x}_n - {\boldsymbol \mu} \right\|}^2}} \right\} \right. \left. {(2\pi )^{ - q/2}}\exp \left\{ { - \frac{1}{2} {\boldsymbol x}_n^{\rm T} {\boldsymbol x}_n} \right\}\right\} Lc​=∑n=1N​ln{p(tn​,xn​)}=∑n=1N​ln{(2πσ2)−d/2exp{−2σ21​∥tn​−Wxn​−μ∥2}(2π)−q/2exp{−21​xnT​xn​}}

在EM算法中，模型的后验均值和协方差为

$⟨{\mathbf{x}}_n⟩={\mathbf{M}}^{-1}{\mathbf{W}}^{\mathrm{T}}\left({\mathbf{t}}_n-\mathbf{\mu }\right)$

$⟨{\mathbf{x}}_n{\mathbf{x}}_n^{\mathrm{T}}⟩={\sigma }^2{\mathbf{M}}^{-1}+⟨{\mathbf{x}}_n⟩{⟨{\mathbf{x}}_n⟩}^{\mathrm{T}}$

混合概率PCA

为了能够对更复杂的数据进行建模，引入了混合概率PCA（MPPCA。通过在K个局部PCA模型预测输出，它为处理非线性和缺失数据提供了更强大的基础。根据概率规则，考虑K个局部PCA模型的混合，而不是用单一模型来表示系统：

$p(\mathbf{t})={\sum }_{i=1}^{K}p(i)p(\mathbf{t}|i)={\sum }_{i=1}^K{\pi }_{i}p(\mathbf{t}|i))$

约束为 ${\pi }_i\ge 0$，$\sum {\pi }_i=1$。$p(i)$ 为选择第$i$个模型的概率，每个 $p(\mathbf{t}|i)$，PCA模型为
$\mathbf{t}={\mathbf{W}}_i\mathbf{x}+{\mathbf{\mu }}_i+{\mathbf{\xi }}_i,i=1,...,K$。 每个局部PCA模型的后验均值和协方差为

$⟨{\mathbf{x}}_n^{(i)}⟩={\mathbf{M}}_i^{-1}{\mathbf{W}}_i^{\mathrm{T}}\left({\mathbf{t}}_n-{\mathbf{\mu }}_i\right)\in {\Re }^q$

$⟨{\mathbf{x}}_n^{(i)}{{\mathbf{x}}_n^{(i)}}^{\mathrm{T}}⟩={\sigma }_i^2{\mathbf{M}}_i^{-1}+⟨{\mathbf{x}}_n^{(i)}⟩{⟨{\mathbf{x}}_n^{(i)}⟩}^{\mathrm{T}}\in {\Re }^{q\times q}$

局部模型数量确定

各局部模型参数受到局部模型数量K影响，本文用Baye Ying-Yang 确定最优K，准则如下：
$K^{\circ }=\arg \underset{i}{\min }H\left(i\right)$
$H\left(i\right)\equiv -\frac{1}{N}\sum _{n=1}^N\sum _{i=1}^{K}p\left(i|{\mathbf{t}}_n,\mathbf{\theta }\right)\ln \left(p\left({\mathbf{t}}_n|i\right)\right)-\sum _{i=1}^K{\pi }_i\ln {\pi }_i$

其中，$\mathbf{\theta }$ 包含模型的所有参数，即 $\mathbf{\theta }={\left\{{\mathbf{W}}_i,{\mathbf{\mu }}_i,{\sigma }_i^2,{\pi }_i\right\}}_{i=1,...,K}$

在MPPCA 中，K 和 $\mathbf{\theta }$ 同时优化，用最大期望方法求解参数。