机器学习-使用梯度上升求解PCA

什么是PCA

PCA 是一种数据处理分析方式，可以用来进行数据降维、降噪，在机器学习领域中，有时候由于样本采集的特征维度太多，造成维度灾难，我们可以使用PCA对样本数据进行降维

关于PCA的原理，这里有一篇介绍的非常详细的文章，简洁易懂又不失专业严谨

文章中求PCA使用的是 协方差矩阵、矩阵对角化等数学知识，但是今天这篇文章是使用梯度上升来求PCA主成分，关于PCA的原理 上面链接中的文章中已经写得很详细，我们直接从定义 效用函数开始

目标

找到一个向量 w1 使得，样本投影到该向量上之后，样本间分布尽可能稀疏，在统计学和概率论中 方差 用于衡量数据的离散程度，所以刚好用来衡量 样本投影后的 分布稀疏程度

定义效用函数

$$Var(X_{pro})=\frac{\sum _{i=1}^m(x_{pro}^i-{\overline{x}}_{pro}{)}^2}{m}$$

优化目标

最大化 Var(x)

注意事项

1. 我们要找的是一个向量 w, 使得样本投影到该方向上后 方差最大，也就是说我们关注的是 w 的方向，与 w的模没有关系，所以为了方便运算，我们只需要找到 w的单位向量就可以了
2. 需要对样本数据进行 demean 归零处理，是样本在各个维度上均值归零，这样做不会影响样本的相对分布，也就不会影响 w 向量的结果，并且对一个样本归零后，将样本投影到新的向量后，投影后的样本 均值仍然为0

样本均值归零后，对应的 效用函数也变成：

$$Var(X_{pro})=\frac{\sum _{i=1}^m{x_{pro}^{(i)}}^2}{m}$$

只要求出 $x_{pro}^{(i)}$ 就可以了

$$x_{pro}=cos\theta \cdot ||\vec{x}||=\frac{\vec{x}\cdot \vec{w}}{||\vec{x}||\cdot ||\vec{w}||}\cdot ||\vec{x}||$$
化简后
$$x_{pro}=\frac{\vec{x}\cdot \vec{w}}{||\vec{w}||}$$

由于 $\vec{w}$ 是单位向量， 模为1 所以

$$x_{pro}^{(i)}=\vec{x}\cdot \vec{w}$$

那么效用函数就变成

$$Var(X_{pro})=\sum _{i=1}^m({\vec{x}}^{(i)}\cdot \vec{w}{)}^2$$
假设样本x 有 n 个特征,我们将向量点乘 展开

$$Var(X_{pro})=\sum _{i=1}^m({\vec{x}}_1^{(i)}\cdot {\vec{w}}_1+{\vec{x}}_2^{(i)}\cdot {\vec{w}}_2...+{\vec{x}}_n^{(i)}\cdot {\vec{w}}_n{)}^2$$

现在我们使用梯度上升来最大化效用函数，需要求w各个维度上的偏导数

我们随机生成100个样本， 两个维度，y 值与x 值存在一定线性关系

import numpy as np

X = np.empty((100,2))

X[:,0] = np.random.uniform(0,100,100)

X[:,1] = X[:,0] * .75 + 3 + np.random.normal(0, 10,100)

import matplotlib.pyplot as plt

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

在这里插入图片描述

对数据进行demean 均值归零

对数据进行demean 均值归零 ，这样向量映射到 新的向量后 均值仍然为0

def demean(X):
    for i in range(0, X.shape[1]):
        X[:,i] -= np.mean(X[:,i])


demean(X) # 对各个维度归零

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

归零化后

效用函数

def effect(X, w):
    return np.sum(X.dot(w)**2)/len(X);

求梯度

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

梯度上升

在梯度上升或下降过程中，由于计算机浮点数精度影响，我们可能永远也无法算出导数为零，所以我们需要定义一个 epsilon 误差，当效用函数的值 与 上一次的值 误差绝对值小于误差，我们我们就认为到达极值

# 初始化一个方向
def gradient_ascent(init_w, X, eta, max_iters = 1000, epsilon = 1e-8):
    count = 0
    effect_array = []
    init_w /= np.linalg.norm(init_w)
    w = init_w
    while count < max_iters:
        effect_array.append(effect(X, w));
        gradient = df(X, w)
        print(gradient)
        w += (eta*gradient)
        w /= np.linalg.norm(w)
        count += 1
        if len(effect_array) > 1 and abs(effect_array[-1] - effect_array[-2]) < epsilon :
            print("上升{}次后到达极值".format(count))
            return w
    
    print("未到达极值")
    return w

# 使用梯度上升求出 w 方向，使得样本投影到w 方向后 样本 方差最大
w = gradient_ascent(np.random.random(X.shape[1]), X, 0.01)

[2044.48163869 1647.38290851]
[2176.24815427 1815.43332636]
[2174.93222982 1819.21344751]
[2174.75206482 1819.44212419]
[2174.73759786 1819.45949649]
[2174.73647111 1819.46084375]
[2174.73638355 1819.4609484 ]
上升7次后到达极值

plt.scatter(X[:,0],X[:,1])
plt.plot([w[0]*-50,w[0]*50],[w[1]*-50,w[1]*50], color ="red")

求第二主成分

将所有样本看成一个向量，首先所有样本减去映射到第一主成分上的向量，然后对剩余分量再次使用梯度上升求 w2 向量

X_pro = X.dot(w).reshape((100,1))

水平堆叠一下，方便使用 numpy的乘法来计算坐标

X_pro = np.hstack([X_pro, X_pro]) 

计算投影到第一主成分后的坐标

X_pro*=w 

plt.scatter(X_pro[:,0], X_pro[:,1])

所有样本减去第一主成分上的分量

X_second = X - X_pro

继续使用梯度上升求第二主成分

w2 = gradient_ascent(np.random.random(X_second.shape[1]), X_second, 0.01)

[-2.98245214  3.56481805]
[-6.78282896  8.10727212]
[-15.25169868  18.22980827]
[-32.5089275   38.85675475]
[-57.58452703  68.82871925]
[-75.03383634  89.68525265]
[-80.63882372  96.38469298]
[-81.86783669  97.8536881 ]
[-82.11065404  98.14391897]
[-82.15759047  98.2000204 ]
[-82.16662452  98.21081847]
[-82.16836191  98.21289511]
[-82.16869598  98.21329441]
[-82.16876021  98.21337119]
[-82.16877256  98.21338595]
[-82.16877494  98.21338879]
[-82.1687754   98.21338934]
[-82.16877548  98.21338944]
[-82.1687755   98.21338946]
[-82.1687755   98.21338946]
上升20次后到达极值

X_pro2 = X_second.dot(w2).reshape((100,1))

X_pro2 = np.hstack([X_pro2, X_pro2])

X_pro2 *= w2

plt.scatter(X[:,0],X[:,1])
plt.plot([w2[0]*-50,w2[0]*50],[w2[1]*-50,w2[1]*50], color ="green")
plt.plot([w[0]*-50,w[0]*50],[w[1]*-50,w[1]*50], color ="red")
plt.scatter(X_pro[:,0], X_pro[:,1],color="red",alpha = .2)
plt.scatter(X_pro2[:,0], X_pro2[:,1], color ="green", alpha = .2)
plt.xlim(-60,60)
plt.ylim(-60,60)

如图中所示 红色线表示 w1 第一主成分， 半透明红点 表示样本投影到该向量上后的 位置，绿色直线表示第二主成分w2 方向，半透明绿点 表示 样本在w2 方向上的投影