（六）PCA与梯度上升法

一，什么是PCA（Principal Component Analysis)

• 一个非监督的机器学习方法

• 主要用于数据的降维

• 通过降维，可以发现更便于人类理解的特征

• 其他应用：可视化；去噪

理解：PCA的主要思想是将n为特征映射到k维上，这k维是全新的正交特征也被称为主成分，是在原有n维特征的基础上重新构造出来的k维特征。PCA的工作就是从原始的空间中顺序地找一组相互正交的坐标轴，新的坐标轴的选择与数据本身是密切相关的。其中，第一个新坐标轴选择是原始数据中方差最大的方向，第二个新坐标轴选取是与第一个坐标轴正交的平面中使得方差最大的，第三个轴是与第1,2个轴正交的平面中方差最大的。依次类推，可以得到n个这样的坐标轴。通过这种方式获得的新的坐标轴，我们发现，大部分方差都包含在前面k个坐标轴中，后面的坐标轴所含的方差几乎为0。于是，我们可以忽略余下的坐标轴，只保留前面k个含有绝大部分方差的坐标轴。事实上，这相当于只保留包含绝大部分方差的维度特征，而忽略包含方差几乎为0的特征维度，实现对数据特征的降维处理。

当样本只有俩个特征时，并且这俩个特征的分布满足某种关系，比如线性关系，那我们就可以把这写样本点近似的映射到一条直线上，也就是把特征1和特征2看作一条直线，就做到把二维的点降维到了一维。并且此时点与点之间的间距更大，更容易区分。

因此问题简化为：

如何找到让这个样本间间距最大的轴（即映射的直线）？

如何定义样本间距？

使用方差Variance，表示样本的疏密！方差最大化的意义是尽最大限度的保留数据样本之间的关系特征。

=》找到一个轴，使得样本空间的所有点映射到这个轴后方差最大

首先要使样本的均值为0，即对所有样本进行demean处理

之后我们想要求一个轴的方向w=(w1,w2,…)，使得所有的样本映射到w以后，

对于二维数据，有

其中X的平均是0，X的投影是俩个向量的点乘

对于多维样本：

一个目标函数的最优化问题，使用梯度上升法解决

二，使用梯度上升法求解PCA问题

注：这样的符号中，右上角的(1)表示行，也就是第几个样本，右下角的1表示样本的特征，也就是列。

虽然上述的推导已经可以交给计算机实现了，但是我们为了简化程序的逻辑，可以进一步转化为：

想清楚那些符号表示矩阵，哪些符号表示向量（小写），以及向量的方向，维度等

三，求数据的主成分PCA

数据生成：

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)

demean数据：

def demean(X):

return X - np.mean(X, axis=0) #减去的是每一个特征的均值，axis=0是列

X_demean = demean(X)

使用梯度上升法：

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)

defdf_debug(w, X, epsilon=0.0001): #验证梯度的准确性

res = np.empty(len(w)) #epsilon取值较小是因为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 direction(w):

return w / np.linalg.norm(w) #norm就是求模

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

initial_w =np.random.random(X.shape[1]) # 注意2：不能用0向量开始

eta =0.001 # 注意3：不能使用StandardScaler标准化数据

因为对于用户给定的数据，如果俩种数据之间有线性关系，当我们归一化数据，也就是把数据的方差定为1的话就破坏了这种线性关系，使得主成分的轴发生改变。

而demean的过程是将坐标轴进行位移，并没有改变数据的分布。

gradient_ascent(df_debug,X_demean, initial_w, eta)

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')

这个轴是计算出的第一个主成分，因此也被称为第一主成分

四，求数据的前n个主成分

第一个主成分的方差最大，第二个主成分次之，以此类推

求出第一主成分后，如何求出下一个主成分？

数据进行改变，将数据在第一个主成分上的分量去掉

映射到绿线，绿线与蓝线垂直！认真理解

当样本中有俩个主成分，当我们去除第一个主成分的分量后，剩余样本均分布在第二主成分上，即剩余样本的特征的方差是0

程序实现：

deffirst_n_components(n, X, eta=0.01, 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) * w

return res

五，高维数据映射为低维数据

Wk是新建立的k维坐标系，Xk是X在这个坐标系下的映射

当我们想要恢复原来的维度时，其实是不能完全复原的，因为在映射到主成分的过程中，我们丢失了一部分的信息。

scikit-learn中的PCA

scikit-learn中的PCA使用的方法不是梯度上升，是另外的数学求解方法

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[:,0],X[:,1], color='b', alpha=0.5)

plt.scatter(X_restore[:,0],X_restore[:,1], color='r', alpha=0.5)

plt.show()

PCA的真正威力：

import numpy as np

importmatplotlib.pyplot as plt

from sklearn importdatasets

digits =datasets.load_digits()

X = digits.data

y =digits.target #创建数据

fromsklearn.model_selection import train_test_split

X_train,X_test,y_train,y_test=train_test_split(X,y,random_state=666)

%%time

#使用k邻近计算预测准确率

from sklearn.neighborsimport KNeighborsClassifier

knn_clf=KNeighborsClassifier()

knn_clf.fit(X_train,y_train)

使用PCA对数据降维：

fromsklearn.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)

#必须用训练X的方法来训练测试数据集

%%time

knn_clf=KNeighborsClassifier()

knn_clf.fit(X_train_reduction,y_train)

knn_clf.score(X_test_reduction,y_test)

识别精度太低！因此降低的维数太小，需要调整参数

PCA提供了一个算法：

pca.explained_variance_ratio_#主成分可解释的方差相应的比例

表示第一个轴可以解释0.14比例原数据的方差，第二个轴同理

PCA的手段是维持的方差越大越好，当这个数太小就表示原始数据的很多信息都丢失了，我们维持方差最大是为了尽可能的保留原始数据间的信息。

可以看到一共64个轴，每个轴对应的重要成分，值越小表示越不重要

plt.plot([i for i inrange(X_train.shape[1])],

[np.sum(pca.explained_variance_ratio_[:i]) for i inrange(X_train.shape[1])])

plt.show()

如果我们希望数据可以保持95%以上的信息，就可以在上图中找到对应主成分数目。这个函数在sklearn中也已经封装好

pca=PCA(0.95) #确定精度范围

pca.fit(X_train)

之后再用KNN训练，计算时间明显减少，并且精度也保持较高水平，只是依然没有使用全样本数据得到的精度高

使用PCA对数据进行降维可视化

pca =PCA(n_components=2) #将数据降到2维

pca.fit(X)

X_reduction =pca.transform(X)

for i in range(10):

plt.scatter(X_reduction[y==i,0],X_reduction[y==i,1], alpha=0.8)

plt.show()

从上图可以看出，有些特征即使降到二维也可以有很高的区分度，比如蓝色和紫色

使用MNIST数据集

import numpy as np

from sklearn.datasetsimport fetch_openml

mnist =fetch_openml('mnist_784')

X,y=mnist['data'],mnist['target']

使用KNN

from sklearn.neighborsimport KNeighborsClassifier

knn_clf=KNeighborsClassifier()

%timeknn_clf.fit(X_train,y_train)

%timeknn_clf.score(X_test,y_test)

这个样本的数据在使用之前并没有对样本进行归一化处理，是因为样本不同维度的数据都表示图像中像素点的灰度值，因此在同一个数量级上

使用PCA进行降维

fromsklearn.decomposition import PCA

pca = PCA(0.90)

pca.fit(X_train)

X_train_reduction =pca.transform(X_train)

X_test_reduction =pca.transform(X_test)

X_train_reduction.shape

从数据可以看出，经过降维后，数据从原来的784降到了87，但仍然保留了原始数据90%的信息

knn_clf =KNeighborsClassifier()

%timeknn_clf.fit(X_train_reduction, y_train)

%timeknn_clf.score(X_test_reduction, y_test)

可以发现，对样本进行降维之后，预测的准确率反而提高了，是因为原始数据中有噪声，通过PCA的方式可以达到降噪的作用

六，使用PCA对数据进行降噪

手写识别的例子

from sklearn importdatasets

digits =datasets.load_digits()

X = digits.data

y = digits.target

noisy_digits= X + np.random.normal(0, 4, size=X.shape)

#没有噪音的数据集加噪声矩阵

example_digits =noisy_digits[y==0,:][:10]

for num in range(1,10):

example_digits = np.vstack([example_digits,noisy_digits[y==num,:][:10]])

defplot_digits(data): #绘制数据图像

fig, axes = plt.subplots(10, 10,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(data[i].reshape(8, 8),

cmap='binary',interpolation='nearest',

clim=(0, 16))

plt.show()

plot_digits(example_digits)

可以看出噪声相对较大

components =pca.transform(example_digits)

filtered_digits= pca.inverse_transform(components) #过滤后的数据

plot_digits(filtered_digits)

处理后的图像更清楚了些

七，人脸识别与特征脸

通过PCA计算出来的前N个主成分就可以称为特征脸，也就是说，形成的k维正交向量组，第一主成分也是这个正交基的第一个轴，这个轴最能反应原始数据的信息，因为这条轴保留了原始数据最多的信息。

绘制一些人脸图像：

import numpy as np

importmatplotlib.pyplot as plt

from sklearn.datasetsimport fetch_lfw_people

faces =fetch_lfw_people()

图像是62乘以47的图像

random_indexes =np.random.permutation(len(faces.data))

#先对样本进行随机排列

X =faces.data[random_indexes]

example_faces =X[:36,:] #要绘制的人脸

example_faces.shape

defplot_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

fromsklearn.decomposition import PCA

pca =PCA(svd_solver='randomized') #使用随机的方式求解PCA

因为数据样本比较多，使用随机的方式可以加快计算速度

pca.fit(X)

pca.components_.shape

plot_faces(pca.components_[:36,:])

绘制出经过PCA降维后的人脸特征，当主成分越多，对应下图越靠后的位置，图像的人脸轮廓越明显