PCA简介
PCA（PrincipalComponents Analysis）即主成分分析，是一个非监督的机器学习算法，它是最常用的降维方法之一，通过正交变换将一组可能存在相关性的变量数据转换为一组线性不相关的变量，转换后的变量被称为主成分。
例如下图所示，样本有2个特征，现在对该样本进行降维处理。首先考虑直接选择特征1或者特征2降维，经过降维后的样本由2维降到1维，如图所示。
可以看出剔除特征2降维比剔除特征1降维的样本间的间距更大，即样本可区分度更大。那是否还有其他的映射方式，使得映射后样本的间距更大，事实上还可以选择某个轴线，例如下图所示，样本映射到该轴线之后，有更大的间距。
PCA降维的思想就是寻找某个轴线，使得样本映射到该轴线上后使得样本区分度更大。衡量可区分度的指标即方差（事实上，方差背后的意义就是数据的稀疏程度），现在的问题就是如何求得该轴线，使得方差的值最大。
求解思路
用方差来定义样本的间距，方差越大表示样本分布越稀疏，方差越小表示样本分布越密集，方差的公式如下。
在求解最大方差前，为了方便计算，可以先对样本进行demean（去均值）处理，即减去每个特征的均值，这种处理方式不会改变样本的相对分布（效果就像坐标轴进行了移动）。去均值后，样本x每个特征维度上的均值都是0，方差的公式转换成图示的公式。
对于2个维度的样本，我们想要求一个轴的方向 w = (w1, w2) 使得我们所有的样本，映射到w以后方差最大，有：
如下图所示，原始样本和映射后的样本有如下的空间几何关系。w为单位向量，所以w的模为1。
样本x为已知，现在的目标转换成如下:
梯度上升发求解
公式推导
代码实现
初始样本准备
import numpy as np
import matplotlib.pyplot as plt
X = np.empty((200, 2))
X[:,0] = np.random.uniform(0., 100., size=200)
X[:,1] = 0.75 * X[:,0] + 1.5 + np.random.normal(0, 10., size=200)
plt.scatter(X[:,0], X[:,1])
plt.show()
开始对该样本进行降维处理
def demean(X):
# axis=0 是 矩阵X - X每一列的均值（即axis=0） 矩阵运算
return X - np.mean(X, axis=0)
# pca 目标函数推导公式
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 direction(w):
# 由于w不一定是一个单位向量，把它转化成一个单位向量， w/w.模
return w / np.linalg.norm(w)
def gradient_ascent(df, X, initial_w, eta, n_iters = 1e4, epsilon=1e-8):
w = direction(initial_w)
cur_iter = 0
while cur_iter < n_iters:
gradient = df(w, X)
last_w = w
w = w + eta * gradient
w = direction(w) # 注意1：每次求一个单位方向，方便能得到一个合适的值方便计算
if(abs(f(w, X) - f(last_w, X)) < epsilon):
break
cur_iter += 1
return w
对样本X进行降维
X_mean = demean(X) #去均值处理
initial_w = np.random.random(X.shape[1]) # 注意：不能用0向量开始
eta = 0.001 # 初始化步长
w = gradient_ascent(df_math, X_mean, initial_w, eta) #降维并得出第一个方向，即第一个主成分，使得在该方向上进行映射可以有最大的间距（方差）
plt.scatter(X_mean[:,0], X_mean[:,1])
plt.plot([0, w[0]*30], [0, w[1]*30], color='r') # 将w表示的方向进行绘制，为了更清楚地展示，这里适当扩大w的表示范围（这里扩大了30倍），但是w0和w1比例要保持一致。
plt.show()
对样本X降维效果图.png
这里求的是第一主成分，在更高维度的空间，可能还需要继续求出其他的主成分，即降更多的维度。
求前n主成分
对于只有2维特征的样本X，我们求解出了第一主成分，即w=（w1,w2），为了求解第二主成分，这里需要将样本X在第一主成分上的分量去除掉，这里使用的方法即空间几何的向量减法，得到的结果即下图中的绿线部分。
在得到绿线部分的分量后，再对该分量重新求第一主成分，以此类推，可以得到前n个主成分。
求前n主成分代码实现
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 direction(w):
return w / np.linalg.norm(w)
def first_component(X, initial_w, eta, n_iters = 1e4, epsilon=1e-8):
w = direction(initial_w)
cur_iter = 0
while cur_iter < n_iters:
gradient = df(w, X)
last_w = w
w = w + eta * gradient
w = direction(w)
if(abs(f(w, X) - f(last_w, X)) < epsilon):
break
cur_iter += 1
return w
求解第2个主成分，首先求去除掉第一个主成分后样本的值
X = demean(X)
initial_w = np.random.random(X.shape[1])
eta = 0.01
w = first_component(X, initial_w, eta) #这里求解出第一个主成分
X2 = np.empty(X.shape) # 初始1个矩阵，和X的维度保持一致，用来存储去除掉第一个主成分分量后的样本
for i in range(len(X)):
X2[i] = X[i] - X[i].dot(w) * w # 得到每个样本在去除掉第一个主成分分量后的值
#将该分量绘制出来
plt.scatter(X2[:,0], X2[:,1])
plt.show()
对剩下的分量求第一主成分
w2 = first_component(X2, initial_w, eta)
# 第一主成分和第二主成分结果相乘接近0，因为第一主成分和第二主成分垂直，所以矩阵相乘为0
w.dot(w2)
代码封装
#Author:Sunshine丶天
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 fit(self, X, eta=0.01, n_iters=1e4):
"""获得数据集X的前n个主成分"""
assert self.n_components <= X.shape[1], \
"n_components must not 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 direction(w):
return w / np.linalg.norm(w)
def first_component(X, initial_w, eta=0.01, n_iters=1e4, epsilon=1e-8):
w = direction(initial_w)
cur_iter = 0
while cur_iter < n_iters:
gradient = df(w, X)
last_w = w
w = w + eta * gradient
w = direction(w)
if (abs(f(w, X) - f(last_w, X)) < epsilon):
break
cur_iter += 1
return w
X_pca = demean(X)
self.components_ = np.empty(shape=(self.n_components, X.shape[1]))
for i in range(self.n_components):
initial_w = np.random.random(X_pca.shape[1])
w = first_component(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_)
def __repr__(self):
return "PCA(n_components=%d)" % self.n_components
sklearn中的PCA
#Author:Sunshine丶天
from sklearn.decomposition import PCA
from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.neighbors import KNeighborsClassifier
digits = datasets.load_digits() # 加载波士顿房价数据
x = digits.data
y = digits.target
# 切分训练集和测试集
x_train,x_test,y_train,y_test = train_test_split(x,y,random_state = 666)
# 使用knn进行预测。
knn = KNeighborsClassifier()
knn.fit(x_train,y_train)
# 预测准确率
print(knn.score(x_test,y_test))
print('未降维之前的维数：', x.shape[1]) # 结果为64
pca = PCA(0.95) # 传入的参数为0.95，代表保留的95%的信息
pca.fit(x_train)
print(pca.n_components_) # 降维之后结果为28，由原来的64位降到28纬，保留95%的信息
x_train_reduction = pca.transform(x_train) #训练数据集降维
x_test_reduction = pca .transform(x_test) #测试数据集也必须降维
knn = KNeighborsClassifier()
knn.fit(x_train_reduction,y_train)
print(knn.score(x_test_reduction,y_test))
