主成分分析Principle Component Analysis
主成分分析（PCA）在统计学领域应用非常广泛，同时也是很重要的非监督机器学习算法，PCA主要用于数据的降维。在机器学习中，降维是很重要的一步预处理操作，通过降维，可以发现便于人类理解的特征和提取数据集主要特征。这样一来可以减少要处理的数据量，同时又不破坏数据整体特征，提高了算法的效率。PCA广泛应用于可视化和去噪过程。
以一个简单的数据降维例子介绍PCA：
数据特征
如图所示是某个数据集的两个特征，现在思考一个问题：能不能将两个特征压缩为一个特征呢？一个简单的处理方法是只取特征1或者特征2，这样就达到了降维的目的：
特征降维.png
而且显然取特征1会比特征2有更好的区分度（样本间距更大）。不过，还有没有更好的降维办法呢？考虑这样一条直线：
better_降维
此时将特征投影到红色直线上进行降维，显然和原始的特征更加接近，同时也更符合特征原始分布。那么如何找到这样一条让样本间距最大的轴呢？
首先给出样本间距定义，考虑到统计学中方差Variance是衡量样本分布离散程度的量，使用方差作为样本间距：
于是问题转化为找到一个轴，使得样本空间所有点映射到这个轴后，方差最大。
注：这里的是映射到新轴后的。
在更高维空间，可以类似推广为寻找超平面，PCA具体步骤如下：
- 将样本均值归0（demean），于是：
- 求一个轴的方向w，使得所有样本特征映射到w方向后，有：
达到最大。由于有demean处理，即：
对上式进行向量运算，于是即求方向向量w，使得
达到最大值。
至此，问题已转化为求一个目标函数的最优化问题，将使用梯度上升法解决此问题。
梯度上升法求解PCA问题
推导过程和线性回归的梯度下降推导过程类似，不再给出过程，这里直接给出向量化后的梯度表达式：
接下来用梯度上升法在模拟的二维特征数据集上求解主成分：
'''模拟数据集'''
import numpy as np
import matplotlib.pyplot as plt
X = np.empty(shape=(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)
'''demean之后的X绘图'''
plt.scatter(X_demean[:,0],X_demean[:,1])
plt.show()
demean之后的图示：
demean后图示
demean之后的数据分布没有改变，只是坐标轴的位置移动了。
梯度上升法：
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):
'''方向向量模为1,epsilon取小些比较好'''
'''为验证梯度正确与否，调试梯度'''
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):
'''将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
这里为了验证梯度向量化是否正确，同样进行了梯度调试函数的编写。另外，在用梯度上升法求解PCA问题时需要注意以下问题：
- 传入的初始化w不能为0向量，因为0本身是一个极值点。
- 每次对w更新之后都要对w进行单位化，这对求解带来便利。
- 梯度上升求解PCA问题前不能对数据进行标准化，标准化会改变数据方差。
求解过程
首先初始化w的值，然后对梯度进行调试，接着用向量化的梯度表达式求解，由结果看出调试求解和向量化后的梯度上升求解结果是一样的，这说明我们的梯度表达式是正确的。
然后将求解结果直线绘制出来：
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(shape=(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)
w2
result
这个结果和实际数据集的斜率0.75对应的(4,3)单位向量是一致的。
plt.scatter(X2_demean[:,0],X2_demean[:,1])
plt.plot([0,w2[0]*30],[0,w2[1]*30],color='r')
plt.show()
ex_result
对于PCA的目标函数，和线性回归一样也是有数学解的，当然可以用数学解直接求解。也可以用随机梯度法，小批量梯度法。
在上面的PCA求解例子中，使用的是二维特征。二维的映射到1维，求得的一个主成分也叫第一主成分，但是如果是1000维总不太可能只映射到1维，可能要10个或更多维度。所以除了第一主成分，还要求第二第三...主成分。这在下篇会讲解。
网友评论