Python机器学习算法实现
Author:louwill
Machine Learning Lab
作为一种常见的多元统计分析方法,主成分分析法(Principal Component Analysis,PCA)也是一种经典的无监督学习算法。PCA通过正交变换将一组由线性相关变量表示的数据转换为少数几个由线性无关变量表示的数据,这几个线性无关的变量就是主成分。PCA通过将高维数据维度减少到少数几个维度,本质上属于一种数据降维方法,也可以用来探索数据的内在结构。
PCA原理与推导
PCA的基本想法如下:首先对需要降维的数据各个变量标准化(规范化)为均值为0、方差为1的数据集。然后对标准化后的数据进行正交变换,将原来的数据转换为若干个线性无关的向量表示的新数据。这些新向量表示的数据不仅要求相互线性无关,而且需要是所包含的信息量最大。
PCA使用方差来衡量新变量的信息量大小,按照方差大小排序依次为第一主成分、第二主成分等。下面对PCA原理进行简单推导。假设原始数据为 维随机变量 ,其均值向量为 :
其协方差矩阵为
:
现假设由
维随机变量
到
维随机变量
的线性变换:
其中
。
然后可得变换后的随机变量 的统计特征:
当上述线性变换满足如下条件时,变换后的 分别为 的第一主成分、第二主成分以及第 主成分。
变换的系数向量 为单位向量,有 。
变换后的变量 与 线性无关,即 。
变量 是 所有线性变换中方差最大的, 是与 不相关的所有线性变换中方差最大的。
以上条件给出了求解PCA主成分的方法。根据约束条件和优化目标,我们可以使用拉格朗日乘子法来求出主成分。先来求解第一主成分。
根据前述条件,我们可以形式化第一主成分的数学表达为:
定义拉格朗日函数:
将上述拉格朗日函数对
求导并令为0可得:
所以 是 的特征值, 是对应的单位特征向量。假设 是 的最大特征值 对应的单位特征向量,则 和 均为上述优化问题的最优解。所以 为第一主成分,其方差为对应协方差矩阵的最大特征值:
同样方法也可求第二主成分。但约束条件中需要加上与第一主成分不相关的约束,具体推导这里略过。按上述方法可一致计算到第 个主成分,且第 个主成分的方差等于 的第 个特征值:
假设原始数据为 行 列的数据 ,梳理PCA的计算流程如下:
对 行 列的数据按照列进行均值为0方差为1的标准化;
计算标准化后的 的协方差矩阵 ;
计算协方差矩阵的特征值和对应的特征向量;
将特征向量按对应特征值大小从上到下按行排列成矩阵,取前 行构成的矩阵 ;
即为PCA降维后的 维数据。
PCA的numpy实现
虽然sklearn中提供了PCA降维的API,但其背后算法是用SVD来实现的。numpy模块下提供了强大的矩阵运算函数,下面我们用numpy来实现一个PCA算法。
import numpy as np
class PCA():
# 计算协方差矩阵
def calculate_covariance_matrix(self, X):
m = X.shape[0]
# 数据标准化
X = X - np.mean(X, axis=0)
return 1 / m * np.matmul(X.T, X)
def pca(self, X, n_components):
# 计算协方差矩阵
covariance_matrix = self.calculate_covariance_matrix(X)
# 计算协方差矩阵的特征值和对应特征向量
eigenvalues, eigenvectors = np.linalg.eig(covariance_matrix)
# 对特征值排序
idx = eigenvalues.argsort()[::-1]
# 取最大的前n_component组
eigenvectors = eigenvectors[:, idx]
eigenvectors = eigenvectors[:, :n_components]
# Y=PX转换
return np.matmul(X, eigenvectors)
由上述代码可以看到,基于numpy我们仅需要数行核心代码即可实现一个PCA降维算法。下面以sklearn digits数据集为例看一下降维效果:
from sklearn import datasets
import matplotlib.pyplot as plt
import matplotlib.cm as cmx
import matplotlib.colors as colors
# 导入sklearn数据集
data = datasets.load_digits()
X = data.data
y = data.target
# 将数据降维到2个主成分
X_trans = PCA().pca(X, 2)
x1 = X_trans[:, 0]
x2 = X_trans[:, 1]
# 绘图展示
cmap = plt.get_cmap('viridis')
colors = [cmap(i) for i in np.linspace(0, 1, len(np.unique(y)))]
class_distr = []
# 绘制不同类别分别
for i, l in enumerate(np.unique(y)):
_x1 = x1[y == l]
_x2 = x2[y == l]
_y = y[y == l]
class_distr.append(plt.scatter(_x1, _x2, color=colors[i]))
# 图例
plt.legend(class_distr, y, loc=1)
# 坐标轴
plt.suptitle("PCA Dimensionality Reduction")
plt.title("Digit Dataset")
plt.xlabel('Principal Component 1')
plt.ylabel('Principal Component 2')
plt.show();
参考资料:
统计学习方法第二版
多元统计分析 何晓群
https://github.com/RRdmlearning/Machine-Learning-From-Scratch
往期精彩:
数学推导+纯Python实现机器学习算法28:奇异值分解SVD
数学推导+纯Python实现机器学习算法17:XGBoost
数学推导+纯Python实现机器学习算法13:Lasso回归
数学推导+纯Python实现机器学习算法10:线性不可分支持向量机
一个算法工程师的成长之路
长按二维码.关注机器学习实验室
喜欢您就点个在看!