一、目标函数的梯度求解公式
- PCA 降维的具体实现，转变为：
- 方案：梯度上升法优化效用函数，找到其最大值时对应的主成分 w ；
- 效用函数中，向量 w 是变量；
- 在最终要求取降维后的数据集时，w 是参数；
1）推导梯度求解公式
- 变形一
- 变形二
- 变形三：向量化处理
- 最终的梯度求解公式：▽f = 2 / m * XT . (X . dot(w) )
二、代码实现（以二维降一维为例）
-
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)
-
2）查看数据分布
plt.scatter(X[:, 0], X[:, 1]) plt.show() -
3）对原始数据集 demean 处理：得到新的数据集，数据集的每一种特征的均值都为 0
def demean(X): return X - np.mean(X, axis=0) X_demean = demean(X)
- 数据集的每一个特征，都减去该列特征的均值，使得新的数据集的每一个特征的均值为 0；
- np.mean(X, axis=0)：得到矩阵 X 每一列的均值，结果为一个列向量；
-
4）使用梯度上升法求解主成分
-
求当前参数 w 对应的目标函数值（按推导公式求解）
def f(w, X): return np.sum((X.dot(w)**2)) / len(X)
-
求当前参数 w 对应的梯度值（按推导公式求解）
def df_math(w, X): return X.T.dot(X.dot(w))*2 / len(X)
-
求当前参数 w 对应的梯度值（按调试公式求解）
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
-
将向量 w 转化为单位向量
def direction(w): return w / np.linalg.norm(w)
# 因为推导公式时人为的将 w 向量的模设为 1，如果不将 w 向量转化为单位向量，梯度上升法的搜索过程会不顺畅；
# np.linalg.norm(向量)：求向量的模
# 向量 / 向量的模 == 单位向量 -
梯度上升法的优化过程
def gradient_ascent(df, X, initial_w, eta, n_iters=10**4, epsilon=10**-8): # 将初始化的向量 initial_w 转化为单位向量 w = direction(initial_w) cur_iter = 0 while cur_iter < n_iters: gradient = df(w, X) last_w = w w += eta * gradient # 注意1：将每一步优化后的向量 w 也转化为单位向量 w = direction(w) if (abs(f(w, X) - f(last_w, X)) < epsilon): break cur_iter += 1 return w
# 注意1：每一步优化后的向量 w 都要转化为单位向量
-
5）求解数据的第一主成分
-
初始化
initial_w = np.random.random(X.shape[1]) eta = 0.001# 注意2：初始化 w 不能为 0 向量，因为在梯度公式中，若 w 为 0 ，无论什么数据，计算的结果都是没有任何方向的 0；
-
求解并绘制出第一主成分
w = gradient_ascent(df_math, X_demean, initial_w, eta) plt.scatter(X_demean[:,0], X_demean[:,1]) # 此处绘制直线，对应的两个点是（0, 0）、（w[0]*30, w[1]*30）,构成了红色的直线 # 由于 w[0]、w[1]数值太小，为了让直线看的更清晰，将该点扩大 30 倍 # 这也是绘制向量的方法 plt.plot([0, w[0]*30], [0, w[1]*30], color='r') plt.show()
# 注意3：不能用 StandardScaler 标准化数据
# 原因：由于PCA的过程本身就是求一个轴，使得所有的样本映射到轴上之后方差最大，但是，如果对数据标准化后，则样本的方差就变为 1 了，就不存在方差最大值了；
# 其实 demean 的过程就是对数据标准化处理的一部分过程，只是没让数据的标准差为 1；
# 在梯度下降法求解线性回归问题时需要对数据做归一化处理；
6）求前 n 个主成分
-
求下一个主成分的思路
- 思路：去掉当前数据集在上一个主成分方向上的分量，得到一个新的数据集，求新的数据集的主成分；
- 操作：将当前数据集的每一个样本，减去该样本在上一个主成分上的分量
- 步骤：
# 1）第一步：求新的数据集 X2 = np.empty(X.shape) for i in range(len(X)): X2[i] = X[i] - X[i].dot(w) * w # for 循环向量化：X2 = X - X.dot(w).reshape(-1, 1) * w # X[i].dot(w)：表示样本 X[i] 在上一主成分 w 方向上的膜 # X[i].dot(w) * w：表示样本 X[i] 在上一主成分 w 上的分量（为向量） # X2[i] = X[i] - X[i].dot(w) * w：表示样本 X[i] 减去其在上一主成分上的分量后的新的样本向量 # 2）第二步：初始化数据，求主成分 initial_w = np.random.random(X.shape[1]) eta = 0.01 w2 = first_component(X2, initial_w, eta)
-
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) # demean函数，将原始数据集 demean 处理 def demean(X): return X - np.mean(X, axis=0) # 求当前变量 w 对应的目标函数值 def f(w, X): return np.sum((X.dot(w)**2)) / len(X) # 求当前变量 w 对应的梯度的值 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=10**4, epsilon=10**-8): w = direction(initial_w) cur_iter = 0 while cur_iter < n_iters: gradient = df(w, X) last_w = w w += eta * gradient w = direction(w) if (abs(f(w, X) - f(last_w, X)) < epsilon): break cur_iter += 1 return w # 求取 X 的前 n 个主成分 # 求 X 的前 n 个主成分时，只需要对 X 做一次 demean 操作 def first_n_components(n, X, eta=0.01, n_iters=10**4, epsilon=10**-8): X_pca = X.copy() X_pca = demean(X_pca) res = [] # 循环 n 次，每次求出一个主成分，并将 n 个主成分放在列表 res 中返回 for i in range(n): # 初始化搜索点：initial_w initial_w = np.random.random(X_pca.shape[1]) # 求出主成分 w w = first_component(X_pca, initial_w, eta) res.append(w) # 获取下一次主成分的数据集：X_pca，用于计算下一个主成分 X_pca = X_pca - X_pca.dot(w).reshape(-1,1) * w return res # 求数据集 X 的前 2 个主成分 first_n_components(2, X) # 输出：[array([0.7993679 , 0.60084187]), array([-0.60084187, 0.7993679 ])]