数据降维2：PCA算法的实现及使用

最新推荐文章于 2022-11-28 10:57:57 发布

木东居士

最新推荐文章于 2022-11-28 10:57:57 发布

阅读量785

点赞数

0x01 PCA算法梯度求解

我们在上篇文章《数据降维1：主成分分析法思想及原理》的最后已经通过推导，将数据映射转化为求目标函数的最优化问题：

求，使得最大

对于最优化问题，除了求出严格的数据解以外，还可以使用搜索策略求极值。

在求极值的问题中，有梯度上升和梯度下降两个最优化方法。梯度上升用于求最大值，梯度下降用于求最小值。

1.1 梯度上升&梯度下降

我们在《还不了解梯度下降法？看完这篇就懂了！》这篇文章中已经详细地介绍了梯度下降的思想及推导问题。已知：使损失函数值变小的迭代公式：

我们现在要求极大值，则使用梯度上升法。梯度的方向就是函数值在该点上升最快的方向，顺着这个梯度方向轴，就可以找到极大值。即将负号变为正号：

1.2 求梯度

现在回到对PCA算法求解的问题上来，对于 PAC 的目标函数：

求使得最大，可以使用梯度上升法来解决。

首先要求梯度。在上式中，是未知数是非监督学习提供的已知样本信息。因此对中的每一个维度进行求导：

对上式进行合并，写成两个向量点乘的形式。更进一步对表达式进行向量化处理

得到梯度为：

有了梯度，就可以使用梯度上升法求最大值了。

0x02 求解第一主成分代码实现

我们已经知道了PCA算法的目标函数（方差），以及如何使其最大（梯度上升）。下面使用代码来实现公式及算法流程。

下面我们以二维数据为例，将其映射到一维，即求出一系列样本点的第一主成分。

2.1 数据准备

首先准备数据：

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)


plt.scatter(X[:,0],X[:,1])
plt.show()

2.2 函数实现

接下来对数据进行主成分分析，第一步是均值归零，定义相应的函数。

将样本进行均值归0（demean），即所有样本将去样本的均值。样本的分布没有改变，只是将坐标轴进行了移动。

其方法是：X表示一个矩阵，每一行代表一个样本，每一个样本的每一个特征都减去这个特征的均值。实现如下：

def demean(X):
    # axis=0按列计算均值，即每个属性的均值，1则是计算行的均值
    return (X - np.mean(X, axis=0))


X_demean = demean(X)
# 注意看数据分布没变，但是坐标已经以原点为中心了
plt.scatter(X_demean[:, 0], X_demean[:, 1])
plt.show()

接下来就是对目标（方差）函数和梯度（导数）函数的定义。

首先定义目标函数：

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

在这里，有一处需要注意的地方：

对于要求的轴，向量，实际上一个单位向量，即要让向量的模为1

因此我们使用np中的线性代数库linalg里面的norm()方法，求第二范数，也就是求算术平方根

def direction(w):
    return w / np.linalg.norm(w)

然后就实现梯度上升代码的流程了：gradient_ascent()方法为梯度上升法的过程，在n_iters次循环中，每次获得新的，如果新的和旧的对应的函数的值的差值满足精度epsilon，则表示找到了合适的 $w$ ：

此处要注意，对于每一次计算向量之前都要进行单位化计算，使其模为1。

# 梯度上升法代码
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_math(w,X)
        last_w = w
        w = last_w + eta * gradient
        w = direction(w)    # 将w转换成单位向量
        if (abs(f(w,X) - f(last_w,X)) < epsilon):
            break
        cur_iter += 1
    return w

接着使用梯度上升法，先设定一个初始的和学习率：

这里面还需要注意一点：线性回归中，通常将特征系数θ的值设为全部为0的向量。但在主成分分析中w的初始值不能为0！！！ 这是因为如果将w=0带入梯度求导的公式中，也是每次求梯度得到的都是没有任何方向的0。所以要设置一组随机数。

initial_w = np.random.random(X.shape[1])
eta = 0.001

然后执行gradient_ascent()方法，计算梯度上升的结果，即要求的轴。

w = gradient_ascent(df_debug, X_demean, initial_w, eta)
# 输出
array([0.76567331, 0.64322965])

2.3 结果可视化

进行可视化，轴对应的方向，即将样本映射到该轴上的方差最大，这个轴就是一个主成分（第一主成分）：

plt.scatter(X_demean[:,0],X_demean[:,1])
plt.plot([0,w[0]*30],[0,w[1]*30], color='red')
plt.show()

这样，就求出了二维数据的1个主成分，二维数据映射到一维就足够了。

但是如果是高维数据，可能就要映射到多维上。即求出第n主成分。

0x03 求前n个主成分

3.1 思想

在实际的降维过程可能会涉及到数据在多个维度的降维，这就需要依次求解多个主成分。

求解第一个主成分后，假设得到映射的轴为所表示的向量，如果此时需要求解第二个主成分怎么做?

答案是：需要先将数据集在第一个主成分上的分量去掉，然后在没有第一个主成分的基础上再寻找第二个主成分。

如上图所示，样本在第一主成分上的分量为，也就是图中的蓝色向量，其模长是向量和向量的乘积，又因为是单位向量，向量的模长乘以方向上的单位向量就可以得到这个向量。

求下一个主成分就是将数据在第一主成分上的分量去掉，再对新的数据求第一主成分。

那么如何去掉第一主成分呢？用样本减去分量，得到结果的几何意义就是一条与第一主成分垂直的一个轴。这个轴就是样本数据去除第一主成分分量后得到的结果，即图中的绿色部分：。

然后在新的样本中求第一主成分，得到的就是第二主成分。循环往复就可以求第n主成分。

3.2 求第二主成分的实现

首先，我们要将数据集在第一个主成分上的分量去掉。即用样本减去样本在第一主成分分量上的映射，得到。

这里可以用向量化的思想：

矩阵（的矩阵）与向量（的矩阵）点乘，得到的是一个的向量，这个元素表示中的每一个样本映射到向量方向上的模长
然后对向量reshape，使其成真正的的矩阵
再乘以，得到的矩阵，也就是把矩阵中每个样本的每个维度在方向上的分量。即矩阵
用样本减去样本在第一主成分上的映射

X_new = X - X.dot(w).reshape(-1,1) * w

是什么样的，可以可视化看一下：

plt.scatter(X_new[:,0], X_new[:,1])
plt.show()

所有数据把第一主成分的分量去掉，得到的就是与第一主成分的分量垂直的向量。

由于整个数据只有两个维度，第一个维度去掉之后，剩下的第二维度就是所有的内容了。因此样本在第二维度上的分布完全在一条直线上，不再有其他方向上分布的方差。

如何第二主成分的轴？就是将去掉第一主成分的分量的带入到gradient_ascent()函数中，即可以得到新的轴。

w_new = gradient_ascent(df_math, X_new, initial_w, eta)
# 输出：
array([-0.64320916,  0.76569052])

和彼此之间互相垂直，点乘为0。

3.3 求前n主成分

下面接可以整合代码。first_n_component()方法用于求的前n个主成分。

首先进行数据归零操作，使用列表存储前n个主成分的方向向量。在for循环中每次求主成分时都设定一个随机的（initial_w），使用梯度上升法得到一个主成分后就去掉在这个主成分上的分量X_pca = X_pca - X_pca.dot(w).reshape(-1, 1) * w：

def first_n_component(n, X, eta=0.001, 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 = gradient_ascent(df_math, X_pca, initial_w, eta)
        res.append(w)
        X_pca = X_pca - X_pca.dot(w).reshape(-1, 1) * w   
    return res