PCA简化数据

1. 降维

在已标注与未标注的数据上都有降维技术。这里我们将主要关注未标注数据上的降维技术，该技术同时也可以应用于已标注的数据

• 主成分分析（Principal Component Analysis，PCA）：
  将数据从原来的坐标系转移到新的坐标系，新坐标系的选择由数据本身决定，新坐标系的第一个坐标轴是原始数据中方差最大的方向，新坐标系的第二个坐标轴和第一个坐标轴正交、并且具有最大方差。该过程一直重复，次数为原始数据中维度。大部分方差都包含在前面几个新坐标轴中，因此可以忽略剩下的坐标轴，即对数据进行了降维处理。
• 因子分析（Factor Analysis）：
  假设观察数据的生成中有一些观察不到的隐变量，即观察数据是由这些隐变量和某些噪声的线性组合，那么隐变量的数据可能比观察数据的数目少，找到隐变量就可以实现数据的降维。
• 独立成分分析（Independent Component Analysis，ICA）：
  假设数据从N个数据源生成，类似因子分析，假设这些数据源之间在统计上相互独立，如果数据源数目少于观察数据数目，就实现降维过程。
  -独立成分分析（Independent Component Analysis，ICA）：
  ICA假设数据是从N个数据源生成的，这一点和因子分析有些类似。假设数据为多个数据源的混合观察结果，这些数据源之间在统计上是相互独立的，而在PCA中只假设数据是不相关的。同因子分析一样，如果数据源的数目少于观察数据的数目，则可以实现降维过程。

在上述3种降维技术中，PCA的应用目前最为广泛。

2. PCA

2.1 理论知识

在正交属性空间中的样本点，可以用一个超平面对所有的样本进行一个合适的表达，而这个超平面应该具有如下的性质：

• 最近重构性：样本点到这个超平面的距离足够近
• 最大可分性：样本点到这个超平面的投影尽可能分开

基于最近重构性和最大可分性，能分别得到主成分分析的两种推到：
（1）最近重构性

• 首先对数据进行中心化，即$\sum_ix_i=0$
• 再假定投影变换后得到的坐标为${w_1,w_2,......w_d}$，其中$w_i$是标准正交基向量，即$||w_i||_2=1, w_i^Tw_j=0(i\neq j)$
• 这时如果丢弃新坐标系中的部分坐标，即把维度降低到${d}'<d$，就达到了降维的目的，此时样本在低维中的投影是$z_i=(z_{i1},z_{i2}.......z_{id'})$，其中$z_{ij}=w_j^Tx_i$是$x_i$在低维坐标系下第$j$维的坐标，现在基于投影的点重构原样本点得到：$\hat{{x_i}}=\sum_{j=1}^{d'}z_{ij}w_{j\cdot}$

则此时对整个训练集要计算原样本点$x_i$与基于投影得到的样本点$\hat{{x_i}}$之间的距离：

$$\sum_{i=1}^m\left \| \sum_{j=1}^{d}'z_{ij}w_j-x_i \right \|_2^2=\sum_{i=1}^mz_i^Tz_i-2\sum_{i=1}^mz_i^TW^Tx_i+const$$
$$=>-tr\left ( W^T\left (\sum_{i=1}^mx_ix_i^T \right ) W\right )$$
所以根据最近重构性，上面的式子应被最小化，同时 $w_j$是标准正交基， $\sum_{i=1}^mx_ix_i^T$是协方差矩阵，有
$$\begin{matrix} \min_w & -tr(W^TXX^TW)\\ s.t &W^TW=1 \end{matrix}$$
这就是主成分分析的优化目标。

（2）最大可分性

我们知道样本点$x_i$在新空间中超平面的投影是$W^Tx_i$，此时，要想让所有的样本的投影尽可能的分开，就要使投影后样本点的方差最大化，如下图：

投影后样本点的方差是$\sum_iW^Tx_ix_i^TW$,于是优化目标为：

$$\begin{matrix} \max_w & tr(W^TXX^TW)\\ s.t &W^TW=1 \end{matrix}$$
可以看出两个优化目标都要使用拉格朗日乘子法可得：
$$XX^TW=\lambda W$$
其中 $\lambda$是特征值，于是可知只需对协方差矩阵 $XX^T$进行特征值分解，将求得的特征值排序，然后取前n个特征值对应的特征向量构成 $W=(w_1,w_2,.......w_n)$，这就是主成分分析的解。

注意：
PCA也可看作是逐一选取方差最大的方向，即先对协方差矩阵$\sum _i x_ix_i^T$做特征值分解，取最大特征值对应的特征向量$w_1$，再对$\sum _i x_ix_i^T-\lambda_1w_1w_1^T$做特征分解，取最大特征值对应的特征向量$w_2$，。。。理论可以证明上述逐一选取方差最大方向的做法和直接选取最大n个特征值等价。

2.2 移动坐标轴

• 面对二维图上的大量数据点，如果要求我们画出一条直线，这条线要尽可能覆盖这些点，那么最长的线可能是哪条？那肯定是B这条直线最长，在PCA中，我们对数据的坐标进行了旋转，该旋转的过程取决于数据的本身。B坐标轴旋转到覆盖数据的最大的方差位置，数据的最大方差给出了数据的最重要的信息。在选择了覆盖数据最大差异性的坐标轴之后，我们选择了第二条坐标轴。也就是正交的那条坐标轴。当然，在二维平面下，垂直和正交是一回事。利用PCA我们将数据坐标轴旋转至数据角度上的那些最重要的方向。
• 如果实现了坐标轴的旋转，那么接下来就是降维。其实坐标轴的旋转并没有减少数据的维度。通过PCA进行降维处理，我们就可以同时获得SVM和决策树的优点：一方面，得到了和决策树一样简单的分类器，同时分类间隔和SVM—样好。

示例：

上图中的数据来自于未经PCA处理的数据经PCA转换之后绘制而成的，如果仅使用原始数据，那么这里的间隔会比决策树的间隔更大。另外，由于只需要考虑一维信息，因此数据就可以通过比SVM简单得多的很容易采用的规则进行区分，这其实就是PCA的主要作用。

前面提到的第一个主成分就是从数据差异性最大（即方差最大）的方向提取出来的，第二个主成分则来自于数据差异性次大的方向，并且该方向与第一个主成分方向正交。通过数据集的协方差矩阵及其特征值分析，我们就可以求得这些主成分的值。一旦得到了协方差矩阵的特征向量，我们就可以保留最大的外个值。这些特征向量也给出了N个最重要特征的真实结构。我们可以通过将数据乘上这N个特征向量而将它转换到新的空间。

2.3 PCA

将数据转换成前N个主成分的伪码大致如下：

去除平均值
计算协方差矩阵
计算协方差矩阵的特征值和特征向量
将特征值从大到小排序
保留最上面的N个特征向量
将数据转换到上述N个特征向量构建的新空间中

提前的数据处理：中心化和标准化

PCA代码：

# -*- coding: utf-8 -*-
"""
Spyder Editor
"""
from numpy import *
import matplotlib
import matplotlib.pyplot as plt

def loadDataSet(fileName, delim='\t'):
    fr = open(fileName)
    stringArr = [line.strip().split(delim) for line in fr.readlines()]
    datArr = [map(float,line) for line in stringArr] 
    return mat(datArr) # 返回矩阵

# PCA算法
def pca(dataMat, topNfeat=9999999):  # topNfeat是应用的99999999个特征
    meanVals = mean(dataMat, axis=0) # 求平均值(按每一列)
    meanRemoved = dataMat - meanVals   
    covMat = cov(meanRemoved, rowvar=0)  # 计算协方差矩阵
    eigVals,eigVects = linalg.eig(mat(covMat)) # 得到特征值和特征向量
    #print 'eigVals,eigVects:',eigVals,eigVects
    eigValInd = argsort(eigVals)       # 特征值从小到大排序时的下标值
    eigValInd = eigValInd[:-(topNfeat+1):-1]  # 得到下标值的排序，倒着排
    #print type(eigValInd)
    redEigVects = eigVects[:,eigValInd]  # 得到topNfeat个最大的特征向量，构成对数据进行转换的矩阵
    lowDDataMat = meanRemoved * redEigVects # 该矩阵则利用N个特征将原始数据转换到新空间中,低维度的矩阵
    reconMat = (lowDDataMat * redEigVects.T) + meanVals # 降维后的矩阵
    return lowDDataMat, reconMat  # 返回低维矩阵和降维后的矩阵


# 绘图
def plt_dat(dataMat,reconMat):
    fig = plt.figure()
    ax = fig.add_subplot(111)
    ax.scatter(dataMat[:,0].flatten().A[0], dataMat[:,1].flatten().A[0], marker='^', s=90)
    ax.scatter(reconMat[:,0].flatten().A[0], reconMat[:,1].flatten().A[0], marker='o', s=20,c='red')
    plt.show()

# 主函数

dataMat=loadDataSet('testSet.txt')
lowDMat,reconMat=pca(dataMat, 1)
print shape(lowDMat)
plt_dat(dataMat,reconMat)
lowDMat,reconMat=pca(dataMat, 2)
plt_dat(dataMat,reconMat)

用到的基本语法：
map（）
mean（）
cov（）
flatten（）
numpy矩阵运算

运行结果：

可以看出lowDMat包含了降维之后的矩阵，这里是个一维矩阵。这里图示1是降维后的数据和原始数据的对比。图示2是改变参数后的重构的数据和原始数据的对比，可以看出已经几乎完全重合了。

注意的几点：

• 实践中通过对 X <script type="math/tex" id="MathJax-Element-30">X</script>进行奇异值分解来代替协方差局针的特征值分解。
• 低维空间的维数通常由用户指定，也可以通过交叉验证来选取较好的维数值
• PCA仅需保留W和样本的均值向量即可通过简单的向量减法和矩阵-向量乘法将新样本投影至低维空间中。而其中保存均值向量是为了对新样本同样进行中心化。

3. 利用PCA对半导体制造数据降维

• 数据集来自UCI机器学习数据库，包含590个特征，其中几乎所有样本都存在特征缺失，用NaN标识，对于这些缺失值，我们用一些简单的方法处理，比如用0去替换NaN，但是由于并不知道这些值的意义，所以该方法不合适。下面通过replaceNanWithMean将缺失的NaN数据用其他样本的相同特征值平均值填充。

3.1 示例代码：

# 缺失值替换
def replaceNanWithMean(): 
    datMat = loadDataSet('secom.data', ' ')
    numFeat = shape(datMat)[1]
    for i in range(numFeat):
        meanVal = mean(datMat[nonzero(~isnan(datMat[:,i].A))[0],i]) # 计算所有非NaN色平均值
        datMat[nonzero(isnan(datMat[:,i].A))[0],i] = meanVal # 将所有的NaN置为平均值
    return datMat

dataMat=replaceNanWithMean()
meanVals = mean(dataMat, axis=0) # 求平均值(按每一列)
meanRemoved = dataMat - meanVals   
covMat = cov(meanRemoved, rowvar=0)  # 计算协方差矩阵
eigVals,eigVects = linalg.eig(mat(covMat)) # 得到特征值和特征向量
print eigVals

nonzero()

运行结果：

[ 53415197.85687514  21746671.9046592    8248376.61529073 ...,         0.
         0.                 0.        ]

• 从特征值可以看出，其中很多值都是0，这就是说这些特征都是其他特征的副本，可以通过其他特征来表示，本身没有提供额外的信息。
• 再看下部分数值的数量级，最前面的值比较大，而后面的特征值变得非常小，这就是说只有部分的特征是主要特征
• 最后发现还有部分特征值是负值，他们主要源自数值误差应四舍五入。

3.2 主成分数目和总方差的关系

可以看出前六个主成分就覆盖了数据96.8%的方差，而前20个主成分覆盖了99.3%的方差，这就说明保留前6个主成分而去除后584个主成分，我们就可以实现大概100：1的压缩***，另外舍弃了噪声的主成分，将后面的主成分去除便使得数据更加干净。*

4. 笔记

（1）eigValInd = argsort(eigVals)使用

In [7]: a
Out[7]: [1, 2, 3, 4, 5, 6, 7, 8]

In [8]: argsort(a)
Out[8]: array([0, 1, 2, 3, 4, 5, 6, 7], dtype=int64)

In [9]: a=[3, 4, 5, 6, 7,1, 2,]  # 下标值

In [19]: bb=argsort(a)

In [21]: bb
Out[21]: array([5, 6, 0, 1, 2, 3, 4], dtype=int64)


In [23]: bb[:-3:-1]
Out[23]: array([4, 3], dtype=int64) # 得到了从最后一个元素-1到-2，包括-3(下标)

In [41]: mm=[[1,2,3,4,5],[6,7,8,9,10]]

In [42]: mm=mat(mm)

In [43]: mm
Out[43]: 
matrix([[ 1,  2,  3,  4,  5],
        [ 6,  7,  8,  9, 10]])

In [44]: bb
Out[44]: array([5, 6, 0, 1, 2, 3, 4], dtype=int64)

In [45]: bb=bb[:-3:-1]

In [46]: bb
Out[46]: array([4, 3], dtype=int64)

In [47]: mm[:,bb]
Out[47]: 
matrix([[ 5,  4],
        [10,  9]])