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参数给出了要降到的维度),然后能清楚看到三个鸢尾花的类别,效果很好。
PCA的另类用法:
在著名的神经网络论文AlexNet中曾提出了一种使用PCA进行数据增强的算法,称之为PCA jitter。所谓数据增强,就是通过对训练数据人为的加噪声,提高模型鲁棒性与预测能力的方法。从观感来讲,PCA jitter近似对光效的变化,对比直接对RGB通道进行加噪声会获得更少的颜色的失真。
而这个方法的主要流程就是:对于一张图片,我们把每个像素看做一个数据,其有三个特征(即R,G,B三通道颜色),把每个特征标准化(均值为0,方差为1,即对每个数据减去均值后除以标准差即可),然后计算出这个矩阵的协方差矩阵,把协方差矩阵对角化,依据协方差矩阵对角化的结果对原图进行抖动,具体原理可以参考论文,这里给出一种实现:
def solve(val,fea): ran=np.random.normal(0,50,3) tval=val*ran D=tval@fea re=np.zeros((1080,1080,3)) for i in range(0,1080): for j in range(0,1080): for k in range(0,3): re[i][j][k]=np.clip(img[i][j][k]+D[k],0,255) plt.imshow(re/255) l=[] for i in img: for j in i: l.append(j) data=(np.array(l)).T data=data.astype(np.float64) for i in range(0,3): ave=data[i,:].mean() std=data[i,:].std() for j in range(0,len(data[i,:])): data[i][j]=(data[i][j]-ave)/std Q=data@data.T/len(data[i,:]) val,fea=np.linalg.eig(Q) solve(val,fea)
这里的抖动是从均值为0,方差为50的正态分布里抽样得到的,这个方差可以自己设定,根据效果选取。实际上抖动的过程就是如果设特征值为$\alpha=(\lambda_{1},\lambda_{2},\lambda_{3})$,特征向量构成的矩阵为$P$,随机数为$r$,那么我们计算$\beta=r\alpha P$得到一个$1*3$的向量即为RGB三通道各自的抖动值。
(上述过程可能有误,欢迎大佬指出,不胜感激!)
而实现结果:如果原图是这样的
那么一个随机的情况是这样的:
类似于加了清晨朦胧的光效?可以看到还是有一定效果的。