主成分分析的课堂小结(1)

这篇课堂笔记来源于视频：
https://www.bilibili.com/video/av22634854/?spm_id_from=333.338.recommend_report.4

PCA 可以降维，去噪，有利于可视化。

那么PCA的原理是什么呢？
举个简单的例子，我们有几个存在两个特征的样本点，如下：

我们希望减少到一个特征，那怎么做呢？显然我们最先想到的是，直接只保留其中一个特征就可以啦。那么分别保留特征一和特征二，可以有如下的结果：

那么现在大家判断下，哪个结果更好些呢？
显然是右侧的方案更好一些，因为点与点之间的距离更大些，也就是说，可区分性更强一些。

那么，有没有更好的方案了呢？
看以下这个维度，是不是间距更大了呢。

所以现在问题来了，我们怎么才能找到这个维度使得样本间的可区分度越大呢？
所谓的样本间距离最大，也就是方差更大，那么我们的任务就是：

首先，我们需要将样本的均值归为零，这会带来什么变化呢，请看下图：

第二步，求一个单位方向向量w，使得映射到此方向上后的样本点间的方差最大：

由以上可得，我们把问题归纳成了求一个目标函数最大值的问题，那么我们可以用梯度上升法来求解：

X代表样本空间，W代表想求的方向，m代表样本的数量，n代表样本的特征数目，

下面用代码实现上述过程：

首先设置一些二维样本，两个特征间有线性关系，使得之后有一个更明显的隆维效果

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()

生成的样本图：

下一步将样本的均值变为零：
np.mean(x, axis=0)的意思是对矩阵x的每一列求均值即每个维度的特征求均值。

def demean(x):
    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)

梯度上升法：
direction是为了将w变为单位方向向量。df为求导后的函数，eta为学习率，n_iters为上升次数，epsilon为设定的最少下降程度，一旦小于这个值就停止上升。

def direction(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)
        if(abs(f(w,x)-f(last_w,x))<epsilon):
            break
        cur_iter += 1
    return w

随机初始化initial_w,注意不能初始化为0,因为当w＝0时，是目标函数的极小值，梯度为0：

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

下面展示w的位置：

w = gradient_ascent(df_math,x_demean,initial_w,eta)
plt.scatter(x_demean[:,0],x_demean[:,1])
plt.plot([0,w[0]*30],[0,w[1]*30],color='r')
plt.show()

其实以上所求的的w表示第一个主成分，那么，怎么求其它成分即前n个主成分呢？

下面我们用代码来实现一下：
x[i].dot(w)表示x[i]在w方向上的模，再乘以w就表示，在w方向上的向量。
w为单位方向。
x［i］表示第i个样本。

x2 = np.empty(x.shape)
for i in range(len(x)):
    x2[i] = x[i] - x[i].dot(w)*w

之后就可重复之前的过程，求x2的主成分。