机器学习基础-5.PCA和梯度上升

一、PCA

1.PCA思想

PCA（principle component analysis），即主成分分析法，是一个非监督的机器学习算法，主要用于对数据的降维，通过降维可以发现更便于人理解的特征，加快对样本有价值信息的处理速度，此外还可以应用于可视化（降到二维）和去噪。

例如下图所示，样本有2个特征，现在对该样本进行降维处理。首先考虑直接选择特征1或者特征2降维，经过降维后的样本由2维降到1维，如图所示。

选择特征1降维。

选择特征2降维。

  

可以发现，相对于选择特征2，选择特征1降维后，样本间的距离更大，即样本的可区分度更高。那是否还有其他的映射方式，使得映射后样本的间距更大，事实上还可以选择某个轴线，例如下图所示，样本映射到该轴线之后，有更大的间距。

PCA降维的思想就是寻找某个轴线，使得样本映射到该轴线上后，能够有最大的可区分度。衡量可区分度的指标即方差（事实上，方差背后的意义就是数据的稀疏程度），现在的问题就是如何求得该轴线，使得方差的值最大。

2.求解思路

用方差来定义样本的间距，方差越大表示样本分布越稀疏，方差越小表示样本分布越密集，方差的公式如下。

在求解最大方差前，为了方便计算，可以先对样本进行demean（去均值）处理，即减去每个特征的均值，这种处理方式不会改变样本的相对分布（效果就像坐标轴进行了移动）。去均值后，样本x每个特征维度上的均值都是0，方差的公式转换成图示的公式。在这里，xi代表已经经过映射后的某样本。

对于只有2个维度的样本，现在的目标就是：求一个轴的方向w=（w1,w2），使得映射到w方向后，方差最大。目标函数表示如下。表示经过映射后的样本。

如下图所示，原始样本和映射后的样本有如下的空间几何关系。w为单位向量，所以w的模为1。

经过什么的推导，最终求解的公式可以转换成如下形式。

样本x为已知，现在的目标转换成如下。

这就转换成了求目标函数的最优化问题，可以用梯度上升法进行求解。

二、梯度上升法求解

1.公式求解

最终公式可以化简为如下形式。

梯度上升法的思路和梯度下降法的思路完全一样，只是前者求极大值，后者求极小值。可以参考梯度下降法博客内容。

2.代码实现

import numpy as np  
import matplotlib.pyplot as plt  
x = np.empty((100,2))   #x为100个样本，2个特征值的矩阵
x[:,0]= np.random.uniform(0.,100.,size = 100) #样本x的第1个特征值为均值为0，方差为100的100个值
x[:,1]= 0.75*x[:,0]+3.+np.random.normal(0,10.,size=100) #样本x的第2个特征值和第1个特征值保持一定的线性关系，这里引入噪声。这里之所以使特征1和特征2有一定的线性关系，仅仅是为了降维后，能够有较好的展示。

plt.scatter(x[:,0],x[:,1]) #将特征1和特征2绘制出来
plt.show()

初始样本准备完毕，现在就开始对该样本进行降维处理。

def demean(x):
    return x-np.mean(x,axis=0)  #首先对每个特征值进行去均值，axis=0代表按列求均值
def f (w,x):
    return np.sum((x.dot(w)**2))/len(x) #即方差公式，len(x)即样本的个数
def df_math(w,x):
    return x.T.dot(x.dot(w))*2./len(x)  #即通过数学推导得出的求梯度公式
def direstion(w):
    return w /np.linalg.norm(w)    #将w进行归一化，因为w仅代表方向，而每次求解后，w模不一定为1
def gradietn_ascent(df,x,initial_w,eta,n_iters = 1e4,epsilon=1e-8): #梯度上升法，传入初始w，迭代次数和停止迭代的最小差值
    w = direstion(initial_w) #先将传入的w进行归一化
    cur_iter =0
    while cur_iter <n_iters: #这里的思想和梯度下降法完全一致，只是此时，eta为正，即朝着函数值增大方向
        gradient = df(w,x)
        last_w = w
        w = w+eta*gradient
        w = direstion(w)
        if(abs(f(w,x)-f(last_w,x))<epsilon): #如果目标函数的值差距小于某个值，则表示差不多已经到达极值点，此时停止迭代
            break
    return w  

编写完核心代码后，现在就开始对样本x进行降维。

x_mean = demean(x)  #首先去均值处理
initial_w = np.random.random(x.shape[1])  #初始化不能为0向量，因为0方向的向量并不代表某个具体方向，相当于没初始化w，此外w的维度和x特征维度保持一致，由上面推导的公式可以得出该结论。
eta = 0.01 #初始化步长
w=gradietn_ascent(df_math,x_mean,initial_w,eta) #降维并得出第一个方向，即第一个主成分，使得在该方向上进行映射可以有最大的间距（方差）
plt.scatter(x_mean[:,0],x_mean[:,1]) #首先将原始数据绘制
plt.plot([-w[0]*60,w[0]*60],[-w[1]*60,w[1]*60],color = 'r') #将w表示的方向进行绘制，为了更清楚地展示，这里适当扩大w的表示范围，但是w0和w1比例要保持一致。 
plt.show()

这里求得的仅仅是1个主成分，在更高维度的空间，可能还需要继续求出其他的主成分，即降更多的维度。

三、求前n个主成分

1.求解思路

在上面，对于只有2维特征的样本x，我们求解出了第一个主成分，即w=（w1,w2），为了求解第二个主成分，这里需要将样本x在第一个主成分上的分量去除掉，这里使用的方法即空间几何的向量减法，得到的结果即下图中的绿线部分。

在得到绿线部分的分量后，再对该分量重新求第一主成分，以此类推，可以得到前n个主成分。

2.代码实现

首先给出求第一主成分的核心代码，这个和上面的核心代码完全一样，只是函数名字稍微有点变化。

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 direstion(w):
    return w /np.linalg.norm(w)  #将w进行归一化，因为w仅代表方向
def first_conponent(df,x,initial_w,eta,n_iters = 1e4,epsilon=1e-8):
    w = direstion(initial_w)
    cur_iter =0
    while cur_iter <n_iters:
        gradient = df(w,x)
        last_w = w
        w = w+eta*gradient
        w = direstion(w)
        if(abs(f(w,x)-f(last_w,x))<epsilon):
            break
    return w 

接下来求解第2个主成分，首先求去除掉第一个主成分后样本的值。

x_mean = demean(x)
initial_w = np.random.random(x.shape[1])    
eta = 0.01
w = first_conponent(df,x_mean,initial_w,eta) #这里求解出第一个主成分
x2 = np.empty(x.shape) #初始1个矩阵，和x的维度保持一致，用来存储去除掉第一个主成分分量后的样本
for i in range(len(x)): 
    x2[i] =x[i]-x[i].dot(w)*w #得到每个样本在去除掉第一个主成分分量后的值
plt.scatter(x2[:,0],x2[:,1])  #将该分量绘制出来
plt.show()

接下来，对剩下的分量求第一主成分。

w2 = first_conponet(df,x2,initial_w,eta)
w2.dot(w)  #验证，结果相乘接近0，因为第一主成分和第二主成分垂直，所以矩阵相乘为0

3.封装代码

将上面的代码全部进行封装。

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 direstion(w):
    return w /np.linalg.norm(w)  
def first_conponent(x,initial_w,eta,n_iters = 1e4,epsilon=1e-8):
    w = direstion(initial_w)
    cur_iter =0
    while cur_iter <n_iters:
        gradient = df(w,x)
        last_w = w
        w = w+eta*gradient
        w = direstion(w)
        if(abs(f(w,x)-f(last_w,x))<epsilon):
            break
    return w       
def first_n_componets(n,x,eta=0.01,n_iters =1e4,epsilon=1e-8):#这里将求前n个主成分进行封装
    x_pca = x.copy() #首先将原始样本拷贝给x_pca
    x_pca = demean(x_pca)
    res =[]          #用来存n个主成分
    for i in range(n):
        initial_w = np.random.random(x_pca.shape[1])
        w = first_conponent(x_pca,initial_w,eta)
        res.append(w)
        
        x_pca = x_pca-x_pca.dot(w).reshape(-1,1)*w #之前用的是for循环，实际上可以进行向量化，向量化后可以优化计算速度
    return res

使用代码，这个时候可以传入参数，指明需要求几个主成分，即保留多少个特征。

first_n_componets(2,x)

上面仅仅使用2个特征值的样本，接下来使用真实的数据集进行降维，并求原始样本向n个主成分降维的结果。

四、高维数据向低维数据转化

1.求解思路

现在有原始样本x，并求得了前k个主成分，形成wk矩阵，如下图所示。现在需要将x映射到新的轴上，形成只有k个特征的新样本xk，用该样本来描述原来的样本x，求解公式如下。

需要注意的是，映射后数据已经丢失了一部分信息，所以做逆运算时，并不能回到最初x的数据。

2.代码实现

在这里创建PCA.py文件（放在ml文件夹下），写入如下代码。

import numpy as np

class PCA:

    def __init__(self, n_components):
        """初始化PCA"""
        assert n_components >= 1, "n_components must be valid"
        self.n_components = n_components #传入需要保留多少
        self.components_ = None

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

    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) #这个数据集和上面使用的是一样的数据
from ml.PCA import PCA
pca = PCA(n_components=1)
pca.fit(x)
x_reduction = pca.transform(x)
x_restore = pca.inverse_transform(x_reduction) #即将一维的数据重新映射回二维，仅仅是展示用
plt.scatter(x[:,0],x[:,1],color='b',alpha=0.5)
plt.scatter(x_restore[:,0],x_restore[:,1],color='r',alpha=0.5)
plt.show() #可以发现重新映射回原来矩阵形式时，只剩下主成分方向上的数据，restore的功能只是在高维的空间表达低维的数据

3.sklearn中的PCA

from sklearn.decomposition import PCA
pca = PCA(n_components=1)
pca.fit(x)
pca.components_

可以发现求得的方向和之前的自己编写的相反。需要注意的是，sklearn降维并不是用梯度上升的计算方法，且其内部有一定的优化。自己编写的代码只是为了更好地理解PCA原理。

4.对手写识别数据进行降维

import numpy as np  
import matplotlib.pyplot as plt  
from sklearn import datasets
digits = datasets.load_digits()
x = digits.data
y = digits.target
from sklearn.model_selection import train_test_split #这里依然需要区分测试数据集和测试数据集，不理解的可以参考knn博客
x_train,x_test,y_train,y_test = train_test_split(x,y,random_state = 666)
from sklearn.neighbors import KNeighborsClassifier #因此目前为止，仅仅使用过knn算法进行分类，所以使用knn进行预测。
knn = KNeighborsClassifier()
knn.fit(x_train,y_train)
knn.score(x_test,y_test)

接下来使用降维的方式，再次用knn算法进行分类，并求分值。

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)#测试数据集也必须降维
knn = KNeighborsClassifier()
knn.fit(x_train_reduction,y_train)
knn.score(x_test_reduction,y_test)

可以发现，经过降维后，因为维度就剩下2个，得出的分值很低，为了提高分值，应该保留更多的特征值。但是具体保留多少特征合适？有什么衡量指标可以作为参考？

sklearn-PCA中有个变量为pca.explained_variance_ratio_，可以用来衡量映射后的数据保留了原始数据的多少信息（方差）。对于上面求得的结果，有如下的值。

这表示：经过映射后的数据，只保留原始数据大约百分之（14+13=27）的方差，丢掉了原来约73的方差，这个损失是比较大的，所以得出来的分数很低。

现在的问题是，应该保留多少的维度使得信息丢失在允许法范围内。假设选取的主成分的个数和原始样本特征维度一样，由于explained_variance_ratio_是按照大小顺序排列，这时候观察explained_variance_ratio_发现，越到后面，值越小，甚至几乎可以忽略。

pca = PCA(n_components=x_train.shape[1]) #保留所有的维度
pca.fit(x_train)
pca.explained_variance_ratio_

现在将这些值从前往后依次进行累加，并在坐标轴上展现。

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

观察该图可以发现，当横轴为0时，纵轴也为0（丢失了所有信息），当横轴为特征值大小时，纵轴为1（保留了全部信息）。因此，比如想要保留95%的信息，就可以通过该图知道应该保留多少个特征维度。

pca = PCA(0.95) #传入的参数为0.95，代表保留的百分之95的信息
pca.fit(x_train)
pca.n_components_  #结果为28

x_train_reduction = pca.transform(x_train)
x_test_reduction = pca .transform(x_test)
knn = KNeighborsClassifier()
knn.fit(x_train_reduction,y_train)
knn.score(x_test_reduction,y_test)

值明显得到了提高。而且如果进行时间测试，会发现时间缩短了很多。这就是pca的作用。

之前已经说过，pca经过降维可以用于可视化。将高维度难以想象的数据，降为2维可以在平面上展示。手写识别数据降为2维的图像如下。

pca = PCA(n_components=2)
pca.fit(x)
x_reduction = pca.transform(x)
for i in range(10):
    plt.scatter(x_reduction[y==i,0],x_reduction[y==i,1],alpha=0.8) #将每种数字进行统一绘制，并用不同颜色展示
plt.show() #可以发现即使是2维，也有较好的区分度

5.MNIST数据集降维

import numpy as np
from sklearn.datasets import fetch_mldata
mnist = fetch_mldata("MNIST original") #第一次加载会需要一些时间
x,y = mnist['data'],mnist['target']
x_train = np.array(x[:60000],dtype=float)
y_train = np.array(y[:60000],dtype=float)
x_test = np.array(x[:60000],dtype=float)
y_test = np.array(y[60000:],dtype=float)

不使用pca时，直接经过knn会需要很长的时间。

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

使用pca后发现，时间提高很多。奇怪的是，虽然特征减低了，但是分值却提高了，这就涉及到pca的另一个重要的作用-降噪，即在降维的过程中也将噪声减小了。

from sklearn.decomposition import PCA
pca = PCA(0.9)
pca.fit(x_train)
x_train_reduction = pca.transform(x_train)
x_test_reduction = pca.transform(x_test)    #测试数据集也必须降维
knme knn.fit(x_train_reduction,y_train)
%time knn.score(x_test_reduction,y_test)

6.降噪

首先看之前的一个数据。

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,5.,size=100)

from sklearn.decomposition import PCA
pca = PCA(n_components=1)
pca.fit(x)
x_reduction = pca.transform(x)
x_restore = pca.inverse_transform(x_reduction)
plt.scatter(x[:,0],x[:,1],color='r')
plt.scatter(x_restore[:,0],x_restore[:,1])
plt.show() 

原来的数据可能是因为测量的过程中出现了问题，所以发生了一定的抖动（即产生了噪声），经过降维后可以发现最终的数据在一条直线上，也就是将噪声进行了去除，因此说pca有一定的降噪能力。下面使用手写识别例子继续探讨降维中的去燥功能。

from sklearn import datasets

digits = datasets.load_digits()
x = digits.data
y = digits.target
noisy_digits = x + np.random.normal(0,4,size=x.shape) #这里加入随机噪声

example_digits = noisy_digits[y==0,:][:10] 
for num in range(1,10):
    x_num = noisy_digits[y==num,:][:10]
    example_digits = np.vstack([example_digits,x_num])#将全部样本按照0-10的顺序叠起来，共100个样本
    
def plot_digits(data): #这个函数当做模块使用即可，不必深究
    fig,axes = plt.subplots(10,10,figsize=(10,10),subplot_kw={'xticks':[],'yticks':[]},gridspec_kw=dict(hspace=0.1,wspace=0.1))
    for i,ax in enumerate(axes.flat):
        ax.imshow(data[i].reshape(8,8),cmap='binary',interpolation='nearest',clim=(0,16))
    plt.show()
plot_digits(example_digits) #调用函数

可以看到图像十分模糊，现在通过pca进行降噪处理。

pca = PCA(0.5) #因为本身噪声比较大，所以这里仅保留原始数据0.5的信息
pca.fit(noisy_digits)
pca.n_components_ #使用了12个特征

components = pca.transform(example_digits)
filtered_digits = pca.inverse_transform(components)
plot_digits(filtered_digits)

发现图像的清晰度有所提高。