近期在完成信息论的作业,发现网上的资料大多是直观解释,对其中的数学原理介绍甚少,并且只介绍了向量降维,而没有介绍向量重构的问题(重构指的是:根据降维后的低维向量来恢复原始向量),因此在这里做一个总结,配合另一篇博客看效果可能会更好参考博客。
简介
PCA降维就是要把
m
m
m维空间的n个样本点
x
i
,
i
=
1
⋯
n
x_i,i=1\cdots n
xi,i=1⋯n映射成
l
l
l维的低维空间中的向量
y
i
,
i
=
1
⋯
n
y_i,i=1\cdots n
yi,i=1⋯n,其中
l
≪
m
l\ll m
l≪m。这个映射不是随意的,而是要确保利用低维的向量
y
i
y_i
yi重构出的
x
i
′
x_i^{'}
xi′与原向量
x
i
x_i
xi的误差最小化。这个过程可以描述成
y
l
×
1
=
W
l
×
m
x
m
×
1
y_{l\times1}=W_{l\times m}x_{m\times1}
yl×1=Wl×mxm×1
于是PCA的任务就是找到这样一个W矩阵。
算法数学推导
我们考虑一个二维空间中的数据压缩成一维的问题:
给定五个样本点
x
i
,
i
=
1
⋯
5
x_i,i=1\cdots5
xi,i=1⋯5,每个样本点都是一个二维的列向量,把它们拼成一个矩阵
X
=
(
x
1
,
x
2
,
x
3
,
x
4
,
x
5
)
=
[
1
1
2
4
2
1
3
3
4
4
]
X=(x_1,x_2,x_3,x_4,x_5)=\begin{bmatrix}1&1&2&4&2\\1&3&3&4&4\end{bmatrix}
X=(x1,x2,x3,x4,x5)=[1113234424]。
这个例子和文章开头提到的那篇博客中的例子是一致的,可以参考那篇博客的图片。
我们知道每个向量是通过一组基和在这组基下的坐标来描述的,例如
x
1
=
[
1
,
1
]
T
x_1=[1,1]^T
x1=[1,1]T是该向量在单位正交基
ξ
1
=
[
0
,
1
]
T
,
ξ
2
=
[
1
,
0
]
T
\xi_1=[0,1]^T,\xi_2=[1,0]^T
ξ1=[0,1]T,ξ2=[1,0]T下的坐标表示,坐标实际上就是向量向各个基上的投影值。同样的,如果我们要将一个二维向量压缩成一维向量,那么只需要找到一条直线,直线的单位方向向量
w
w
w作为基,然后用向量向这条直线的投影值就可以描述压缩后的一维向量。
一条直线可以用它经过的点
μ
\mu
μ和单位方向向量
w
w
w来描述,即
x
=
μ
+
α
w
x=\mu+\alpha w
x=μ+αw,(这里我们用的
μ
\mu
μ是上述五个样本点的均值
[
2
,
3
]
T
[2,3]^T
[2,3]T)这里的
α
\alpha
α就可以理解为坐标,
w
w
w是基。那么我们要寻找的二维向量
x
i
x_i
xi经过压缩后的一维向量其实就是
α
i
\alpha_i
αi,现在需要确定
w
w
w,这样才能求出二维向量向直线的投影值。正如前一节所述,这个
w
w
w不是任意的,而是应该确保重构后的误差最小化,它可以描述成如下的优化问题:
min
w
f
(
w
)
=
1
2
∑
i
=
1
n
∥
(
μ
+
α
i
w
)
−
x
i
∥
2
2
\min _{w} f(w)=\frac{1}{2} \sum_{i=1}^n\left\|\left(\mu+\alpha_i w\right)-x_i\right\|_2^2
wminf(w)=21i=1∑n∥(μ+αiw)−xi∥22
先求出
α
i
\alpha_i
αi的取值,
∂
f
∂
α
i
=
w
T
(
μ
+
α
i
w
−
x
i
)
=
0
\frac{\partial f}{\partial \alpha_i}=w^T(\mu+\alpha_iw-x_i)=0
∂αi∂f=wT(μ+αiw−xi)=0
由于
w
w
w是单位向量,即
w
T
w
=
1
w^Tw=1
wTw=1,由上式可得:
α
i
=
w
T
(
x
i
−
μ
)
\alpha_i=w^T(x_i-\mu)
αi=wT(xi−μ)
这个公式就给出了
α
i
\alpha_i
αi的求解方法。下面继续推导确定
w
w
w的过程,定义散布矩阵如下:
S
=
∑
i
=
1
n
(
x
i
−
μ
)
(
x
i
−
μ
)
T
S=\sum_{i=1}^n(x_i-\mu)(x_i-\mu)^T
S=i=1∑n(xi−μ)(xi−μ)T
对于上面的例子,我们把
X
X
X中的每个样本
x
i
x_i
xi都减去均值
μ
\mu
μ得到一个新的矩阵记为
X
~
=
[
−
1
−
1
0
2
0
−
2
0
0
1
1
]
\tilde{X}=\begin{bmatrix}-1&-1&0&2&0\\-2&0&0&1&1\end{bmatrix}
X~=[−1−2−10002101],那么上面的散步矩阵其实可以简单地记为
S
=
X
~
X
~
T
S=\tilde{X}\tilde{X}^T
S=X~X~T
说明它是一个对称矩阵。
将上面求出的
α
i
\alpha_i
αi的表达式代入到
f
(
w
)
f(w)
f(w)中,得到:
f
(
w
)
=
1
2
∑
i
=
1
n
∥
α
i
w
−
(
x
i
−
μ
)
∥
2
2
=
1
2
(
∑
i
=
1
n
α
i
2
∥
w
∥
2
2
−
2
∑
i
=
1
n
α
i
w
T
(
x
i
−
μ
)
+
∑
i
=
1
n
∥
x
i
−
μ
∥
2
2
)
=
−
1
2
∑
i
=
1
n
α
i
2
+
1
2
∑
i
=
1
n
∥
x
i
−
μ
∥
2
2
=
−
1
2
∑
i
=
1
n
[
w
T
(
x
i
−
μ
)
]
2
+
1
2
∑
i
=
1
n
∥
x
i
−
μ
∥
2
2
=
−
1
2
∑
i
=
1
n
w
T
(
x
i
−
μ
)
(
x
i
−
μ
)
T
w
+
1
2
∑
i
=
1
n
∥
x
i
−
μ
∥
2
2
=
−
1
2
w
T
(
∑
i
=
1
n
(
x
i
−
μ
)
(
x
i
−
μ
)
T
)
w
+
1
2
∑
i
=
1
n
∥
x
i
−
μ
∥
2
2
=
−
1
2
w
T
S
w
+
1
2
∑
i
=
1
n
∥
x
i
−
μ
∥
2
2
\begin{aligned} f(\boldsymbol{w}) & =\frac{1}{2} \sum_{i=1}^n\left\|\alpha_i \boldsymbol{w}-\left(\boldsymbol{x}_i-\boldsymbol{\mu}\right)\right\|_2^2 \\ & =\frac{1}{2}\left(\sum_{i=1}^n \alpha_i^2\|\boldsymbol{w}\|_2^2-2 \sum_{i=1}^n \alpha_i \boldsymbol{w}^{\mathrm{T}}\left(\boldsymbol{x}_i-\boldsymbol{\mu}\right)+\sum_{i=1}^n\left\|\boldsymbol{x}_i-\boldsymbol{\mu}\right\|_2^2\right) \\ & =-\frac{1}{2} \sum_{i=1}^n \alpha_i^2+\frac{1}{2} \sum_{i=1}^n\left\|\boldsymbol{x}_i-\boldsymbol{\mu}\right\|_2^2 \\ & =-\frac{1}{2} \sum_{i=1}^n\left[\boldsymbol{w}^{\mathrm{T}}\left(\boldsymbol{x}_i-\boldsymbol{\mu}\right)\right]^2+\frac{1}{2} \sum_{i=1}^n\left\|\boldsymbol{x}_i-\boldsymbol{\mu}\right\|_2^2 \\ & =-\frac{1}{2} \sum_{i=1}^n \boldsymbol{w}^{\mathrm{T}}\left(\boldsymbol{x}_i-\boldsymbol{\mu}\right)\left(\boldsymbol{x}_i-\boldsymbol{\mu}\right)^{\mathrm{T}} \boldsymbol{w}+\frac{1}{2} \sum_{i=1}^n\left\|\boldsymbol{x}_i-\boldsymbol{\mu}\right\|_2^2 \\ & =-\frac{1}{2} \boldsymbol{w}^{\mathrm{T}}\left(\sum_{i=1}^n\left(\boldsymbol{x}_i-\boldsymbol{\mu}\right)\left(\boldsymbol{x}_i-\boldsymbol{\mu}\right)^{\mathrm{T}}\right) \boldsymbol{w}+\frac{1}{2} \sum_{i=1}^n\left\|\boldsymbol{x}_i-\boldsymbol{\mu}\right\|_2^2 \\ & =-\frac{1}{2} \boldsymbol{w}^{\mathrm{T}} S \boldsymbol{w}+\frac{1}{2} \sum_{i=1}^n\left\|\boldsymbol{x}_i-\boldsymbol{\mu}\right\|_2^2 \end{aligned}
f(w)=21i=1∑n∥αiw−(xi−μ)∥22=21(i=1∑nαi2∥w∥22−2i=1∑nαiwT(xi−μ)+i=1∑n∥xi−μ∥22)=−21i=1∑nαi2+21i=1∑n∥xi−μ∥22=−21i=1∑n[wT(xi−μ)]2+21i=1∑n∥xi−μ∥22=−21i=1∑nwT(xi−μ)(xi−μ)Tw+21i=1∑n∥xi−μ∥22=−21wT(i=1∑n(xi−μ)(xi−μ)T)w+21i=1∑n∥xi−μ∥22=−21wTSw+21i=1∑n∥xi−μ∥22
上式的第二项与
w
w
w无关,因此要极小化
f
(
w
)
f(w)
f(w),只要使第一项极小化,于是优化问题转化为
min
−
1
2
w
T
S
w
s
.
t
.
w
T
w
=
1
\begin{aligned}\min &-\frac{1}{2}w^TSw\\ {\rm s.t.\ } &w^Tw=1\end{aligned}
mins.t. −21wTSwwTw=1
这个优化问题可以用拉格朗日乘子法求解,令拉格朗日函数为
L
(
w
,
λ
)
=
−
1
2
w
T
S
w
+
λ
2
(
w
T
w
−
1
)
L(w,\lambda)=-\frac{1}{2}w^TSw+\frac{\lambda}{2}(w^Tw-1)
L(w,λ)=−21wTSw+2λ(wTw−1)
令
∂
L
∂
w
=
−
S
w
+
λ
w
=
0
\frac{\partial L}{\partial w}=-Sw+\lambda w=0
∂w∂L=−Sw+λw=0
从而
S
w
=
λ
w
Sw=\lambda w
Sw=λw
到这里,结果已经逐渐清晰了,我们要求的
w
w
w正是矩阵
S
S
S的特征向量。稍作变形:
w
T
S
w
=
λ
w
T
w
=
λ
w^TSw=\lambda w^Tw=\lambda
wTSw=λwTw=λ
我们要最小化
−
1
2
w
T
S
w
-\frac{1}{2}w^TSw
−21wTSw,就是要最大化
w
T
S
w
w^TSw
wTSw,则
w
w
w应该是
S
S
S的最大特征值
λ
max
\lambda_{\max}
λmax对应的特征向量。
算法总结
至此,我们可以总结一下二维向量压缩成一维的PCA的方法:
(1)求矩阵
S
=
∑
i
=
1
n
(
x
i
−
μ
)
(
x
i
−
μ
)
T
=
X
~
X
~
T
S=\sum_{i=1}^n(x_i-\mu)(x_i-\mu)^T=\tilde{X}\tilde{X}^T
S=i=1∑n(xi−μ)(xi−μ)T=X~X~T
(2)求
S
S
S最大的特征值对应的特征向量,即
w
w
w
(3)求
α
i
=
w
T
(
x
i
−
μ
)
\alpha_i=w^T(x_i-\mu)
αi=wT(xi−μ)
于是 X = ( x 1 , x 2 , x 3 , x 4 , x 5 ) X=(x_1,x_2,x_3,x_4,x_5) X=(x1,x2,x3,x4,x5)经过压缩之后得到的结果就是 Y = ( α 1 , α 2 , α 3 , α 4 , α 5 ) Y=(\alpha_1,\alpha_2,\alpha_3,\alpha_4,\alpha_5) Y=(α1,α2,α3,α4,α5)。
投影到方向向量
w
w
w所对应的直线之后,
w
w
w成了唯一的一个基,于是一维空间中的样本
x
i
′
x_i^{'}
xi′可以由基向量
w
w
w表示:
x
i
′
=
μ
+
α
i
w
x_i^{'}=\mu+\alpha_iw
xi′=μ+αiw
在原来的2维空间中,我们用基的系数来表示样本
x
i
x_i
xi,而在1维空间中,同样以基
w
w
w的系数
α
i
\alpha_i
αi来表示一维向量,它被称为主成分。
将二维向量压缩成一维向量 α i \alpha_i αi有时候只是为了减少传输时的数据量,一维向量是无法直接使用的,需要根据一维向量重构出原来的二维向量。如何重构呢?其实上面关于的 x i ′ x_i^{'} xi′的公式已经给出了答案。
二维样本PCA降维的例子
还是上面的例子,我们按照这个流程走一遍:
μ
=
[
2
,
3
]
T
\mu=[2,3]^T
μ=[2,3]T
X
~
=
X
−
μ
=
[
−
1
−
1
0
2
0
−
2
0
0
1
1
]
\tilde{X}=X-\mu=\begin{bmatrix}-1&-1&0&2&0\\-2&0&0&1&1\end{bmatrix}
X~=X−μ=[−1−2−10002101]
S
=
X
~
X
~
T
=
[
−
1
−
1
0
2
0
−
2
0
0
1
1
]
[
−
1
−
2
−
1
0
0
0
2
1
0
1
]
=
[
6
4
4
6
]
S=\tilde{X}\tilde{X}^T=\begin{bmatrix}-1&-1&0&2&0\\-2&0&0&1&1\end{bmatrix}\begin{bmatrix}-1&-2\\-1&0\\0&0\\2&1\\0&1\end{bmatrix}=\begin{bmatrix}6&4\\4&6\end{bmatrix}
S=X~X~T=[−1−2−10002101]
−1−1020−20011
=[6446]
S的最大特征值为10,对应特征向量
w
=
[
1
2
,
1
2
]
T
w=[\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}}]^T
w=[21,21]T
于是降维后的表示为:
Y
=
(
α
1
,
α
2
,
α
3
,
α
4
,
α
5
)
=
w
T
X
~
=
[
1
2
,
1
2
]
[
−
1
−
1
0
2
0
−
2
0
0
1
1
]
=
[
−
3
/
2
,
−
1
/
2
,
0
,
3
/
2
,
−
1
/
2
]
Y=(\alpha_1,\alpha_2,\alpha_3,\alpha_4,\alpha_5)=w^T\tilde{X}=[\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}}]\begin{bmatrix}-1&-1&0&2&0\\-2&0&0&1&1\end{bmatrix}=[-3/\sqrt{2},-1/\sqrt{2},0,3/\sqrt{2},-1/\sqrt{2}]
Y=(α1,α2,α3,α4,α5)=wTX~=[21,21][−1−2−10002101]=[−3/2,−1/2,0,3/2,−1/2]
要重构第一个样本
x
1
′
x_1{'}
x1′,方法是:
x
1
′
=
μ
+
α
1
w
=
[
2
,
3
]
T
+
−
3
2
[
1
2
,
1
2
]
T
=
[
1
2
,
3
2
]
T
x_1^{'}=\mu+\alpha_1w=[2,3]^T+\frac{-3}{\sqrt{2}}[\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}}]^T=[\frac{1}{2},\frac{3}{2}]^T
x1′=μ+α1w=[2,3]T+2−3[21,21]T=[21,23]T
当然,与原本的
x
1
=
[
1
,
1
]
T
x_1=[1,1]^T
x1=[1,1]T还是有一些误差的。
更高维的情况
前面介绍的是二维降为一维的情况,更一般地,对于
m
m
m维向量
x
i
x_i
xi如果要降维为
l
l
l维的
y
i
y_i
yi,算法也是类似的,不加证明地给出以下步骤:
(1)求矩阵
S
=
∑
i
=
1
n
(
x
i
−
μ
)
(
x
i
−
μ
)
T
=
X
~
X
~
T
S=\sum_{i=1}^n(x_i-\mu)(x_i-\mu)^T=\tilde{X}\tilde{X}^T
S=i=1∑n(xi−μ)(xi−μ)T=X~X~T
(2)求
S
S
S的所有特征值,从大到小排列,选取前
l
l
l个特征值所对应的特征向量,即
W
=
(
w
1
,
w
2
,
⋯
,
w
l
)
W=(w_1,w_2,\cdots,w_l)
W=(w1,w2,⋯,wl)。
(3)求各个样本点
x
i
x_i
xi对应于基
w
1
,
w
2
,
⋯
,
w
l
w_1,w_2,\cdots,w_l
w1,w2,⋯,wl的系数,即主成分,
α
i
,
k
=
w
k
T
(
x
i
−
μ
)
,
k
=
1
⋯
l
\alpha_{i,k}=w_k^T(x_i-\mu),k=1\cdots l
αi,k=wkT(xi−μ),k=1⋯l,得到低维的表示
y
i
=
(
α
i
,
1
,
α
i
,
2
,
⋯
α
i
,
l
)
T
y_i=(\alpha_{i,1},\alpha_{i,2},\cdots \alpha_{i,l})^T
yi=(αi,1,αi,2,⋯αi,l)T
这个过程可以写成矩阵的形式:
Y
=
W
T
X
~
Y=W^T\tilde{X}
Y=WTX~
Y
=
(
y
1
,
y
2
,
⋯
,
y
n
)
Y=(y_1,y_2,\cdots,y_n)
Y=(y1,y2,⋯,yn)是压缩后的样本点组成的矩阵。
(4) 原向量的重构:
x
i
′
=
μ
+
∑
k
=
1
L
α
i
,
k
w
k
x_i^{'}=\mu+\sum_{k=1}^L\alpha_{i,k}w_k
xi′=μ+k=1∑Lαi,kwk
写成矩阵的形式为:
X
′
=
μ
+
W
Y
X^{'}=\mu+WY
X′=μ+WY
PCA用于图像处理
原图:
PCA处理后:(选40个特征向量)