3

终于明白PCA降维的数学原理了

 2 years ago
source link: https://weisenhui.top/posts/22264.html
Go to the source link to view the article. You can view the picture content, updated content and better typesetting reading experience. If the link is broken, please click the button below to view the snapshot at that time.
neoserver,ios ssh client

一、相关资料

本文主要参考了PCA的数学原理 ,除此之外加上了自己的一些理解和笔记,以及PCA代码实战。

在机器学习领域中,我们对原始数据提取特征后的维度可能是很高的,比如10万维,在这个高维空间中,包含了很多的冗余和噪声。我们希望通过降维的方式寻找数据内部的特性,从而提高特征的表达能力,降低训练的复杂度。

PCA(Principal Component Analysis)作为降维中最经典的方法,至今已有100多年的历史,PCA属于一种线性、非监督、全局的降维算法,是在面试中经常被问到的问题。

机器学习中经常要处理成千上万甚至几十万维的数据在这种情况下,机器学习的资源消耗是不可接受的,因此我们必须对数据进行降维。

但是降维必然意味着信息的丢失,不过鉴于实际数据本身常常存在的相关性,我们可以想办法在降维的同时将信息的损失尽量降低。

例如淘宝店铺的数据,从经验我们可以知道,“浏览量”和“访客数”往往具有较强的相关关系,而“下单数”和“成交数”也具有较强的相关关系。这里我们非正式的使用“相关关系”这个词,可以直观理解为“当某一天这个店铺的浏览量较高(或较低)时,我们应该很大程度上认为这天的访客数也较高(或较低)”。

这种情况表明,如果我们删除浏览量或访客数其中一个指标,我们应该期待并不会丢失太多信息。因此我们可以删除一个,以降低机器学习算法的复杂度。

上面给出的是降维的朴素思想描述,可以有助于直观理解降维的动机和可行性,但并不具有操作指导意义。

  • 我们是单纯地删除几列?
  • 还是通过某些变换将原始数据变成更少的列,来使得丢失的信心最小呢?

三、PCA数学原理(最大方差理论,最大可分性)

设向量B的模为1,则A与B的内积值等于A向B所在直线投影的矢量长度
A⋅B=|A||B|cos(θ)=|A|cos(θ)

要准确描述向量,首先要确定一组基,然后给出在基所在的各个直线上的投影值(该向量与基进行点积),就可以了

一个二维向量可以对应二维笛卡尔直角坐标系中从原点出发的一个有向线段

例如,在二维空间中向量3,2,他在e1=(1,0)和e2=(0,1)基下的坐标为3,2。实际上我们可以选择任意两组线性无关的向量作为基。α1=(1/√2,1/√2)和α2=(−1/√2,1/√2)也可以成为一组基,一般来说,我们希望基的模是1,因为从内积的意义可以看到,如果基的模是1,那么就可以方便的用向量 点积 基而直接获得其在新基上的坐标了!(==这个结论在之后会反复用到==)

如何获得3,2在新基上的坐标?

向量3,2在这组新基下的坐标为((3,2)⋅α1,(3,2)⋅α2)=(5√2,−1√2)

通常来说,基并不要求是正交的,但是呢,通常来说正交性质比较好,所以一般使用的基都是正交的。

3.2 基变换

标准正交基:基两两之间点积为0,并且每个基的范数为1

问题:设向量v=(5,1,3,4,0,..,3)假设(e1,e2,…,en)是n维空间中的自然基(也是标准正交基),(α1,α2,…,αn)是n维空间中的另一组标准正交基,则n维向量v在(α1,α2,…,αn)这组基下的坐标是多少?

根据前面的介绍向量v在新基下的坐标为
(α1⋅v,α2⋅v,…,αn⋅v)=(3,7,18,..,0,4)

那如果我现在有很多个向量 v1,v2,..vm,那我就可以转换成矩阵的形式来表示这m个向量在新基下的坐标

(α1 α2 ⋮ αn)(v1v2⋯vm)=(α1v1α1v2⋯α1vm α2v1α2v2⋯α2vm ⋮⋮⋱⋮ αnv1αnv2⋯αnvm)

显然选择不同的基可以对同样一组数据(m个向量)可以给出不同的表示

3.3 如何将数据从高维降到低维,同时信息损失更少?

上述的操作没有实战性,因为我们没法衡量要删哪几列,而PCA给了一种具体的衡量标准。

3.4 均值、方差、协方差

3.5 PCA降维的优化目标

简单来说,PCA希望找到这样一组基(α1,α2,…,αn),使得:

原始的m条数据在这组基下的m个坐标(m个新数据)

  • m个新数据,同一个维度上的m个值方差尽量大

    • 因为我们我们希望投影到每个基上的投影值尽可能分散(熵越大所含信息越多),而这种分散程度,可以用数学上的方差来表述。
    • 寻找这样的基,使得所有数据变换为这个基上的坐标表示后,方差值最大
  • m个新数据,不同维度之间互不线性相关,即两个维度的协方差为0(协方差为 0 时,两个变量线性不相关,仍可能有高阶相关)

    • 从直观上说,让两个字段尽可能表示更多的原始信息,我们是不希望它们之间存在(线性)相关性的,因为相关性意味着两个字段不是独立的,必然存在重复表示的信息。

至此,我们得到了降维问题的优化目标:将一组N维向量降为K维(K大于0,小于N),其目标是选择K个单位(模为1)正交基,使得原始数据变换到这组基上后,各字段两两间协方差为0,而字段的方差则尽可能大(在正交的约束下,取最大的K个方差)。

用数学公式来表示上面的优化目标:

解释下为什么D=1m−1Y⋅YT是m个新数据的协方差矩阵?

我们在做PCA前,会先进行去中心化(平移一下,让每个维度的m个数据的平均值为0),好处是表示方差和协方差的时候比较简单。

3.6 最终的优化目标

换句话说,现在优化目标变成了寻找一个矩阵α,满足αCαT是一个对角矩阵(其实Cn×n=1m−1vvT,显然C是个实对称矩阵,一定有n个特征值),并且对角元素按从大到小依次排列,那么α的前K行就是要寻找的基,用α的前K行组成的矩阵乘以V就使得V从10000维降到了K维并满足上述优化条件。

至于这个矩阵α怎么求,在数学上叫做 特征值分解,也就是常看到的去求C的特征值、特征向量的做法。由于这个比较容易,就不再解释了,可参考线性代数书籍关于“实对称矩阵对角化”的内容。

四、PCA算法总结

PCA算法的两种实现方法:

特征值和特征向量是针对方阵才有的,而对任意形状的矩阵都可以做奇异值分解。

4.1 基于特征值分解协方差矩阵实现PCA算法

输入:数据集Xn维,需要降到k维

  1. 零均值化,每一位特征减去自己的平均值
  2. 计算协方差矩阵Cov
  3. 计算协方差矩阵的特征值和对应的特征向量
  4. 对特征值按大小进行排序,选择其中最大的k个。然后将其对应的特征向量分别作为行向量组成特征向量矩阵P
  5. 将数据转换到由k个特征向量构成的新空间中,即Y=PX

4.2 基于SVD分解协方差矩阵实现PCA算法!

输入:数据集Xn维,需要降到k维

  1. 零均值化,每一位特征减去自己的平均值
  2. 计算协方差矩阵Cov
  3. 通过SVD(奇异值分解)计算协方差矩阵的特征值和对应的特征向量
  4. 对特征值按大小进行排序,选择其中最大的k个。然后将其对应的特征向量分别作为列向量组成特征向量矩阵P
  5. 将数据转换到由k个特征向量构成的新空间中,即Y=PX

4.3 PCA 特征值分解与SVD分解对比

这一小节主要参考:https://zhuanlan.zhihu.com/p/77151308

而实际上 Sklearn 的 PCA 就是用 SVD 进行求解的,原因有以下几点:

  1. 当样本维度很高时,协方差矩阵计算太慢;
  2. 方阵特征值分解计算效率不高;
  3. SVD 除了特征值分解这种求解方式外,还有更高效更准球的迭代求解方式,避免了ATA的计算;
  4. 其实 PCA 与 SVD 的右奇异向量的压缩效果相同。

五、进一步讨论

根据上面对PCA的数学原理的解释,我们可以了解到一些PCA的能力和限制。PCA本质上是将方差最大的方向作为主要特征,并且在各个正交方向上将数据“离相关”,也就是让它们在不同正交方向上没有相关性。

因此,PCA也存在一些限制,例如它可以很好的解除线性相关,但是对于高阶相关性就没有办法了,对于存在高阶相关性的数据,可以考虑Kernel PCA,通过Kernel函数将非线性相关转为线性相关,关于这点就不展开讨论了。另外,PCA假设数据各主特征是分布在正交方向上,如果在非正交方向上存在几个方差较大的方向,PCA的效果就大打折扣了。

最后需要说明的是,PCA是一种无参数技术,也就是说面对同样的数据,如果不考虑清洗,谁来做结果都一样,没有主观参数的介入,所以PCA便于通用实现,但是本身无法个性化的优化。

希望这篇文章能帮助朋友们了解PCA的数学理论基础和实现原理,借此了解PCA的适用场景和限制,从而更好的使用这个算法。

六、实战演练

6.1 代码:Python手写PCA(特征值分解)

X降维到Y, 再从Y重构X, 有信息损失

from sklearn import datasets
import numpy as np
import matplotlib.pyplot as plt


class myPCA():
    def __init__(self, k):
        self.k = k

    # 数据去均值(中心化)
    def stand_data(self, data):
        mean_vector = np.mean(data, axis=0)
        return mean_vector, data - mean_vector

    # 计算协方差矩阵
    def getCovMat(self, standData):
        # rowvar=0表示数据的每一列代表一个维度
        return np.cov(standData, rowvar=0)

    # 计算协方差矩阵的特征值和特征向量
    def getFValueAndFVector(self, covMat):
        fValue, fVector = np.linalg.eig(covMat)
        return fValue, fVector

    # 得到特征向量矩阵
    def getVectorMatrix(self, fValue, fVector):
        # 从大到小排序,并返回排序后的原索引值
        fValueSort = np.argsort(-fValue)
        fValueTopN = fValueSort[:self.k]
        return fVector[:, fValueTopN]

    # 得到降维后的数据
    def getResult(self, data, vectorMat):
        return np.dot(data, vectorMat)


if __name__ == "__main__":
    k = 1

    # -----------------------------
    # # 加载鸢尾花数据集中的特征作为PCA的原始数据集
    # iris = datasets.load_iris()
    # X = iris.data
    # iris_y = iris.target

    # 没有去中心化的原始数据X
    X = np.array([[-1, -1],
                  [-1, 1],
                  [0, 1],
                  [2, 2],
                  [0, 2]])  # X: mxn, m条n维的数据

    # -----------------------------
    pca = myPCA(k)

    print("原始数据:\n%s" % X)
    (mean_vector, standdata) = pca.stand_data(X)
    print("均值向量为:%s \n标准化数据:\n%s" % (mean_vector, standdata))
    cov_mat = pca.getCovMat(standdata)
    print("协方差矩阵:\n%s" % cov_mat)
    fvalue, fvector = pca.getFValueAndFVector(cov_mat)
    print("特征值为:%s\n特征向量为:\n%s" % (fvalue, fvector))
    fvectormat = pca.getVectorMatrix(fvalue, fvector)
    print("最终需要的特征向量:\n%s" % fvector)
    Y = pca.getResult(standdata, fvectormat)
    print("降维后的数据:\n%s" % Y)

    # X降维成Y后,对数据进行重构得到X',通常是有信息损失的,即X不等于X'
    print("最终重构结果为:\n{}".format(np.mat(Y) * fvectormat.T + mean_vector))

    # # 注意令k=2
    # plt.scatter(Y[:, 0], Y[:, 1], c=iris_y)
    # plt.show()
原始数据:
[[-1 -1]
 [-1  1]
 [ 0  1]
 [ 2  2]
 [ 0  2]]
均值向量为:[0. 1.] 
标准化数据:
[[-1. -2.]
 [-1.  0.]
 [ 0.  0.]
 [ 2.  1.]
 [ 0.  1.]]
协方差矩阵:
[[1.5 1. ]
 [1.  1.5]]
特征值为:[2.5 0.5]
特征向量为:
[[ 0.70710678 -0.70710678]
 [ 0.70710678  0.70710678]]
最终需要的特征向量:
[[ 0.70710678 -0.70710678]
 [ 0.70710678  0.70710678]]
降维后的数据:
[[-2.12132034]
 [-0.70710678]
 [ 0.        ]
 [ 2.12132034]
 [ 0.70710678]]
最终重构结果为:
[[-1.5 -0.5]
 [-0.5  0.5]
 [ 0.   1. ]
 [ 1.5  2.5]
 [ 0.5  1.5]]

由于基每个维度同时乘以-1,仍然是-1,所以最后的结果可能不唯一(某些轴调个方向),但描述数据的性质是一致的

6.2 代码:sklearn PCA(SVD分解)

sklearn中的PCA使用的是SVD分解

import numpy as np
from sklearn.decomposition import PCA

# 2维(高维)数据降到2维(低维)
k = 2
pca = PCA(n_components=k)

# m = 5, n = 2, k = 1
X = [[-1, -2],
     [-1, 0],
     [0, 0],
     [2, 1],
     [0, 1]]
X = np.array(X)  # X: mxn, m条n维的数据

Y = pca.fit_transform(X)  # Y_T=alpha矩阵乘以X: mxk
print("在基alpha下m条数据的新坐标:\n", Y)

print("实对称矩阵C的前k个大的特征值:\n", pca.explained_variance_)

# PCA选出来的k组基(也是C的特征向量) alpha_1, ... alpha_k
base_vector_alpha = pca.components_  # 选前k个基 alpha_1,alpha_2,...alpha_k, ... alpha_n
print("PCA选出来的k组基(kxn):\n", base_vector_alpha)  # kxn


print("--------实际上sklearn PCA用的是SVD分解--------")
print("奇异值:\n", pca.singular_values_)
在基alpha下m条数据的新坐标:
 [[ 2.12132034  0.70710678]
 [ 0.70710678 -0.70710678]
 [-0.          0.        ]
 [-2.12132034  0.70710678]
 [-0.70710678 -0.70710678]]
实对称矩阵C的前k个大的特征值:
 [2.5 0.5]
PCA选出来的k组基(kxn):
 [[-0.70710678 -0.70710678]
 [ 0.70710678 -0.70710678]]
--------实际上sklearn PCA用的是SVD分解--------
奇异值:
 [3.16227766 1.41421356]

另外如果对iris数据集进行降维150x4->150x2,最终可视化效果如下图所示:

注意,算协方差的时候是样本协方差所以除以的是m-1不是m。

实际上,sklearn PCA用的是SVD分解,不是特征值分解,具体请看:https://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html

7.1 PCA的最小平方误差理论(最近重构性)!

除了基于最大可分性的角度去理解PCA,还可以用最小平方误差理论。我们可以将其转换为线型回归问题,其目标是求解一个线性函数使得对应直线能够更好地拟合样本点集合。这就使得我们的优化目标从方差最大转化为平方误差最小,因为映射距离越短,丢失的信息也会越小。

7.2 零均值化(训练集和测试集)

当对训练集进行 PCA 降维时,也需要对验证集、测试集执行同样的降维。而对验证集、测试集执行零均值化操作时,均值必须从训练集计算而来,不能使用验证集或者测试集的中心向量。

其原因也很简单,因为我们的训练集时可观测到的数据,测试集不可观测所以不会知道其均值,而验证集再大部分情况下是在处理完数据后再从训练集中分离出来,一般不会单独处理。如果真的是单独处理了,不能独自求均值的原因是和测试集一样。

我之前参加过一个竞赛,里面也用到了PCA,不过他不是对训练集做的,而是对整个物品集做的,所以不需要考虑训练集和测试集的问题。


About Joyk


Aggregate valuable and interesting links.
Joyk means Joy of geeK