降维PCA

##1. 维度灾难

###1.1 The curse of dimensionality

Richard E.Bellman(发明动态规划的美国应用数学家) 在１９６１年提出了这个术语：

The “Curse of dimensionality”, is a term coined by Bellman to describe the problem caused by the exponential increase in volume associated with adding extra dimensions to a (mathematical) space.

###1.2 用３类模式识别问题举例

一个简单的方法是：

将特征空间划分为小bins
计算每个bins里边每个样本对于每一类的ratio
对于一个新的样本，找到它所在的bin然后将它划分到该bin里最主要的类里
这个方法类似于k-nn算法，和桶算法类似

比如我们用单特征来划分三个类的９个实例：
　
在一维我们会注意到类之间会有太多重叠,于是我们就决定引入第二个特征试图增加可分性．

这时，在二维空间，如果

我们选择保留样本密度，那么样本点数量就从９激增至9*3=27
我们选择保留样本数量，那么２维的散点图就十分稀疏
该怎么抉择呢？再增加一个feature？

在３维空间，事情变得更糟糕：

bins的数量增加到２７
维持样本密度不变，样本点的数量就增加到了８１
维持样本点个数不变，那么3D散点图几乎就是空的

很明显，我们用相同的bins来划分样本空间的办法是十分低效率的

1.3 我们如何打败维度灾难？

引入先验知识？
提高目标函数的光滑性（例如光滑核方法），可使问题的推论性增强，甚至变为解析问题．然而，这需要大量的训练样本和训练时间．
减少维度！（是否和前面加入特征矛盾）

在实践中，维度灾难意味着，给定一个样本大小，存在一个维数的上限，超过了这个上线，分类器的性能就会不升反降．
pca4
很多时候，在低维空间做更准确的映射比在高维空间直接映射有更好的分类性能．

1.4 维度灾难更加深刻的含义

保持样本密度需要指数级增长的样本数量
- 给定一个密度为N的样本和D维，需要的总样本数为$N^D$
密度估计的目标函数复杂性也随着维度增长呈指数增长

定义在高维空间的函数比定义在低维空间的函数复杂的多，并且那些复杂性是难以辨明的复杂．这意味着，为了更好的学习它，越复杂的函数需要越密集的样本点 –by Friedman( famous for his friedman test)

如果样本的分布不是高斯分布情况会更糟糕
-　因为在教科书里只有大量的１维分布函数，到了高维情况就只剩下了高斯分布，而且在维度很高的情况下，高斯密度函数只能在简化的形势下处理．
人脑也对４维以上的特征失去了辨识能力（但是人脑为什么可以辨识复杂的模式呢）

##2. 降维

2.1 特征选择　vs. 特征提取

Feature extraction:从现有的特征组合中生成一个特征子集
Feature selection:选择包含全部特征信息的子集（新的特征由原来的一维或者多维特征映射获得）

2.1.1 特征抽取的定义

给定一个特征空间　$x_i \in R^N$ 找出一个映射：　$y=f(x):R^N \rightarrow R^M$ 其中　$M \lt N$，这样，变换后的特征向量　$y_i \in R^M$　保存了（大多数）原来特征空间　$R^N$内的信息（或结构）
一个最优的映射的引入$y=f(x)$不会增加分类错误（就是对同一个算法说和原来的分类正确率一样）

2.2 最优映射 y=f(x)　

通常，一个最优映射　y=f(x)是一个非线性函数，但是，没有一个系统和成熟的方法进行非线性变换，非线性变换的数学理论还很不成熟．（对特定特征子集的选择还要结合具体问题具体分析）．因此，特征抽取常常被限定为了线性变换　$y=Wx$.

这就是说，y 是　x　的一个线性投影
注意：当映射一个非线性函数时，减少后的特征空间叫做manifold(流形)．参见Nonlinear dimensionality reduction

2.3 信号表示 vs. 分类

选择特征抽取映射$y=f(x)$的办法是找到一个目标函数，对其进行最大化或者最小化．按照目标函数的条件，特征提取技术被分为了两类：

Signal representation: 目标是在更低的维度对样本信息精确的表示
classification:目标是在低维上增强样本的可分类性．

pca7
　　　
在线性特征提取的领域内，常用的技术有两个：

Pricipal Components Analysis主成分分析：用于信号表示
Linear Discriminant Analysis线性判别分析(fisher 1936)：用于样本分类

##3. Principal Component Analysis(PCA)

3.1 PCA的目标

PCA的目标是在降维的同时尽可能多的保留高维空间的随机性（方差）

3.2 PCA 的推导

假设　x 是一个　N-维随机向量，表示为一组正交基向量[$\phi_1|\phi_2|….|\phiN$]的线性组合
$$
x = \sum{i=1}^N y_i \phi_i \quad where \quad \phi_i\cdot \phi_j =
\begin{cases}
0, &\text{i $\neq$ j}\[2ex]
1, &\text{ i = j }
\end{cases}
$$

假设我们只选择了基向量中的M(Ｍ＜N)维来表示x,我们可以用预先选择的常量$bj$来替换$[y{M+1},…,yN]^T$
$$
\hat x (M) = \sum{i=1}^M y_i\phii + \sum{i=M+1}^N b_i \phii
$$
两种表示形式的方差为：
$$
\Delta x(M) = x- \hat x(M) = \sum{i=1}^N y_i \phii - \left \lbrace \sum{i=1}^M y_i\phii + \sum{i=M+1}^N b_i\phii \right\rbrace = \sum{i=M+1}^N (y_i-b_i)\phi_i
$$
我们可以用均方差来表示这个误差的大小．　　　
　　　
我们的目标是找到使均方误差最小的参数：基向量$\phi_i$和常量$b_i$
$$
\overline \epsilon^2(M) = E [|\Delta x(M)|^2] = E \left [ \sum{i=M+1}^N \sum{j=M+1}^N (y_i-b_i)(y_j-b_j) \phi_i^T\phij \right] = \sum{i=M+1}^N E \left [ (y_i-b_i)^2 \right]
$$
这里之所以能抵消是因为，当 $i \neq j $ 时， $ \phi_i \cdot \phi_j = 0 $
当 $ i = j$ 时，$ \phi_i \cdot \phi_j = 1 $

问题现在转化成了求 $b_j$ 和 $\phi_i$,我们可以用对目标函数求偏导数并且令偏导数为0的方法来求$b_i$.

现在,把$b_i = E [y_i] $ 代入上边的公式,均方误差可以变换成方差的正交化形式:

其中 $\sum_X$ 是 x 的协方差矩阵.

现在,引入拉格朗日因子:

对每个 $\phi_i$ 求偏导:

于是,$\phi_i$和$\lambda_i$就是协方差矩阵 $\sum_X$ 的特征值和特征向量

这样,我们就能把均方差的和表示为:

为了最小化这个式子,$\lambda_i$ 必须是最小的特征值

因此,用最小方差和来表示 x ,我们将会选择最大的特征值$\lambda_i$对应的特征向量$\phi_i$

3.3 总结

把随机向量 X 投影到其协方差矩阵$\Sigma_X$的最大特征值$\lambda_i$对应的特征向量$\phi_i$上,我们就可以得到N维随机向量 $x \in \mathscr R^N$ 的一个最优(有误差的最小方差和)近似: M 个独立向量的线型组合.

注意:

由于PCA使用了协方差矩阵的特征向量,它能够找到单峰正态分布下数据的独立分布
- 对于非高斯分布或者多峰正态分布数据,PCA只是简单的去相关
PCA的主要限制是它没有考虑类别的可分性,因为它没有考虑类标签
- PCA只是简单的将坐标转向了有最大方差的方向
- 而有最大方差的方向并不保证有良好的可分类特征
  例如:
  
  用PCA的话,会降维到Y轴上,因为Y轴有最大的方差.然而降维到y轴之后,样本几乎不可分,因为都混杂在一起了.

3.4 举例

计算下列2维数据集的 pricipal components:
$ X = (x_1,x_2) = \lbrace (1,2),(3,3),(3,5),(5,4),(5,6),(6,5),(8,7),(9,8) \rbrace$
求解过程
- 协方差矩阵是:
  $$
  \Sigma_X = \begin{bmatrix} 6.25& 4.25\ 4.25 & 3.5 \end{bmatrix}
  $$
- 求特征值:
- 特征向量:

4.实现PCA算法

4.1 数据预处理

在应用PCA算法之前,我们常常要对数据进行预处理,因为PCA算法对原始数据的可靠性,一致性,以及规范性十分敏感,如果不做相应的处理,算法的效果就会大打折扣.
通常我们希望所有的特征 $x_1,x_2,…,x_n$都有相似的取值范围(并且均值接近0).在部分应用中,我们都有必要对每个特征做预处理,即通过估算每个特征$x_j$的均值和方差,而后将其取值范围规整化为零均值和单位方差.

通常,预处理主要有2步:

feature scaling/mean normalization 均值规整化为 0
白化,一些文献中也叫 sphering ,白化的目的就是降低输入的冗余性,通过白化使得学习算法的输入具有一下性质:
- (1) 特征之间相关性较低
- (2) 所有特征具有相同的方差(分布的离散程度相同)

4.1.1 均值规整化

假设我们有训练集: $ x^{(1)},x^{(2)},…x^{(m)}$
均值规整化过程:

求每个特征的均值 $\mui = \frac1m \sum{i=1}^m x_j^{(i)}$
把每个特征的值$x_j^{(i)}$ 替换为 $x_j-\mu_j$
如果不同的特征之间的取值范围不一样(例如: $x_1$ = size of house(大多取值 1000多) $x_2$ = number of bedrooms(取值不超过10) )，就需要把每个特征的刻度转换成和其他特征差不多的取值范围:
$$ x_j^{(i)} = \frac{x_j^{(i)}=\mu_j}{S_j} $$
这里 $S_j$是某种对特征j的取值范围的度量.(常常是最大值减最小值或者是标准差)

4.1.2 白化

消除输入特征之间的相关性:

将原始样本值映射到新基上的过程 $ x{rot}^{(i)} = U^T x^{(i)} $ 实际上已经消除了输入特征 $ x^{(i)}$ 之间的相关性.因为对于映射后的样本 $ x{rot} $ 来说,它的协方差矩阵 $\Sigma_{rot}$是一个对角矩阵,而且对角线上的元素值为特征值,由于非对角线元素都为0,所以特征之间是不相关的.

使所有特征具有单位方差:

为了使所有特征具有单位方差,我们可以直接使用 $ \frac1{ \sqrt{ \lambdai } }$作为缩放因子来缩放每个特征 $ x{rot,i} $.具体的,我们定义白化后的数据 $x_{PCAwhite} \in \mathscr R^N $ 如下:

$$
x{PCAwhite,i} = \frac{x{rot,i}}{ \sqrt{ \lambdai } }
$$
这样一来,白化后的数据$x{PCAwhite} $ 的协方差矩阵就变成了 单位矩阵.我们说, $x{PCAwhite} $ 是数据经过 PCA白化后的版本,$x{PCAwhite} $ 中不同的特征之间部相关并且具有单位方差.

正则化

在实践中,白化过程中,有些特征值 $ \lambda_i $ 在数值上接近于0,这样在缩放步骤时,我们除以 $ \sqrt{ \lambdai } $将导致除以一个接近0的值;这可能使数据上溢或者造成数值不稳定.因而在实践中,我们使用少量的正则化实现这个缩放过程,即在取平方根的倒数之前给特征值加上一个很小的常数 $ \epsilon $ :
$$
x{PCAwhite,i} = \frac{x_{rot,i}}{ \sqrt{ \lambda_i } + \epsilon }
$$
当 $ x $ 在区间 $ [-1, 1] $ 上时,一般取值为 $ \epsilon \approx 10^{-5} $

4.2 PCA 算法

首先,计算 协方差矩阵
$$
\Sigma = \frac1m \sum_{i=1}^n (x^{(i)})(x^{(i)})^T
$$
其中 $\Sigma$ 是 n x n 大小的协方差矩阵 , 因为 $ x^{(i)}$ 是 n x 1 (一个样本), $(x^{(i)})^T$ 是 1 x n .

计算协方差矩阵 $ \Sigma $ 的特征向量:

1
2
3

[U,S,V] = svd(Sigma);
//或者
[U,S,V] = eig(Sigma);

关于这里为什么要用 svd函数,因为 PCA 根据应用领域的不同,还有很多别名:
线性代数中叫 :SVD(singular value decomposition of X ) EVD(eigenvlaue decomposition of $X^TX$) ，KLT(信号处理),POD(机械工程).

这样,我们得到了 [U,S,V],其中 U 就是把特征向量按照对应的特征值大小从大到小排列所得到的一个 n x n 特征向量矩阵,这就是一个变换后的正交基
$$
\begin{align}
U =
\begin{bmatrix}
| & | & & | \
u_1 & u_2 & \cdots & u_n \
| & | & & |
\end{bmatrix}
\end{align}
$$
我们将要使用这个矩阵来降维,要将原来的数据 $x \in R^n$ 映射到 $z \in R^k (k<n)$,需要使用前 K 个最大的特征值对应的特征向量(因特征值大的更能代表数据特征)组生成一组新的正交基:

$$
\begin{align}
U_{reduce} =
\begin{bmatrix}
| & | & & | \
u_1 & u_2 & \cdots & u_k \
| & | & & |
\end{bmatrix}
\end{align}
$$

1	Ureduce = U(:,1:k);

这样以来,我们就可以在新基上用投影的方式表示原来的数据了:
$$
\begin{align}
Z =
\begin{bmatrix}
| & | & & | \
u_1 & u_2 & \cdots & u_k \
| & | & & |
\end{bmatrix}^T \cdot X =
\begin{bmatrix}

& (u^1)^T & - \
& (u^2)^T & - \
& \vdots & \
& (u^k)^T & -
\end{bmatrix} \cdot X =
\begin{bmatrix}
z^1 \
z^2 \
\vdots \
z^k
\end{bmatrix}
\end{align}
$$
这样,
$$
\begin{align}
\begin{bmatrix}
x^1 \
x^2 \
\vdots \
x^n
\end{bmatrix} \Longrightarrow
\begin{bmatrix}
z^1 \
z^2 \
\vdots \
z^k
\end{bmatrix}
\qquad\color{lime}{(k<n)}\quad\color{red}{高维数据被投影到低维的正交基向量组上,实现了降维}
\end{align}
$$
1
z = Ureduce'*x;

4.3 python 代码实现

数据下载地址

代码github地址

#coding=utf-8
import numpy as np
__author__ = 'anboqing'
"""
feature scaling and mean normalization
零均值化
"""
def zeroMean(dataMat):
    print("feature scaling ... ")
    meanVal = np.mean(dataMat,axis=0) # axis = 0 means that calculate mean value by column(every feature)
    print("mean normalization ...")
    newVal = dataMat - meanVal
    return newVal,meanVal
"""
    求预处理之后的样本矩阵的协方差矩阵
"""
def covarianceMat(normal_data):
    '''
    计算协方差矩阵
    :param normal_data: 零均值规范化过的数据
    :return: 协方差矩阵
    '''
    print('calculate the covariance matrix ... ')
    covMat = np.cov(normal_data,rowvar=0)
    """
    numpy中的cov函数用于求协方差矩阵,
    参数rowvar很重要,
    若rowvar=0,说明传入的数据一行代表一个样本,
    若rowvar!=0,说明一列代表一个样本
    """
    return covMat
def eigenValAndVector(covMat):
    """
    求协方差矩阵的特征值和特征向
    :param covMat: 协方差矩阵
    :return: 协方差矩阵的特征值,特征向量
    """
    print('calculate eigen value and eigen vector of covariance matrix...')
    eigVals,eigVects = np.linalg.eig(np.mat(covMat))
    return eigVals,eigVects
def reduceDimensionality(eigVals,eigVects,dataMat,k=1):
    """
    降维:保留主要成分
    用前k个特征向量组成新的正交基,把原始数据映射到新的正交基上
    :param eigVals: 原始数据的协方差矩阵的特征值
    :param eigVects: 原始数据的协方差矩阵的特诊向量
    :param dataMat: 原始数据
    :param k: 要降维到多少维
    :return: 降维后的数据,投影的正交基
    """
    print('reduce the dimension ... ')
    eigValIndices = np.argsort(eigVals) # 对特征值从小到大排序
    n_eigValIndices = eigValIndices[-1:-(k+1):-1] # 最大的k个特征值的下标
    data_mat_trans = eigVects[:,n_eigValIndices] # 取出前k个最大的特征向量
    low_dimensional_data = dataMat * data_mat_trans # 将原始数据投影到新的向量空间的正交基上
    return low_dimensional_data,data_mat_trans
def load_libsvm_data(file_path):
    fin = open(file_path)
    print('load data : '+file_path)
    data_mat=[]
    label_vec = []
    for strline in fin.readlines():
        line_arr = strline.strip().split()
        label_vec.append(int(line_arr[0]))
        line_arr = line_arr[1:]
        row_arr = []
        idx = 1
        for pair in line_arr:
            str_indice,str_val=pair.split(":")
            indice = int(str_indice)
            fval = float(str_val)
            #print indice,fval
            if indice == idx:
                row_arr.append(fval)
            else:
                row_arr.append(0.0)
            idx=idx+1
        data_mat.append(row_arr)
    return data_mat,label_vec
def choosePCASize(eig_vals,percentage=0.99):
    '''
    计算主成分个数的函数
    :param eig_vals: 特征值数组
    :param percentage: 要保留的方差的百分比
    :return:应该降到的维数
    '''
    asc_eigs = np.sort(eig_vals) # 升序
    desc_eigs = asc_eigs[-1::-1] # 翻转成从大到小
    eig_sum = sum(desc_eigs)
    tmpSum =0
    num = 0
    for i in desc_eigs:
        tmpSum += i
        num+=1
        if tmpSum >= eig_sum * percentage:
            return num
if __name__ == "__main__":
    file_path = 'madelon'
    # 读取livsvm格式的数据
    data_mat,label_vec = load_libsvm_data(file_path)
    # 转换成 numpy matrix 结构
    dataMat = np.mat(data_mat)
    # 数据零均值化
    new_dat,mean_dat = zeroMean(dataMat)
    # 计算协方差矩阵
    cov_mat = covarianceMat(new_dat)
    #print np.shape(cov_mat)
    # 计算协方差矩阵的特征值和特征向量
    eigen_values,eigen_vectors = eigenValAndVector(cov_mat)
    # 选择要降到的维数(保存99的方差)
    num_ = choosePCASize(eigen_values)
    # 试着把这个数据集降维为num_2维
    low_dim_mat,orth_base = reduceDimensionality(eigen_values,eigen_vectors,new_dat,num_2)
    print np.shape(low_dim_mat)
    print np.shape(orth_base)

4.4 选择主成分的个数

在应用PCA的时候,对于一个1000维的数据,我们怎么知道要降维到多少维才是合理的?也就是要在保留最多信息的同时取出最多的噪声.Andrew Ng的机器学习课讲的方法是:

算出 Average squared projection error: $$\frac1m \sum{i=1}^m \left\lVert x^{(i)} - x{approx}^{(i)} \right\rVert^2$$
算出 Total variation in the data: $$\frac1m \sum_{i=1}^m \left\lVert x^{(i)} \right\rVert^2$$
然后,选择满足下列条件的 k :
$$ \frac{\frac1m \sum{i=1}^m \left\lVert x^{(i)} - x{approx}^{(i)} \right\rVert^2}{\frac1m \sum{i=1}^m \left\lVert x^{(i)} \right\rVert^2} \le 0.01 \qquad(1\%)$$
也就是说: 保留了 99 % 的方差,其中 0.01可以根据应用的不同来自己设定
上式也可以改成:
$$
1 - \frac{\sum{i=1}^k \lambdai}{\sum{i=1}^n \lambda_i} \le 0.01
$$
其中$\lambda_i$ 是协方差矩阵的特征值

自动选择维数的算法:

Automatic choice of dimensionality for PCA
Thomas P. Minka

欢迎访问我的博客

Prev Home Next

降维PCA

1.3 我们如何打败维度灾难？

1.4 维度灾难更加深刻的含义

2.1 特征选择 vs. 特征提取

2.1.1 特征抽取的定义

2.2 最优映射 y=f(x)

2.3 信号表示 vs. 分类

3.1 PCA的目标

3.2 PCA 的推导

3.3 总结

3.4 举例

4.实现PCA算法

4.1 数据预处理

4.1.1 均值规整化

4.1.2 白化

4.2 PCA 算法

4.3 python 代码实现

数据下载地址

代码github地址

4.4 选择主成分的个数

自动选择维数的算法:

Comments

2.1 特征选择　vs. 特征提取

2.2 最优映射 y=f(x)　