基于自适应Lasso阈值的EWMA协方差矩阵控制图

An Adaptive Lasso Thresholding-Based Covariance Matrix EWMA Control Chart

摘要: 许多工业多元质量控制应用中通常过程失控时偏移只发生在协方差矩阵中的小部分元素中，即协方差矩阵的偏移具有稀疏性。本文针对这一现实提出一种基于自适应lasso阈值的EWMA (E_ALT)控制图，并通过蒙特卡洛模拟得到该控制图在不同偏移量下的平均运行链长指标，实验结果证明其在中小偏移下监控效率较高。之后进一步对控制图中包含的参数进行优化，以提高控制图性能。

Abstract: In many industrial multivariate quality control applications, when the process is out of control, the offset usually occurs only in a small part of the covariance matrix, which calls that the offset of the covariance matrix is sparse. In view of this reality, this paper proposes an EWMA(E_ALT) control chart based on adaptive lasso threshold, and obtains the average running chain length index of the control chart at different offsets through Monte Carlo simulation, and the experimental results prove that its monitoring efficiency is high under medium and small offsets. The parameters contained in the chart are then further optimized to improve chart performance.

1. 引言

统计质量控制目前广泛用于生产制造过程和服务产品的质量改进中，其中控制图是最强大的技术，其目的是监控一个正在进行的过程中相互关联的特征情况。在相关工业制造业和服务业中，大多数产品的质量需要多个相关变量的联合结果来描述，随着数据采集系统和计算技术的进步，多元统计过程控制图在监控和改进制造过程中可以发挥更大作用 [1] 。

随着SPC理论体系的成熟，该技术也被广泛应用于信号处理、网络安全、智能制造等领域，在相关应用中往往会出现监控变量特征数远大于获得子组样本数的情况，此时许多传统监控方法由于矩阵的奇异性或接近奇异而无法使用。许多学者针对该问题提出改进方案，如研究在基于单个观测值的情况下设计用于监控协方差矩阵的多元控制图 [2] [3] [4] ，并以实例验证控制图在实际应用中有很好的表现；Li和Wang [5] 通过使用惩罚似然函数得到的估计值替换未知的协方差矩阵后进行假设检验设计一个新的多元协方差矩阵控制图。此外，基于对一些工程和操作原理的理解，在许多工业多元质量控制应用中通常发生失控时偏移只发生在协方差矩阵中的小部分元素中，即协方差矩阵的偏移具有稀疏性。对此Maboudou-Tchao和Diawara [6] 利用Lasso方法求解协方差矩阵逆矩阵的稀疏估计过程设计多元控制图进行协方差矩阵的监控，该控制图被证明在仅方差偏移和方差与相关性同时偏移情况下效果较好；最近，Galal和Jinho [7] 提出了可以使用阈值法估计稀疏协方差矩阵作为正则化估计方法的潜在替代方案，设计出通过自适应lasso阈值估计稀疏矩阵的Shewhart控制图，并验证了该控制图在某些偏移情况下表现优于基于似然比的控制图，然而当监控过程偏移量较小时，其效果并不是十分理想，为此本文将原有的Shewhart控制图推广到EWMA控制图，并重新对各类偏移情况进行充分的模拟研究。

2. 基于自适应Lasso阈值的EWMA协方差控制图

设过程观测值 $X={\left(X_1,X_2,\cdots ,X_p\right)}^{\prime }$ 为p维随机向量，独立同分布于多元正态分布，其均值向量和协

方差矩阵分别为 $E\left(X\right)={\mu }_X={\left({\mu }_{X_1},{\mu }_{X_2},\cdots ,{\mu }_{X_p}\right)}^{\prime }$ 和 $\mathrm{cov}\left(X\right)={\Sigma }_X=E\left[\left(X-{\mu }_X\right){\left(X-{\mu }_X\right)}^{\prime }\right]$ ，其中 ${\Sigma }_X$ 为

$p\times p$ 维对称正定矩阵。假定过程受控时过程均值向量为 ${\mu }_0$ ，协方差矩阵为 ${\Sigma }_0$ ，且总体参数已知。对于监控阶段，以 $X_t$ 表示在多元过程中随机抽取的第 $t\left(t=1,2,3,\cdots \right)$ 个子组，每个子组包含n个独立的样品，则 $X_{ti}={\left(X_{ti1},X_{ti2},\cdots ,X_{tip}\right)}^{\prime }$ 表示第t个子组的第 $i\left(i=1,2,3,\cdots ,n\right)$ 观测值，其中 $X_{tij}$ 表示t时刻第i个观测值的第 $j\left(j=1,2,3,\cdots ,p\right)$ 个指标值。

首先以样本协方差矩阵 $S_t=\frac{1}{n-1}{\displaystyle {\sum }_{i=1}^n\left(X_{ti}-{\bar{X}}_t\right)}{\left(X_{ti}-{\bar{X}}_t\right)}^{\prime }$ 估计过程协方差矩阵，基于稀疏性假设，过程失控时差异矩阵 $\Sigma -{\Sigma }_0$ 中仅有少数元素不为0。记差异矩阵 $D_t={\left(d_{ijt}\right)}_{1\le i,j\le p}=S_t-{\Sigma }_0$ ，则当过程受控且样本量足够大时 $D_t$ 趋于0矩阵。为此我们采取自适应阈值法对差异矩阵进行压缩。通常统计质量控制在监控过程时假设受控协方差矩阵已知或通过阶段I估计得到，记为 ${\Sigma }_0={\left({\sigma }_{ij}\right)}_{p\times p}$ ，若样本量为n，正则化参数 $\phi$ ，可得到阈值函数

$g_{ij}=\phi \sqrt{\left({\sigma }_{ii}{\sigma }_{jj}+{\sigma }_{ij}^2\right)n^{-1}\log p}$ (1)

$G={\left(g_{ij}\right)}_{p\times p}$ 中每个元素为估计差异矩阵 $D_t$ 对应元的阈值，若 $d_{ijt}$ 的绝对值小于对应阈值，则将其压缩为0。采用连续变换将 $d_{ijt}$ 调整为

${\tilde{d}}_{ijt}=d_{ijt}\cdot \max \left(0,1-{\left|\frac{g_{ij}}{d_{ijt}}\right|}^{\eta }\right)$ (2)

其中 $\eta$ 为压缩参数。

同时对矩阵 ${\tilde{D}}_t={\left({\tilde{d}}_{ijt}\right)}_{{}_{1\le i,j\le p}}$ 中各元素进行标准化处理，得到 ${\tilde{d}}_{ijt}^{*}={\tilde{d}}_{ijt}/\sqrt{n{\theta }_{ij}}$ ，其中 ${\theta }_{ij}$ 为元素 ${\tilde{d}}_{ijt}$ 的方差。

检验一个矩阵是否为0矩阵的常用方法是计算其Frobenius范数并检验其是否为0，进而可得到矩阵 ${\tilde{D}}_t^{*}$ 的F范数为

$F_t={\Vert {\tilde{D}}_t^{*}\Vert }_F^2={\displaystyle \underset{i=1}{\overset{p}{\sum }}{\displaystyle \underset{j=1}{\overset{p}{\sum }}{\left({\tilde{d}}_{ijt}^{*}\right)}^2}}=\frac{\text{1}}{n}{\displaystyle \underset{i=1}{\overset{p}{\sum }}{\displaystyle \underset{j=1}{\overset{p}{\sum }}\frac{{\tilde{d}}_{ijt}^2}{{\sigma }_{ii}{\sigma }_{jj}+{\sigma }_{ij}^2}}}$ (3)

对其进行指数加权平均，得到E_ALT控制图的监控统计量。记 $ME{F}_1=F_1$ ， $t>1$ 时

$ME{F}_t=\left(1-\lambda \right)ME{F}_{t-1}+\lambda F_t$ (4)

若给定参数 $\lambda$ ( $0<\lambda <1$ )，给定无偏移协方差矩阵 ${\Sigma }_0$ ，给定无偏移时平均链长 $L_0$ ，进一步确定控制图上控制限。利用 $ME{F}_t$ 的经验分布函数求控制限

$F_N\left(u\right)=\frac{诸ME{F}_t中不超过u的个数}{N}$ ，

$H=F_N^{-1}\left(1-\alpha \right)=F_N^{-1}\left(1-1/L_0\right)$ (5)

记m为 $\left(1-\alpha \right)N$ 四舍五入之后的整数，将 $ME{F}_1,ME{F}_2,\cdots ,ME{F}_N$ 由小到大排序，则

$H=ME{F}_{\left(m\right)}$ (6)

3. 控制图的性能比较

在统计质量控制领域，比较控制图的性能优劣通常使用平均运行链长(Average Run Length，记为ARL)作为指标，即从监控开始至其发出报警信号为止抽取的平均样本组数。我们希望受控时的ARL(记为ARL₀)越大越好，而失控时的ARL(记为ARL₁)越小越好。通常比较控制图好坏的具体方法是：固定二者的受控ARL₀，之后比较失控时的ARL₁，ARL₁越小则表示该控制图灵敏度更高。

本文主要研究高维情况下的控制图性能，故取样本容量 $n=15$ ，变量维数 $p=15$ ，在固定受控平均运行链长为200的前提下，通过随机模拟得到各控制图在不同平滑系数以及不同偏移情况下的失控平均运行链长，从而进行性能比较。事实上对于协方差矩阵的监控，我们无法罗列所有失控情况，因为在p维的协方差阵中包含任意 $p\left(p+1\right)/2$ 种偏移。因此本文不仅基于偏移量大小比较控制图性能，同时也会考虑偏移位置，即按照仅有对角线元素元素偏移，仅有非对角线元素发生偏移，以及对角线元素和非对角线元素共同发生偏移进行控制图性能比较。具体偏移情况分类如下：

1) 仅有对角线元素发生偏移

$O{C}_1:{\sigma }_{ii}^{OC}:=1+\delta ,1\le i\le p$

$O{C}_2:{\sigma }_{11}^{OC}:=1+\delta$

2) 仅有非对角线元素发生偏移

$O{C}_3:{\sigma }_{ij}^{OC}:=\delta ,1\le i,j\le k\le p,i\ne j$

$O{C}_4:{\sigma }_{ij}^{OC}:=\delta ,1\le i,j\le 2\le p,i\ne j$

3) 对角线元素和非对角线元素同时偏移

$O{C}_5:{\sigma }_{ij}^{OC}:={\sigma }_{ij}^{IC}+\delta ,1\le i,j\le k\le p$

$O{C}_6:{\sigma }_{ij}^{OC}:=\delta ,{\sigma }_{ii}^{OC}:=1+{\delta }^2,1\le i,j\le 2\le p$

模拟过程中对第三种偏移情况的k取5，对第5种偏移情况的k取3。

当用自适应lasso阈值估计稀疏矩阵时，存在正则化参数 $\phi$ 和压缩参数 $\eta$ ，本文初步选取通过交叉验证法得到的最优参数 $\phi =1$ ， $\eta =4$ ，固定参数的情况下，研究偏移量 $\delta$ 从0到1时的各控制图性能，为了减少模拟过程中的随机性，对每一个偏移量进行10,000次模拟，得到评价指标平均运行链长，具体结果见图1~3。

图1中为只考虑对角线元素发生变化时各控制图的ARL₁。首先由图1分析可知这两类控制图均随着偏移量 $\delta$ 的增大灵敏度有所提升。在 ${\Sigma }_{OC1}$ 偏移情况下，随着 $\lambda$ 的减小，E_ALT控制图也更为灵敏。当所有对角线元素都发生偏移时，E_ALT控制图性能与ALT_norm控制图较为相似；而对于 ${\Sigma }_{OC2}$ 偏移情况，E_ALT具有一致优越性，对于小的偏移效果更为明显，且 $\lambda$ 越小，控制图灵敏性越高。

Figure 1. ARL comparison of E_ALT and ALT_norm control chart when shift in diagonal elements (ARL₀ = 200, $n=15$ , $p=15$ )

图1. 对角线元素发生偏移时E_ALT和ALT_norm控制图的ARL比较(ARL₀ = 200, $n=15$ , $p=15$ )

图2为仅有非对角线元素发生偏移时各控制图的性能表现，分析可知当协方差矩阵中对角线元素不变，部分非对角线元素发生偏移时，偏移量 $\delta <0.4$ 时，ALT_norm控制图性能较好，而当偏移量较大时，E_ALT控制图比ALT_norm控制图更灵敏，且调节参数越小E_ALT控制图的性能越好。而对于协方差阵中仅有一组非对角线元素发生偏移，当 $\delta <1$ 时E_ALT控制图ARL₁小于其他控制图的ARL₁，即E_ALT对于中小偏移更灵敏。

Figure 2. ARL comparison of E_ALT and ALT_norm control chart when shift in off-diagonal elements (ARL₀ = 200, $n=15$ , $p=15$ )

图2. 仅有非对角线元素发生偏移时E_ALT和ALT_norm控制图的ARL比较(ARL₀ = 200, $n=15$ , $p=15$ )

图3中则为考虑变量对角线元素和非对角线元素均发生偏移的情况，对于样本对角线元素和非对角线元素均发生偏移的情况，由图3可知当偏移量 $\delta <0.4$ 时，E_ALT控制图更加灵敏，而 $\delta \ge 0.6$ 时，较大的调节参数 $\lambda$ 下E_ALT控制图的ARL₁与ALT_norm控制图的ARL₁较为接近，也就是说在选取合适的 $\lambda$ 值后，对于固定的 $\delta$ 值E_ALT控制图具有一致的优越性。

Figure 3. ARL comparison of E_ALT and ALT_norm control chart when shift in diagonal and off-diagonal elements (ARL₀ = 200, $n=15$ , $p=15$ )

图3. 对角线元素及非对角线元素同时发生偏移时E_ALT和ALT_norm控制图的ARL比较(ARL₀ = 200, $n=15$ , $p=15$ )

4. 参数优化

在使用自适应Lasso阈值法估计矩阵时，存在正则化参数 $\phi$ 与压缩参数 $\eta$ ，正则化参数 $\phi$ 表示由(2)式求出阈值对应方差的正则化水平，压缩参数 $\eta$ 则决定着压缩过程中的压缩程度。前文使用通过交叉验证法求得的参数进行控制图设计，然而从多元统计过程控制角度来看，参数的选择也与实际偏移情况有关，因此本文进一步按照偏移量情况，对各参数进行讨论。具体结果见表1和表2。

表1是对正则化参数分析结果，加粗的部分是ARL1最小的情况。对于第一种偏移情况，当偏移量较小，小于0.4时，正则化参数取0.5控制图具有最佳的ARL1，而偏移量达到0.6之后，正则化参数的取值对控制图性能基本没有影响。对于第二、第四及第六种偏移情况，偏移量低于0.2时，越小的正则化参数使得控制图效果越好，而对于较大偏移量则正则化参数越大越好，因为这三种偏移情况相较于其他情况稀疏性更强。对于第三、第五种失控情况，正则化参数取较小值时控制图整体表现更佳。而表2展示了压缩参数变化时的对应指标，正如上文所说，若偏移量较小，则更大的压缩参数对应的控制图效果更好。综上，若某个系统可从历史数据中获得一些相关情况，则可以按照分析预先设定参数，实现更高效的监控。

Table 1. ARL performance of E_ALT according to φ (ARL0 = 200, n = 15 , p = 15 , λ = 0.5 , η = 4 )

表1. 关于正则化参数 $\phi$ 的E_ALT控制图ARL表现(ARL₀ = 200, $n=15$ , $p=15$ , $\lambda =0.5$ , $\eta =4$ )

Table 2. ARL performance of E_ALT according to η (ARL0 = 200, n = 15 , p = 15 , λ = 0.5 , φ = 1 )

表2. 关于压缩参数 $\eta$ 的E_ALT控制图ARL表现(ARL₀ = 200, $n=15$ , $p=15$ , $\lambda =0.5$ , $\phi =1$ )

5. 结束语

本文将使用自适应Lasso阈值法估计稀疏矩阵的思想与EWMA控制图相结合，提出了一种用于监控协方差矩阵的E_ALT控制图，来监控高维过程中协方差矩阵偏移情况。通过蒙特卡洛模拟得出该控制图的ARL指标，数值模拟结果表明过程失控时，E_ALT控制图对监控过程中可能出现的不同偏移量具有一致优越性，尤其在中小偏移量下优势更为明显。同时对控制图中涉及的相关参数进行模拟分析，在通过历史数据得到相关信息时可以对参数进行调整以提高监控效率。