机器学习 -- 主成分分析

主成分分析

Principal Component Analysis PCA

一个非监督的机器学习算法
主要用于数据的降维
通过降维，可以发现更便于人类理解的特征
其他应用，可视化，去噪

原理：

进行降维，保留特征1

进行降维，保留特征2
上面哪种方案更好？可以看出保留特征1 的点与点之间的间距更大，拥有更好的可区分度。这种方案比较好
还可以有更好的方案吗？

将点映射到这条直线上

这种方案是不是和原来的更接近，也更有区分度
所以，变成：如何找到这个让样本间间距最大的轴？
如何定义样本间间距？使用方差（Variance）
Var(x)=1m∑i=1m(xi−x‾)2Var(x) = \frac{1}{m} \sum_{i=1}^m(x_i -\overline{x})^2Var(x)=m1∑i=1m(xi−x)2
找到一个轴，使得样本空间的所有点映射到这个轴后，方差最大

第一步：将样本的均值归为0 （demean处理）

样本就变成如下图：

将所有样本进行demean处理
我们想要求一个轴的方向 w = (w1,w2),
使得我们所有的样本，映射到w以后，有：
Var(Xprroject)=1m∑i=1m(Xprojecti−X‾project)2Var(X_{prroject})=\frac{1}{m} \sum_{i=1}^m(X^i_{project}- \overline{X}_{project})^2Var(Xprroject)=m1∑i=1m(Xprojecti−Xproject)2 最大
因为 X‾project\overline{X}_{project}Xproject = 0
所以： Var(Xprroject)=1m∑i=1m∣∣Xprojecti∣∣2Var(X_{prroject})=\frac{1}{m} \sum_{i=1}^m ||X^i_{project}||^2Var(Xprroject)=m1∑i=1m∣∣Xprojecti∣∣2

将一个向量映射到另一个向量上，就是两个向量 ⋅\cdot⋅ 点乘的定义，就是两个向量的相乘然后乘以两个向量之间的夹角
Xi⋅w=∣∣Xi∣∣⋅∣∣w∣∣⋅cosθX^i \cdot w = ||X^i|| \cdot ||w|| \cdot cos\thetaXi⋅w=∣∣Xi∣∣⋅∣∣w∣∣⋅cosθ
www可以用方向向量来表示就可以，所以可以使其为 www=1
Xi⋅w=∣∣Xi∣∣⋅cosθX^i \cdot w = ||X^i|| \cdot cos\thetaXi⋅w=∣∣Xi∣∣⋅cosθ
↓\downarrow↓
Xi⋅w=∣∣Xprojecti∣∣X^i \cdot w = ||X^i_{project}||Xi⋅w=∣∣Xprojecti∣∣
所以 Var(Xprroject)=1m∑i=1m∣∣Xi⋅w∣∣2Var(X_{prroject})=\frac{1}{m} \sum_{i=1}^m ||X^i \cdot w||^2Var(Xprroject)=m1∑i=1m∣∣Xi⋅w∣∣2 ,
↓\downarrow↓
目标：求 www，使得 Var(Xprroject)=1m∑i=1m(Xi⋅w)2Var(X_{prroject})=\frac{1}{m} \sum_{i=1}^m (X^i \cdot w)^2Var(Xprroject)=m1∑i=1m(Xi⋅w)2 最大

Var(Xprroject)=1m∑i=1m(X1iw1+X2iw2+...+Xniwn)2Var(X_{prroject})=\frac{1}{m} \sum_{i=1}^m (X^i_1w_1+X^i_2w_2+...+X^i_nw_n)^2Var(Xprroject)=m1∑i=1m(X1iw1+X2iw2+...+Xniwn)2
Var(Xprroject)=1m∑i=1m(∑j=1nXjiwj)2Var(X_{prroject})=\frac{1}{m} \sum_{i=1}^m (\sum_{j=1}^n X^i_jw_j)^2Var(Xprroject)=m1∑i=1m(∑j=1nXjiwj)2
一个目标函数的最优化问题，使用梯度上升法解决

主成分分析和线性回归的区别

这是主成分分析，是多个特征之间的，是垂直于这个轴的
线性回归是一个标记结果和特征之间的，是垂直于特征的，

使用梯度上升法来解决PCA问题

目标：求w,使得 f(X)=1m∑i=1m(X1iw1+X2iw2+...+Xniwn)2f(X) =\frac{1}{m} \sum_{i=1}^m (X^i_1w_1+X^i_2w_2+...+X^i_nw_n)^2f(X)=m1∑i=1m(X1iw1+X2iw2+...+Xniwn)2 最大

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

import numpy as np
import matplotlib.pyplot as pltX = 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)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])
plt.show()def f(w,f):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的模等于1 ，所以此处需要将 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) # 因为在推到公式时候需要 w的模等于1 ，所以此处需要将 w进行处理cur_iter=0while cur_iter < n_iters:gradient = df(w,X)last_w = ww = w + eta * gradient w = direction(w)  # 使用pca 求解梯度法需要注意的，需要将w转为单位向量。if (abs(f(w,X) - f(last_w,X)) < epsilon):breakcur_iter +=1return winitial_w = np.random.random(X.shape[1])   # 初始向量不能为0， 不能从0 开始eta = 0.001
# 注意不能使用 standardScaler 标准化数据，因为 standardScaler 本来就是将数据标准化，使方差为1，
# 而 pca 求的就是将方差最大值。 demeanw = 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')
plt.show()

注意事项：

先对数据进行 demean 处理
不要对数据进行 standardScaler 标准化处理，因为pca求得就是最大方差，standardScaler 就是使得方差变成1
在给 w初始化值的时候，不能为0向量，可以进行非0随机化处理
在梯度上升计算中，要是w变成w的单位向量

求前n个主成分

上面代码将两个特征轴变成了一个特征轴，
求出第一主成分后，如何求出下一个主成分呢？
将数据进行改变，将数据在第一个主成分上的分量去掉

上图中的 X‘iX^{`i}X‘i就是新的数据。
在新的数据上求第一主成分

# 根据上面的运算 继续
X2 = np.empty(X.shape)
for i in range(len(X)):X2[i] = X[i] - X[i].dot(w)*w  # 将样本在第一主成分上相应的分量去掉的结果# 上面公式 可以进行向量化
X2= X - X.dot(w).reshape(-1,1) * wplt.scatter(X2[:,0],X2[:,1])
plt.show()w2 = gradient_ascent(df, X2, initial_w, eta)
w2w.dot(w2)

前n个主成分计算

def first_component(X, initial_w, eta, n_iters=1e4, epsilon=1e-8):w = direction(initial_w) # 因为在推到公式时候需要 w的模等于1 ，所以此处需要将 w进行处理cur_iter=0while cur_iter < n_iters:gradient = df(w,X)last_w = ww = w + eta * gradient w = direction(w)  # 使用pca 求解梯度法需要注意的，需要将w转为单位向量。if (abs(f(w,X) - f(last_w,X)) < epsilon):breakcur_iter +=1return wdef first_n_components(n, X,n_iters=1e4, epsilon=1e-8):X_pca = X.copy()X_pca = demean(X_pca)res = []for i in range(n):initial_w = np.random.random(X_pca.shape[1])w = first_component(X_pca, initial_w, eta)res.append(w)X_pca = X_pca - X_pca.dot(w).reshape(-1,1) * wreturn res

对数据进行降维

高维数据向低维数据映射

降维：

低维映射会高维但是生成的数据和原来的高维数据不一样，在降维过程中会丢失一些数据，所以无法恢复回来

import numpy as npclass PCA:def __init__(self, n_components):assert n_components >=1, 'n_components must be valid'self.n_components = n_componentsself.components = Nonedef fit(self, X, eta=0.01, n_iters=1e4):assert self.n_components <= X.shape[1], 'n_components must not be greater than the feature number of X'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 demean(X):return X- np.mean(X,axis=0)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= 0while cur_iter < n_iters:gradient = df(w,X)last_w = ww = w + eta * gradientw = direction(w)if abs(f(w,X) - f(last_w, X)) < epsilon:breakcur_iter +=1return wX_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,:] =wX_pca = X_pca - X_pca.dot(w).reshape(-1,1)*wreturn self# 降维  会丢失一些信息def transform(self,X):assert X.shape[1] == self.components_.shape[1]return X.dot(self.components_.T)# 升维  只是在高维空间里表达低维的信息def inverse_transform(self,X):assert X.shape[1] == self.components_.shape[0]return X.dot(self.components_)def __repr__(self):return f'PCA(n_components={self.n_components})'

scikit-learn 中 PCA

from sklearn.decomposition import PCA
pca =PCA(n_components=1)
pca.fit(X)
pca.components_
X_reduction = pca.transform(X)
X_restore  = pca.inverse_transform(X_reduction)

使用digits数据

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasetsdigits = datasets.load_digits()
X = digits.data
y = digits.targetfrom sklearn.model_selection import train_test_split
X_train, X_test, y_train,y_test = train_test_split(X,y,random_state=666)%%time
from sklearn.neighbors import KNeighborsClassifier
knn_clf = KNeighborsClassifier()
knn_clf.fit(X_train, y_train)  # Wall time: 2 ms
# KNeighborsClassifier(algorithm='auto', leaf_size=30, metric='minkowski',metric_params=None, n_jobs=None, n_neighbors=5, p=2,weights='uniform')knn_clf.score(X_test, y_test)  # 0.9866666666666667from sklearn.decomposition import PCA
pca = PCA(n_components=2)
pca.fit(X_train)
X_train_reduction = pca.transform(X_train)  # 多数据进行降维处理
X_test_reduction = pca.transform(X_test)
%%time
knn_clf = KNeighborsClassifier()  # 再使用knn 处理
knn_clf.fit(X_train_reduction, y_train)# Wall time: 1 ms
# KNeighborsClassifier(algorithm='auto', leaf_size=30, metric='minkowski',metric_params=None, n_jobs=None, n_neighbors=5, p=2,weights='uniform')knn_clf.score(X_test_reduction,y_test) # 0.6066666666666667# 上面因为将64维降成2维，很多数据丢失，所以准确率低
# pca 为我们提供了一个很好的接口， pca = PCA(0.95)  # 这里的0.95表示我们如果想要将64维数据压缩成保留全部数据的95%的数据，pca内部会为我们计算出合适的 维度
pca.fit(X_train)
pca.n_components_   # 28  pca 计算出 保留95%数据时，应该选取维度是28
# 此时转换 数据
X_train_reduction = pca.transform(X_train)
X_test_reduction = pca.transform(X_test)
# 将数据进行knn操作%%time  # 2ms
knn_clf = KNeighborsClassifier()
knn_clf.fit(X_train_reduction, y_train)knn_clf.score(X_test_reduction,y_test)  # 0.98   维度增加，时间提升，但是准确率也提升了

有时数据由很多维降到2维也是很有效的

pca =PCA(n_components=2)
pca.fit(X)
X_reduction = pca.transform(X)
X_reduction.shapefor i in range(10):plt.scatter(X_reduction[y==i,0], X_reduction[y==i,1], alpha=0.8)
plt.show()

可以看出启动蓝色和紫色的数据区分很明显。不需要高维数据也可以区分。

我们可以通过 pca.explained_variance_ratio_ 来查看多维数据每一维所占的总数据比例：

pca = PCA(n_components=64)
pca.fit(X_train)
pca.explained_variance_ratio_
#array([1.45668166e-01, 1.37354688e-01, 1.17777287e-01, 8.49968861e-02,5.86018996e-02, 5.11542945e-02, 4.26605279e-02, 3.60119663e-02,3.41105814e-02, 3.05407804e-02, 2.42337671e-02, 2.28700570e-02,1.80304649e-02, 1.79346003e-02, 1.45798298e-02, 1.42044841e-02,1.29961033e-02, 1.26617002e-02, 1.01728635e-02, 9.09314698e-03,8.85220461e-03, 7.73828332e-03, 7.60516219e-03, 7.11864860e-03,6.85977267e-03, 5.76411920e-03, 5.71688020e-03, 5.08255707e-03,4.89020776e-03, 4.34888085e-03, 3.72917505e-03, 3.57755036e-03,3.26989470e-03, 3.14917937e-03, 3.09269839e-03, 2.87619649e-03,2.50362666e-03, 2.25417403e-03, 2.20030857e-03, 1.98028746e-03,1.88195578e-03, 1.52769283e-03, 1.42823692e-03, 1.38003340e-03,1.17572392e-03, 1.07377463e-03, 9.55152460e-04, 9.00017642e-04,5.79162563e-04, 3.82793717e-04, 2.38328586e-04, 8.40132221e-05,5.60545588e-05, 5.48538930e-05, 1.08077650e-05, 4.01354717e-06,1.23186515e-06, 1.05783059e-06, 6.06659094e-07, 5.86686040e-07,7.44075955e-34, 7.44075955e-34, 7.44075955e-34, 7.15189459e-34])# 画图看下各个维度的数据占比
plt.plot([i for i in range(X_train.shape[1])], [np.sum(pca.explained_variance_ratio_[:i+1]) for i in range(X_train.shape[1])])
plt.show()

MNIST 数据集处理

降维能很大的节省处理时间

import numpy as np
from sklearn.datasets  import fetch_mldatamnist = fetch_mldata("MNIST original")
X,y = mnist['data'], mnist['target']
## 这个数据中已经处理好。前60000个是训练数据，后面10000个是测试数据
X_train = np.array(X[:60000],dtype=float)
y_train = np.array(y[:60000],dtype=float)
X_test = np.array(X[60000:],dtype=float)
y_test = np.array(y[60000:],dtype=float)# 直接使用 knn
from sklearn.neighbors import KNeighborsClassifier
knn_clf = KNeighborsClassifier()
%time knn_clf.fit(X_train, y_train)  # 13s
%time knn_clf.score(X_test, y_test) # 8min 59s   0.9688# 使用pca降维， 然后使用knn
# 时间明显降低
from sklearn.decomposition import PCA
pca = PCA(0.9)
pca.fit(X_train)
X_trian_reduction = pca.transform(X_train)knn_clf = KNeighborsClassifier()
%time knn_clf.fit(X_trian_reduction, y_train)  # 291 ms%time X_test_reduction = pca.transform(X_test) #68.1 ms%time knn_clf.score(X_test_reduction, y_test)  # 1min 4s  0.9728

pca 对数据的降噪处理

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,5,size=100)
plt.scatter(X[:,0],X[:,1])
plt.show()

from sklearn.decomposition import PCA
pca =PCA(n_components=1)
pca.fit(X)
X_reduction =pca.transform(X)  # 先降维
X_restore = pca.inverse_transform(X_reduction) # 再还原  plt.scatter(X_restore[:,0],X_restore[:,1])  # 此时可以看到变成一条直线，进行了降噪
plt.show()

人脸识别和特征脸

import numpy as np
import matplotlib.pyplot as plt# 获取数据
from sklearn.datasets import fetch_lfw_people
faces = fetch_lfw_people()faces.keys()
faces.data.shapefaces.images.shaperandom_indexes = np.random.permutation(len(faces.data))
X = faces.data[random_indexes]example_faces =X[:36,:]  # 获取36个图片数据
example_faces.shapedef plot_faces(faces):fig, axes = plt.subplots(6,6,figsize=(10,10),subplot_kw={'xticks':[],'yticks':[]},gridspec_kw=dict(hspace=0.1,wspace=0.1))for i ,ax in enumerate(axes.flat):ax.imshow(faces[i].reshape(62,47),cmap='bone')plt.show()
plot_faces(example_faces)  # 显示图片数据%%time   # 查看数据的特征的重要性
from sklearn.decomposition import PCA
pca = PCA(svd_solver='randomized')
pca.fit(X)# 查看特征 数据 ，这些特征就是特征脸
pca.components_.shape
plot_faces(pca.components_[:36,:])