本文深入介绍了主成分分析PCA的概念，通过实例展示了PCA如何降维并保留样本间最大方差。PCA通过寻找样本数据的最大方差轴，实现数据的线性变换，降低数据复杂性。文中还探讨了PCA与线性回归的区别，并提供了使用梯度上升法实现PCA的Python代码。此外，还介绍了PCA在高维数据降维和可视化中的应用，以及在手写数字识别任务中的降噪效果。
PCA简介
PCA（Principal Component Analysis）：也是一个梯度分析的应用，不仅是机器学习的算法，也是统计学的经典算法
1.1 举个栗子
例如下面一个两个特征的一个训练集，我们可以选择一个特征，扔掉一个特征
下图分别是扔掉了特征一和特征二的两种方案，很明显右边这种的效果会更好一些，因为访问二扔掉特征二以后，点之间的分布情况更接近与原图，但是这不是更好的
我们希望有一根直线，是斜着的，我们希望将所有的点都映射到这条直线上，那么这个时候我们就成功的将二维降到了一维，与此同时，这些点更加趋近与原来的点的分布情况，换句话说，点和点之间的距离比无论是映射到x还是映射到y周，他们之间的区分度都更加的大，也就更加容易区分
1.2 总结
那么如何找到这个让样本间间距最大的轴? 如何定义样本间间距? 事实上有一个指标可以之间定义样本间的距离，就是方差（Variance）（方差：描述样本整体之间的疏密的一个指标，方差越大，代表样本之间越稀疏，方差越小，代表样本之间越紧密）
方差
第一步: 将样例的均值归为0 (demean)（归0：所有样本都减去他们的均值），使得均值为0，这样可以简化方差的公式
1.3 推导
进行均值归0操作以后，就是下面的式子
注：|Xproject|的平均值也是一个向量
X(i)映射到w的距离实际上就是X(i)与w的点乘（蓝色的线），根据定义推导，其值实际上就是Xproject
此时我们的目标函数就可以化简成
这是一个目标函数的最优化问题，使用梯度上升法解决。 当然我们也可以之间使用数学原理推导出结果，这里我们主要关注使用搜索的策略来求解主成分分析法，这样我们对梯度上升发和梯度下降法也可以有一个更深刻的认识
1.4 与线性回归的区别
1.主成分分析法的两个轴都是特征，线性回归y轴是目标结果值 2.主成分分析法的点是垂直于方差轴直线的，线性回归的点事垂直于x轴的
使用梯度上升法解决PCA问题
- 1.注意上面式子里的每一个**(X1(i)·w1+X2(i)·w2+…Xn(i)·wn)**都是一个X(i)和w的点乘，所以式子可以进一步化解，
- 2.化简过后可以进行向量化，即每一个∑(X(i)·w1)·X1(i) 可以看成是**(X·w)这个向量的转置（本来是个行向量，转置后是1行m列的列向量）与X**这个矩阵（m行n列）做点乘等到的其中一项的相乘相加的结果
- 3.最后根据转置法则 ((AB)T=BTAT)转换成最后的结果
使用梯度上升法实现PCA
import numpy as np
import matplotlib.pyplot as plt
X = np.empty((100, 2))
X[:,0] = np.random.uniform(0., 100., size=100) # 0~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()
np.mean(X_demean[:,0])
1.3500311979441904e-15
np.mean(X_demean[:,1])
-7.6738615462090825e-15
梯度上升法
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 direction(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
初始化w，相当于梯度下降时候设置theta的值，但是初始位置不能是0，因为将0带入df_math结果还是零， X.shape[1]样本的个数， 列数
initial_w = np.random.random(X.shape[1]) ## 注意2：不能用0向量开始
initial_w
array([ 0.37061708, 0.28515471])
eta = 0.001
## 注意3： 不能使用StandardScaler标准化数据， 标准化的方差就变成1了， demean已经将均值归0了， 做了一半的工作
gradient_ascent(df_debug, X_demean, initial_w, eta)
array([ 0.78121351, 0.62426392])
gradient_ascent(df_math, X_demean, initial_w, eta)
array([ 0.78121351, 0.62426392])
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')
plt.show()
使用极端数据集测试
X2 = np.empty((100, 2))
X2[:,0] = np.random.uniform(0., 100., size=100)
X2[:,1] = 0.75 * X2[:,0] + 3.
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_demean[:,0], X2_demean[:,1])
plt.plot([0, w2[0]*30], [0, w2[1]*30], color='r')
plt.show()
求数据的前N个主成分
求出第一主成分以后，如何求出下一个主成分? 1.数据进行改变，将数据在第一个主成分上的分量去掉
X(i)·w = ||Xproject(i)|| 即X(i)映射到w上的值，那么||Xproject(i)||（大小） ·w（方向）就是X(i)在w上的分向量记为Xproject(i)= ||Xproject(i)|| ·w X(i)-Xproject(i)就可以实现将X样本在Xproject相应上的分量去掉，相减之后的集合意义就是讲X样本映射到了Xproject向量相垂直的一个轴上，记为X`(i) = Xproject(i)
2.在新的数据上求第一主成分 得到的X` 是X中的所有样本都去除了第一主成分上的分量得到的结果，要求第二主成分，只要在新的数据上，重新求一下第一主成分
获得前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 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
initial_w = np.random.random(X.shape[1])
eta = 0.01
w = first_component(X, initial_w, eta)
w
array([ 0.78307181, 0.62193129])
X2 = np.empty(X.shape)
for i in range(len(X)):
X2[i] = X[i] - X[i].dot(w) * w
plt.scatter(X2[:,0], X2[:,1])
plt.show()
## 向量化，X.dot(w)为m*1的向量，reshape后变成了1*m的列向量，
# 再乘以w（方向）就是X的每一个值在w上 的分量矩阵
X2 = X - X.dot(w).reshape(-1, 1) * w
## 相减得到的样本分布几乎垂直于原来的样本分布
plt.scatter(X2[:,0], X2[:,1])
plt.show()
w2 = first_component(X2, initial_w, eta)
w2
array([-0.6219279 , 0.78307451])
## 因为w和w2都是单位向量，所以他们两个点乘得到的结果就是他们夹角的cos值，
## 又因为w和w2应该是互相垂直的，所以他们夹角的cos值等于0
w.dot(w2)
4.3300901158005445e-06
增加一个参数n，求出样本X的前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 = []
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
first_n_components(2, X)
[array([ 0.78307192, 0.62193116]), array([-0.62192683, 0.78307536])]
高维数据向低维数据进行映射
对于一个数据集X来说，这个X有m行n列，代表有m个样本n个特征，通过我们前面学习的主成分分析法，假设我们已经求出了针对这个数据来说的前k个主成分，每一个主成分对应一个单位方向，用W来表示,W也是一个矩阵，他有k行，代表我们求出的前K个主成分，每一行有n列，代表每一个主成分的坐标轴应该是有n个元素的。这是因为我们的主成分分析法主要就是将数据从一个坐标系转化成了另外一个坐标系，原来这个坐标系有n个维度，现在这个坐标系也应该有n个维度，只不过对于转化的坐标系来说，我们取出来前k个，这k个方向更加重要。
如何将我们的样本X从n维转化成k维呢，回忆们之前学到的，对于一个X样本，与一个W进行点乘，其实就是讲一个样本映射到了w这个坐标轴，得到的模，如果讲这一个样本和这k个w分别做点乘，得到的就是这一个样本，在这k个方向上做映射后每一个方向上的大小，这k个元素合在一起，就代表这一个样本映射到新的k个轴所代表的坐标系上相应的这个样本的大小
X1分别乘以W1到Wn，得到的k个数组成的向量，就是样本1映射到Wk这个坐标系上得到的k维的向量，由于k<n，所以我们就完成了一个样本从n维到k维的映射，这个过程依次类推从样本1到样本m都这么做，我们就将m个样本都从N维映射到了k维-----其实我们就是做了一个乘法X·WT(为什么是转置呢，因为我们是拿X的每一行去和W的每一行做点乘的，但是矩阵乘法规定是拿X的每一行和W的每一列做乘法)
我们得到新的降维后的矩阵Xk以后，是可以通过和Wk想乘回复回来的，但是由于我们在降维的过程中丢失了一部分信息，这时及时回复回来也和原来的矩阵不一样了，但是这个从数据角度成立的
从高维数据向低维数据的映射
class PCA:
def __init__(self,n_components):
"""初始化PCA"""
assert n_components>=1, "n_components must be vaild"
# 前面一直说的k，想要多少个主成分
self.n_components = n_components
# 每个主成分是谁，就是上面式子的Wk
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 be greater then 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 direction(w):
"""计算单位向量"""
return w / np.linalg.norm(w)
def first_componet( X,inital_w,eta,n_iters = 1e4,epsilon=1e-8):
w = direction(inital_w)
cur_iter = 0
while cur_iter < n_iters:
gradient = df(w,X)
last_w = w
w = w + eta * gradient
## 注意1：每次求单位向量
w = direction(w)
if abs(f(w,X)-f(last_w,X)) < epsilon:
break
cur_iter = 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_componet(X_pca,initial_w,eta)
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
PCA 降维的基本原理:找到另外一个坐标系，这个坐标系每一个轴依次可以表达原来的样本他们的重要程度，也就是称为所有的主成分，我们取得前k个最重要的主成分，就可以将所有的样本映射到这k个轴上，获得一个低维的数据信息
测试封装的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, 10., size=100)
from playML.PCA import PCA
pca = PCA(n_components=2)
pca.fit(X)
PCA(n_components=2)
pca.components_
array([[ 0.76676948, 0.64192256],
[-0.64191827, 0.76677307]])
pca = PCA(n_components=1)
pca.fit(X)
PCA(n_components=1)
X_reduction = pca.transform(X)
X_reduction.shape
(100, 1)
X_restore = pca.inverse_transform(X_reduction)
X_restore.shape
(100, 2)
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()
scikit-learn中的PCA
from sklearn.decomposition import PCA
pca = PCA(n_components=1)
pca.fit(X)
PCA(copy=True, iterated_power='auto', n_components=1, random_state=None,
svd_solver='auto', tol=0.0, whiten=False)
pca.components_
array([[-0.76676934, -0.64192272]])
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()
digits测试scikit-learn中的PCA
import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets
digits = datasets.load_digits()
X = digits.data
y = digits.target
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=666)
X_train.shape
(1347, 64)
%%time
from sklearn.neighbors import KNeighborsClassifier
knn_clf = KNeighborsClassifier()
knn_clf.fit(X_train, y_train)
CPU times: user 19.9 ms, sys: 7.47 ms, total: 27.4 ms
Wall time: 64.5 ms
knn_clf.score(X_test, y_test)
0.98666666666666669
from 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_clf.fit(X_train_reduction, y_train)
CPU times: user 2.13 ms, sys: 767 µs, total: 2.9 ms
Wall time: 2.93 ms
knn_clf.score(X_test_reduction, y_test)
0.60666666666666669
主成分所解释的方差
explained_variance_ratio_-解释方差的比例
0.14566817 代表第一个主成分可以解释14%的原数据
0.13735469 代表第二个主成分可以解释13%的原数据
两个主成分加起来可以解释百分之27的原数据，而其他的信息丢失了
可以使用explained_variance_ratio_这个参数来查看每个主成分所解释的原数据，来判断要取多少个主成分
pca.explained_variance_ratio_
array([ 0.14566817, 0.13735469])
pca.explained_variance_
array([ 175.77007821, 165.73864334])
from sklearn.decomposition import PCA
pca = PCA(n_components=X_train.shape[1])
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])
绘制曲线观察取前i个主成分的时候，所能解释的原数据比例
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()
sklearn中的PCA算法支持传入一个小于1的数来表示我们希望能解释多少比例的主成分
pca = PCA(0.95)
pca.fit(X_train)
PCA(copy=True, iterated_power='auto', n_components=0.95, random_state=None,
svd_solver='auto', tol=0.0, whiten=False)
pca.n_components_
28
X_train_reduction = pca.transform(X_train)
X_test_reduction = pca.transform(X_test)
%%time
knn_clf = KNeighborsClassifier()
knn_clf.fit(X_train_reduction, y_train)
CPU times: user 4.21 ms, sys: 1.28 ms, total: 5.49 ms
Wall time: 19.7 ms
knn_clf.score(X_test_reduction, y_test)
0.97999999999999998
使用PCA对数据进行降维可视化
PCA降维到两维的意义-可以方便可视化展示，帮助人们理解
下图每个颜色代表一个数字在降维到二维空间中的分布情况
仔细观察后可以发现，很多数字的区分还是比较明细的
比如如果只是区分蓝色的数字和紫色的数字，那么使用二个维度就足够了
pca = PCA(n_components=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.datasets import fetch_mldata
mnist = fetch_mldata('MNIST original')
mnist
{'COL_NAMES': ['label', 'data'],
'DESCR': 'mldata.org dataset: mnist-original',
'data': array([[0, 0, 0, ..., 0, 0, 0],
[0, 0, 0, ..., 0, 0, 0],
[0, 0, 0, ..., 0, 0, 0],
...,
[0, 0, 0, ..., 0, 0, 0],
[0, 0, 0, ..., 0, 0, 0],
[0, 0, 0, ..., 0, 0, 0]], dtype=uint8),
'target': array([ 0., 0., 0., ..., 9., 9., 9.])}
X, y = mnist['data'], mnist['target']
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)
X_train.shape
(60000, 784)
y_train.shape
(60000,)
X_test.shape
(10000, 784)
y_test.shape
(10000,)
使用kNN
sklearn 封装的KNeighborsClassifier，在fit过程中如果数据集较大，会以树结构的过程进行存储，以加快knn的预测过程，但是会导致fit过程变慢
没有进行数据归一化，是因为这里的每个维度都标示的是每个像素点的亮度，他们的尺度是相同的，这个时候比较两个样本之间的距离是有意义的
from sklearn.neighbors import KNeighborsClassifier
knn_clf = KNeighborsClassifier()
%time knn_clf.fit(X_train, y_train)
CPU times: user 57.6 s, sys: 681 ms, total: 58.3 s
Wall time: 59.4 s
KNeighborsClassifier(algorithm='auto', leaf_size=30, metric='minkowski',
metric_params=None, n_jobs=1, n_neighbors=5, p=2,
weights='uniform')
%time knn_clf.score(X_test, y_test)
CPU times: user 14min 20s, sys: 4.3 s, total: 14min 24s
Wall time: 14min 29s
0.96879999999999999
PCA进行降维
from sklearn.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
(60000, 87)
knn_clf = KNeighborsClassifier()
%time knn_clf.fit(X_train_reduction, y_train)
CPU times: user 588 ms, sys: 5.23 ms, total: 593 ms
Wall time: 593 ms
KNeighborsClassifier(algorithm='auto', leaf_size=30, metric='minkowski',
metric_params=None, n_jobs=1, n_neighbors=5, p=2,
weights='uniform')
%time knn_clf.score(X_test_reduction, y_test)
CPU times: user 1min 55s, sys: 346 ms, total: 1min 56s
Wall time: 1min 56s
0.9728
使用PCA进行降维后的数据集进行训练，不光时间变短了，准确度也变高了
这是因为PCA的过程中，不仅仅是进行了降维，还在降维的过程中将数据包含的噪音给消除了
这使得我们可以更加好的，更加准确的拿到我们数据集对应的特征，从而使得准确率大大提高
降维的过程可以理解成是去噪。
PCA降噪
手写识别的例子
from sklearn import datasets
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]])
example_digits.shape
(100, 64)
def plot_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)
pca = PCA(0.5).fit(noisy_digits)
pca.n_components_
12
components = pca.transform(example_digits)
filtered_digits = pca.inverse_transform(components)
plot_digits(filtered_digits)
400
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
