PCA主成份分析

原理

PCA主成分分析是一种非监督机器学习算法，主要用于数据的降维，相较于knn线性回归这两个监督学习的算法，其理论部分要稍显复杂一些，因此篇幅较长，若有不准确之处欢迎指正，同时后期发现需要删改或是增加的内容，我也会及时更新
首先我们来看一个具有两个特征的数据，其进行可视化后如图

原始两特征
我们会想，如果我们删掉其中一个特征，那么剩下另一个特征能不能表征这个原本具有两个特征的数据呢？
分别保留特征2,1
我们会认为右侧那个靠着特征1的情形要更加合适一些，因为其数据更加分散，有了区分度，那么我们能不能找到这样一条直线，使得数据点映射到他上面，具有“最好”的区分度呢？
某一条斜线
我们会觉得，原来的分布变化不是很大，但是数据点都集中在了这样一条直线上，于是我们实现了从二维到一维的这样一个过程。
或许你也发现了，这实际上就是一个坐标变换的过程
但是如何去选择这样的一条直线呢？

• 我们使用方差来代表样本间间距，方差大小代表出样本之间的紧密稀疏关系

  
  方差
• 于是问题变成了：找到这样一个轴/直线，使得样本空间所有点映射到这个轴之后，方差最大（此时最稀疏）
  等价问题
  具体有如下步骤
  1、将样本的均值归零（demean），这时坐标轴和方差公式变换成
  均值归零
  2、求出一个轴的方向w(w1,w2,...)使得样本映射到w方向后方差最大
  根据初中数学知识可知（其中X_project为均值归零后的X） 映射到w方向后的x的长度求法
• 于是问题变成：求该方差的目标函数的最优化问题，即求最大值，我们使用梯度上升法解决
  目标变成
  下面开始进行梯度上升法的过程

• 我们对其求导得

  
  求导得
  

  我们进行化简，同时观察右侧的矩阵和向量

  
  化简得
  我们观察到，左式等于，右侧上面的行向量和右侧下面矩阵的点乘然后将结果进行转置所得到的列向量，于是我们进一步化简得
  最终化简
  

  其中XT（Xw）是（Xw）T·X的转置结果

基本实现

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)#使用自己编造的数据

我们可以先看看这堆数据的散点图

原始数据

def demean(x):
    return x-np.mean(x,axis=0)#均值归零的函数
x_demean=demean(x)
#下面开始梯度上升法
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)#上面所求出的f导数方程函数
def direction(w):
    return w/np.linalg.norm(w)      #将向量转换为单位向量
def gradient_ascent(df,x,intial_w,eta,n_itera=1e4,epsilon=1e-8):
    w=direction(intial_w)
    cur_iter=0                       #记录当前循环次数
    while cur_iter<n_itera:
        gradient=df(w,x)
        last_w=w                      #记录上次的w
        w=w+eta*gradient
        w=direction(w)               #w要注意变成单位向量
        if(abs(f(w,x)-f(last_w,x))<epsilon):
            break
        cur_iter += 1
    return w                         #返回单位方向向量
intial_w=np.random.random(x.shape[1])#初始不能为零，因为0处为极小值
eta=0.001
gradient_ascent(df_math,x_demean,intial_w,eta)

最后得到结果

w向量结果

求出前N个主成分

正如我们前面所讨论的，PCA的实质就是改变了坐标系，一个原来是n维的数据，我们无论怎样去变化坐标系，他还是n维的
但是我们使用PCA很大一个目的就是降维，因为在变化坐标系后，我们可以发现，n维坐标系里有k维是很重要的，如果他们足以描述这组数据90%以上的话，我们何不妨就采用这k维数据呢，要知道，维数了下降可以很大程度上减小计算负担！

• 之前我们只在二维数据中求了一个主成分，那么我们现在开始求第二个主成分（当然了，二维最多也就两个成分）
  实际上原理类似于高中物理中力的分解，前面我们已经将X分解出了一个w向量，另一个于是就很容易得出了
  1、将数据进行改变，将数据在第一个主成分上的分量去掉。
  去除第一主成分分量
  2、去除第一主成分分量后，对剩下的分量执行主成分分析即得第二主成分分量
  代码前部分和上面基本相同

import numpy as np
#import PCA
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)#使用自己编造的数据
def demean(x):
    return x-np.mean(x,axis=0)#均值归零的函数
x=demean(x)
#下面开始梯度上升法
def f(w,x):
    return np.sum((x.dot(w)**2))/len(x)#求出方差的函数
def df(w,x):
    return x.T.dot(x.dot(w))*2/len(x)#上面所求出的f导数方程结果
def direction(w):
    return w/np.linalg.norm(w)      #将向量转换为单位向量
def first_component(df,x,intial_w,eta,n_itera=1e4,epsilon=1e-8):
    w=direction(intial_w)
    cur_iter=0                       #记录当前循环次数
    while cur_iter<n_itera:
        gradient=df(w,x)
        last_w=w                      #记录上次的w
        w=w+eta*gradient
        w=direction(w)               #w要注意变成单位向量
        if(abs(f(w,x)-f(last_w,x))<epsilon):
            break
        cur_iter += 1
    return w
intial_w=np.random.random(x.shape[1])#初始不能为零，因为0处为极小值
eta=0.001
w=first_component(df,x,intial_w,eta)
x2=np.empty(x.shape)#先创建一个和原来大小相同的空矩阵
for i in range(len(x)):
    x2[i]=x[i]-x[i].dot(w)*w
w2=first_component(df,x2,intial_w,eta)
plt.scatter(x2[:,0],x2[:,1])

分别得到w，w2,如下

w
w2
w2散点图

• 下面我们写一个可以求前n个主成分的函数

def first_n_components(n,x,eta=0.01,n_itera=1e4,epsilon=1e-8):
    x_pca=x.copy()
    x_pca=demean(x_pca)
    res=[]
    for i in range(n):
        intial_w=np.random.random(x_pca.shape[1])
        w=first_component(df,x_pca,intial_w,eta)
        res.append(w)
        x_pca=x_pca-x_pca.dot(w).reshape(-1,1)*w
    return res

利用该函数计算前面的例子

高维数据向低维数据转换

• 前面我们提到了，n维的X，在进行完坐标变换后仍为n维，但是我们只取了k维，那么我们如何将mn的n维X矩阵转换成k维呢？
  之前我们将一个样本与W（kn）每行相乘，得到的应该是样本在k个维度上的分量
  用矩阵乘法表示如下
  矩阵乘法
  当然同样的，我们可以反向得到一个Xm,只不过已经损失的信息不会恢复
  反向得到m*n
  这个实例和上面基本相同，区别只是增加了transform部分

class PCA:

    def __init__(self, n_components):
        """初始化PCA"""
        assert n_components >= 1, "n_components must be valid"
        self.n_components = n_components   #为传入的k值
        self.components_ = None           #为最后返回的w的矩阵

    def fit(self, X, eta=0.01, n_iters=1e4):
        """获得数据集X的前n个主成分"""
        assert self.n_components <= X.shape[1], \
            "n_components must not be greater than the feature number of X"

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

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

        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=0.01, 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

        X_pca = demean(X)
        self.components_ = np.empty(shape=(self.n_components, X.shape[1]))#k*n
        for i in range(self.n_components):
            initial_w = np.random.random(X_pca.shape[1])
            w = first_component(X_pca, initial_w, eta, n_iters)
            self.components_[i,:] = w

            X_pca = X_pca - X_pca.dot(w).reshape(-1, 1) * w

        return self

    def transform(self, X):
        """将给定的X，映射到各个主成分分量中"""
        assert X.shape[1] == self.components_.shape[1]

        return X.dot(self.components_.T)#m*n   n*k

    def inverse_transform(self, X):
        """将给定的X，反向映射回原来的特征空间"""
        assert X.shape[1] == self.components_.shape[0]

        return X.dot(self.components_)

    def __repr__(self):
        return "PCA(n_components=%d)" % self.n_components
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)#使用自己编造的数据
pca=PCA(n_components=2)
pca.fit(x)
pca.components_

结果和前面的例子一样

然后我们再进行转换

调用类中的transform函数得

至于inverse_transform也已经给出，这里不做演示了

使用sklearn中的PCA

我们使用sklearn进行前面例子的计算如下

import numpy as np
#import PCA
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)#使用自己编造的数据
from sklearn.decomposition import PCA
pca=PCA(n_components=1)
pca.fit(x)
pca.components_

得到结果

结果

然后进行转换如下

这里将两个转换函数都调用了
我们最后可以得到散点图

下面使用sklearn中的数据集进行分析

(同时我们也说过，pca很大的一个用处就是降维，所以我们在此比较比较降维的优势所在，同时提出一些新特点)

import numpy as np
from sklearn import datasets
import matplotlib.pyplot as plt
digits=datasets.load_digits()
x=digits.data
y=digits.target

from sklearn.model_selection import train_test_split
x_train,x_test,y_train,y_test=train_test_split(x,y,test_size=0.2,random_state=666)
#实际上这里size默认是0.2

我们可以看到x_train的大小，其有64个特征，即为64维

x_train.shape

%%time
from sklearn.neighbors import KNeighborsClassifier
knn=KNeighborsClassifier()
knn.fit(x_train,y_train)
knn.score(x_test,y_test)

可以看到，使用全部维度数据进行fit计算耗时和准确率如下

全维度下耗时和准确率

from sklearn.decomposition import PCA
pca=PCA(n_components=2)
pca.fit(x_train)
x_train_reduction=pca.transform(x_train)
x_test_reduction=pca.transform(x_test)
%%time
knn.fit(x_train_reduction,y_train)
knn.score(x_test_reduction,y_test)

我们先极端的试探一下降至2维，这时确实耗时大大缩减，但是准确率也大打折扣，这显然不是我们所乐意见到的

降维至2维后耗时和准确率

• 那么如何能把控降维的“度”呢，在仅损失较小准确率的基础上大幅缩短时间？
  这时我们需要用到pca中的explained_variance_ratio_函数（解释方差比例），他会返回对应轴能解释方差的所占比例

  
  2维情况下各占比例
  

  我们还可以通过可视化看到随着维数变化，解释效果占比的变化

plt.plot([i for i in range(x_train.shape[1])],
         [np.sum(pca.explained_variance_ratio_[:i+1]) for i in range(x_train.shape[1])])
plt.show

变化曲线
现在我们回到之前的问题，如何来保证其准确率呢？

pca=PCA(0.95)
pca.fit(x_train)

对调用PCA时，我们不再传k值了，而是传准确率，这里他便自动给我们找到相应准确率下最小的维度值进行使用
我们查看n_components_发现

维度选择值
结果

而且结果也是我们所期望的，缩短了大幅时间后准确率下降很小！

总结

PCA算法到这里就完全结束了，他是我们第一个接触的非监督学习算法，同时他的降维作用用处广泛，所以我后面还会聊一聊几个实例，限于篇幅，实在是不想还在同一篇文章里聊了，后续如果有需要改正或者增删的地方我会及时来更新的！！！