记录个人对PCA算法的理解
文章结构
- 我们为什么要使用PCA
- PCA算法理解与数学推导
- 个人PCA算法的简单实现
- 总结
我们为什么要使用PCA
在机器学习领域,一直都有一个特别让人头痛的问题,那就是"维度灾难"(curse of dimensionality)。当一个数据集的特征过于多的时候,那么这个数据集的维度就会特别的大。过大的维度会给我们的模型带来非常大的压力,整个模型的训练需要极大运算资源(computationally expensive)。所以在实际的生产中,我们总是在思考如何将高维数据使用低维数据进行表示。从而在机器学习的领域中,就出现了甚多的降维技术(dimensionality reduction)。PCA(Principal components analysis)就是其中的一种。PCA的主要作用,简单的理解就是提取原数据集中的特征,力争用最少的维度表达最多的信息。 如何理解这句话呢?拿动物识别来做个’栗子’,人可以从高清的照片和相对模糊的照片中识别出不同的动物,虽然我们无法从模糊的照片中看到动物的细节,但是模糊照片中的颜色和轮廓就已经给我们提供了足够多的信息,而高清图片对于识别出动物没有实质性的帮助,它只是提供了更多的细节。同样的对于机器学习模型也是如此。除了减少模型的负担之外,PCA还有两个重要动机:数据集中那些只包含很少信息的成分,往往可以视为噪声,所以丢弃这些成分,可以起到去噪的作用。其二,部分信息的舍弃可使得整个采样密度上升。
PCA算法的理解和数学推导
注:PCA的定义有两种,分别是maximum variance formulation 和 minimum-error formulation,接下来的推导将围绕着前者展开。
那么PCA到底是如何来完成所谓降维的呢?如何来保留最多的信息呢?
首先我想指出,为什么降维之后信息会丢失??是什么原因导致的信息丢失?
在这里举一个将二维数据投射至一维的栗子,如下图所示
从上图这个简单的例子中我们可以观察到在笛卡尔坐标系中的两个红色点和两个绿色点投射至一维之后分别只剩下一个了。不难看出投射导致的信息丢失主要是因为数据点之间的重叠和覆盖。所以为了避免这样的覆盖发生我们需要做的就是在投射的时候尽可能将点投射的分散一点,数据点不会相互的重叠和覆盖或者说大部分的数据点没有相互覆盖,这样可以保证原始数据的信息的最大存留。那么在数学上表示离散数据离散程度的概念是方差,它阐述了每一个数据点离中心的距离。换而言之,新维度上的数据点的方差应该被最大化。
在开始寻找主成分之前(principle components)之前,我们需要对整个数据做中心化(均值化),通俗的理解就是将每一个数据点的中心都变成坐标原点。做过中心化之后,数据集中每一个维度的均值应该0,表示成:
1 m ∑ m X m = 0 \frac{1} {m}\sum_{m}X_{m}=0 m1∑mXm=0
接下来根据数据创建一个协方差矩阵,它可以反映出每一个维度的线性相关性。
协方差矩阵可以表示成
C
O
V
=
1
n
∑
(
(
X
i
j
−
E
(
x
j
)
(
X
i
k
−
E
(
X
k
)
)
COV=\frac{1}{n}\sum((X_{ij}-E(x_{j})(X_{ik}-E(X_{k}))
COV=n1∑((Xij−E(xj)(Xik−E(Xk))
X i j X_{ij} Xij和 X i k X_{ik} Xik表示在两个维度下的单个样本,E表示的是当前维度的期望(平均值)。在中心化之后每一个维度的均值为0, 所以上述公式可以被简化成以下式子,用S来表示:
S = 1 n ∑ X i j X i k S=\frac{1}{n}\sum X_{ij}X_{ik} S=n1∑XijXik
接下来我们需要投射样本到一个低维空间上,我们需要找到一个投射的方向,我们将这个方向表示成
U
T
U^{T}
UT,所以投射之后的样本可以记为:
x
p
r
o
j
e
c
t
e
d
=
U
j
T
x
i
x_{projected}=U^{T}_{j}x_{i}
xprojected=UjTxi
这里的 x i x_{i} xi是原始样本,上述式子表达的是在低维的坐标系下,第 j 维的坐标上的数据。因为投射前每一个维度的均值是0,投射之后也应该是0,所以投射之后的每一个维度的均值表示成:
1 n ∑ i = 1 n U T x p r o j e c t e d ( i ) = 0 \frac{1}{n}\sum_{i=1}^{n}U^{T}x_{projected}^{(i)}=0 n1∑i=1nUTxprojected(i)=0
那么同理,投射之后的数据方差也可以用方差公式求得:
v
a
r
i
a
n
c
e
p
r
o
j
e
c
t
e
d
=
1
n
∑
i
=
1
n
(
U
T
x
p
r
o
j
e
c
t
e
d
(
i
)
−
E
(
x
i
)
)
2
variance_{projected} = \frac{1}{n}\sum_{i=1}^{n}(U^{T}x_{projected}^{(i)}-E(x_{i}))^{2}
varianceprojected=n1∑i=1n(UTxprojected(i)−E(xi))2
因为投射之后的数据之间的均值任是0,即E(x)是0,所以
v a r i a n c e p r o j e c t e d = 1 n ∑ i = 1 n ( U T x p r o j e c t e d ( i ) ) 2 variance_{projected} = \frac{1}{n}\sum_{i=1}^{n} (U^{T}x_{projected}^{(i)})^{2} varianceprojected=n1∑i=1n(UTxprojected(i))2
尝试将上述式子展开得到:
v
a
r
i
a
n
c
e
p
r
o
j
e
c
t
e
d
=
1
n
∑
i
=
1
n
(
U
T
x
p
r
o
j
e
c
t
e
d
⋅
x
p
r
o
j
e
c
t
e
d
T
U
)
variance_{projected} = \frac{1}{n}\sum_{i=1}^{n} (U^{T}x_{projected} \cdot x_{projected}^{T}U)
varianceprojected=n1∑i=1n(UTxprojected⋅xprojectedTU)
因为我们只对U的方向感兴趣,如果我们将U提到式子的最前面:
v a r i a n c e p r o j e c t e d = U T ⋅ U 1 n ∑ i = 1 n ( x p r o j e c t e d ⋅ x p r o j e c t e d T ) variance_{projected} = U^{T} \cdot U\frac{1}{n}\sum_{i=1}^{n} (x_{projected} \cdot x_{projected}^{T}) varianceprojected=UT⋅Un1∑i=1n(xprojected⋅xprojectedT)
我们会发现新的这个式子中,
∑
i
=
1
n
(
x
p
r
o
j
e
c
t
e
d
⋅
x
p
r
o
j
e
c
t
e
d
T
)
\sum_{i=1}^{n} (x_{projected} \cdot x_{projected}^{T})
∑i=1n(xprojected⋅xprojectedT),这一块其实就是一个我们之前推导的当每个维度均值为0的时候,协方差矩阵的表达式,所以将上述式子再一次简化,我们最后可以得到投射之后数据的方差表示为:
v
a
r
i
a
n
c
e
p
r
o
j
e
c
t
e
d
=
U
T
⋅
S
⋅
U
variance_{projected}=U^{T}\cdot S \cdot U
varianceprojected=UT⋅S⋅U
我们的目的就是要让投射的数据点尽可能的分散,从而最大限度避免数据点重合儿导致数据丢失,所以上述式子的结果应该被最大化!很容易看出,在最大化上述式子的时候,我们应该要避免U的无穷大,所有我们可以将约束设置成 U T ⋅ U = 1 U^{T}\cdot{U}=1 UT⋅U=1。
将目标函数拉格朗日化:
(
U
,
λ
)
=
U
T
⋅
S
⋅
U
+
λ
i
(
1
−
U
T
⋅
U
)
(U,\lambda)=U^{T}\cdot S \cdot U +\lambda_{i}(1-U^{T}\cdot{U})
(U,λ)=UT⋅S⋅U+λi(1−UT⋅U)
将上述式展开可得:
L
(
U
,
λ
)
=
U
T
⋅
S
⋅
U
+
λ
i
−
λ
U
T
⋅
U
L(U,\lambda)=U^{T}\cdot S \cdot U +\lambda_{i}-\lambda U^{T}\cdot{U}
L(U,λ)=UT⋅S⋅U+λi−λUT⋅U
针对求U的偏导,并将其偏导设为0,整理可得:
S
U
=
λ
U
SU=\lambda U
SU=λU,若左乘
U
T
U^{T}
UT,我们可以的得到投射之后的方差可以表示为:
U T S U = λ U^{T}SU=\lambda UTSU=λ
这里我们不难看出,这里的方差的大小等于S的特征值 λ \lambda λ的,也就是说最大的特征值 λ \lambda λ应该对应着最大的特征向量 U U U,那么这个最大的特征向量U就是第一个主成分(first principle component),它可以反映较多的原始数据的信息。我们可以根据实际需求,选择最大的几个principle components,将它们与原数据做点积,即可实现高维数据的降维。
pca算法的简单实现
这里记录一下自己实现PCA算法的基本思路。
整个PCA算法的程序结构我设计了如下几个板块:
My_pca 类方法结构一览:
- 这里有两个必须的参数,一个需要降维的训练数据(x_tr)和目标维度(n_component)
- 在egin_pairs类方法中,主要实现获取根据训练数据得到的协方差矩阵的特征值和特征向量的配对和排序。
- 在‘get_principle_components’的函数中,我们将特征值和特征向量进行配对并排序。
- dimensionality_reduce中将多维度的数据降维成指定维度数据
- project函数将所有类函数逻辑连接起来。
在__init__中,只是简单的声明了需要的两个关键变量
def __init__(self,x_tr,n_component):
self.X = x_tr #训练数据
self.n_component = n_component#目标维度
之后求得目标数据的协方差矩阵, 解出其特征值和特征向量。即得到主成分,并将它们从大到小依次排列。
def eigen_pairs(self,x):
#find cov
covoriance=np.cov(x.T)
#eigen_value and eigen_vector
eig_vals,eig_vec = np.linalg.eig(covoriance)
#bind eigenvalues and eigenvectors
eigen_pairs= [(np.abs(eig_vals[i]), eig_vec[:, i])
for i in range(len(eig_vals))]
#sort pairs descendingly
eigen_pairs.sort(reverse=True)
return eigen_pairs
然后根据用户维度的输入,从大到小提取主成分。如果用户输入的维度大于了现有数据维度,应该抛出异常,提示输入有误。
def get_principle_components(self,eig_pairs):
if len(eig_pairs)<self.n_component:
raise Exception('n_component must be smaller than the length of eigen pairs ')
else:
list_prin = []
for i in range(self.n_component):
list_prin.append(eig_pairs[i][1])
principle_components=np.array(list_prin)
return principle_components
整理好principle components之后, 我们用它与原数据点乘即可得到降维数据。
def dimensionaliy_reduce(self, principle_components):
return self.X.dot(principle_components.T)
之后将我们的上述逻辑用一个类方法连接起来即可:
def project(self):
eign_pairs_ = self.eign_pairs(self.X)
pcs = self.get_principle_components(eign_pairs_)
X_pca = self.dimensionaliy_reduce(pcs)
return X_pca
接下来用My_pca和sklearn提供的PCA对同一数据行进降维,并在SVM模型上进行比较
首先了解数据集,wine.data有三个labels,每一条数据有13个特征,一共有177个样本。
实验设计:
- 导入,分割并标准化数据
- 用标准化之后的时间,找到效果最好的SVM的核,以及参数。
- 分别使用My_pca 和 Sklearn_pca 对数据进行降维。
- 分别将原数据和降维之后的数据代入模型预测。
- 并在训练集上做测试。
数据的导入和处理:
X_data = pd.read_csv('wine.data')
x = X_data.iloc[:, 1:]
y = X_data.iloc[:, 0] #数据的第一列是label
#在训练集上测试
x_train, _, y_train, _ = train_test_split(x, y, test_size=0.15,stratify=y, random_state=777)
scale = StandardScaler() #标准化
x_scaled_train = scale.fit_transform(x_train)
接下来试用网格搜索寻找最佳的SVM参数,并定义模型。
param_grid = [
{'kernel': ['linear'], 'C': [1, 5, 10, 50]},
{'kernel': ['rbf'], 'C': [1, 5, 10, 50], 'gamma': [0.0001, 0.0005, 0.001, 0.005, 0.01, 0.1]},
{'kernel': ['poly'], 'C': [1, 5, 10, 50], 'degree':[2,3,4], 'gamma': ['auto']}
]
grid = GridSearchCV(SVC(class_weight='balanced'), param_grid,cv=5)
grid.fit(x_scaled_train,y_train)
best_param = grid.best_params_
model = SVC(kernel=best_param['kernel'],gamma=best_param['gamma'],C=best_param['C'])
用原始数据和降维之后的数据分别训练model并查看结果。
My_pca对数据进行降维:
Mypca_train= My_pca(x_scaled_train,2).project()
使用sklearn提供的PCA对数据进行降维:
sk_pca = PCA(n_components=2)
x_train_sk_pca = sk_pca.fit_transform(x_scaled_train)
将他们分代入SVM模型中进行训练并在训练集上进行测试:
#----------------------------------------fit svm with My_pca----------------------------------------------
print('--------'*5,'My_pca','--------'*5)
model.fit(x_pca_train,y_train)
predicts = model.predict(x_pca_train)
print(classification_report(y_train,predicts))
#----------------------------------------fit svm with sklearn_pca----------------------------------------------
print('--------'*5,'sklean_pca','--------'*5)
model.fit(x_train_sk_pca,y_train)
predicts = model.predict(x_train_sk_pca)
print(classification_report(y_train,predicts))
最后的到效果如下:
总结:
PC是十分常用的降维技术,主要是通过对协方差矩阵的SVD来获取所谓的主成分。在拉格朗日化投射之后数据的方差时,不要忘记限制
U
U
U的大小,不限制会使它无限增长。