机器学习之主成分分析

机器学习 同时被 2 个专栏收录
3 篇文章 0 订阅

一.绪言

前段时间开始入门机器学习,《机器学习实战》(Peter Harrington 著)差不多过了一遍,第一遍主要关注算法的实现和应用。现在把注意力转移到算法的原理,这一次争取把算法总结得尽量全面,便于以后查阅。

这是本系列的第一篇文章——主成分分析。

二.什么是主成分分析

主成分分析是姜维,哦不,降维的一种。

降维就是将维度降低。在机器学习中,降维常常用来做数据的预处理。为什么要对数据进行降维了?那来从数据本身说起。

  • 大数据时代,数据冗余,维度高。例如个人用户信息,存储了身份证,同时也存储了生日,就造成了冗余。
  • 数据维度有相关性。例如,人脸头像具有对称性,去掉一般的像素点也是没有太大问题的。
  • 数据有噪声。噪声对学习会产生干扰,去掉噪声可以提高算法的精度。

那如何进行降维了?难道要随机去掉一些维度吗?答案是否定的。直接去掉维度会导致数据信息的大量确实。主成分分析(PCA)技术可以尽量保证数据信息少量减小的情况下,进行维度的缩减。

主成分分析的原理是将数据从原来的坐标系转化到新的坐标系,新坐标系的选择是由数据本身决定的。第一个新坐标轴选择的是原始数据中方差最大的方向,第二个新坐标轴的选择和第一个坐标轴正交且具有最大方差的方向。重复选择新坐标轴之后,最终大部分方差都包含在最前面的几个新坐标轴中。因此,可以忽略余下的坐标轴,达到降维的目的。

举个栗子:

在上图中,在第一次移动坐标轴时,先找到了数据的最大方差方向,即此时的 u u u方向,接着第二次移动坐标轴,找到的次大方差方向且与 u u u轴正交的 v v v方向。

这里,第一个主成分就是从数据差异性最大的(即方差最大)的方向,即 u u u。第二个主成分则来自于数据差异性次大的方向,并且与第一个主成分方向正交,即 v v v

我们可以只用 u u u方向来研究原来的数据,就达到了降维的目的,而且最大得保留了原数据的差异性。

三.数学原理

总体主成分

x = ( x 1 , x 2 , … , x m ) T x=(x_{1},x_{2},\dots,x_{m})^{T} x=(x1,x2,,xm)T是m维随机变量,其均值向量是 μ \mu μ

μ = E ( x ) = x = ( μ 1 , μ 2 , … , μ m ) T \mu=E(x)=x=(\mu_{1},\mu_{2},\dots,\mu_{m})^{T} μ=E(x)=x=(μ1,μ2,,μm)T

协方差矩阵是 Σ \Sigma Σ

Σ = c o v ( x , x ) = E [ ( x − μ ) ( x − μ ) T ] \Sigma=cov(x,x)=E[(x-\mu)(x-\mu)^{T}] Σ=cov(x,x)=E[(xμ)(xμ)T]

考虑由m维随机变量 x x x到m维随机变量 y = ( y 1 , y 2 , … , y m ) T y=(y_{1},y_{2},\dots,y_{m})^{T} y=(y1,y2,,ym)T的线性变换
y i = α i T x = α 1 i x 1 + α 2 i x 2 + ⋯ + α m i x m y_{i}=\alpha_{i}^{T}x=\alpha_{1i}x_{1}+\alpha_{2i}x_{2}+\dots+\alpha_{mi}x_{m} yi=αiTx=α1ix1+α2ix2++αmixm
其中 α i T = ( α 1 i , α 2 i , … , α m i ) \alpha_{i}^{T}=(\alpha_{1i},\alpha_{2i},\dots,\alpha_{mi}) αiT=(α1i,α2i,,αmi)

于是y有下列性质:

E ( y i ) = α i T μ E(y_{i})=\alpha_{i}^{T}\mu E(yi)=αiTμ
v a r ( y i ) = α i T Σ α i var(y_{i})=\alpha_{i}^{T}\Sigma\alpha_{i} var(yi)=αiTΣαi
c o v ( y i , y j ) = α i T Σ α j cov(y_{i},y_{j})=\alpha_{i}^{T}\Sigma\alpha_{j} cov(yi,yj)=αiTΣαj

总体主成分的定义

  1. α 1 , α 2 , … , α m \alpha_{1},\alpha_{2},\dots,\alpha_{m} α1,α2,,αm是一组标准正交集。
  2. y i y_{i} yi y j y_{j} yj互不相关,即 c o v ( y i , y j ) = 0 ( i ≠ j ) cov(y_{i},y_{j})=0(i\neq j) cov(yi,yj)=0(i=j)
  3. y 1 y_{1} y1 x x x的所有线性变换中方差最大的, y 2 y_{2} y2是与 y 1 y_{1} y1不相关的所有线性变换中方差最大的。依次类推。

从定义中可以知道,求 y 1 y_{1} y1的方法是在 α 1 T α 1 = 1 \alpha_{1}^{T}\alpha_{1}=1 α1Tα1=1的条件下让 α 1 T x \alpha_{1}^{T}x α1Tx的方差最大。

这是一个极值问题,用拉格朗日乘子法可以求解。

求解的结果包含在下面的定理中

定理1
x x x是m维随机变量, Σ \Sigma Σ x x x的协方差矩阵, Σ \Sigma Σ的特征值分别是 λ 1 ≥ λ 2 ≥ ⋯ ≥ λ m ≥ 0 \lambda_{1}\geq\lambda_{2}\geq\dots\geq\lambda_{m}\geq0 λ1λ2λm0,特征值对应的单位特征向量分别是 α 1 ≥ α 2 ≥ ⋯ ≥ α m ≥ 0 \alpha_{1}\geq\alpha_{2}\geq\dots\geq\alpha_{m}\geq0 α1α2αm0 ,则 x x x的第 k k k主成分是
y k = α k T x = α 1 k x i + α 2 k x 2 + ⋯ + α m k x m , k = 1 , 2 , … , m y_{k}=\alpha_{k}^{T}x=\alpha_{1k}x_{i}+\alpha_{2k}x_{2}+\dots+\alpha_{mk}x_{m},k =1,2,\dots,m yk=αkTx=α1kxi+α2kx2++αmkxm,k=1,2,,m
x x x的第 k k k主成分方差是
v a r ( y k ) = α k T Σ α k = λ k var(y_{k})=\alpha_{k}^{T}\Sigma\alpha_{k}=\lambda_{k} var(yk)=αkTΣαk=λk
证明过程暂时省略,立个flag,以后会补充的。(可以参考《统计学习方法第二版》P301)

样本主成分

假设对m维随机变量 x = ( x 1 , x 2 , … , x m ) T x=(x_{1},x_{2},\dots,x_{m})^{T} x=(x1,x2,,xm)T进行n次独立观测, z 1 , z 2 , … , z n z_{1},z_{2},\dots,z_{n} z1,z2,,zn表示观测样本,其中 y j = ( z 1 j , z 2 j , … , z m j ) T y_{j}=(z_{1j},z_{2j},\dots,z_{mj})^{T} yj=(z1j,z2j,,zmj)T表示的 j j j个观测样本, z i j z_{ij} zij表示 j j j个观测样本的第 i i i个变量 j = 1 , 2 , … , n j=1, 2, \dots,n j=1,2,,n
观测数据用样本矩阵 X X X表示
X = [ z 1 , z 2 , … , z n ] X=[z_{1},z_{2},\dots,z_{n}] X=[z1,z2,,zn]

样本主成分与总体主成分的区别在于:样本主成分用样本的均值向量 z ‾ \overline{z} z代替总体主成分的 μ \mu μ,用协方差矩阵 S S S来代替总体主成分的 Σ \Sigma Σ
其中:
z ‾ = 1 n ∑ j = 1 n z j = ( z ‾ 1 , z ‾ 2 , … , z ‾ m ) \overline{z}=\frac{1}{n}\sum_{j=1}^{n}z_{j}=(\overline{z}_{1},\overline{z}_{2},\dots,\overline{z}_{m}) z=n1j=1nzj=(z1,z2,,zm)
S = [ s i j ] S=[s_{ij}] S=[sij]
s i j = 1 n − 1 ∑ k = 1 n ( z i k − z ‾ i ) ( z j k − z ‾ j ) , i , j = 1 , 2 , … , m s_{ij}=\frac{1}{n-1}\sum_{k=1}^{n}(z_{ik}-\overline{z}_{i})(z_{jk}-\overline{z}_{j}),i,j=1, 2, \dots, m sij=n11k=1n(zikzi)(zjkzj),i,j=1,2,,m

考虑由m维随机变量 x x x到m维随机变量 y = ( y 1 , y 2 , … , y m ) T y=(y_{1},y_{2},\dots,y_{m})^{T} y=(y1,y2,,ym)T的线性变换
y = A T x y=A^{T}x y=ATx
其中 y i = α i T x y_{i}={\alpha}^{T}_{i}x yi=αiTx
y i y_{i} yi的样本均值 y ‾ i = 1 n ∑ j = 1 n α i T z j = α i T z ‾ \overline{y}_{i}=\frac{1}{n}\sum_{j=1}^{n}\alpha^{T}_{i}z_{j}=\alpha^{T}_{i}\overline{z} yi=n1j=1nαiTzj=αiTz
y i y_{i} yi的样本方差为 v a r ( y i ) = 1 n − 1 ∑ j = 1 n ( α i T z j − α i T z ‾ ) 2 = α i T S α i var(y_{i})=\frac{1}{n-1}\sum_{j=1}^{n}(\alpha^{T}_{i}z_{j}-\alpha^{T}_{i}\overline{z})^{2}=\alpha^{T}_{i}S\alpha_{i} var(yi)=n11j=1n(αiTzjαiTz)2=αiTSαi
y i y_{i} yi y j y_{j} yj的样本协方差为 c o v ( y i , y j ) = 1 n − 1 ∑ k = 1 n ( α i T z k − α i T z ‾ ) ( α j T z k − α j T z ‾ ) = α i T S α j cov(y_{i},y_{j})=\frac{1}{n-1}\sum_{k=1}^{n}(\alpha^{T}_{i}z_{k}-\alpha^{T}_{i}\overline{z})(\alpha^{T}_{j}z_{k}-\alpha^{T}_{j}\overline{z})=\alpha^{T}_{i}S\alpha_{j} cov(yi,yj)=n11k=1n(αiTzkαiTz)(αjTzkαjTz)=αiTSαj

样本主成分的定义

  1. α 1 , α 2 , … , α m \alpha_{1},\alpha_{2},\dots,\alpha_{m} α1,α2,,αm是一组标准正交集。
  2. y i y_{i} yi y j y_{j} yj互不相关,即 c o v ( y i , y j ) = 0 ( i ≠ j ) cov(y_{i},y_{j})=0(i\neq j) cov(yi,yj)=0(i=j)
  3. y 1 y_{1} y1 x x x的所有线性变换中方差最大的, y 2 y_{2} y2是与 y 1 y_{1} y1不相关的所有线性变换中方差最大的。依次类推。

求解 y 1 , y 2 , … , y m y_{1},y_{2},\dots,y_{m} y1,y2,,ym的过程与总体主成分的过程几乎一样,只是用样本的均值向量 z ‾ \overline{z} z代替总体主成分的 μ \mu μ,用协方差矩阵 S S S来代替总体主成分的 Σ \Sigma Σ

思维导图

这里放上《统计学习方法》第十六章主成分分析的思维导图,这一章还有一些内容这里为涉及到,日后如果需要,再扩展一下。
在这里插入图片描述

四.算法步骤

实际问题中,遇到的往往是对一些观测到的样本进行降维,因此使用的是样本主成分分析。
具体步骤如下:

  1. 去除平均值
  2. 计算协方差矩阵
  3. 计算协方差矩阵的特征值和特征向量
  4. 将特征值从大到小排列
  5. 保留前N个特征值对应的特征向量(N可以自己选取)
  6. 将数据转换到上述N各特征向量构成的空间中

五.代码实现

PCA过程实现:第一个主成分从数据差异性最大的方向获取,可以通过数据集的协方差矩阵 Convariance和特征值分析求得。下面pca函数的流程为,首先去除平均值,之后计算协方差矩阵cov,计算协方差矩阵的特征值和特征向量linalg.eig,将特征值从大到小排序,保留对应的最上面的N个特征向量,最后将数据转换到上述N个特征向量构建的新空间中。

import numpy as np
import matplotlib.pyplot as plt

"""
函数说明:解析文本数据

Parameters:
    filename - 文件名
    delim - 每一行不同特征数据之间的分隔方式,默认是tab键‘\t’

Returns:
    j将float型数据值列表转化为矩阵返回

"""


def loadDataSet(filename, delim='\t'):
    fr = open(filename)
    stringArr = [line.strip().split(delim) for line in fr.readlines()]
    datArr = [list(map(float, line)) for line in stringArr]
    return np.mat(datArr)


"""
函数说明:PCA特征维度压缩函数

Parameters:
    dataMat - 数据集数据
    topNfeat - 需要保留的特征维度,即要压缩成的维度数,默认4096

Returns:
    lowDDataMat - 压缩后的数据矩阵
    reconMat - 压缩后的数据矩阵反构出原始数据矩阵
"""


def pca(dataMat, topNfeat=4096):
    # 求矩阵每一列的均值
    meanVals = np.mean(dataMat, axis=0)
    # 数据矩阵每一列特征减去该列特征均值
    meanRemoved = dataMat - meanVals
    # 计算协方差矩阵,处以n-1是为了得到协方差的无偏估计
    # cov(x, 0) = cov(x)除数是n-1(n为样本个数)
    # cov(x, 1)除数是n
    covMat = np.cov(meanRemoved, rowvar=0)
    # 计算协方差矩阵的特征值及对应的特征向量
    # 均保存在相应的矩阵中
    eigVals, eigVects = np.linalg.eig(np.mat(covMat))
    # sort():对特征值矩阵排序(由小到大)
    # argsort():对特征矩阵进行由小到大排序,返回对应排序后的索引
    eigValInd = np.argsort(eigVals)
    # 从排序后的矩阵最后一个开始自下而上选取最大的N个特征值,返回其对应的索引
    eigValInd = eigValInd[: -(topNfeat + 1): -1]
    # 将特征值最大的N个特征值对应索引的特征向量提取出来,组成压缩矩阵
    redEigVects = eigVects[:, eigValInd]
    # 将去除均值后的矩阵*压缩矩阵,转换到新的空间,使维度降低为N
    lowDDataMat = meanRemoved * redEigVects
    # 利用降维后的矩阵反构出原数据矩阵(用作测试,可跟未压缩的原矩阵比对)
    # 此处用转置和逆的结果一样redEigVects.I
    reconMat = (lowDDataMat * redEigVects.T) + meanVals
    print(reconMat)
    # 返回压缩后的数据矩阵及该矩阵反构出原始数据矩阵
    return lowDDataMat, reconMat


if __name__ == '__main__':
    dataMat = loadDataSet('testSet.txt')
    lowDmat, reconMat = pca(dataMat, 1)
    print(np.shape(lowDmat))
    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=90, c='red')
    plt.show()

结果:
在这里插入图片描述
代码中的一点解释:

# 将去除均值后的矩阵*压缩矩阵,转换到新的空间,使维度降低为N
lowDDataMat = meanRemoved * redEigVects
# 利用降维后的矩阵反构出原数据矩阵(用作测试,可跟未压缩的原矩阵比对)
# 此处用转置和逆的结果一样redEigVects.I
 reconMat = (lowDDataMat * redEigVects.T) + meanVals

第一步,每个数据乘以特征向量,相当于往特征方向投影。所以这一步得到的结果是每条数据在以选出的特征向量方向构建的坐标系下的坐标值。
第二步,相当于求压缩后的数据在原坐标系下的坐标值。

一个例子

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

数据集来自UCI机器学习数据库,包含590个特征,其中几乎所有样本都存在特征缺失,用NaN表示,通过replaceNanWithMean将缺失的NaN数据用其他样本的相同特征值平均值填充。
这个栗子中将590维的数据压缩为6维。

# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt

"""
函数说明:解析文本数据

Parameters:
    filename - 文件名
    delim - 每一行不同特征数据之间的分隔方式,默认是tab键‘\t’
    
Returns:
    j将float型数据值列表转化为矩阵返回
"""
def loadDataSet(filename, delim='\t'):
    fr = open(filename)
    stringArr = [line.strip().split(delim) for line in fr.readlines()]
    datArr = [list(map(float, line)) for line in stringArr]
    return np.mat(datArr)


"""
函数说明:PCA特征维度压缩函数

Parameters:
    dataMat - 数据集数据
    topNfeat - 需要保留的特征维度,即要压缩成的维度数,默认4096
    
Returns:
    lowDDataMat - 压缩后的数据矩阵
    reconMat - 压缩后的数据矩阵反构出原始数据矩阵
"""
def pca(dataMat, topNfeat=4096):
    # 求矩阵每一列的均值
    meanVals = np.mean(dataMat, axis=0)
    # 数据矩阵每一列特征减去该列特征均值
    meanRemoved = dataMat - meanVals
    # 计算协方差矩阵,处以n-1是为了得到协方差的无偏估计
    # cov(x, 0) = cov(x)除数是n-1(n为样本个数)
    # cov(x, 1)除数是n
    covMat = np.cov(meanRemoved, rowvar=0)
    # 计算协方差矩阵的特征值及对应的特征向量
    # 均保存在相应的矩阵中
    eigVals, eigVects = np.linalg.eig(np.mat(covMat))
    # sort():对特征值矩阵排序(由小到大)
    # argsort():对特征矩阵进行由小到大排序,返回对应排序后的索引
    eigValInd = np.argsort(eigVals)
    # 从排序后的矩阵最后一个开始自下而上选取最大的N个特征值,返回其对应的索引
    eigValInd = eigValInd[: -(topNfeat+1): -1]
    # 将特征值最大的N个特征值对应索引的特征向量提取出来,组成压缩矩阵
    redEigVects = eigVects[:, eigValInd]
    # 将去除均值后的矩阵*压缩矩阵,转换到新的空间,使维度降低为N
    lowDDataMat = meanRemoved * redEigVects
    # 利用降维后的矩阵反构出原数据矩阵(用作测试,可跟未压缩的原矩阵比对)
    # 此处用转置和逆的结果一样redEigVects.I
    reconMat = (lowDDataMat * redEigVects.T) + meanVals
    print(reconMat)
    # 返回压缩后的数据矩阵及该矩阵反构出原始数据矩阵
    return lowDDataMat, reconMat


"""
函数说明:缺失值处理函数

Parameters:
    None
    
Returns:
    datMat - 处理后的数据集
"""
def replaceNaNWithMean():
    # 解析数据
    datMat = loadDataSet('secom.data', ' ')
    # 获取特征维度
    numFeat = np.shape(datMat)[1]
    for i in range(numFeat):
        # 利用该维度所有非NaN特征求取均值
        meanVal = np.mean(datMat[np.nonzero(~np.isnan(datMat[:, i].A))[0], i])
        # 若均值也是NaN则用0代替
        if (np.isnan(meanVal)):
            meanVal = 0
        # 将该维度中所有NaN特征全部用均值替换
        datMat[np.nonzero(np.isnan(datMat[:, i].A))[0], i] = meanVal
    return datMat


if __name__ == '__main__':
    dataMat = replaceNaNWithMean()
    # lowDmat, reconMat = pca(dataMat, topNfeat=20)
    fig = plt.figure()
    ax = fig.add_subplot(111)
    # ax.scatter(reconMat[:, 0].flatten().A[0], reconMat[:, 1].flatten().A[0], marker='o', s=90, c='red')
    meanVals = np.mean(dataMat, axis=0)
    meanRemoved = dataMat - meanVals
    covMat = np.cov(meanRemoved, rowvar=0)
    eigVals, eigVects = np.linalg.eig(np.mat(covMat))
    # 特征数
    i = 6
    ax.scatter(range(i), eigVals[:i], marker='^', s=50, c='red')
    ax.plot(range(i), eigVals[:i])
    plt.show()
    lowDmat, reconMat = pca(dataMat, topNfeat=i)
    # 提取的6个特征
    print(lowDmat)

结果:
在这里插入图片描述

六.算法评价

优点:降低数据的复杂性,识别最重要的多个特征
缺点:不一定需要,且可能损失有用信息

参考资料

1.《统计学习方法》第二版 李航 著
2.《机器学习实战》Peter Harrington 著
3.机器学习实战之主成分分析(PCA)
4.机器学习笔记(Chapter 13 - PCA简化)
5.机器学习实战——PCA(主成分分析)

  • 1
    点赞
  • 0
    评论
  • 4
    收藏
  • 一键三连
    一键三连
  • 扫一扫,分享海报

©️2021 CSDN 皮肤主题: 大白 设计师:CSDN官方博客 返回首页
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值