机器学习系列笔记五:主成分分析PCA[上]

文章目录

机器学习系列笔记五:主成分分析PCA[上]IntroductionPCA的原理特征降维的引入PCA的求解数学推导 梯度上升法解决主成成分分析问题代码实现梯度上升法求解主成分demean梯度上升法 求解其他主成分求解前n个主成分对求解n个主成分方法的最终封装 高维数据向低维数据映射

Introduction

主成分分析是一个非监督的机器学习算法：

主要用于数据的降维，通过降维，可以发现更便于人类理解的特征，比如在人脸识别中，就必须通过降维对数据进行处理。其他应用：可视化、去噪

定义： PCA（Principal Component Analysis）是一种常用的数据分析方法。PCA通过线性变换将原始数据变换为一组各维度线性无关的表示，可用于提取数据的主要特征分量，常用于高维数据的降维。PCA的思想是将n维特征映射到k维空间上k<n，这k维特征是全新的正交特征，是重新构造出来的k维特征，而不是简单地从n维特征中去除其余n−k维特征。

主成分分析法本质上是将特征从一组坐标系转移到了另一组坐标系

其实学过信号与系统的话，PCA和所谓的傅里叶变换，拉普拉斯变换，Z变换的作用和原理是一样的

PCA的原理

特征降维的引入

假设把右上角所表示的是特征1和特征2所组成的散点，X轴为特征1，y轴为特征2，然后通过对X轴Y轴的映射可以变为下面的两种形式，这就是一个特别简单的降维方法，对于图2，我们只需要去掉特征2，对于图3，我们是去掉特征1的考量。

然后根据特征点的密度可以推算各个特征的可区分度，故根据上图所示方案2的效果显然要比方案3的要好。

基于此，就有了下面的方案3：

用这个方式我们将二维的特征映射到了一维，只不过这个"轴"不再是传统的XY轴，利用这种方式得到的映射结果，更加趋近于原始特征的分布。

因此，问题就转换为了如何找到让这个样本(特征)间间距最大的轴

PCA的求解数学推导

首先我们需要定义样本间间距：

使用方差的形式: $$Var(x)=\frac{1}{m}\sum _{i=1}^m(x_i-\bar{x}{)}^2$$ 问题就变成了找到一个轴，使得样本空间的所有点映射到这个轴后，方差最大

从而通过主成分分析使得样本空间降维的步骤如下：

step.1 ：将样例的均值归为0(demean处理)

step.2：求一个轴的方向$w=(w_1,w_2)$,使得我们所有的样本在映射到$w$ 之后,有：

$Var\left(X_{project}\right)=\frac{1}{m}{\sum }_{i=1}^m||X_{project}^{\left(i\right)}-{\bar{X}}_{project}|{|}^2$ 最大，而$\bar{X}=0$,故我们的目标是： $$max\{Var\left(X_{project}\right)\}=max\{\frac{1}{m}\sum _{i=1}^m||X_{project}^{\left(i\right)}|{|}^2\}$$ 其中，$||X_{project}^{(i)}||$可以根据向量乘法的定义，用如下定义所表示： $$||X_{project}^{(i)}||=X^{(i)}\cdot w$$ 目标就变为了： $$\underset{w}{\max }\left\{Var\left(X_{project}\right)\right\}=\underset{w}{\max }\left\{\frac{1}{m}\sum _{i=1}^m{\left(X^{\left(i\right)}\cdot w\right)}^2\right\}$$

显然主成分分析就变为了一个目标函数的最优化问题，由于是使得目标函数最大化，所以可以采用梯度上升法解决该优化问题。

梯度上升法解决主成成分分析问题

根据上面的分析，我们可以将目标函数定义为： $$f(X,w)=\frac{1}{m}\sum _{i=1}^m(X_1^i{w}_1+X_2^i{w}_2+...+X_n^i{w}_n{)}^2$$ 求$w$,使得目标函数$f(X,w)$ 最大。

注意在目标函数中X是已知量，而w才是未知量

$$\nabla f=\left[\begin{array}{c}\frac{\partial f}{\partial w_1}\\ \frac{\partial f}{\partial w_2}\\ ⋮\\ \frac{\partial f}{\partial w_n}\end{array}\right]=\frac{2}{m}{\left[\begin{array}{c}{\sum }_{i=1}^m\left(X^{\left(i\right)}w\right)X_1^{\left(i\right)}\\ {\sum }_{i=1}^m\left(X^{\left(i\right)}w\right)X_2^{\left(i\right)}\\ ⋮\\ {\sum }_{i=1}^m\left(X^{\left(i\right)}w\right)X_m^{\left(i\right)}\end{array}\right]}_{m\times 1}$$

向量化： $$\nabla f=\frac{2}{m}{\left[\begin{array}{c}{\sum }_{i=1}^m\left(X^{\left(i\right)}w\right)X_1^{\left(i\right)}\\ {\sum }_{i=1}^m\left(X^{\left(i\right)}w\right)X_2^{\left(i\right)}\\ ⋮\\ {\sum }_{i=1}^m\left(X^{\left(i\right)}w\right)X_n^{\left(i\right)}\end{array}\right]}_{n\times 1}=>\nabla f=\frac{2}{m}\left(X^{\left(1\right)}w,X^{\left(2\right)}w,...,X^{\left(m\right)}w\right)\cdot \left[\begin{array}{ccc}{X}_1^{\left(1\right)}& X_2^{\left(1\right)}\cdots & X_n^{\left(1\right)}\\ X_1^{\left(2\right)}& X_2^{\left(2\right)}\cdots & X_n^{\left(2\right)}\\ ⋮& \ddots & ⋮\\ X_1^{\left(m\right)}& X_2^{\left(m\right)}\cdots & X_n^{\left(m\right)}\end{array}\right]$$ $$=\frac{2}{m}\cdot {\left(Xw\right)}^T\cdot X=\frac{2}{m}\cdot X^T(Xw)$$

对其中的维度做一个说明：

m表示样本数: $X=(X^{(1)},X^{(2)},...,X^{(m)}{)}^T$

n表示每个样本所包含的特征个数: $X^{(i)}=(X_1^{(i)},X_2^{(i)},...,X_n^{(i)})$

w表示最佳映射方向，因与特征个数有关，所以$w=(w_1,w_2,...,w_n)$

有了上面的推导和结论，我们很容易就能利用梯度上升法把主成分求解出来。

代码实现梯度上升法求解主成分

import numpy as np import matplotlib.pyplot as plt X = np.empty((100,2)) X[:,0] = np.random.uniform(0.,100.,size=100) # 特征1 X[:,1] = 0.75 * X[:,0]+3.+np.random.normal(0,10,size=100) # 特征2 plt.scatter(X[:,0],X[:,1]) plt.show()

demean

def demean(X): return X - np.mean(X,axis=0) # 每一列特征分别减去其对应的列平均值 X_demean = demean(X) plt.scatter(X_demean[:,0],X_demean[:,1]) # 对比线 x_0 = X_demean[:,0] y_0 = np.zeros(100) x_1 = np.zeros(100) y_1 = X_demean[:,1] plt.plot(x_0,y_0,color='r') plt.plot(x_1,y_1,color='r') plt.axis([-50,50,-50,50]) plt.show()

可以看到demean后的特征分布均值都在0的位置

np.mean(X_demean[:,0]) -1.3073986337985844e-14 np.mean(X_demean[:,1]) -1.1510792319313624e-14

梯度上升法

def f(w,X): """目标函数""" return np.sum((X.dot(w))**2)/len(X) def df_math(w,X): """目标函数的梯度""" return X.T.dot(X.dot(w)) * 2 / len(X) def df_debug(w,X,epsilon = 0.0001): """测试梯度上升法的样例函数""" res = np.empty(len(w)) for i in range(len(w)): w_1 = w.copy() w_1[i] += epsilon w_2 = w.copy() w_2[i] -= epsilon res[i] = (f(w_1,X) - f(w_2,X)) / (2 * epsilon) return res def directionize(w): """单位化向量""" return w/np.linalg.norm(w) def gradient_ascent(df,X,initial_w,eta,n_iters = 1e4,epsilon = 1e-8): """梯度上升法""" w = directionize(initial_w) cur_iter = 0 while cur_iter < n_iters: gradient = df(w,X) last_w = w w = w + eta*gradient w = directionize(w) # 注意：每次求一个单位方向 if (abs(f(w,X)-f(last_w,X))<epsilon): break cur_iter += 1 return w

测试：

initial_w = np.random.random(X.shape[1]) # 注意：不能用0向量开始 eta = 0.001 initial_w array([0.4547083 , 0.32103896]) # 注意：不能使用StandardScaler标准化数据，标准化后的数据方差是为一的，这与我们想通过最大方差找到映射轴的思想是冲突的 gradient_ascent(df_debug,X_demean,initial_w,eta) array([0.79226227, 0.61018071]) gradient_ascent(df_math,X_demean,initial_w,eta) array([0.79226227, 0.61018071])

可视化结果

w = gradient_ascent(df_math,X_demean,initial_w,eta) plt.scatter(X_demean[:,0],X_demean[:,1]) plt.plot([0, w[0]*30],[0, w[1]*30], color ='r',label='the first Principal Component') plt.legend() plt.show()

我们看看比较极端的情况下PCA的表示

X2 = np.empty((100,2)) X2[:,0] = np.random.uniform(0.,100.,size=100) # 特征1 X2[:,1] = 0.75 * X2[:,0]+3. # 特征2，斜率是3/4 plt.scatter(X2[:,0],X2[:,1]) plt.show()

X2_demean = demean(X2) w2 = gradient_ascent(df_math,X2_demean,initial_w,eta) plt.scatter(X2[:,0],X2[:,1]) plt.plot([0, w2[0]*30],[0, w2[1]*30], color ='r',label='the first Principal Component') plt.legend() plt.show()

求解其他主成分

求出了第一主成分以后，如何求解下一个主成分：

把数据$X^{(i)}$进行改变，将数据的第一个主成分上的分量去掉(向量减法)后得到$X^{\prime (i)}$ $$X^{\prime (i)}=X^{(i)}-X_{project}^{(i)}$$

$$X_{project}^{(i)}=X^{(i)}.w\ast w$$

$\cdot$ 表示矩阵内积.dot* 表示矩阵内各对应位置相乘

得到了新的数据$X^{\prime (i)}$ ,对于这个新的数据重新求解第一主成分作为原始数据的第二主成分

由此类推，可以求解第三、第三…主成分

求解前n个主成分

import numpy as np import matplotlib.pyplot as plt X = np.empty((100,2)) X[:,0] = np.random.uniform(0,100,size=100) X[:,1] = 0.75 * X[:,0] + 3 + np.random.normal(0,10,size=100) def demean(X): return X - np.mean(X,axis=0) X = demean(X) plt.scatter(X[:,0],X[:,1]) plt.show()

def f(w,X): """目标函数""" return np.sum((X.dot(w))**2)/len(X) def df(w,X): """目标函数的梯度""" return X.T.dot(X.dot(w)) * 2 / len(X) def directionize(w): """单位化向量""" return w/np.linalg.norm(w) def first_componet(X,initial_w,eta,n_iters = 1e4,epsilon = 1e-8): """梯度上升法求解给定X的第一主成分""" w = directionize(initial_w) cur_iter = 0 while cur_iter < n_iters: gradient = df(w,X) last_w = w w = w + eta*gradient w = directionize(w) # 注意：每次求一个单位方向 if (abs(f(w,X)-f(last_w,X))<epsilon): break cur_iter += 1 return w initial_w = np.random.random(X.shape[1]) eta = 0.01 w = first_componet(X,initial_w,eta) w array([0.78343282, 0.62147648])

去除原始数据的第一主成分分量，得到新的数据X2

X2 = np.empty(X.shape) # for i in range(len(X)): # X2[i] = X[i] - X[i].dot(w) * w # 将上式向量化 ,X为(m,n),w为(n,1),X.dot(w)为(m,1),X.dot(w).reshape(-1,1)*w用到了numpy的广播机制 X2 = X-X.dot(w).reshape(-1,1)*w plt.scatter(X2[:,0],X2[:,1]) plt.show()

w2 = first_componet(X2,initial_w,eta) w2 array([-0.62147116, 0.78343704]) w.dot(w2) # 趋紧于0，说明两个向量垂直 6.792886741657789e-06 plt.scatter(X2[:,0],X2[:,1]) plt.plot(w2,color='r')

对求解n个主成分方法的最终封装

def first_n_components(n,X,eta=0.01,n_iters =1e4,epsilon=1e-8): X_pca = X.copy() X_pca = demean(X_pca) res = [] # 通过n次循环拿到n个主成分 for i in range(n): initial_w = np.random.random(X_pca.shape[1]) w = first_componet(X_pca,initial_w,eta) res.append(w) X_pca = X_pca-X_pca.dot(w).reshape(-1,1)*w return res first_n_components(2,X) [array([0.78343286, 0.62147643]), array([-0.62147234, 0.78343611])]

高维数据向低维数据映射

将上图所示的功能封装到一个类中：

import numpy as np class PCA: def __init__(self, n_components): """初始化PCA""" assert n_components >= 1, \ "n_components must be valid" self.n_components = n_components self.components_ = None def __repr__(self): return "PCA(n_components=%d)" % self.n_components def fit(self, X, eta=0.01, n_iters=1e4): """根据训练集X获得n_components个主成分""" assert X.shape[1] >= self.n_components, \ "n_components must be greater than the feature number of X" def demean(X): return X - np.mean(X, axis=0) def f(w, X): """目标函数""" return np.sum((X.dot(w)) ** 2) / len(X) def df(w, X): """目标函数的梯度""" return X.T.dot(X.dot(w)) * 2 / len(X) def directionize(w): """单位化向量""" return w / np.linalg.norm(w) def first_componet(X, initial_w, eta, n_iters=1e4, epsilon=1e-8): """梯度上升法求解给定X的第一主成分""" w = directionize(initial_w) cur_iter = 0 while cur_iter < n_iters: gradient = df(w, X) last_w = w w = w + eta * gradient w = directionize(w) # 注意：每次求一个单位方向 if (abs(f(w, X) - f(last_w, X)) < epsilon): break cur_iter += 1 return w X_pca = demean(X) # components_相当于W_k矩阵：W(k,n) self.components_ = np.empty(shape=(self.n_components, X.shape[1])) res = [] # 通过n次循环拿到n个主成分 for i in range(self.n_components): initial_w = np.random.random(X_pca.shape[1]) w = first_componet(X_pca, initial_w, eta, n_iters) self.components_[i, :] = w X_pca = X_pca - X_pca.dot(w).reshape(-1, 1) * w return self def transform(self, X): """将给定的X,映射到各个主成分分量中""" assert X.shape[1] == self.components_.shape[1] return X.dot(self.components_.T) def inverse_transform(self, X): """将给定的低维特征空间X转换回高维空间中""" assert X.shape[1] == self.components_.shape[0] return X.dot(self.components_)

测试

import numpy as np import matplotlib.pyplot as plt from relatedCode.PCA import PCA if __name__ == '__main__': # 准备数据 X = np.empty((100, 2)) X[:, 0] = np.random.uniform(0, 100, size=100) X[:, 1] = 0.75 * X[:, 0] + 3 + np.random.normal(0, 10, size=100) # 寻找两个主成分 n_components = 2 pca = PCA(n_components) pca.fit(X) print("the priciple {} component:{}".format(n_components, pca.components_)) # 寻找1个主成分 pca2 = PCA(1) pca2.fit(X) print("the priciple {} component:{}".format(1, pca2.components_)) ## 根据寻找到的1个主成分对数据进行降维 X_reduction = pca2.transform(X) print("after Dimensionality Reduction,X_reduction.shape=", X_reduction.shape) ## 根据降维后的矩阵恢复为高维的矩阵 X_restore = pca2.inverse_transform(X_reduction) print("after Dimensionality Restore,X_restore.shape=", X_restore.shape) # 可视化 plt.scatter(X[:, 0], X[:, 1], color='b', alpha=0.5) plt.scatter(X_restore[:, 0], X_restore[:, 1], color='r', alpha=0.5) plt.legend(['original', 'after PCA process']) plt.savefig("PCA_contrast_fig") plt.show()

结果：

C:\Users\ASUS\Anaconda3\envs\tf2\python.exe C:/Users/ASUS/Desktop/思维导图/机器学习编程实践/day04/test/PCATest.py the priciple 2 component:[[ 0.7653877 0.64356948] [-0.64356549 0.76539105]] the priciple 1 component:[[0.76538767 0.64356951]] after Dimensionality Reduction,X_reduction.shape= (100, 1) after Dimensionality Restore,X_restore.shape= (100, 2) Process finished with exit code 0