机器学习实战教程(⑤):使用PCA实战人脸降维

1.引言

在互联网大数据场景下,我们经常需要面对高维数据,在对这些数据做分析和可视化的时候,我们通常会面对「高维」这个障碍。在数据挖掘和建模的过程中,高维数据也同样带来大的计算量,占据更多的资源,而且许多变量之间可能存在相关性,从而增加了分析与建模的复杂性。

我们希望找到一种方法,在对数据完成降维「压缩」的同时,尽量减少信息损失。由于各变量之间存在一定的相关关系,因此可以考虑将关系紧密的变量变成尽可能少的新变量,使这些新变量是两两不相关的,那么就可以用较少的综合指标分别代表存在于各个变量中的各类信息。机器学习中的降维算法就是这样的一类算法。

主成分分析(Principal Components Analysis,简称PCA)是最重要的数据降维方法之一。在数据压缩消除冗余和数据噪音消除等领域都有广泛的应用。本篇我们来展开讲解一下这个算法。

2.相关概念

协方差矩阵

协方差(Covariance)目的是度量两个变量(只能为两个)线性相关的程度。
在这里插入图片描述
cov=0为可以说明两个变量线性无关,但不能证明两个变量相互独立,当cov>0时,二者呈正相关,cov<0时,二者呈负相关。

  • 协方差矩阵可以处理多维度问题。
  • 协方差矩阵是一个对称的矩阵,而且对角线是各个维度上的方差。
  • 协方差矩阵计算的是不同维度之间的协方差,而不是不同样本之间的。
  • 样本矩阵中若每行是一个样本,则每列为一个维度。
  • 假设数据是3维的,那么对应协方差矩阵为:
    在这里插入图片描述
    这里简要概括一下协方差矩阵是怎么求得的,假设一个数据集有3维特征、每个特征有m个变量,这个数据集对应的数据矩阵如下:
    在这里插入图片描述若假设他们的均值都为0,可以得到下面等式:
    在这里插入图片描述
    可以看到对角线上为每个特征方差,其余位置为两个特征之间的协方差,\frac{1}{m}XX^T
    求得的就为协方差矩阵。
    推导:
    在这里插入图片描述
    如果列是特征,公式为:
    在这里插入图片描述
'''
假设列是矩阵特征,代数里面是行表示特征
[
 [x1,y2]
 [x2,y2]
]
求协方差是
[
[cov(x,x),cov(x,y)],
[cov(y,x),cov(y,y)],
]

'''
pc=np.array([[-1,4],
            [-2,8],
            [-7,2]
            ]);
mean_pa=np.mean(pc,axis=0)
print("均值",mean_pa)
pc_zero=pc-mean_pa
print(pc_zero)
print(pc_zero.T.dot(pc_zero)/(len(pc_zero)-1))  #注意样本的话是n-1啊,全部数据集是n,否则和np.cov对不上
print("conv",np.cov(pc,rowvar=False)) #rowvar=False表示列是特征,默认行是特征

结果为:

均值 [-3.33333333  4.66666667]
[[ 2.33333333 -0.66666667]
 [ 1.33333333  3.33333333]
 [-3.66666667 -2.66666667]]
[[10.33333333  6.33333333]
 [ 6.33333333  9.33333333]]
conv [[10.33333333  6.33333333]
 [ 6.33333333  9.33333333]]

协方差矩阵一定是方形矩阵,第一矩阵列数n(特征数)决定了是(n,n)方形矩阵。

矩阵的行列式

行列式在数学中,是一个函数,其定义域为det的矩阵A,取值为一个标量,写作det(A)或|A|。无论是在线性代数、多项式理论,还是在微积分学中(比如说换元积分法中),行列式作为基本的数学工具,都有着重要的应用。

行列式等于零可以得出结论

  1. A的行向量线性相关;
  2. A的列向量线性相关;
  3. 方程组Ax=0有非零解;
  4. A的秩小于n。(n是A的阶数)
  5. A不可逆

比如,如果行列式=0,说明x1y2-x2y1=0, x1=y1/y2 * x2 说明行线性相关。

[x1,x2      
y1,y2]

[1,5
3,15] 该矩阵就是行列式=0例子。

矩阵的行列式是一个可以从方形矩阵(方阵)计算出来的特别的数。
在这里插入图片描述
这矩阵的行列式是(待会儿会解释计算方法):

3×6 − 8×4 = 18 − 32 = −14

行列式告诉我们矩阵的一些特性,这些特性对解线性方程组很有用,也可以帮我们找逆矩阵,并且在微积分及其他领域都很有用.
2×2 矩阵
2×2 矩阵 (2行和2列):
在这里插入图片描述
行列式是:

|A| = ad - bc
“A 的行列式等于 a 乘 d 减 b 乘 c”

3×3 矩阵
在这里插入图片描述
行列式是:

|A| = a(ei - fh) - b(di - fg) + c(dh - eg)
“A 的行列式等于。。。。。。”

乍看很复杂,但这是有规律的:

在这里插入图片描述
求 3×3 矩阵的行列式:

  • 把 a 乘以不在 a 的行或列上的 2×2 矩阵的行列式。
  • 以 b 和 c 也做相同的计算
  • 把结果加在一起,不过 b 前面有个负号!
    公式是(记着两边的垂直线 || 代表 “的行列式”):
    在这里插入图片描述
    更多维计算参考:
    shuxuele
    github

numpy求行列式

#行列式一定是个方形矩阵
p_pet=np.array([[-1,4],
            [-2,8]
            ]);
print(np.linalg.det(p_pet))

输出:0

特征向量和特征值

得到了数据矩阵的协方差矩阵,下面就应该求协方差矩阵的特征值和特征向量,先了解一下这两个概念,如果一个向量v是矩阵A的特征向量,那么一定存在下列等式:

在这里插入图片描述

其中A可视为数据矩阵对应的协方差矩阵,是特征向量v的特征值。数据矩阵的主成分就是由其对应的协方差矩阵的特征向量,按照对应的特征值由大到小排序得到的。最大的特征值对应的特征向量就为第一主成分,第二大的特征值对应的特征向量就为第二主成分,依次类推,如果由n维映射至k维就截取至第k主成分。

求矩阵特征值的例子
在这里插入图片描述

如果入=4
在这里插入图片描述
numpy

#求特征值和特征向量
p_eig=np.array([[1.2,0.8],
            [0.8,1.2]
            ]);
eigenvalue, featurevector =(np.linalg.eig(p_eig))
print(eigenvalue)
# 这里特征向量是进行单位化(除以所有元素的平方和的开方)的形式
# 比如入=0.4时特征向量是[-1,1],单位化[-1/开根(1**2+1**2)=-0.70710678 , -1/开根(1**2+1**2)=0.70710678]
# 比如入=2时特征向量是[1,1],单位化[1/开根(1**2+1**2)=0.70710678 , 1/开根(1**2+1**2)=0.70710678]
# 所以注意特征值2是第一列,对应的是特征向量第一列的值作为向量
print(featurevector)
#组合特征值和特征向量
eig=[(eigenvalue[i],featurevector[:,i])for i in range(len(eigenvalue))]
print(eig)    

输出

[2.  0.4]
[[ 0.70710678 -0.70710678]
 [ 0.70710678  0.70710678]]
[(2.0, array([0.70710678, 0.70710678])), (0.3999999999999997, array([-0.70710678,  0.70710678]))]

可以理解为特征向量就是原始特征的一个基础向量,最后生成一个和原始数据相同维度,行是降维维度值的矩阵,比如一个50行60列(特征)的矩阵需要降维到30个特征,特征向量会是一个(30,60)的矩阵

3. 降维实现

通过上述部分总结一下PCA降维操作的步骤:

  1. 去均值化
  2. 依据数据矩阵计算协方差矩阵
  3. 计算协方差矩阵的特征值和特征向量
  4. 将特征值从大到小排序
  5. 保留前k个特征值对应的特征向量
  6. 将原始数据的n维映射至k维中

公式手推

原始数据集矩阵,每行代表一个特征:
在这里插入图片描述
对每个特征去均值化:
在这里插入图片描述

计算对应的协方差矩阵:
在这里插入图片描述

依据协方差矩阵计算特征值和特征向量,套入公式:
在这里插入图片描述

拆开计算如下:
在这里插入图片描述

可以求得两特征值:
在这里插入图片描述

在这里插入图片描述
时,对应向量应该满足如下等式:
在这里插入图片描述

对应的特征向量可取为:
在这里插入图片描述

同理当在这里插入图片描述
时,对应特征向量可取为:
在这里插入图片描述

这里我就不对两个特征向量进行标准化处理了,直接合并两个特征向量可以得到矩阵P:
在这里插入图片描述

选取大的特征值对应的特征向量乘以原数据矩阵后可得到降维后的矩阵A:
在这里插入图片描述

综上步骤就是通过PCA手推公式实现二维降为一维的操作。

手推转自:https://juejin.cn/post/6844904177571758088

numpy实现降维

输出结果,查看gitlab

#%%
import numpy as np;
import matplotlib.pyplot as plot;
# 数学中行是表示特征
X=np.array([
    [1,1,2,4,2],
    [1,3,3,4,4]
])
# 转置为列为特征
X=X.T
print(X)
#%%

#均值归零
x_demean=np.mean(X,axis=0)
X_Zero=X-x_demean
print("均值:",x_demean)
print(X_Zero)

#%%

#计算对应的协方差矩阵,手撕或者使用np.conv
con=X_Zero.T.dot(X_Zero)/(len(X_Zero)-1)
con1=np.cov(X_Zero,rowvar=False)
print("协方差:",con,"\n",con1)
#%%

#依据协方差矩阵计算特征值和特征向量,转换参考 矩阵计算.ipynb
f_value,f_vector=np.linalg.eig(con)
eig=[(f_value[i],f_vector[:,i])for i in range(len(f_value))]
print("特征值-特征向量",eig)
#%%
#获取最大值的索引
max_f_value_index=np.argsort(f_value)[::-1][0]
print(np.array([eig[max_f_value_index][1]]).dot(X_Zero.T))
print(X_Zero.dot(np.array([eig[max_f_value_index][1]]).T))


代码部分是公式的套用,每一步后都有注释,不再过多解释。可以看到得到的结果和上面手推公式得到的有些出入,上文曾提过特征向量是可以随意缩放的,这也是导致两个结果不同的原因,eig方法计算的特征向量是归一化后的。

sklean实现降维

#%%

from sklearn.decomposition import PCA
import numpy as np

X = [[1, 1], [1, 3], [2, 3], [4, 4], [2, 4]]
X = np.array(X)
pca = PCA(n_components=1)
PCA_mat = pca.fit_transform(X)
print(PCA_mat)

这里只说一下参数n_components,如果输入的是整数,代表数据集需要映射的维数,比如输入3代表最后要映射至3维;如果输入的是小数,则代表映射的维数为原数据维数的占比,比如输入0.3,如果原数据20维,就将其映射至6维。

4. pca人脸数据降维

fetch_lfw_people人脸识别数据集是由n个人不同时间、不同角度、不同表情等图像组成的数据集;
从我这统计目录结构有5760个人
在这里插入图片描述
这是小布什部分图(530个)530
下载并抓取图库

#读取人脸数据
from sklearn.decomposition import PCA
from sklearn.datasets import fetch_lfw_people
import matplotlib.pyplot as plt

faces=fetch_lfw_people(data_home="d:/test/face",min_faces_per_person=60)
print("图片数据维度:",faces.images.shape) #数据维度:1348张照片,每张照片是一个62*47=2914的矩阵
print("图片二维数据维度:",faces.data.shape)
X=faces.data #sklearn降维算法只接受二维特征矩阵,把数据换成特征矩阵,维度是1348*2914,
X = faces.data

min_faces_per_person 提取的数据集将仅保留具有至少min_faces_per_person不同图片的人的图片
比如小布什的图片超过了60就加载出来,比如Abdullah只有3张该用户就不会加载。

如果数据403无法下载请百度,下载其他用户上传到类似百度盘数据解压到data_home指定的目录即可
输出

图片数据维度: (1348, 62, 47)
图片二维数据维度: (1348, 2914)

绘制原始图片,绘制前32个

# 在matplotlib中,整个图像为一个Figure对象。在Figure对象中可以包含一个或者多个Axes对象
# figsize代表画布的大小 3行8列表示子图axes大小
fig, axes = plt.subplots(3,8 #创建一个画布有3*8个子图
                         ,figsize = (8,4) #创建一个大小为8*4的黄布
                         ,subplot_kw = {"xticks":[],"yticks":[]} # 每个子图都不显示坐标轴
                        )
for i,ax in enumerate(axes.flat):
    ax.imshow(faces.images[i,:,:], cmap = "gray")

在这里插入图片描述
pca降维到150维

pca=PCA(n_components=150)
V1 = pca.fit_transform(X)
x_inv = pca.inverse_transform(V1)
print("逆转升维维度:",x_inv.shape)
V = pca.components_
print("降维后特征向量:",V.shape)
print("降维后数据:",V1.shape)

显示特征向量

fig, axes = plt.subplots(3,8 #创建一个画布有3*8个子图
                         ,figsize = (8,4) #创建一个大小为8*4的黄布
                         ,subplot_kw = {"xticks":[],"yticks":[]} # 每个子图都不显示坐标轴
                        )
for i,ax in enumerate(axes.flat):
    ax.imshow(V[i,:].reshape(62,47), cmap = "gray")

在这里插入图片描述
显示降维数据

fig, axes = plt.subplots(3,8 #创建一个画布有3*8个子图
                         ,figsize = (8,4) #创建一个大小为8*4的黄布
                         ,subplot_kw = {"xticks":[],"yticks":[]} # 每个子图都不显示坐标轴
                        )
for i,ax in enumerate(axes.flat):
    ax.imshow(x_inv[i].reshape(62, 47), cmap='binary_r') 

在这里插入图片描述

4. pca+knn识别手写数据

MNIST是一个手写体数字的图片数据集,该数据集来由美国国家标准与技术研究所(National Institute of Standards and Technology (NIST))发起整理,一共统计了来自250个不同的人手写数字图片,其中50%是高中生,50%来自人口普查局的工作人员。该数据集的收集目的是希望通过算法,实现对手写数字的识别。
sklearn.datasets中提供了fetch_openml的方法抓取:https://www.openml.org/search?type=data&sort=runs&status=active
免费的数据,其中mnist_784就是手写数据。
在这里插入图片描述

knn预测

抓取数据集

#%%

from sklearn.decomposition import PCA
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
from sklearn.model_selection import train_test_split
from sklearn.datasets import fetch_openml
## https://www.openml.org/可以搜索到对应数据集
mnist = fetch_openml(data_home="d:/test/face",name='mnist_784')
X, y = mnist['data'], mnist['target']
X_train, X_test, y_train, y_test=train_test_split(X,y,train_size=0.9)
#共有7万张图片,每张图片有784个特征 784开方是:28*28。
print(X.shape, y.shape)

输出:(70000, 784) (70000,)
绘制25张图片

fig, axes = plt.subplots(5,5 #创建一个画布有3*8个子图
                         ,figsize = (8,4) #创建一个大小为8*4的黄布
                         ,subplot_kw = {"xticks":[],"yticks":[]} # 每个子图都不显示坐标轴
                        )
for i,ax in enumerate(axes.flat):
    ax.imshow(X[i].reshape(28, 28), cmap = "gray")

在这里插入图片描述
使用knn训练后进行预测

#%%

from sklearn.neighbors import KNeighborsClassifier

knn=KNeighborsClassifier()
#导入time模块 训练数据将近一分钟左右
%time knn.fit(X_train,y_train)

#%%
#获取第几个模型的测试数据,用来预测
preindex=101;
plt.imshow(X_test[preindex,].reshape(28, 28), cmap = "gray");
%time print("预测的数字:",knn.predict(X_test[preindex:preindex+1,]))
plt.show()

注意由于knn计算量大,fit使用时间为:
Wall time: 1min 2s
计算准确率

#%%
#通过测试数据获取该模型的得分。
%time print(knn.score(X_test,y_test))

输出:

0.9738571428571429   准确率
Wall time: 9min 9s 用时  

用于score是使用剩余的测试数据来测试准确性,用时很长

pca降维

使用pca从784降维成100,会发现训练和求score时间大幅下降
完整代码参考:

from sklearn.neighbors import KNeighborsClassifier
from sklearn.decomposition import PCA
#保留多少个主成分维度,如果是数组是保留多少个 如果是比例 用0-1的数字,比如0.9保留90%主成分
pca = PCA(n_components=100)
#注意要transform多个 一定要调用fit方法,而不是调用两次fit_transform否则导致两次的维度不一致,fit会根据数据行算出特征的。
pca.fit(X_train,y_train)
PCA_trainmat = pca.transform(X_train)
PCA_testmat = pca.transform(X_test)
print(PCA_trainmat.shape,PCA_testmat.shape)
knn1=KNeighborsClassifier()
%time knn1.fit(PCA_trainmat,y_train)
preindex1=101;
x_inv = pca.inverse_transform(PCA_testmat) 
plt.imshow(x_inv[preindex1,].reshape(28, 28), cmap = "gray");
plt.show()

%time print("预测的数字:",knn1.predict(PCA_testmat[preindex1:preindex1+1,]))
#%%

%time print(knn1.score(PCA_testmat,y_test))

在这里插入图片描述

  • 2
    点赞
  • 26
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值