Python教程

python机器学习——PCA降维算法

本文主要是介绍python机器学习——PCA降维算法,对大家解决编程问题具有一定的参考价值,需要的程序猿们随着小编来一起学习吧!

背景与原理:

PCA(主成分分析)是将一个数据的特征数量减少的同时尽可能保留最多信息的方法。所谓降维,就是在说对于一个$n$维数据集,其可以看做一个$n$维空间中的点集(或者向量集),而我们要把这个向量集投影到一个$k<n$维空间中,这样当然会导致信息损失,但是如果这个$k$维空间的基底选取的足够好,那么我们可以在投影过程中尽可能多地保留原数据集的信息。

数据降维的目的在于使得数据更直观、更易读、降低算法的计算开销、去除噪声。

接下来我们讨论下如何选取$k$维子空间:

假设原数据集有$m$条数据,每条数据有$n$维,那么可以将其拼成一个$n*m$的矩阵M,而我们想投影到的$k$维空间的一个单位正交基底为$(p_{1},...,p_{k})$,那么我们想把这$m$维向量投影到这个空间中实际上就是进行一次矩阵乘法$\begin{pmatrix} p_{1}\\p_{2}\\...\\p_{k} \end{pmatrix} M$

这个道理是简单易懂的。对于一个向量$\alpha$,其在另一个向量$\beta$方向上的投影是$\dfrac{\alpha \cdot \beta }{|\beta|}$(高中数学)

如果$|\beta|$是一个单位向量,那么这个投影即为$\alpha \cdot \beta=\beta^{T} \alpha$,于是投影到$k$个单位向量为基底的空间中的情况即如上述所示。

因此我们要找到这$k$个单位向量作为基底,然后拼出$P=\begin{pmatrix} p_{1}\\p_{2}\\...\\p_{k} \end{pmatrix}$即可。

那么怎么找呢?我们考虑降维之后我们需要什么,由于我们在降维之后要尽可能多地保留原始信息,因此降维之后的数据要提供最大的信息量,那么这个信息量在这里可以用数据的方差来反映,方差越大,数据的离散程度越高,那么数据的自身特征保留的就越好(个人理解:PCA降维的目的在于突出数据的个体特征,减少信息损失,而如果降维之后数据离散程度低,意味着这些数据全都堆在一起,那数据的特征体现的就不明显了——所有数据全都差不多,这样个体信息保留的就不好了)。

对于一个特征,在$m$组数据中的方差为$\sigma^{2}=\dfrac{1}{m}\sum_{i=1}^{m}(x_{i}-\overline{x})^{2}$,为了便于讨论,我们对所有特征零均值化(即把每个$x_{i}$预先减去$\overline{x}$),这样一个特征的方差即为$\sigma^{2}=\dfrac{1}{m}\sum_{i=1}^{m}x_{i}^{2}$

但是降维过程中只考虑方差是不够的——如果我们发现两个特征之间有很强的线性相关性,那么这两个特征其实差别就不大了,我们当然不需要同时保留这两个特征,因此我们还希望降维之后任意两个特征的协方差($cov(a,b)=\dfrac{1}{m}\sum_{i=1}^{m}a_{i}b_{i}$,因为已经进行过零均值化了)为零,也就是在说我们选取的子空间的基底一定是正交的。

那么现在的问题就转化成了:对于一个$n$维$m$组数据的$n*m$数据矩阵$X$,我们希望将其投影到$n$维空间的一个$k$维子空间中,因此我们要找到$k$个单位正交基$(p_{1},...,p_{k})$,而如果这$k$个单位正交基构成的矩阵$P=\begin{pmatrix} p_{1}\\p_{2}\\...\\p_{k} \end{pmatrix}$,那么投影过程即为$Y=PX$,$Y$即为降维后所得的$k$维数据集

而结合上述讨论,我们希望$Y$各个特征的方差最大,同时$Y$的两特征的协方差为零,这怎么操作呢?

对于一个$n$维有$m$组数据的$n*m$矩阵$X$,我们考察$C=\dfrac{1}{m}XX^{T}$,那么我们看到如果$X=\begin{pmatrix} x_{11} & x_{12} & ... & x_{1m}\\...\\x_{n1} & x_{n2} &... & x_{nm}\end{pmatrix}$,我们有:

$XX^{T}=\begin{pmatrix} \sum_{i=1}^{m}x_{1i}^{2} & \sum_{i=1}^{m} x_{1i}x_{2i} &...& \sum_{i=1}^{m}x_{1i}x_{ni}\\...\\ \sum_{i=1}^{m}x_{ni}x_{1i} & \sum_{i=1}^{m} x_{ni}x_{2i} &...& \sum_{i=1}^{m}x_{ni}^{2}\end{pmatrix}$

我们称$\dfrac{1}{m}XX^{T}$为协方差矩阵,因为我们看到按照我们上面的解释,这个矩阵是一个实对称矩阵,其主对角线上的元素是一个特征维度的方差,而其余位置上的元素是两个对应特征的协方差!

那么我们的目的是要最大化主对角线上的元素,同时让其余位置上的元素为$0$,那么我们进行的不就是实对称矩阵的正交相似对角化嘛!

形式化地解释一下:我们设$Y$的协方差矩阵为$D$,那么我们希望$D$是一个对角矩阵,同时$D$的主对角线上的元素要尽可能大,那么我们有:

$D=\dfrac{1}{m}YY^{T}=\dfrac{1}{m}(PX)(X^{T}P^{T})=P(\dfrac{1}{m}XX^{T})P^{T}$

那么我们实际进行的不就是把$C=\dfrac{1}{m}XX^{T}$这个协方差矩阵正交相似对角化嘛!

至于我们希望主对角线元素尽可能大,那我们就选取$C$的前$k$大的特征值组成$D$就好了嘛,而此时的$P$就对应于前$k$大的特征值对应的$k$个正交的特征向量构成的矩阵。

那么我们的算法步骤如下:

对于$n$行$m$列的矩阵$X$,我们解释成其有$m$组数据,每组数据有$n$个特征,现在我们欲将其变成$k*m$的矩阵$Y$,表示降维后每组数据只有$k$个特征。

(1)零均值化:对$X$的每个元素,减去自己所在的均值(即我们是逐特征操作,一行对应于同一个特征,不要搞错这一点)

(2)计算协方差矩阵$C=\dfrac{1}{\textbf{m}}XX^{T}$

(3)对协方差矩阵对角化$C=P\Sigma P^{T}$,找到其单位正交的特征向量$e_{1},...,e_{n}$

(4)选取最大的$k$个特征值对应的特征向量$e_{i_{1}},...,e_{i_{k}}$,拼成一个变换矩阵$P_{k}=\begin{pmatrix} e_{i_{1}}\\e_{i_{2}}\\...\\e_{i_{k}} \end{pmatrix}$

(5)降维后的数据即为$Y=P_{k}X$

(6)如果希望根据降维后的数据集$Y$近似还原原数据集$\hat{X}$,我们有$\hat{X}=P_{k}^{T}Y$(这里的逻辑是如果我们不降维,那么$P_{k}=P$就是一个正交矩阵,那么$P^{T}P=I$,相当于此时数据集没有损失,那么类比这个过程就能导出近似还原方法$\hat{X}=P_{k}^{T}Y$)

代码实现:

import numpy as np
from sympy.matrices import Matrix,GramSchmidt
np.random.seed(1)

x = 7*np.random.rand(100)
y = 0.5*x + 1 + 3*np.random.rand(100)

X = np.hstack([x.reshape(100, 1), x.reshape(100, 1), y.reshape(100, 1), x.reshape(100, 1)])

def centerData(X):
    X = X.copy()
    X -= np.mean(X, axis=0)
    return X

X = centerData(X)print(X[7][2])
C= (np.transpose(X)@X)/100
val,fea=np.linalg.eig(C)
dic=dict()
for i in range(0,4):
    dic[val[i]]=fea[:,i]
val=abs(np.sort(-val))
P=np.vstack([dic[val[0]],dic[val[1]]])
Y=X@P.T
reconstruct_X=Y@P
print(reconstruct_X[7][2])

值得注意的问题是这段代码中生成的数据每组数据是一个行向量,每列对应于一个特征,因此所有的计算和上面的推导都构成一个转置。

此外这里使用了numpy里面的linalg.eig方法用来求一个实对称矩阵的特征值和特征向量,返回的val是特征值,fea是特征向量,要特别说明的是不出意外的情况下这里的val都是按从大到小排序的,而fea实际上是一个矩阵,这个矩阵每个列向量对应于一个特征值,因此一定要注意选取的方法。上面代码使用前两个特征向量作为主成分恢复原矩阵,可以看到恢复效果还是不错的。

import numpy as np
from sklearn.decomposition import PCA
import matplotlib.pyplot as plt
from sklearn.datasets import load_iris

data=load_iris()
X=data.data
Y=data.target

pca=PCA(n_components=2)

X_2d=pca.fit_transform(X)

plt.scatter(X_2d[:,0],X_2d[:,1],c=Y)
plt.show()

当然,PCA也可以直接使用sklearn里面的包,上述代码加载了经典的鸢尾花数据集,然后进行PCA降维(降成二维,这个n_components参数给出了要降到的维度),然后能清楚看到三个鸢尾花的类别,效果很好。

这篇关于python机器学习——PCA降维算法的文章就介绍到这儿,希望我们推荐的文章对大家有所帮助,也希望大家多多支持为之网!