数据降维之PCA方法与python实现

基本数学知识

• 向量内积：通过内积可以判断两个数据信息之间的相关性
  设两个向量：$\overrightarrow{a_1}=\left(x_1,y_1\right),\overrightarrow{a_2}=\left(x_2,y_2\right)$，夹角为$\theta$
  有$\vec{a}\cdot \vec{b}=|\vec{a}|\ast |\vec{b}|\ast \cos \theta$
• 矩阵转置：矩阵A的行和列互相交换所产生的矩阵称为A的转置矩阵
  $A=\left[\begin{array}{ll}{a}_{11}& a_{12}\\ a_{21}& a_{22}\end{array}\right]$ ，$A^T=\left[\begin{array}{ll}{a}_{11}& a_{21}\\ a_{12}& a_{22}\end{array}\right]$
• 对称矩阵：$A=A^T$
• 矩阵的逆：设A是数域上的一个n阶方阵，若在相同的数域上存在另一个n阶方阵B，使得AB=BA=E，那么称B为A的逆矩阵，记作$A^{-1}$，而A被称为可逆矩阵或非奇异矩阵。如果A不存在逆矩阵，那么A称为奇异矩阵。
• 特征值、特征向量：定义为：$Ax=\lambda x$，$x$为特征向量，$\lambda$为特征值
  $|\lambda E-A|=0$、$(\lambda E-A)x=0$
• 特征值分解：矩阵A的n个特征值${\lambda }_1\le {\lambda }_2\le ...\le {\lambda }_n$,以及这n个特征值所对应的特征向量 $w_1,w_2,...w_n$，如果这n个特征向量线性无关，那么矩阵A就可以用下式的特征分解表示：$A=W\Sigma W^{-1}$
  W是这n个特征向量所张成的n×n维矩阵，并对n个特征向量标准化，而Σ为这n个特征值为主对角线的n×n维矩阵。若A为实对称矩阵，另有
  $A=W\Sigma W^T$

SVD

• 特征值分解只可以在方阵中进行处理，不能应用的非方阵中 .特征值分解可以利用到降维的处理中，提升使用效率,但是大多数矩阵都是非方阵形式，所以不可能进行特征值分解：$A=W\Sigma W^T$,这里应运而生SVD（奇异值）分解 ,可以解决非方阵使用特征值分解的方法，这种方法被称作奇异值分解。
• 奇异值分解（SVD, singular value decomposition），一种对矩阵较好的分解形式

算法理论

• 矩阵的奇异值分解是指，将一个非零的m×n实矩阵A，$A\in R_{m\times n}$，表示为以下三个实矩阵乘积形式 的运算，即进行矩阵的因子分解：$A=U\Sigma V^T$
• 其中U是m阶正交矩阵，V是n阶正交矩阵，Σ是由降序排列的非负的对角元素组成的m×n矩形对角 矩阵，满足：
  $$\begin{array}{l}U{U}^T=I\\ V{V}^T=I\\ \Sigma =\mathrm{diag}\left({\sigma }_1,{\sigma }_2,\cdots ,{\sigma }_p\right)\\ {\sigma }_1\ge {\sigma }_2\ge \cdots \ge {\sigma }_p\ge 0\\ p=\min (m,n)\end{array}$$
• $U\Sigma V^T$ 称为矩阵A的奇异值分解，${\sigma }_i$称为矩阵A的奇异值，U的列向量称为左奇异向量，V的列向量称 为右奇异向量。
• 注意奇异值分解不要求矩阵A是方阵。
• 正常求上面的U,V,Σ不便于求，我们可以利用如下性质：
  $A{A}^T=U\Sigma V^{T}V{\Sigma }^T{U}^T=U\Sigma {\Sigma }^T{U}^T$
  $A^{T}A=V{\Sigma }^T{U}^{T}U\Sigma V^T=V{\Sigma }^T\Sigma V^T$
• 这里$\Sigma {\Sigma }^T$与${\Sigma }^T\Sigma$在矩阵的角度上来讲，它们是不相等的，因为它们的维数不同$\Sigma {\Sigma }^T$∈R_{m×m}，而 ${\Sigma }^T\Sigma$∈R_{n×n}，但是它们在主对角线的奇异值是相等的
  $$\Sigma {\Sigma }^T={\left(\begin{array}{cccc}{\sigma }_1^2& 0& 0& 0\\ 0& {\sigma }_2^2& 0& ⋮\\ 0& 0& \ddots & 0\\ 0& \dots & 0& {\sigma }_m^2\end{array}\right)}_{m\times m}$$
  $${\Sigma }^T\Sigma ={\left(\begin{array}{cccc}{\sigma }_1^2& 0& 0& 0\\ 0& {\sigma }_2^2& 0& ⋮\\ 0& 0& \ddots & 0\\ 0& \dots & 0& {\sigma }_n^2\end{array}\right)}_{n\times n}$$
• 奇异值分解是一种矩阵因子分解方法，是线性代数概念，但在统计学习中被广泛使用，成为其重要 工具。
• 应用：主成分分析、潜在语义分析
• 矩阵的奇异值分解一定存在，但不唯一。奇异值分解可以看做矩阵数据压缩的一种方法。

SVD小例子

实现一张图的(R,G,B)通道简单分解，相当于打码

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib as mpl
from PIL import Image
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import Binarizer
from sklearn.decomposition import PCA,TruncatedSVD

class SVD_():
    def __init__(self):
        pass

    def restore1(self,u,sigma,v,k):
        #refactoring image
        a = np.dot(u[:,:k],np.diag(sigma[:k])).dot(v[:k,:])
        a[a<0] = 0
        a[a>255] = 255
        return np.rint(a).astype('uint8')

    def main_(self):
        A = Image.open(r'data/godness.jpg')
        a = np.array(A)
        u_r,sigma_r,v_r = np.linalg.svd(a[:,:,0])
        u_g,sigma_g,v_g = np.linalg.svd(a[:,:,1])
        u_b,sigma_b,v_b = np.linalg.svd(a[:,:,2])
        plt.figure(facecolor='w',figsize=(10,10))
        #Number of singular values in sequence from 1 to 20
        K = 20
        for k in range(1,K+1):
            R = self.restore1(u_r,sigma_r,v_r,k)
            G = self.restore1(u_g,sigma_g,v_g,k)
            B = self.restore1(u_b,sigma_b,v_b,k)
            I = np.stack((R,G,B),axis=2)
            plt.subplot(5,4,k)
            plt.imshow(I)
            plt.axis('off')
            plt.title('奇异值个数：%d'%k)
        plt.suptitle('SVD与图像分解',fontsize=20)
        plt.tight_layout(0.1,rect=(0,0,1,0.92))
        plt.show()

if __name__=='__main__':
	svd = SVD_()
	svd.main_()

PCA

• PCA是一种无监督学习的降维技术
• 信噪比：信噪比就是信号与噪声的方差比，在信号处理中认为信号具有较大方差，而噪声的方差应尽量小，因此，信噪比越大越好
• 做法：将所有的样本点向一条直线进行投影；让投影后样本的方差达到最大，即让投影后的样本尽量分散

算法理论

计算过程

这里先举个小例子理解PCA的计算过程：

x	2.5	0.5	2.2	1.9	3.1	2.3	2	1	1.5	1.1
y	2.4	0.7	2.9	2.2	3.0	2.7	1.6	1.1	1.6	0.9

• 第一步分别求$x$、$y$的平均值，然后对于所有样例都减去对应均值（这里只减去了均值，还应该对数据除以方差，进行归一化处理），得到：

x	0.69	-1.31	0.39	0.09	1.29	0.49	0.19	-0.81	-0.31	-0.71
y	0.49	-1.21	0.99	0.29	1.09	0.79	-0.31	-0.81	-0.31	-1.01

• 第二步，求特征协方差矩阵：
  $$cov=\left(\begin{array}{cc}cov(x,x)& cov(x,y)\\ cov(y,x)& cov(y,y)\end{array}\right)=\left(\begin{array}{cc}0.616555556& 0.615444444\\ 0.615444444& 0.716555556\end{array}\right)$$
• 第三步，对协方差矩阵求特征值与特征向量：
  $$eig\_values=\left(\begin{array}{c}0.0491\\ 1.284\end{array}\right)$$
  $$eig\_vectors=\left(\begin{array}{cc}-0.7352& -0.67787\\ 0.67789& -0.7352\end{array}\right)$$
• 第四步，将特征值从大到小顺序排序，选择其中最大的$k$个，然后将其$k$个特征向量分别作为列向量组成特征向量矩阵。
  这里的特征只有两个，选取最大的那个，即1.284，对应特征向量为$(-0.67787,-0.7352{)}^T$
• 第五步，将样本点投影到选取的特征向量上，假设样本个数为$m$，特征数为$n$，减去均值后的样本矩阵为$W_{m\times n}$，协方差矩阵$co{v}_{n\times n}$，选取的$k$个特征向量组成的矩阵为$E{V}_{n\times k}$，则投影后的数据为：
  $$dat{a}_{m\times k}=W_{m\times n}\cdot E{V}_{n\times k}$$
  所以，这里投影后的数据计算是：
  $$dat{a}_{10\times 1}=W_{10\times 2}\cdot E{V}_{2\times 1}$$
  计算得到：
  $data=(-0.82797,1.77758,-0.992197,\dots ,0.438046,1.22382056)$

这样，原始数据就从2维变成了1维。

理论基础

投影

红色点表示样例$x^{(i)}$，蓝色点表示$x^{(i)}$在$u$上的投影，$u$ 是直线的斜率也是直线的方向向量，而且是单位向量。蓝色点是$x^{(i)}$在 $u$ 上的投影点，离原点的距离是$<x^{(i)},u>$(即$x^{(i{)}^T}u$或者$u^T{x}^{(i)})$，由于这些样本点（样例）的每一维特征均值都为 0，因此投影到 u 上的样本点（只有一个到原点的距离值）的均值仍然是 0。

最大方差理论

PCA的思想是要让投影后的点方差越大越好

• 若样本没有进行中心化，$$\underset{u}{\max }\sum _i{\left[u^T\left(x_i-\bar{x}\right)\right]}^2$$
• 若样本进行中心化后，$$\underset{u}{\max }\sum _i{\left[u^T\left(x_i\right)\right]}^2$$

这里以中心化后的式子进行推导：
$$\frac{1}{m}\sum _{i=1}^m{\left(x^{(i)T}u\right)}^2=\frac{1}{m}\sum _{i=1}^m{u}^T{x}^{(i)}{x}^{(i{)}^T}u=u^T\left(\frac{1}{m}\sum _{i=1}^m{x}^{(i)}{x}^{(i)T}\right)u$$
括号中间部分是协方差矩阵，这里用$\lambda$表示$\frac{1}{m}{\sum }_{i=1}^m{\left(x^{(i)T}u\right)}^2$，用$\Sigma$表示$\left(\frac{1}{m}{\sum }_{i=1}^m{x}^{(i)}{x}^{(i)T}\right)$，得到：$$\lambda =u^T\Sigma u$$
由于$u$是单位向量，则有$u^{T}u=1$，$u\lambda =\lambda u=u{u}^T\Sigma u=\Sigma u$，即$\Sigma u=\lambda u$

• 这样就得到$\lambda$是$\Sigma$的特征值，$u$是特征向量。因此，我们只需要对协方差矩阵进行特征值分解，得到的前 k 大特征值对应的特征向量就是最佳的 k 维新特征，而且这 k 维新特征是正交的。得到前 k 个 u 以后，样例通过变换可以得到新的样本。

PCA实现

import numpy as np
import pandas as pd
from sklearn.decomposition import PCA

def pca(X,k):
    m_samples,n_feayures = X.shape
    #Centralization De-mean Mean for zero
    mean = np.mean(X,axis=0)
    normX = X - mean
    #transpose a matrix
    cov_mat = np.dot(np.transpose(normX),normX)
    vals,vecs = np.linalg.eig(cov_mat)
    eig_pairs = [(np.abs(vals[i]),vecs[i]) for i in range(n_feayures)]
    #Eigenvalues from large to small
    eig_pairs = eig_pairs.sort(reverse=True)
    #Get the number of dimensions
    feature = np.array(eig_pairs[0][k])
    #Data restoration
    data = np.dot(normX,np.transpose(feature))
    #sklearn pca
    p = PCA(n_components=k)
    a = p.fit_transform(X)
    #sklearn的pca使用
    pca_sklearn = PCA(n_components=K)
    x_ = pca.fit_transform(X)
  
if __name__=='__main__':
    X = np.array([[-1,1],[-2,-1],[-3,-2],[1,1],[2,1],[3,2]])
    pca(X,1)