主成分分析PCA+编程实现

主成分分析

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

对数据进行简化的一系列原因：

• 使得数据集更易使用
• 降低很多算法的计算开销
• 去除噪声
• 使得结果易懂

什么是主成分分析？

在PCA中，数据从原来的坐标系转换到了新的坐标系，新的坐标系的选择是由数据本身决定的。第一个新坐标轴选择的是原始数据中方差最大的方向，第二个新坐标轴的选择和第一个坐标轴正交且具有最大方差的方向。该过程一直重复，重复次数为原始数据中特征的数目，我们会发现，大部分方差都包含在最前面的几个新坐标轴中。因此，可以忽略余下的坐标轴，及对数据进行了降维处理。

PCA优缺点

• 优点：降低数据的复杂性，识别最重要的多个特征
• 缺点：不一定需要，且可能损失有用信息
• 使用数据类型：数值型数据

举个例子理解一下，二维空间上有一些样本点如下，将二维降维至一维：

原始样本点（二维）	映射到特征1上	映射到特征2上

从图中可以看到，将样本点映射到特征1这种方案感觉会好一些，因为相较另一种方案，样本与样本之间的间距较大一些，也就是说点和点之间的区分度更大一些。那有没有更好地映射方案呢？来看看下面这中映射成一维的方案：

原始样本点（二维）	找一根直线	将所有的点都映射到这根直线上

在这种映射方案下，我们可以看到它将所有的点都映射到直线上，虽然线是斜的，我们也可以把它看成是一个轴，一个维度，可以看到所有的点在映射之后，可以发现它们与原来的样本点的分布并没有更大的差距，并且样本之间的距离（可区分度）也比无论映射到x轴还是y轴都要大。

那么我们需要思考的就是如何找到这个让样本间间距最大的轴？如何定义样本间的间距？通常情况下，我们使用方差（Variance），因为方差可以作为反映样本的疏密的一个指标。

这样上面的降维问题，就转化成为了找到一个轴，使得样本空间的所有点映射到这个轴后，方差最大的问题。

为了解决这个问题，在进行主成分分析之前，要做一些工作：

第一步：将样例的均值归为0（demean），所有样本减去整体的均值，可以看到处理之后，样本分布没有改变，只是坐标轴变了。

原始样本分布	demean之后

demean之后样本均值为0，方差公式就可以做出一些改变，更加方便

第二步，求这个轴的方向w=(w1,w2)(两维的表达)，使得我们所有的样本映射到w维度上方差最大，即

如何求	图解
w为要映射的轴的方向，蓝色的点为原始的样本点，映射到w上，其实就是蓝色的那跟线段，其实也就是点乘的定义

由这个问题，扩展到N维,主成分分析的目标就是求w,使得映射到此维度后，方差最大

插曲，光看图我们会感觉到pca和线性回归有些相似，但需要注意他们之间的差别

PCA	线性回归
横轴轴是特征 推导式子不同 虽然都是样本与一条直线的关系，但pca是距离时点垂直于直线	横轴特征，纵轴输出标记 推导式子不同 点垂直于特征

梯度上升法解决主成分分析问题

1.生成数据，进行demean

import numpy as np
import matplotlib.pyplot as plt
X=np.empty(shape=(100,2))
#从一个均匀分布[low,high)中随机采样，注意定义域是左闭右开，即包含low，不包含high.
X[:,0]=np.random.uniform(1,100,size=100)
X[:,1]=0.75*X[:,0]+3+np.random.normal(0,10,size=100)
plt.scatter(X[:,0],X[:,1])
plt.title('X')
plt.show()

def demean(X):
    #求每列的均值，也是每个特征维度上的均值
    return X-np.mean(X,axis=0)
X_demean=demean(X)
plt.scatter(X_demean[:,0],X_demean[:,1])
plt.title('demean(X)')
plt.show()

2.梯度上升法

def f(w,X):
    '''
    目标函数
    :param w:映射向量
    :param X: 样本数据
    :return:
    '''
    return np.sum(np.dot(X,w)**2)/len(X)

def df(w,X):
    '''
    目标函数求梯度
    :param w:
    :param X:
    :return:
    '''
    return 2/len(X)*np.dot(X.T,np.dot(X,w))

def direction(w):
    '''
    将向量w转化为单位向量
    :param w:
    :return:
    '''
    return w/np.linalg.norm(w)

def gradient_ascent(df,X,init_w,eta,max_iters=1e4,epsilon=1e-8):
    '''
    梯度上升
    :param df: 
    :param X: 
    :param init_w: 
    :param eta: 步长学习率
    :param max_iters: 最大迭代步数
    :param epsilon: 精度
    :return: 
    '''
    cur_iter=0
    w=direction(init_w)
    while cur_iter<max_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

init_w=np.random.random(X.shape[1])
eta=0.01
w=gradient_ascent(df,X,init_w,eta)

plt.scatter(X_demean[:,0],X_demean[:,1])
plt.plot([-50,w[0]*80],[-50,w[1]*80],c='r')
plt.title('demean(X)')
plt.show()

求出第一主成分以后，如何求出下一个主成分？

将数据进行改变，将数据在第一个主成分上的分量去掉，即下图中的X',然后在新的数据上求第一主成分

# X2=np.empty(X.shape)
# for i in range(len(X)):
#     X2[i]=X[i]-np.dot(X[i],w)*w
X2=X-np.dot(X,w).reshape(-1,1)*w
plt.scatter(X2[:,0],X2[:,1])
plt.title('X2')
plt.show()
w2=gradient_ascent(df,X2,init_w,eta)
print(w2)
#垂直应该等于0
print(w.dot(w2))

将求前n个主成分，封装成函数

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

print(first_n_components(2,X))

[array([ 0.75823566,  0.65198059]), array([ 0.65198379, -0.7582329 ])]

高维数据向低维数据映射

样本数据（m行n列）	X的前k个主成分（k行，n列）

将高维数据向低维数据映射，其实就是求

反之，从低维映射到高维也是可以的

import numpy as np
class PCA:
    def __init__(self,n_components):
        '''
        :param n_components: 指定几个主成分
        '''
        self.n_components=n_components
        self.components_=None

    def __repr__(self):
        return 'pca%d个主成分'%self.n_components

    def fit(self,X,eta=0.01,n_iter=1e4):
        '''根据训练数据集X获得数据的均值和方差'''
        def demean(X):
            # 求每列的均值，也是每个特征维度上的均值
            return X - np.mean(X, axis=0)

        def f(w, X):
            '''
            目标函数
            :param w:映射向量
            :param X: 样本数据
            :return:
            '''
            return np.sum(X.dot(w) ** 2) / len(X)

        def df(w, X):
            '''
            目标函数求梯度
            :param w:
            :param X:
            :return:
            '''
            return 2 / len(X) * np.dot(X.T, np.dot(X, w))

        def direction(w):
            '''
            将向量w转化为单位向量
            :param w:
            :return:
            '''
            return w / np.linalg.norm(w)

        def gradient_ascent(df, X, init_w, eta=eta, max_iters=n_iter, epsilon=1e-8):
            '''
            梯度上升
            :param df:
            :param X:
            :param init_w:
            :param eta: 步长学习率
            :param max_iters: 最大迭代步数
            :param epsilon: 精度
            :return:
            '''
            w = direction(init_w)
            cur_iter = 0

            while cur_iter < max_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):
            init_w=np.random.random(x_pca.shape[1])
            w=gradient_ascent(df,x_pca,init_w)
            self.components_[i,:]=w
            x_pca=x_pca-np.dot(x_pca,w).reshape(-1,1)*w
        return self

    def transform(self,X):
        '''将给定的X,映射到各个主成分分量中'''
        return np.dot(X,np.transpose(self.components_))

    def inverse_transform(self,X):
        '''将给定的X，反向映射回原来的特征空间'''
        return np.dot(X,self.components_)

import matplotlib.pyplot as plt
X=np.empty(shape=(100,2))
#从一个均匀分布[low,high)中随机采样，注意定义域是左闭右开，即包含low，不包含high.
X[:,0]=np.random.uniform(1,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)
print(np.transpose(pca.components_))
#降维
pca=PCA(n_components=1)
pca.fit(X)
x_reduction=pca.transform(X)
print(x_reduction.shape)

#恢复原始维度
x_restore=pca.inverse_transform(x_reduction)
print(x_restore.shape)
plt.scatter(X[:,0],X[:,1],c='b')
plt.scatter(x_restore[:,0],x_restore[:,1],c='g')
plt.show()

SKlearn中的PCA

import numpy as np
import matplotlib.pyplot as plt
X=np.empty(shape=(100,2))
#从一个均匀分布[low,high)中随机采样，注意定义域是左闭右开，即包含low，不包含high.
X[:,0]=np.random.uniform(1,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)
print(pca)
print(pca.components_)
x_reduction=pca.transform(X)
print(x_reduction.shape)
x_restore=pca.inverse_transform(x_reduction)
plt.scatter(X[:,0],X[:,1],c='b')
plt.scatter(x_restore[:,0],x_restore[:,1],c='g')
plt.show()

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
x_train,x_test,y_train,y_test=train_test_split(x,y,random_state=666)
print(x_train.shape)

from sklearn.neighbors import KNeighborsClassifier
knn_clf=KNeighborsClassifier()
knn_clf.fit(x_train,y_train)
print(knn_clf.score(x_test,y_test))

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_clf=KNeighborsClassifier()
knn_clf.fit(x_train_reduction,y_train)
print(knn_clf.score(x_test_reduction,y_test))
#解释各轴方差相应的比例
print(pca.explained_variance_ratio_)

pca=PCA(n_components=64)
pca.fit(x_train)
print(pca.explained_variance_ratio_)
plt.plot([i for i in range(64)],[np.sum(pca.explained_variance_ratio_[:i+1]) for i in range(64)])
plt.show()

pca=PCA(0.95)
pca.fit(x_train)
print(pca.n_components_)
x_train_reduction=pca.transform(x_train)
x_test_reduction=pca.transform(x_test)
knn_clf=KNeighborsClassifier()
knn_clf.fit(x_train_reduction,y_train)
print(knn_clf.score(x_test_reduction,y_test))

#降到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.5)
plt.show()

PCA+MNIST

import numpy as np
from sklearn.datasets import fetch_mldata
mnist=fetch_mldata("MNIST original",data_home='./mnist')
print(mnist)
x=mnist['data']
y=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)

from sklearn.neighbors import KNeighborsClassifier
knn_clf=KNeighborsClassifier()
knn_clf.fit(x_train,y_train)
print(knn_clf.score(x_test,y_test))


#pca降维
from sklearn.decomposition import PCA
pca=PCA(0.9)
pca.fit(x_train)
x_train_reduction=pca.transform(x_train)
knn_clf.fit(x_train_reduction,y_train)
x_test_reduction=pca.transform(x_test)
print(knn_clf.score(x_test_reduction,y_test))

{'COL_NAMES': ['label', 'data'], 'data': array([[0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       ..., 
       [0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0]], dtype=uint8), 'target': array([ 0.,  0.,  0., ...,  9.,  9.,  9.]), 'DESCR': 'mldata.org dataset: mnist-original'}
0.9688
0.9728

可以看到在降维之后，不仅节省了计算开销，正确率也上升了，一箭双雕。

使用PCA对数据进行降噪

import numpy as np
from sklearn import datasets
import matplotlib.pyplot as plt
digits=datasets.load_digits()
x=digits.data
y=digits.target
#生成噪声数据
noise_digits=x+np.random.normal(0,4,size=x.shape)
#y=0的数字10个
example_digits=noise_digits[y==0,:][:10]
print(example_digits.shape)
#将另外的1-9标签各类图片10张，压入example_digits中
for num in range(1,10):
    x_num=noise_digits[y==num,:][:10]
    example_digits=np.vstack([example_digits,x_num])
print(example_digits.shape)
#画图
def plot_digits(data,title):
    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.title(title)
    plt.show()
plot_digits(example_digits,'降噪前')
#因为噪声较多，只选取50%的信息
from sklearn.decomposition import PCA
pca=PCA(0.5)
pca.fit(noise_digits)
print(pca.n_components_)
components=pca.transform(example_digits)
filtered_digits=pca.inverse_transform(components)
plot_digits(filtered_digits,'降噪后')

降噪前	降噪后

特征脸

import numpy as np
from sklearn.datasets import fetch_lfw_people
import matplotlib.pyplot as plt
faces=fetch_lfw_people()


#对所有的样本进行一个随机的排列
random_indexes=np.random.permutation(len(faces.data))
x=faces.data[random_indexes]

example_faces=x[:36,:]
print(example_faces.shape)
def plot_faces(data):
    fig,axes=plt.subplots(6,6,figsize=(10,10),
                         subplot_kw={'xticks':[],'yticks':[]},gridspec_kw=dict(hspace=0.1,wspace=0.1))
    for i,ax in enumerate(axes.flat):
        #图片w*h=62,47
        ax.imshow(data[i].reshape(62,47))
    plt.show()
plot_faces(example_faces)

from sklearn.decomposition import PCA
pca=PCA(svd_solver='randomized')
pca.fit(x)
plot_faces(pca.components_[:36,:])

降维前	降维后