主要内容

PCA简介
使用梯度上升法解决PCA问题

1.PCA简介

PCA（Principal Component Analysis : 主成分分析）：也是一个梯度分析的应用，不仅是机器学习的算法，也是统计学的经典算法

一个非监督的机器学习算法
主要用于数据的降维
通过降维，可以发现更便于人类理解的特征
其他应用：可视化，去噪。

1.1 PCA举例

例如下面一个两个特征的一个训练集，我们可以选择一个特征，扔掉一个特征：

下图分别是扔掉了特征一和特征二的两种方案，很明显右边这种的效果会更好一些，因为访问二扔掉特征二以后，点之间的分布情况更接近与原图，但是这不是更好的

我们希望有一根直线，是斜着的，我们希望将所有的点都映射到这条直线上，那么这个时候我们就成功的将二维降到了一维，与此同时，这些点更加趋近与原来的点的分布情况。
换句话说，点和点之间的距离比无论是映射到x还是映射到y周，他们之间的区分度都更加的大，也就更加容易区分。

1.2 分析

那么如何找到这个让样本间间距最大的轴? 如何定义样本间间距? 事实上有一个指标可以之间定义样本间的距离，就是方差（Variance）（方差：描述样本整体之间的疏密的一个指标，方差越大，代表样本之间越稀疏，方差越小，代表样本之间越紧密）。

方差公式：

为了简化方差计算，我们可以将样例的均值归为0 (demean)，使得均值为0，这样可以简化方差的公式：

对于右上角的(w1,w2),我们将其理解为一个二维向量，因为是在二维空间中，这个维度是由样本数据为维度决定的，当样本的维度为多维时，这个w也是多维向量，类似于线性回归中的theta向量，但是theta值是直线x的系数，w向量就只是一个向量。

向量知识可以参考本节后的数学知识--向量。

1.3 推导

我们想要求一个轴的方向w = (w1,w2),是的我们所有的样本映射到w后，有：

其中X(project)为样本映射到w向量后对应的特征值，也可以理解为一个向量：

进行方差归一化后：

X(i)映射到w的距离实际上就是X(i)与w的点乘（棕色的线），根据定义推导，其值实际上就是Xproject，其中w看作单位向量。

那么：

此时我们的目标函数就可以化简为：

其中X(i)为样本，w为我们要找的向量，这个向量要做到上式结果最大。
我们的样本X的特征维度时几维，我们要找的w就是几维。

1.4 与线性回归的区别

线性回归：

PCA降维

1.主成分分析法（PCA）的两个轴都是特征，线性回归y轴是目标结果值。

主成分分析法的点是垂直于方差轴直线的，目标是让这些点投影到这个线上后的间距尽可能的大。
而线性回归的点是垂直于x轴的，目标是让这些点特征值对应值与该线预测值之间结果尽可能的小。

2. 使用梯度上升法解决PCA问题

2.1 公式推导--梯度值

对于梯度值得计算，与J函数一样，对f函数对每一个w求偏导组成得向量：

注意上面式子里的每一个(X1(i)·w1+X2(i)·w2+......Xn(i)·wn)都是一个X(i)和w的点乘，所以式子可以进一步化解：

化简过后可以进行向量化，也就是去掉∑(m)，其每一项拆开开都是m个标量相加。每一项的前半部分都一样(向量)，后半部分组合起来的值为X矩阵的转置，那么这样的话就变成了X矩阵的转置与向量的点乘：(X (*) w )*X(T)。
其中：

X：m行n列
w：n行1列

这样计算出来的结果为n行1列。

2.2使用梯度上升法代码实现PCA

1.创建数据--画图

import numpy as np
import matplotlib.pyplot as plt

X = np.empty((100, 2))
# 将X的第0列设置为0-100的随机小数
X[:,0] = np.random.uniform(0., 100., size=100)
# 第2列数据 = 0.75 * X[:,0] + 3 + noise
X[:,1] = 0.75 * X[:,0] + 3. + np.random.normal(0, 10., size=100)

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

2.demean(归0)

上述公式使用的前提有两个：

要对数据进行归0处理。
w必须为单位向量。

那么这里我们一一处理：

1.归0

假设X为m行n列的矩阵：
np.mean(X, axis=0)：对X中数据求均值。

axis = 0，对每一列求均值，返回一行n列向量。
axis = 1，对每一行求均值，返回m行一列向量。

归0操作就是每一列数据减去该列的平均值组成的新的向量。
demean：

def demean(X):
    return X - np.mean(X, axis=0)
X_mean = demean(X)

re = np.mean(X_mean[:,0])
print(re)       # 2.5437429940211587e-14

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

2.转化单位向量函数:

np.linalg.norm(w)，w向量的模

# 单位向量w
def dirextion(w):
    return w/np.linalg.norm(w)

3.梯度上升法

1.f函数
f函数，也就是效用函数。

f():

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

2.df函数

df函数分为df_math()和df_debug():

df_math()按照公式计算：

df_debug()与之前的梯度下降法调试中的一样：

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

def df_debug(X,w,epsilon=0.0001):
    gradient = np.empty(len(w))
    
    for i in range(len(w)):
        w_1 = w.copy()
        w_1[i] += epsilon
        w_2 = w.copy()
        w_2[1] -=epsilon
        
        gradient[i] = (f(X,w_1) - f(X,w_2)) /(2 *epsilon)
        
    return gradient

3.梯度上升函数

该函数的逻辑与梯度下降相似，只是在每次循环中，都要调用单位向量函数将w转为单位向量进行计算。

def gradient_ascent(df, X, initial_w, eta, n_iters = 1e4, epsilon=1e-8):
    w = dirextion(initial_w)
    i_iter = 0

    while i_iter < n_iters:
        gradient = df(X,w)
        last = w
        w = w + eta * gradient
        w = dirextion(w)

        if abs(f(X,w) - f(X,last)) < epsilon:
            break
        i_iter +=1

    return w

4. 测试：

# 设置参数
init_w = np.random.random(X.shape[1])
eta = 0.001

# df_debug
re = gradient_ascent(df_debug,X_mean,init_w,eta)
print(re)    # [0.74936349 0.66215886]

# df_math
re = gradient_ascent(df_math,X_mean,init_w,eta)
print(re)      # [0.78531688 0.61909401]

5.画图

将所有样本点描在图中，并将w向量放大100倍画在图中：

plt.scatter(X_mean[:,0],X_mean[:,1])
plt.plot([0,re[0]*100],[0,re[1]*100],color='r')
plt.show()

6.使用极端数据集测试

1.生成数据集：

数据集的第二列不加噪音，也就是这些样本点全部在一条直线上。

import numpy as np
import matplotlib.pyplot as plt


X2 = np.empty((100, 2))
X2[:,0] = np.random.uniform(0., 100., size=100)
X2[:,1] = 0.75 * X2[:,0] + 3.

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

2.数据处理,测试数据，画图：

from PCA.pca_01 import demean,df_math,gradient_ascent

# 数据集归0
X2_mean = demean(X2)

# 设置参数
init_w = np.random.random(X2.shape[1])
eta = 0.001
# 梯度上升
re = gradient_ascent(df_math,X2_mean,init_w,eta)
print(re)

# 画图对比
plt.scatter(X2_mean[:,0],X2_mean[:,1])
plt.plot([0,re[0]*10],[0,re[1]*10],color='r')
plt.show()