机器学习 -- PCA（Ⅱ 梯度上升法解决主成分分析问题）

一. 原理

经过上一节，我们有如下结论。

1. 目标：求w，使得f(X)最大。其中：

$f(X)= \frac{1}{m} \sum _{i=1}^m(X^{(i)}_1w_1 + X^{(i)}_2w_2 + ... + X^{(i)}_nw_n)^2$

2. 梯度可表示为：

即：

对其进行向量化处理：

也可写成：

故将梯度向量化结果为：

二. 实现

1. 生成数据集

（1）导入需要的模块和包

import numpy as np
import matplotlib.pyplot as plt

（2）生成虚拟的测试用例。生成的X有两个特征，特征1为1~100之间随机分布的100个样本，特征2和相应的特征1有基本的线性关系0.75 * 特征1 + 3 + 噪音。（线性关系时降维效果比较明显）

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)

查看图像：

2. demean过程

（1）demean函数：

def demean(X):# np.mean(X, axis=0)：对X矩阵在行方向上求得均值，最终结果为每一列的均值。# 若X为m * n，则操作后得到1 * n的向量return X - np.mean(X, axis=0)

（2）测试一下函数：

# 对X矩阵数据均值归零操作
X_demean = demean(X)

查看图像，发现具体绘制数据情况和原先X大多一致。观察坐标轴，中心位置大概为(0, 0)。

验证：查看X和demean后的X在两个维度上的平均值，发现demean后的平均值都接近于0。

np.mean(X_demean[:, 0]), np.mean(X_demean[:, 1])

3. 使用梯度上升法求解数据的主成分

① 相关函数的实现：

（1）根据目标函数：

$f(X)= \frac{1}{m} \sum _{i=1}^m(X^{(i)}_1w_1 + X^{(i)}_2w_2 + ... + X^{(i)}_nw_n)^2$

计算目标函数的方法f(w, X)：

# X为归零化的矩阵
def f(w, X):return np.sum((X.dot(w) ** 2)) / len(X)

（2）根据目标函数对应梯度的计算公式：

$\triangledown f = \frac{2}{m} \cdot X^T(Xw)$

计算目标函数对应的梯度的方法df_math(w, X)：

def df_math(w, X):return X.T.dot(X.dot(w)) * 2. / len(X)

（3）验证梯度计算是否正确的函数df_debug(w, X, epsilon)：

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] += epsilonw_2 = w.copy()w_2[i] -= epsilonres[i] = (f(w_1, X) - f(w_2, X)) / (2 * epsilon)return res

（4）单位向量化函数direction(w)（公式为：向量 / 模长）：

# 计算单位向量：w / |w|
def direction(w):return w / np.linalg.norm(w)

（5）梯度上升函数 gradient_ascent（注意先将w变成单位向量）：

# 梯度上升法过程，代码结构与梯度下降法完全一样
# 注意：w为单位向量
def gradient_ascent(df, X, inital_w, eta, n_iters = 1e4, epsilon=1e-8):w = direction(inital_w)cur_iter = 0while cur_iter < n_iters:gradient = df(w, X)last_w = ww = w + eta * gradient# 注意：每次求单位向量w = direction(w) if (abs(f(w, X) - f(last_w, X)) < epsilon):breakcur_iter += 1return w

注意：对于PCA问题对样本数据标准化归一化过程不能使用StandardScaler！！！

为什么？由于PCA过程本身就是要求一个轴使得所有的样本映射到那轴之后样本的方差最大，但是一旦我们将我们样本数据进行标准化，样本的方差就变成1，方差的最大值就不存在了。这样一来就求不出PCA真正想最大化的结果了。

实际上demean过程已经做了一半，将均值变为0了，少做的一半就是不让我们数据的标准差为1。

② 调用梯度上升法：

（1）初始化w值：不能为零向量。

# 调用梯度上升法
# 初始化w值。注意：不能为零向量，否则带入都为0
init_w = np.random.random(X.shape[1])  # 个数为样本特征数
init_w
# array([0.51062981, 0.24947237])

（2）分别使用测试和计算的方法，得出的w向量相同。

eta = 0.001
# 梯度上升法求得的原来的数据将它映射到此轴上方差最大# 使用测试的方法
gradient_ascent(df_debug, X_demean, init_w, eta)
# array([0.76696199, 0.64169254])# 使用计算的方法
w = gradient_ascent(df_math, X_demean, init_w, eta)
w
# array([0.76696199, 0.64169254])

（3）绘制出图像对w进行观察：

plt.scatter(X_demean[:, 0], X_demean[:, 1])
# 乘以30是以便观察， 因为方向向量非常小
plt.plot([0, w[0] * 30], [0, w[1] * 30], color='r')
plt.show()

绘制出的红线即为样本要映射到的方向向量，使得样本间方差最大。相应的轴即为一个主成分，又称为第一主成分。