7.4.5 鲁棒主成分分析 PCA
根据每个样本点数据
a
i
\mathbf{a}_{i}
ai 在第一主方向
u
1
\mathbf{u}_1
u1 上的投影的方差最大,知样本点在此方向最分散,为第一主方向。我们可以反过来思考,每个样本点数据
a
i
\mathbf{a}_{i}
ai 在某个方向
α
\mathbf{\alpha}
α 上投影,计算其方差,方向不同时,方差是不同的,当方差最大时,此时方向就是第一主方向,所以第一主方向可如下定义:
u
1
=
a
r
g
m
a
x
α
∑
a
i
∈
A
(
a
i
T
α
)
2
=
a
r
g
m
a
x
α
(
A
T
α
)
T
A
T
α
=
=
a
r
g
m
a
x
α
α
T
(
A
A
T
)
α
\mathbf{u}_1= argmax_{\mathbf{\alpha}} \sum_{\mathbf{a}_{i} \in A} (\mathbf{a}^T_{i}\mathbf{\alpha})^2\\ = argmax_{\mathbf{\alpha}} (A^T\mathbf{\alpha})^TA^T\mathbf{\alpha}== argmax_{\mathbf{\alpha}} \mathbf{\alpha}^T(AA^T)\mathbf{\alpha}
u1=argmaxαai∈A∑(aiTα)2=argmaxα(ATα)TATα==argmaxααT(AAT)α
该定义是计算方差,需进行平方运算,在最小二乘法章节指出,平方运算很容易受到异常点的影响,比如数据矩阵 A A A 中存在离群的样本点,即使这样的离群点数目很少,也会导致第一主方向极大地偏离理想方向。如果点云不存在离群点,即样本点都是采样于高斯分布,则上述计算得到的主方向就是最优主方向。
为了提高鲁棒性即抗离群点能力,需要减弱离群点对方差的贡献,采用一次函数就可达到目的,即使用绝对值,而不是平方,则第一主方向可如下定义:
u
1
=
a
r
g
m
a
x
α
∑
a
i
∈
A
∣
a
i
T
α
∣
\mathbf{u}_1 = argmax_{\mathbf{\alpha}} \sum_{\mathbf{a}_{i} \in A} |\mathbf{a}^T_{i}\mathbf{\alpha}|
u1=argmaxαai∈A∑∣aiTα∣
这是一种常用的鲁棒PCA,原理简单易懂,效果也较好,计算也方便。下面来求解 u 1 \mathbf{u}_1 u1 。
J 1 = ∑ a i ∈ A ∣ a i T α ∣ = ∑ a i ∈ A ∣ α T a i ∣ = ∑ a i ∈ A + α T a i − ∑ a i ∈ A − α T a i = α T ( ∑ a i ∈ A + a i − ∑ a i ∈ A − a i ) = 2 α T ∑ a i ∈ A + a i J_1 = \sum_{\mathbf{a}_{i} \in A} |\mathbf{a}^T_{i}\mathbf{\alpha}|= \sum_{\mathbf{a}_{i} \in A} |\mathbf{\alpha}^T\mathbf{a}_{i}|\\ = \sum_{\mathbf{a}_{i} \in A_+} \mathbf{\alpha}^T\mathbf{a}_{i} - \sum_{\mathbf{a}_{i} \in A_-} \mathbf{\alpha}^T\mathbf{a}_{i} \\ = \mathbf{\alpha}^T (\sum_{\mathbf{a}_{i} \in A_+} \mathbf{a}_{i} - \sum_{\mathbf{a}_{i} \in A_-} \mathbf{a}_{i}) \\ = 2\mathbf{\alpha}^T \sum_{\mathbf{a}_{i} \in A_+} \mathbf{a}_{i} J1=ai∈A∑∣aiTα∣=ai∈A∑∣αTai∣=ai∈A+∑αTai−ai∈A−∑αTai=αT(ai∈A+∑ai−ai∈A−∑ai)=2αTai∈A+∑ai
其中集合 A + = { a i ∈ A : α T a i ≥ 0 } A_+ = \{\mathbf{a}_{i} \in A: \mathbf{\alpha}^T\mathbf{a}_{i} \ge 0\} A+={ai∈A:αTai≥0} 即与向量 α \mathbf{\alpha} α 夹角小于等于 90 90 90 度的样本集 和 A − = { a i ∈ A : α T a i < 0 } A_- = \{\mathbf{a}_{i} \in A: \mathbf{\alpha}^T\mathbf{a}_{i} < 0\} A−={ai∈A:αTai<0} 。
因为样本集进行了中心化,所以 ∑ a i ∈ A a i = 0 \sum_{\mathbf{a}_{i} \in A} \mathbf{a}_{i} = \mathbf{0} ∑ai∈Aai=0 ,因此 ∑ a i ∈ A + a i = − ∑ a i ∈ A − a i \sum_{\mathbf{a}_{i} \in A_+} \mathbf{a}_{i} = - \sum_{\mathbf{a}_{i} \in A_-} \mathbf{a}_{i} ∑ai∈A+ai=−∑ai∈A−ai 。
要最大化
J
1
J_1
J1 ,只需向量
α
\mathbf{\alpha}
α 与向量
∑
a
i
∈
A
+
a
i
\sum_{\mathbf{a}_{i} \in A_+} \mathbf{a}_{i}
∑ai∈A+ai 平行即可,所以
α
=
u
n
i
t
(
∑
a
i
∈
A
+
a
i
)
\mathbf{\alpha} = unit (\sum_{\mathbf{a}_{i} \in A_+} \mathbf{a}_{i})
α=unit(ai∈A+∑ai)
即单位化向量
∑
a
i
∈
A
+
a
i
\sum_{\mathbf{a}_{i} \in A_+} \mathbf{a}_{i}
∑ai∈A+ai 。
所以可以采用迭代法求解最优 u 1 \mathbf{u}_1 u1 。首先随机生成单位向量 α \mathbf{\alpha} α 。
1、根据内积 α T a i \mathbf{\alpha}^T\mathbf{a}_{i} αTai 符号,划分集合 A + A_+ A+ ;
2、令 α = u n i t ( ∑ a i ∈ A + a i ) \mathbf{\alpha} = unit (\sum_{\mathbf{a}_{i} \in A_+} \mathbf{a}_{i}) α=unit(∑ai∈A+ai) 最大化 J 1 J_1 J1 ,得到新的单位向量 α \mathbf{\alpha} α。
重复上面两步,直到收敛即单位向量 α \mathbf{\alpha} α 改变很小,此时 u 1 = α \mathbf{u}_1=\mathbf{\alpha} u1=α 。
可计算第一主方向上方差为 ∥ A T u 1 ∥ 2 \|A^T\mathbf{u}_1\|^2 ∥ATu1∥2 ,注意其不等于矩阵 A A A 奇异值 σ 1 \sigma_1 σ1 。
那如何计算第二主方向呢?
第二主方向是垂直于第一主方向,最大方差方向即
u
2
=
a
r
g
m
a
x
α
∈
⊥
u
1
∑
a
i
∈
A
∣
α
T
a
i
∣
\mathbf{u}_2 = argmax_{\mathbf{\alpha} \in \perp \mathbf{u}_1} \sum_{\mathbf{a}_{i} \in A} |\mathbf{\alpha}^T\mathbf{a}_{i}|
u2=argmaxα∈⊥u1ai∈A∑∣αTai∣
把每个样本点沿第一主方向进行正交分解,即分解为第一主方向分量
a
i
∥
\mathbf{a}^\parallel_{i}
ai∥ 和垂直于第一主方向分量
a
i
⊥
\mathbf{a}^\perp_{i}
ai⊥ ,则
a
i
=
a
i
∥
+
a
i
⊥
a
i
⊥
=
(
E
−
u
1
u
1
T
)
a
i
a
i
∥
=
u
1
u
1
T
a
i
\mathbf{a}_{i} = \mathbf{a}^\parallel_{i} + \mathbf{a}^\perp_{i} \\ \mathbf{a}^\perp_{i} = (E-\mathbf{u}_1\mathbf{u}^T_1)\mathbf{a}_{i} \\ \mathbf{a}^\parallel_{i} = \mathbf{u}_1\mathbf{u}^T_1\mathbf{a}_{i} \\
ai=ai∥+ai⊥ai⊥=(E−u1u1T)aiai∥=u1u1Tai
令
J
=
∑
a
i
∈
A
∣
α
T
a
i
∣
=
∑
a
i
∈
A
∣
(
a
i
∥
+
a
i
⊥
)
(
α
T
∥
+
α
T
⊥
)
∣
=
∑
a
i
∈
A
∣
a
i
∥
α
T
∥
+
a
i
⊥
α
T
⊥
∣
J = \sum_{\mathbf{a}_{i} \in A} |\mathbf{\alpha}^T\mathbf{a}_{i}| \\ = \sum_{\mathbf{a}_{i} \in A} |(\mathbf{a}^{\parallel}_{i} + \mathbf{a}^{\perp}_{i}) (\mathbf{\alpha}^{T\parallel} + \mathbf{\alpha}^{T\perp})|\\ = \sum_{\mathbf{a}_{i} \in A} |\mathbf{a}^{\parallel}_{i}\mathbf{\alpha}^{T\parallel} + \mathbf{a}^{\perp}_{i}\mathbf{\alpha}^{T\perp}| \\
J=ai∈A∑∣αTai∣=ai∈A∑∣(ai∥+ai⊥)(αT∥+αT⊥)∣=ai∈A∑∣ai∥αT∥+ai⊥αT⊥∣
当
α
∈
⊥
u
1
\mathbf{\alpha} \in \perp \mathbf{u}_1
α∈⊥u1 时,即
α
∥
=
0
\mathbf{\alpha}^{\parallel} = \mathbf{0}
α∥=0 ,
J
2
=
J
=
∑
a
i
∈
A
∣
a
i
⊥
α
T
⊥
∣
J_2 = J = \sum_{\mathbf{a}_{i} \in A} |\mathbf{a}^{\perp}_{i}\mathbf{\alpha}^{T\perp}|
J2=J=ai∈A∑∣ai⊥αT⊥∣
要使 J 2 J_2 J2 最大化,同样可以采用使 J 1 J_1 J1 最大化的迭代法,只需保证 α \mathbf{\alpha} α 始终垂直于 u 1 \mathbf{u}_1 u1 ,由于 a i ⊥ \mathbf{a}^{\perp}_{i} ai⊥ 均垂直于 u 1 \mathbf{u}_1 u1 ,故 ∑ a i ∈ A + a i ⊥ \sum_{\mathbf{a}_{i} \in A_+} \mathbf{a}^{\perp}_{i} ∑ai∈A+ai⊥ 垂直于 u 1 \mathbf{u}_1 u1 ,其中集合 A + = { a i ∈ A : a i ⊥ α T ≥ 0 } A_+ = \{\mathbf{a}_{i} \in A: \mathbf{a}^{\perp}_{i}\mathbf{\alpha}^T \ge 0\} A+={ai∈A:ai⊥αT≥0} 。所以只要初始化时 α \mathbf{\alpha} α 垂直于 u 1 \mathbf{u}_1 u1 即可,则令 α = u i n t ( ( E − u 1 u 1 T ) v ) \mathbf{\alpha} = uint((E-\mathbf{u}_1\mathbf{u}^T_1)\mathbf{v}) α=uint((E−u1u1T)v) 其中 v \mathbf{v} v 为单位随机向量。
整理上面结论得到第二主方向的计算过程:
首先样本点投影到垂直于第一主方向的子空间,即 ( E − u 1 u 1 T ) A (E-\mathbf{u}_1\mathbf{u}^T_1)A (E−u1u1T)A ,为了简化记号,投影后的数据矩阵还是记为 A A A ,其次随机生成单位向量 α = u i n t ( ( E − u 1 u 1 T ) v ) \mathbf{\alpha} = uint((E-\mathbf{u}_1\mathbf{u}^T_1)\mathbf{v}) α=uint((E−u1u1T)v) 其中 v \mathbf{v} v 为单位随机向量,最后最大化 J 2 = 2 α T ∑ a i ∈ A + a i J_2 = 2\mathbf{\alpha}^T \sum_{\mathbf{a}_{i} \in A_+} \mathbf{a}_{i} J2=2αT∑ai∈A+ai ,此时 u 2 = α \mathbf{u}_2 = \mathbf{\alpha} u2=α 。
可计算第二主方向上方差为 ∥ A T u 2 ∥ 2 \|A^T\mathbf{u}_2\|^2 ∥ATu2∥2 ,注意其不等于矩阵 A A A 奇异值 σ 2 \sigma_2 σ2 。
第 k + 1 k+1 k+1 个主方向的计算过程如下:
首先样本点投影到垂直于前 k k k 个主方向的子空间,即 A k + 1 = ( E − u k u k T ) A k A_{k+1} = (E-\mathbf{u}_k\mathbf{u}^T_k)A_k Ak+1=(E−ukukT)Ak ,其中 A 0 = A A_0 = A A0=A 为原数据矩阵, u 0 = 0 \mathbf{u}_0 = \mathbf{0} u0=0 为零向量。为了简化记号,投影后的数据矩阵还是记为 A = A k + 1 A=A_{k+1} A=Ak+1 ,其次随机生成单位向量 α = u i n t ( ( E − ( u 1 u 1 T + ⋯ + u k u k T ) ) v ) \mathbf{\alpha} = uint((E-(\mathbf{u}_1\mathbf{u}^T_1+\cdots+\mathbf{u}_k\mathbf{u}^T_k))\mathbf{v}) α=uint((E−(u1u1T+⋯+ukukT))v) 其中 v \mathbf{v} v 为单位随机向量,最后最大化 J k + 1 = 2 α T ∑ a i ∈ A + a i J_{k+1} = 2\mathbf{\alpha}^T \sum_{\mathbf{a}_{i} \in A_+} \mathbf{a}_{i} Jk+1=2αT∑ai∈A+ai ,此时 u k + 1 = α \mathbf{u}_{k+1} = \mathbf{\alpha} uk+1=α 。
据此可以获得全部主方向 u 1 , ⋯ , u r \mathbf{u}_1,\cdots,\mathbf{u}_r u1,⋯,ur ,取前 k k k 个作为最后的主方向。
数据矩阵重构如果采用经典PCA重构方法,由于主方向计算方法不同,即使采用全部主方向,也不能使重构残差为零。
异常点检测
点云位于
m
m
m 维空间的
k
k
k 维子空间,在该子空间中可被超椭球紧致包围,位于该椭球内的点认为是正常点,椭球外的点是异常点,越是远离椭球,异常度越大。据此提出两个定量指标衡量异常程度,一个是显而易见的,异常点位于
k
k
k 维子空间外,其与
k
k
k 维子空间的垂直距离为 OD,距离越大,则越异常,样本点
a
i
\mathbf{a}_{i}
ai 的垂直距离
O
D
i
OD_i
ODi 为
O
D
i
=
∥
(
E
−
∑
j
=
1
k
u
j
u
j
T
)
a
i
∥
OD_i = \|(E-\sum^k_{j=1} \mathbf{u}_j\mathbf{u}^T_j)\mathbf{a}_{i}\|
ODi=∥(E−j=1∑kujujT)ai∥
异常点虽然位于 k k k 维子空间,但远离椭球中心,即异常点投影值 u j T a i \mathbf{u}^T_j\mathbf{a}_{i} ujTai 很大,对每个投影值用该主方向上方差的平方根 σ j \sigma_j σj 归一化,得到投影距离 PD,距离越大,则越异常,样本点 a i \mathbf{a}_{i} ai 的投影距离 P D i PD_i PDi 为
P D i = ∑ j = 1 k ( u j T a i / σ j ) 2 PD_i = \sqrt{\sum^k_{j=1}(\mathbf{u}^T_j\mathbf{a}_{i}/\sigma_j)^2} PDi=j=1∑k(ujTai/σj)2
以三维作为例子直观说明异常点,假设正常点云位于平面上椭圆内,则该平面内远离椭圆的点投影距离 PD 大;该平面外的点距离该平面的距离为垂直距离 OD,距离越大,则越异常。
根据垂直距离和投影距离,可以把样本点分为四类:
1、正常样本点:两个距离都小,在正常范围内;
2、好杠杆点:投影距离大,垂直距离正常;
3、坏杠杆点:投影距离大,垂直距离大;
4、垂直外点:投影距离正常,垂直距离大。
后面三种都是异常点。
这样可以获得诊断图,横坐标为投影距离,纵坐标为垂直距离,绘制每个样本点。聚集在原点附近的点为正常点,远离原点的点为异常点。
当存在特别离群的点或离群点比例较高时,此时可以先去除这些离群点,然后利用剩下的正常点进行常规PCA,可以获得更好的结果。利用此时的主方向和奇异值,重新绘制诊断图。
上面介绍的仅是众多鲁棒主成分方法的一种,鲁棒主成分方法的核心是获得方差 u 1 = a r g m a x α ∑ a i ∈ A ( α T a i ) 2 \mathbf{u}_1 = argmax_{\mathbf{\alpha}} \sum_{\mathbf{a}_{i} \in A} (\mathbf{\alpha}^T\mathbf{a}_{i})^2 u1=argmaxα∑ai∈A(αTai)2 的鲁棒估计,采用绝对值只是一种比较直观的方法。
采用中值估计方差是一个更鲁棒的方法,但寻求最优解很困难。令投影为 z i = α T a i z_i = \mathbf{\alpha}^T\mathbf{a}_{i} zi=αTai ,则基于中值的方差估计公式为
S = 1.483 m e d ( ∣ z i − z ∗ ∣ ) S = 1.483 med(|z_i - z^*|) S=1.483med(∣zi−z∗∣)
其中 z ∗ z^* z∗ 是投影 z i z_i zi 的中值,系数 1.483 是为了使方差估计值与高斯分布一致, S S S 为方差的平方根,称为标准差。具体如何获得最优解,由于很复杂,不打算展开。