机器学习系列（十八）——主成分分析Principle Compo

2019-07-12 本文已影响24人 Ice_spring

主成分分析Principle Component Analysis

主成分分析（PCA）在统计学领域应用非常广泛，同时也是很重要的非监督机器学习算法，PCA主要用于数据的降维。在机器学习中，降维是很重要的一步预处理操作，通过降维，可以发现便于人类理解的特征和提取数据集主要特征。这样一来可以减少要处理的数据量，同时又不破坏数据整体特征，提高了算法的效率。PCA广泛应用于可视化和去噪过程。
以一个简单的数据降维例子介绍PCA：

数据特征

如图所示是某个数据集的两个特征，现在思考一个问题：能不能将两个特征压缩为一个特征呢？一个简单的处理方法是只取特征1或者特征2，这样就达到了降维的目的：

特征降维.png

而且显然取特征1会比特征2有更好的区分度（样本间距更大）。不过，还有没有更好的降维办法呢？考虑这样一条直线：

better_降维

此时将特征投影到红色直线上进行降维，显然和原始的特征更加接近，同时也更符合特征原始分布。那么如何找到这样一条让样本间距最大的轴呢？
首先给出样本间距定义，考虑到统计学中方差Variance是衡量样本分布离散程度的量，使用方差作为样本间距：
$var(x)=\frac{1}{m}\sum_{i=1}^{m}(x_{i}-\overline{x})^{2}$

于是问题转化为找到一个轴，使得样本空间所有点映射到这个轴后，方差最大。
注：这里的 $x_{i},\overline{x}$ 是映射到新轴后的。

在更高维空间，可以类似推广为寻找超平面，PCA具体步骤如下：

将样本均值归0（demean），于是： $var(x)=\frac{1}{m}\sum_{i=1}^{m}x_{i}^{2}$
求一个轴的方向w，使得所有样本特征映射到w方向后，有：
$var(X_{project})=\frac{1}{m}\sum_{i=1}^{m}||X_{project}^{(i)}-\overline{X}_{project}||^{2}$
达到最大。由于有demean处理，即：
$var(X_{project})=\frac{1}{m}\sum_{i=1}^{m}||{X_{project}^{(i)}}||^{2}$

对上式进行向量运算，于是即求方向向量w，使得
$var(X_{project})=\frac{1}{m}\sum_{i=1}^{m}||X^{(i)}.w||^{2}=\frac{1}{m}\sum_{i=1}^{m}(X_{1}^{(i)}w_{1}+X_{2}^{(i)}w_{2}+...+X_{n}^{(i)}w_{n})^{2}$

达到最大值。
至此，问题已转化为求一个目标函数的最优化问题，将使用梯度上升法解决此问题。

梯度上升法求解PCA问题

推导过程和线性回归的梯度下降推导过程类似，不再给出过程，这里直接给出向量化后的梯度表达式：
$\triangledown{f}=\frac{2}{m}X^{T}Xw$

接下来用梯度上升法在模拟的二维特征数据集上求解主成分：

'''模拟数据集'''
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个或更多维度。所以除了第一主成分，还要求第二第三...主成分。这在下篇会讲解。

机器学习系列（十八）——主成分分析Principle Compo

主成分分析Principle Component Analysis

梯度上升法求解PCA问题

猜你喜欢

热点阅读