首页 > 其他分享 >【机器学习-无监督学习】降维与主成分分析

【机器学习-无监督学习】降维与主成分分析

时间:2024-10-01 14:22:49浏览次数:9  
标签:机器 特征向量 boldsymbol 矩阵 学习 降维 PCA mathrm

在这里插入图片描述

【作者主页】Francek Chen
【专栏介绍】 ⌈ ⌈ ⌈Python机器学习 ⌋ ⌋ ⌋ 机器学习是一门人工智能的分支学科,通过算法和模型让计算机从数据中学习,进行模型训练和优化,做出预测、分类和决策支持。Python成为机器学习的首选语言,依赖于强大的开源库如Scikit-learn、TensorFlow和PyTorch。本专栏介绍机器学习的相关算法以及基于Python的算法实现。
【GitCode】专栏资源保存在我的GitCode仓库:https://gitcode.com/Morse_Chen/Python_machine_learning

文章目录


  在上一篇文章聚类中,我们介绍了无监督学习的重要问题之一:聚类问题,并主要讲解了k均值算法。结尾处我们提到,在解决复杂聚类问题时,第一步通常不会直接使用k均值算法,而是会先用其他手段提取数据的有用特征。对于高维复杂数据来说,其不同维度代表的特征可能存在关联,还有可能存在无意义的噪声干扰。因此,无论后续任务是有监督学习还是无监督学习,我们都希望能先从中提取出具有代表性、能最大限度保留数据本身信息的几个特征,从而降低数据维度,简化之后的分析和计算。这一过程通常称为数据降维(dimensionality reduction),同样是无监督学习中的重要问题。本文就来介绍数据降维中最经典的算法——主成分分析(principal component analysis,PCA)。

一、降维

(一)维数灾难

  维数灾难(Curse of Dimensionality)通常是指在涉及到向量的计算的问题中,随着维数的增加,计算量呈指数倍增长的一种现象。在很多机器学习问题中,训练集中的每条数据经常伴随着上千、甚至上万个特征。要处理这所有的特征的话,不仅会让训练非常缓慢,还会极大增加搜寻良好解决方案的困难。这个问题就是我们常说的维数灾难。

在这里插入图片描述

图1 维数灾难

  维数灾难涉及数字分析、抽样、组合、机器学习、数据挖掘和数据库等诸多领域。在机器学习的建模过程中,通常指的是随着特征数量的增多,计算量会变得很大,如特征达到上亿维的话,在进行计算的时候是算不出来的。有的时候,维度太大也会导致机器学习性能的下降,并不是特征维度越大越好,模型的性能会随着特征的增加先上升后下降

(二)降维概述

1. 什么是降维?

  降维(Dimensionality Reduction)是将训练数据中的样本(实例)从高维空间转换到低维空间,该过程与信息论中有损压缩概念密切相关。同时要明白的,不存在完全无损的降维。有很多种算法可以完成对原始数据的降维,在这些方法中,降维是通过对原始数据的线性变换实现的。

2. 为什么要降维

  • 高维数据增加了运算的难度。
  • 高维使得学习算法的泛化能力变弱(例如,在最近邻分类器中,样本复杂度随着维度成指数增长),维度越高,算法的搜索难度和成本就越大。
  • 降维能够增加数据的可读性,利于发掘数据的有意义的结构。

3. 降维的主要作用

(1)减少冗余特征,降低数据维度

  假设我们有两个特征: x 1 x_1 x1​:长度用厘米表示的身高; x 2 x_2 x2​:是用英寸表示的身高。这两个分开的特征 x 1 x_1 x1​和 x 2 x_2 x2​,实际上表示的内容相同,这样其实可以减少数据到一维,只有一个特征表示身高就够了。很多特征具有线性关系,具有线性关系的特征很多都是几余的特征,去掉冗余特征对机器学习的计算结果不会有影响。

(2)数据可视化

  t-distributed Stochastic Neighbor Embedding(t-SNE)将数据点之间的相似度转换为概率。原始空间中的相似度由高斯联合概率表示,嵌入空间的相似度由“学生t分布”表示。虽然Isomap,LLE和variants等数据降维和可视化方法,更适合展开单个连续的低维的manifold。但如果要准确的可视化样本间的相似度关系,如对于下图所示的S曲线(不同颜色的图像表示不同类别的数据),t-SNE表现更好。因为t-SNE主要是关注数据的局部结构。

4. 降维的优缺点

降维的优点:

  • 通过减少特征的维数,数据集存储所需的空间也相应减少,减少了特征维数所需的计算训练时间;
  • 数据集特征的降维有助于快速可视化数据;
  • 通过处理多重共线性消除冗余特征。

降维的缺点:

  • 由于降维可能会丢失一些数据;
  • 在主成分分析(PCA)降维技术中,有时需要考虑多少主成分是难以确定的,往往使用经验法则。

二、主成分与方差

  顾名思义,PCA的含义是将高维数据中的主要成分找出来。设原始数据 D = { x 1 , x 2 , ⋯   , x m } \mathcal D=\{\boldsymbol x_1,\boldsymbol x_2,\cdots,\boldsymbol x_m\} D={x1​,x2​,⋯,xm​},其中 x i ∈ R d \boldsymbol x_i\in\mathbb R^d xi​∈Rd。我们的目标是通过某种变换 f : R d → R k f:\mathbb R^d\rightarrow\mathbb R^k f:Rd→Rk,将数据的维度从 d d d减小到 k k k,且通常 k ≪ d k\ll d k≪d。变换后的数据不同维度之间线性无关,这些维度就称为主成分。也就是说,如果将变换后的数据排成矩阵 Z = ( z 1 , ⋯   , z m ) T \boldsymbol Z=(\boldsymbol z_1,\cdots,\boldsymbol z_m)^{\mathrm T} Z=(z1​,⋯,zm​)T,其中 z i = f ( x i ) \boldsymbol z_i=f(\boldsymbol x_i) zi​=f(xi​),那么 Z \boldsymbol Z Z的列向量是线性无关的。从矩阵的知识可以立即得到 k < m k<m k<m。

  PCA算法不仅希望变换后的数据特征线性无关,更要求这些特征之间相互独立,也即任意两个不同特征之间的协方差为零。因此,PCA算法在计算数据的主成分时,会从第一个主成分开始依次计算,并保证每个主成分与之前的所有主成分都是正交的,直到选取了预先设定的 k k k个主成分为止。关于主成分选取的规则,我们以一个平面上的二元高斯分布为例来具体说明。如图2(a)所示,数据点服从以 ( 3 , 3 ) (3,3) (3,3)为中心、 x 1 x_1 x1​方向方差为1.5、 x 2 x_2 x2​方向方差为0.5、协方差为0.8的二元高斯分布。从图中可以明显看出,当前的 x 1 x_1 x1​和 x 2 x_2 x2​两个特征之间存在关联,并不是独立的。

在这里插入图片描述

图2 一个二元高斯分布和两个维度的投影情况

  图2(b)展示了数据分别向 x 1 x_1 x1​和 x 2 x_2 x2​方向投影的结果。由于 x 1 x_1 x1​方向方差更大,数据在该方向的投影分布也更广。也就是说,数据向 x 1 x_1 x1​方向的投影保留了更多信息, x 1 x_1 x1​是一个比 x 2 x_2 x2​更好的特征。如果我们不再把目光局限于 x 1 x_1 x1​和 x 2 x_2 x2​两个方向,而是考虑平面上的任意一个方向,那么如左图中红色虚线所示,沿直线 x 2 = 0.557 x 1 + 1.330 x_2=0.557x_1+1.330 x2​=0.557x1​+1.330 的方向数据的方差是最大的,我们应当将此方向作为数据的特征之一。既然PCA算法希望选出的主成分保留最多的信息,那么就应当不断选取数据方差最大的方向作为主成分。因此,PCA计算的过程就是依次寻找方差最大方向的过程。

  为了方便计算,我们通常要先对数据进行中心化,把当前每个特征都变为均值为0的分布。用 x i ( j ) x_i^{(j)} xi(j)​表示数据 x i \boldsymbol x_i xi​的第 j j j个维度,那么均值 μ j \mu_j μj​为 μ j = 1 m ∑ i = 1 m x i ( j ) , j = 1 , ⋯   , d \mu_j=\frac{1}{m}\sum_{i=1}^mx_i^{(j)}, \quad j=1,\cdots,d μj​=m1​i=1∑m​xi(j)​,j=1,⋯,d 将数据中心化,得到 x i ( j ) ← x i ( j ) − μ j , i = 1 , ⋯   , m x_i^{(j)}\leftarrow x_i^{(j)}-\mu_j, \quad i=1,\cdots,m xi(j)​←xi(j)​−μj​,i=1,⋯,m

  以图2(a)的高斯分布为例,变换后的图像如图3(b)所示,其中心点已经从 ( 3 , 3 ) (3,3) (3,3)变为了 ( 0 , 0 ) (0,0) (0,0)。接下来在讨论时,我们默认数据都已经被中心化。下面,我们就来具体讲解如何寻找数据中方差最大的方向。

在这里插入图片描述

图3 数据的中心化

三、利用特征分解进行PCA

  为了找到方差最大的方向,我们先来计算样本在某个方向上的投影。设 u \boldsymbol u u为方向向量,满足 ∥ u ∥ = 1 \Vert\boldsymbol u\Vert=1 ∥u∥=1,向量 x \boldsymbol x x在方向 u \boldsymbol u u上的投影为 x T u \boldsymbol x^{\mathrm T}\boldsymbol u xTu。于是所有样本在方向 u \boldsymbol u u上的方差为
σ u = 1 m ∑ i = 1 m ( x i T u ) 2 = 1 m ∑ i = 1 m u T x i x i T u = u T ( 1 m ∑ i = 1 m x i x i T ) u = u T Σ u \begin{aligned} \sigma_{\boldsymbol u} &= \frac{1}{m}\sum_{i=1}^m(\boldsymbol x_i^{\mathrm T}\boldsymbol u)^2 \\ &=\frac{1}{m}\sum_{i=1}^m\boldsymbol u^{\mathrm T}\boldsymbol x_i\boldsymbol x_i^{\mathrm T}\boldsymbol u \\ &= \boldsymbol u^{\mathrm T}\left(\frac{1}{m}\sum_{i=1}^m\boldsymbol x_i\boldsymbol x_i^{\mathrm T}\right)\boldsymbol u \\ &= \boldsymbol u^{\mathrm T}\boldsymbol\Sigma\boldsymbol u \end{aligned} σu​​=m1​i=1∑m​(xiT​u)2=m1​i=1∑m​uTxi​xiT​u=uT(m1​i=1∑m​xi​xiT​)u=uTΣu​

  矩阵 Σ = 1 m ∑ i = 1 m x i x i T ∈ R d × d \boldsymbol\Sigma=\frac{1}{m}\sum\limits_{i=1}^m\boldsymbol x_i\boldsymbol x_i^{\mathrm T} \in \mathbb R^{d\times d} Σ=m1​i=1∑m​xi​xiT​∈Rd×d 称为样本的协方差矩阵。由于 m m m是常数,简单起见,下面省略因子 1 m \frac{1}{m} m1​。要令方差最大,就相当于求解下面的优化问题 max ⁡ u u T Σ u s.t. ∥ u ∥ = 1 \max_{\boldsymbol u}\boldsymbol u^{\mathrm T}\boldsymbol\Sigma\boldsymbol u \quad \text{s.t.}\Vert\boldsymbol u\Vert=1 umax​uTΣus.t.∥u∥=1

  在求解该问题之前,我们先来介绍矩阵的特征值(eigenvalue)与特征分解(eigendecomposition)相关的知识。对于方阵 A ∈ R n × n \boldsymbol A\in\mathbb R^{n\times n} A∈Rn×n,如果向量 ξ ∈ R n \boldsymbol\xi\in\mathbb R^n ξ∈Rn 与实数 λ \lambda λ满足 A ξ = λ ξ \boldsymbol A\boldsymbol\xi=\lambda\boldsymbol\xi Aξ=λξ,就称 λ \lambda λ是矩阵 A \boldsymbol A A的特征值, ξ \boldsymbol\xi ξ为矩阵的特征向量。特征向量 x \boldsymbol x x的任意非零实数倍 r x r\boldsymbol x rx满足 A ( r x ) = r ( A x ) = r ( λ x ) = λ ( r x ) \boldsymbol A(r\boldsymbol x)=r(\boldsymbol A\boldsymbol x)=r(\lambda\boldsymbol x)=\lambda(r\boldsymbol x) A(rx)=r(Ax)=r(λx)=λ(rx) 因而 r x r\boldsymbol x rx也还是特征向量。简单来说,矩阵 A \boldsymbol A A作用到其特征向量 x \boldsymbol x x上的结果等价于对该向量做伸缩变换,伸缩的倍数等于特征值 λ \lambda λ,但不改变向量所在的直线。从这个角度来看, r x r\boldsymbol x rx与 x \boldsymbol x x共线, r x r\boldsymbol x rx自然就是特征向量。因此,我们通常更关心单位特征向量,即长度为1的特征向量。例如,矩阵
A = ( 2 1 0 1 ) \boldsymbol A=\begin{pmatrix}2 & 1 \\ 0 &1\\ \end{pmatrix} A=(20​11​) 的特征值及其对应的单位特征向量有两组,分别为 λ 1 = 1 \lambda_1=1 λ1​=1, ξ 1 = ( − 1 , 1 ) T \xi_1=(-1,1)^{\mathrm T} ξ1​=(−1,1)T 和 λ 2 = 2 \lambda_2=2 λ2​=2, ξ 2 = ( 1 , 0 ) T \xi_2=(1,0)^{\mathrm T} ξ2​=(1,0)T,大家可以自行计算验证。显然,对角矩阵 diag ( d 1 , ⋯   , d n ) \text{diag}(d_1,\cdots,d_n) diag(d1​,⋯,dn​)的特征值就等于其对角线上的元素 d 1 , ⋯   , d n d_1,\cdots,d_n d1​,⋯,dn​。特别地, n n n阶单位阵 I n \boldsymbol I_n In​的特征值是1,且空间中的所有向量都是该特征值对应的特征向量。

  我们可以从线性变换的角度来理解矩阵的特征值和特征向量。一个 n n n阶方阵 A \boldsymbol A A可以看作是 n n n维空间中的变换,将 n n n维向量 x \boldsymbol x x变为另一个 n n n维向量 A x \boldsymbol A\boldsymbol x Ax。由于 A ( x 1 + x 2 ) = A x 1 + A x 2 \boldsymbol A(\boldsymbol x_1+\boldsymbol x_2)=\boldsymbol A\boldsymbol x_1+\boldsymbol A\boldsymbol x_2 A(x1​+x2​)=Ax1​+Ax2​ 以及 A ( r x 1 ) = r A x 1 \boldsymbol A(r\boldsymbol x_1)=r\boldsymbol A\boldsymbol x_1 A(rx1​)=rAx1​ 对任意向量 x 1 \boldsymbol x_1 x1​、 x 2 \boldsymbol x_2 x2​和实数 r r r成立,该变换是一个线性变换。事实上,当坐标轴确定之后,线性变换与矩阵之间有一一对应关系。如果把向量都看作从原点指向向量所表示坐标的有向线段,可以证明,任意一个线性变换总可以分解成绕原点旋转和长度伸缩的组合。因此,矩阵与向量相乘 A x \boldsymbol A\boldsymbol x Ax可以理解为将向量 x \boldsymbol x x旋转某一角度,再将长度变为某一倍数。而矩阵的特征向量,就是在该变换作用下只伸缩不旋转的向量,并且其对应的特征值就是伸缩的倍数。当特征值 λ > 0 \lambda>0 λ>0 时,变换后的向量与原向量方向相同; λ < 0 \lambda<0 λ<0 时方向相反; λ = 0 \lambda=0 λ=0 则表示矩阵会把所有该直线上的向量压缩为零向量。

  从这一角度继续,我们可以引入正定矩阵(positive definite)和半正定矩阵(positive semidefinite)的概念。如果对任意非零向量 x \boldsymbol x x,都有 x T A x ≥ 0 \boldsymbol x^{\mathrm T}\boldsymbol A\boldsymbol x\ge0 xTAx≥0,就称 A \boldsymbol A A为半正定矩阵;如果严格有 x T A x > 0 \boldsymbol x^{\mathrm T}\boldsymbol A\boldsymbol x>0 xTAx>0,就称 A \boldsymbol A A为正定矩阵。如果将 y = A x \boldsymbol y=\boldsymbol A\boldsymbol x y=Ax 看作变换后的向量,那么 x T y ≥ 0 \boldsymbol x^{\mathrm T}y\ge0 xTy≥0 就表示 x \boldsymbol x x与 y \boldsymbol y y之间的夹角小于等于90°。由此,我们可以立即得到半正定矩阵的所有特征值都非负,否则负特征值会使特征向量在变换后反向,与原向量夹角为180°,产生矛盾。

  为什么我们要引入半正定矩阵呢?在上面的优化问题中,我们得到的目标函数是 u T Σ u \boldsymbol u^{\mathrm T}\boldsymbol\Sigma\boldsymbol u uTΣu,其中 Σ \boldsymbol\Sigma Σ是协方差矩阵。事实上,该矩阵一定是半正定矩阵。设 y \boldsymbol y y是任一非零向量,那么 y T Σ y = y T ( ∑ i = 1 m x i x i T ) y = ∑ i = 1 m y T x i x i T y = ∑ i = 1 m ( x i T y ) 2 ≥ 0 \boldsymbol y^{\mathrm T}\boldsymbol\Sigma\boldsymbol y=\boldsymbol y^{\mathrm T}\left(\sum_{i=1}^m\boldsymbol x_i\boldsymbol x_i^{\mathrm T}\right)\boldsymbol y=\sum_{i=1}^m\boldsymbol y^{\mathrm T}\boldsymbol x_i\boldsymbol x_i^{\mathrm T}\boldsymbol y=\sum_{i=1}^m(\boldsymbol x_i^{\mathrm T}\boldsymbol y)^2 \ge 0 yTΣy=yT(i=1∑m​xi​xiT​)y=i=1∑m​yTxi​xiT​y=i=1∑m​(xiT​y)2≥0

  根据定义, Σ \boldsymbol\Sigma Σ是半正定矩阵。因此,我们可以用半正定矩阵的一些性质来帮助我们求解该优化问题。这里,我们不加证明地给出一条重要性质:对于 d d d阶对称半正定矩阵 Σ \boldsymbol\Sigma Σ,总可以找到它的 d d d个单位特征向量 e 1 , ⋯   , e d \boldsymbol e_1,\cdots,\boldsymbol e_d e1​,⋯,ed​,使得这些特征向量是两两正交的,也即对任意 1 ≤ i 1\le i 1≤i, j ≤ d j\le d j≤d,都有 ∥ e i ∥ = 1 , e i T e j = { 1 , i = j 0 , i ≠ j \Vert\boldsymbol e_i\Vert=1, \quad \boldsymbol e_i^{\mathrm T}\boldsymbol e_j=\begin{cases}1, \quad i=j \\ 0, \quad i \ne j\end{cases} ∥ei​∥=1,eiT​ej​={1,i=j0,i=j​ 有兴趣的可以在线性代数的相关资料中找到该性质的证明。利用该性质,记 Q = ( e 1 , ⋯   , e d ) ∈ R d × d \boldsymbol Q=(\boldsymbol e_1,\cdots,\boldsymbol e_d)\in\mathbb R^{d\times d} Q=(e1​,⋯,ed​)∈Rd×d,有 ( Q T Q ) i j = e i T e j = I ( i = j ) (\boldsymbol Q^{\mathrm T}\boldsymbol Q)_{ij}=\boldsymbol e_i^{\mathrm T}\boldsymbol e_j=\mathbb I(i=j) (QTQ)ij​=eiT​ej​=I(i=j) 于是, Q T Q \boldsymbol Q^{\mathrm T}\boldsymbol Q QTQ对角线上的元素全部为1,其他元素全部为0,恰好是单位矩阵 I d \boldsymbol I_d Id​。因此, Q \boldsymbol Q Q的逆矩阵就是其转置, Q − 1 = Q T \boldsymbol Q^{-1}=\boldsymbol Q^{\mathrm T} Q−1=QT。因为组成该矩阵的向量是相互正交的,我们将这样的矩阵称为正交矩阵(orthogonal matrix)。根据逆矩阵的性质,我们有 I d = Q T Q = Q Q T = ( e 1 , ⋯   , e d ) ( e 1 T ⋮ e d T ) = ∑ i = 1 d e i e i T \boldsymbol I_d=\boldsymbol Q^{\mathrm T}\boldsymbol Q=\boldsymbol Q\boldsymbol Q^{\mathrm T}=(\boldsymbol e_1,\cdots,\boldsymbol e_d)\left(\begin{matrix}\boldsymbol e_1^{\mathrm T} \\ \vdots \\ \boldsymbol e_d^{\mathrm T}\end{matrix}\right)=\sum_{i=1}^d\boldsymbol e_i\boldsymbol e_i^{\mathrm T} Id​=QTQ=QQT=(e1​,⋯,ed​) ​e1T​⋮edT​​ ​=i=1∑d​ei​eiT​ 设特征向量 e i \boldsymbol e_i ei​对应的特征值是 λ i \lambda_i λi​,那么矩阵 Σ \boldsymbol\Sigma Σ可以写为
Σ = Σ I d = Σ ∑ i = 1 d e i e i T = ∑ i = 1 d ( Σ e i ) e i T = ∑ i = 1 d λ i e i e i T = ( e 1 , e 2 , ⋯   , e d ) ( λ 1 0 ⋯ 0 0 λ 2 ⋯ 0 ⋮ ⋮ ⋮ 0 0 ⋯ λ d ) ( e 1 T e 2 T ⋮ e d T ) = Q Λ Q T \begin{aligned} \boldsymbol\Sigma &= \boldsymbol\Sigma\boldsymbol I_d = \boldsymbol\Sigma\sum_{i=1}^d\boldsymbol e_i\boldsymbol e_i^{\mathrm T} = \sum_{i=1}^d(\boldsymbol\Sigma\boldsymbol e_i)\boldsymbol e_i^{\mathrm T} = \sum_{i=1}^d\lambda_i\boldsymbol e_i\boldsymbol e_i^{\mathrm T} \\[3ex] &= (\boldsymbol e_1,\boldsymbol e_2,\cdots,\boldsymbol e_d)\begin{pmatrix}\lambda_1 & 0 &\cdots & 0 \\ 0 & \lambda_2 & \cdots & 0 \\ \vdots & \vdots & & \vdots \\ 0 & 0 & \cdots & \lambda_d\end{pmatrix}\begin{pmatrix}\boldsymbol e_1^{\mathrm T} \\ \boldsymbol e_2^{\mathrm T} \\ \vdots \\ \boldsymbol e_d^{\mathrm T}\end{pmatrix} \\[2ex] &= \boldsymbol Q\boldsymbol\Lambda\boldsymbol Q^{\mathrm T} \end{aligned} Σ​=ΣId​=Σi=1∑d​ei​eiT​=i=1∑d​(Σei​)eiT​=i=1∑d​λi​ei​eiT​=(e1​,e2​,⋯,ed​) ​λ1​0⋮0​0λ2​⋮0​⋯⋯⋯​00⋮λd​​ ​e1T​e2T​⋮edT​​ ​=QΛQT​ 其中, Λ = diag ( λ 1 , ⋯   , λ d ) \boldsymbol\Lambda=\text{diag}(\lambda_1,\cdots,\lambda_d) Λ=diag(λ1​,⋯,λd​) 是由特征向量所对应的特征值依次排列而成的对角矩阵。上式表明,一个半正定矩阵可以分解成3个矩阵的乘积,其中 Q \boldsymbol Q Q是其正交的特征向量构成的正交矩阵, Λ \boldsymbol\Lambda Λ是其特征值构成的对角矩阵,这样的分解就称为矩阵的特征分解。

  利用特征分解,我们可以很容易地计算 u T Σ u \boldsymbol u^{\mathrm T}\boldsymbol\Sigma\boldsymbol u uTΣu的值。首先,由于 e 1 , ⋯   , e d \boldsymbol e_1,\cdots,\boldsymbol e_d e1​,⋯,ed​是维空间中的个正交向量, u \boldsymbol u u一定可以唯一表示成这些向量的线性组合,即 u = ∑ i = 1 d α i e i \boldsymbol u=\sum\limits_{i=1}^d\alpha_i\boldsymbol e_i u=i=1∑d​αi​ei​,其中系数 α i \alpha_i αi​等于向量 u \boldsymbol u u在 e i \boldsymbol e_i ei​方向上的投影长度。我们可以将这组向量想象成相互垂直的坐标轴,而 ( α 1 , ⋯   , α d ) (\alpha_1,\cdots,\alpha_d) (α1​,⋯,αd​)就相当于 u \boldsymbol u u在这组坐标轴下的坐标。这样, Q T u \boldsymbol Q^{\mathrm T}\boldsymbol u QTu可以化为
Q T u = Q T ∑ i = 1 d α i e i = ( ∑ i = 1 d α i e 1 T e i ⋮ ∑ i = 1 d α i e d T e i ) = ( α 1 , ⋯   , α d ) = α \begin{aligned} \boldsymbol Q^{\mathrm T}\boldsymbol u=\boldsymbol Q^{\mathrm T}\sum_{i=1}^d\alpha_i\boldsymbol e_i = \begin{pmatrix}\sum\limits_{i=1}^d\alpha_i\boldsymbol e_1^{\mathrm T}\boldsymbol e_i \\ \vdots \\ \sum\limits_{i=1}^d\alpha_i\boldsymbol e_d^{\mathrm T}\boldsymbol e_i\end{pmatrix} = (\alpha_1,\cdots,\alpha_d)=\boldsymbol\alpha \end{aligned} QTu=QTi=1∑d​αi​ei​= ​i=1∑d​αi​e1T​ei​⋮i=1∑d​αi​edT​ei​​ ​=(α1​,⋯,αd​)=α​ 这里,我们用到了 e i T e j = I ( i = j ) \boldsymbol e_i^{\mathrm T}\boldsymbol e_j=\mathbb I(i=j) eiT​ej​=I(i=j) 的性质,把求和中 e i T e j \boldsymbol e_i^{\mathrm T}\boldsymbol e_j eiT​ej​都消去了。接下来, u T Σ u \boldsymbol u^{\mathrm T}\boldsymbol\Sigma\boldsymbol u uTΣu可以计算得: u T u = u T Q Λ Q T u = α T Λ α = ∑ i = 1 d λ i α i 2 \boldsymbol u^{\mathrm T}\boldsymbol u=\boldsymbol u^{\mathrm T}\boldsymbol Q\boldsymbol\Lambda\boldsymbol Q^{\mathrm T}\boldsymbol u=\boldsymbol\alpha^{\mathrm T}\boldsymbol\Lambda\boldsymbol\alpha=\sum_{i=1}^d\lambda_i\alpha_i^2 uTu=uTQΛQTu=αTΛα=i=1∑d​λi​αi2​

  在上面的优化问题 max ⁡ u u T Σ u   s.t. ∥ u ∥ = 1 \max\limits_{\boldsymbol u}\boldsymbol u^{\mathrm T}\boldsymbol\Sigma\boldsymbol u \ \ \text{s.t.}\Vert\boldsymbol u\Vert=1 umax​uTΣu  s.t.∥u∥=1 中,我们还要求 u \boldsymbol u u是方向向量,即 ∥ u ∥ = 1 \Vert\boldsymbol u\Vert=1 ∥u∥=1。这一要求给系数 α i \alpha_i αi​添加了限定条件:
∑ i = 1 d α i 2 = ∑ i = 1 d ( u T e i ) 2 = ∑ i = 1 d u T e i e i T u = u T ( ∑ i = 1 d e i e i T ) u = u T I d u = u T u = ∥ u ∥ 2 = 1 \begin{aligned} \sum_{i=1}^d\alpha_i^2 &= \sum_{i=1}^d(\boldsymbol u^{\mathrm T}\boldsymbol e_i)^2=\sum_{i=1}^d\boldsymbol u^{\mathrm T}\boldsymbol e_i\boldsymbol e_i^{\mathrm T}\boldsymbol u \\ &= \boldsymbol u^{\mathrm T}\left(\sum_{i=1}^d\boldsymbol e_i\boldsymbol e_i^{\mathrm T}\right)\boldsymbol u \\ &= \boldsymbol u^{\mathrm T}\boldsymbol I_d\boldsymbol u \\ &= \boldsymbol u^{\mathrm T}\boldsymbol u = \Vert\boldsymbol u\Vert^2=1 \end{aligned} i=1∑d​αi2​​=i=1∑d​(uTei​)2=i=1∑d​uTei​eiT​u=uT(i=1∑d​ei​eiT​)u=uTId​u=uTu=∥u∥2=1​ 因此,原优化问题等价于 max ⁡ u ∑ i = 1 d λ i α i 2 , s.t. ∑ i = 1 d α i 2 = 1 \max_{\boldsymbol u}\sum_{i=1}^d\lambda_i\alpha_i^2, \quad \text{s.t.} \quad \sum_{i=1}^d\alpha_i^2=1 umax​i=1∑d​λi​αi2​,s.t.i=1∑d​αi2​=1

  由于 Σ \boldsymbol\Sigma Σ是半正定矩阵,其特征值 λ i \lambda_i λi​必定非负,该问题的解就是 Σ \boldsymbol\Sigma Σ的最大特征值 λ max ⁡ \lambda_{\max} λmax​,不妨设其为 λ u \lambda_u λu​。为了使上式取到最大值,应当有 α u = u T e u = 1 \alpha_{\boldsymbol u}=\boldsymbol u^{\mathrm T}\boldsymbol e_{\boldsymbol u}=1 αu​=uTeu​=1,且其他的 α i \alpha_i αi​全部为零。因此,使方差最大的方向就是该特征值 λ u \lambda_{\boldsymbol u} λu​对应的特征向量 e u \boldsymbol e_{\boldsymbol u} eu​的方向。

  上面的计算结果表面,用PCA寻找第一个主成分的过程就是求解PCA最大特征值及其对应特征向量的过程。事实上,由于协方差矩阵 Σ \boldsymbol\Sigma Σ半正定,其所有特征向量正交,恰好也满足我们最开始“每个主成分都与之前的所有主成分正交”的要求。如果排除掉第一主成分,第二主成分就对应 Σ \boldsymbol\Sigma Σ第二大特征值的特征向量,依次类推。因此,如果我们要把 d d d维的数据降到 k k k维,只需要计算 Σ \boldsymbol\Sigma Σ最大的 k k k个特征值对应的特征向量即可。设这 k k k个特征向量依次为 e 1 , ⋯   , e k \boldsymbol e_1,\cdots,\boldsymbol e_k e1​,⋯,ek​,矩阵 W = ( e 1 , ⋯   , e k ) ∈ R d × d \boldsymbol W=(\boldsymbol e_1,\cdots,\boldsymbol e_k)\in\mathbb R^{d\times d} W=(e1​,⋯,ek​)∈Rd×d,原数据矩阵 X = ( x 1 ⋯   , x m ) T ∈ R d × d \boldsymbol X=(\boldsymbol x_1\cdots,\boldsymbol x_m)^{\mathrm T}\in\mathbb R^{d\times d} X=(x1​⋯,xm​)T∈Rd×d,那么降维后的数据为 PCA ( X ) = X W \text{PCA}(\boldsymbol X)=\boldsymbol X\boldsymbol W PCA(X)=XW

四、动手实现PCA算法

  下面,我们在NumPy库中线性代数工具的帮助下实现PCA算法。首先,我们导入数据集PCA_dataset.csv并将其可视化。数据集的每一行包含两个数 x 1 x_1 x1​与 x 2 x_2 x2​,代表平面上点的坐标 ( x 1 , x 2 ) (x_1,x_2) (x1​,x2​)。

import numpy as np
import matplotlib.pyplot as plt

# 导入数据集
data = np.loadtxt('PCA_dataset.csv', delimiter=',')
print('数据集大小:', len(data))

# 可视化
plt.figure()
plt.scatter(data[:, 0], data[:, 1], color='blue', s=10)
plt.axis('square')
plt.ylim(-2, 8)
plt.grid()
plt.xlabel(r'$x_1$')
plt.ylabel(r'$x_2$')
plt.show()

在这里插入图片描述

  接下来,我们按上面讲解的方式实现PCA算法。numpy.linalg中提供了许多线性代数相关的工具,可以帮助我们计算矩阵的特征值与特征向量。

def pca(X, k):
    d, m = X.shape
    if d < k:
        print('k应该小于特征数')
        return X, None

    # 中心化
    X = X - np.mean(X, axis=0)
    # 计算协方差矩阵
    cov = X.T @ X
    # 计算特征值和特征向量
    eig_values, eig_vectors = np.linalg.eig(cov)
    # 获取最大的k个特征值的下标
    idx = np.argsort(-eig_values)[:k]
    # 对应的特征向量
    W = eig_vectors[:, idx]
    # 降维
    X = X @ W
    return X, W

  最后,我们在数据集上测试该PCA函数的效果,并将变换后的数据绘制出来。由于原始数据只有二维,为了演示PCA的效果,我们仍然设置 k = 2 k=2 k=2,不进行降维。但是,从结果中仍然可以看出PCA计算的主成分方向。相比于原始数据,变换后的数据最“长”的方向变成了横轴的方向。

X, W = pca(data, 2)
print('变换矩阵:\n', W)

# 绘图
plt.figure()
plt.scatter(X[:, 0], X[:, 1], color='blue', s=10)
plt.axis('square')
plt.ylim(-5, 5)
plt.grid()
plt.xlabel(r'$x_1$')
plt.ylabel(r'$x_2$')
plt.show()

在这里插入图片描述

五、用sklearn实现PCA算法

  Sklearn库中同样提供了实现好的PCA算法,我们可以直接调用它来完成PCA变换。可以看出,虽然结果图像与我们上面直接实现的版本有180°的旋转,变换矩阵的元素也互为相反数,但PCA本质上只需要找到主成分的方向,因此两者得到的结果是等价的。

使用scikit-learn库中decomposition模块的PCA类可以创建PCA模型,其基本语法格式和参数说明如下。

sklearn.decomposition.PCA(n_components=None, *, copy=True, whiten=False, svd_solver='auto', tol=0.0, iterated_power='auto', n_oversamples=10, power_iteration_normalizer='auto', random_state=None)

参数名称参数说明
n_components接收int或str。表示所要保留的主成分个数n,即保留下来的特征个数n,赋值为int时,表示降维的维度,如n_components=1,将把原始数据降到一个维度。赋值为str时,表示降维的模式,如取值为’mle’时,将自动选取特征个数n,使得满足所要求的方差百分比。默认为None。
copy接收bool。表示是否在运行算法时,将原始训练数据复制一份。若为True,则运行后,原始训练数据的值不会有任何改变,因为是在原始数据的副本上进行运算;若为False,则运行后,原始训练数据的值会发生改变。默认为True。
whiten接收bool。表示是否白化,使得每个特征具有相同的方差。默认为False。
svd_solver接收{‘auto’, ‘full’, ‘arpack’, ‘randomized’},用于奇异值分解的算法。‘auto’会根据输入数据选择最佳的求解器。默认为’auto’。
tol接收float。奇异值的容忍度。在使用’arpack’求解器时使用。默认为0.0。
iterated_power接收int或’auto’。在使用’randomized’求解器时,幂法迭代的次数。默认为’auto’。
n_oversamples接收int。随机SVD的过采样数量。默认值为10。
power_iteration_normalizer接收{‘auto’, ‘LU’, ‘none’}。在幂迭代中使用的归一化方法。默认为’auto’。
random_state接收int、RandomState实例或None。控制估计器的随机性。如果设置为固定整数,则可以确保结果的可重复性。默认为None。
from sklearn.decomposition import PCA

# 中心化
X = data - np.mean(data, axis=0)
pca_res = PCA(n_components=2).fit(X)
W = pca_res.components_.T
print ('sklearn计算的变换矩阵:\n', W)
X_pca = X @ W

# 绘图
plt.figure()
plt.scatter(X_pca[:, 0], X_pca[:, 1], color='blue', s=10)
plt.axis('square')
plt.ylim(-5, 5)
plt.grid()
plt.xlabel(r'$x_1$')
plt.ylabel(r'$x_2$')
plt.show()

在这里插入图片描述

六、拓展:奇异值分解

  本文介绍了数据降维的常用算法之一:PCA算法。数据降维是无监督学习的重要问题,在机器学习中有广泛的应用。由于从现实生活中采集的数据往往非常复杂,包含大量的冗余信息,通常我们必须对其进行降维,选出有用的特征供给后续模型使用。此外,有时我们还希望将高维数据可视化,也需要从数据中挑选2到3个最有价值的维度,将数据投影后绘制出来。除了PCA之外,现在常用的降维算法还有线性判别分析(linear discriminant analysis,LDA)、一致流形逼近与投影(uniform manifold approximation and projection,UMAP)、t分布随机近邻嵌入(t-distributed stochastic neighbor embedding,t-SNE)等。这些算法的特点各不相同,也有不同的适用场景。

  关于PCA的计算方式,由于计算协方差矩阵 Σ = X X T \boldsymbol\Sigma=\boldsymbol X\boldsymbol X^{\mathrm T} Σ=XXT 和特征分解的时间复杂度较高,实践中通常会采用矩阵的奇异值分解(singular value decomposition,SVD)来代替特征分解。上面我们用到的sklearn库中的PCA算法就是以奇异值分解为基础的实现的。相比于特征分解,奇异值分解不需要实际计算出 Σ \boldsymbol\Sigma Σ,并且存在更高效的迭代求解方法。感兴趣的可以查阅相关资料,了解特征分解和奇异值分解的异同。

在这里插入图片描述
1. 奇异值分解的构造

在这里插入图片描述

在这里插入图片描述

在这里插入图片描述
2. 奇异值分解用于数据压缩

在这里插入图片描述

在这里插入图片描述

在这里插入图片描述
3. 奇异值分解的几何解释

在这里插入图片描述

在这里插入图片描述

在这里插入图片描述

4. 奇异值分解——SVD与PCA的关系

  实际上,如果将矩阵 A m × n \boldsymbol A_{m\times n} Am×n​看作是一个样本集合,其中的行看作特征随机变量,列看作每一个样本。当对数据集进行规范化后,矩阵 A T A \boldsymbol A^{\mathrm T}\boldsymbol A ATA就是样本集合的协方差矩阵。这样,SVD分解后的右奇异矩阵就是PCA分析中的特征向量组成的矩阵。

  最后,大家可以在SETOSA网站的PCA算法动态展示平台上观察PCA的结果随数据分布的变化,加深对算法的理解。

:以上文中的数据集及相关资源下载地址:
链接:https://pan.quark.cn/s/15d427a80558
提取码:vYkb

标签:机器,特征向量,boldsymbol,矩阵,学习,降维,PCA,mathrm
From: https://blog.csdn.net/Morse_Chen/article/details/142647797

相关文章

  • Python从0到100(六十一):机器学习实战-实现客户细分
    一、导入数据在此项目中,我们使用UCI机器学习代码库中的数据集。该数据集包含关于来自多种产品类别的各种客户年度消费额(货币单位计价)的数据。该项目的目标之一是准确地描述与批发商进行交易的不同类型的客户之间的差别。这样可以使分销商清晰地了解如何安排送货服务,以便......
  • [rCore学习笔记 028] Rust 中的动态内存分配
    引言想起我们之前在学习C的时候,总是提到malloc,总是提起,使用malloc现场申请的内存是属于堆,而直接定义的变量内存属于栈.还记得当初学习STM32的时候CubeIDE要设置stack和heap的大小.但是我们要记得,这么好用的功能,实际上是操作系统在负重前行.那么为了实现动态内存分配功......
  • 深度学习(计算数据集均值标准差)
      深度学习中有些数据集可能不符合imagenet计算出的均值和标准差,需要根据自己的数据集单独计算。下面这个脚本能够计算当前数据集均值和标准差。 importtorchimportosfromPILimportImagefromtorchvisionimporttransforms#trans=transforms.Compose([#......
  • 【学术】过来人对研究生阶段的学习建议
    【学术】过来人对研究生阶段的学习建议不同阶段的学习模式高中—>本科—>研究生是不同的学习阶段,面对不同的生存学习环境,需要做出调整。诚然,一些科研能力强的研究生会产生一种如鱼得水的感觉——不仅适应良好,且成长很快。两大法则:硕博生和本科生最大的区别在于独立处理......
  • 有监督学习&无监督学习
    有监督学习&无监督学习有监督学习(SupervisedLearning)和无监督学习(UnsupervisedLearning)是机器学习中的两种主要方法,它们在目标、数据使用和应用场景上有显著的区别:目标不同:有监督学习:目的是通过训练数据集(包含输入特征和对应的标签)来学习一个模型,以便对新的、未见......
  • Elasticsearch学习笔记(3)
    RestAPIElasticsearch(ES)官方提供了多种语言的客户端库,用于与Elasticsearch进行交互。这些客户端库的主要功能是帮助开发者更方便地构建和发送DSL(DomainSpecificLanguage)查询语句,并通过HTTP请求与Elasticsearch集群进行通信。官方文档地址:https://www.elastic.co/guide/en/......
  • Python电能质量扰动信号分类(六)基于扰动信号特征提取的超强机器学习识别模型
    创新度高!!!需要发论文的同学即买即用 往期精彩内容:Python-电能质量扰动信号数据介绍与分类-CSDN博客Python电能质量扰动信号分类(一)基于LSTM模型的一维信号分类-CSDN博客Python电能质量扰动信号分类(二)基于CNN模型的一维信号分类-CSDN博客Python电能质量扰动信号分类(......
  • MySql学习笔记:什么是数据库?
    数据库的概念:         数据库(Database),简而言之可视为数字化的文件柜,是一个长期储存在计算机内有组织的、统一管理的数据集合,用于存储和管理大量相关信息。        数据库是一个按数据的结构来存储和管理数据的计算机系统,也就是说,数据库通常有两方面含义:......
  • 基于深度学习的跌倒检测系统:YOLOv5/v6/v7/v8/v10模型实现与UI界面集成、数据集
    1.引言跌倒检测是一个重要的研究领域,尤其在老年人和病人监护中,及时检测并响应跌倒事件可以大大减少伤害和死亡的风险。本博客将介绍如何构建一个基于深度学习的跌倒检测系统,使用YOLOv5进行目标检测,并设计一个用户界面(UI)来实时监控和反馈。本文将详细描述系统的各个组成部分......
  • 【Linux学习】基本指令其二
    前言【Linux学习】基本指令其一第一篇中介绍Linux中最基本的命令,现在继续说明后续指令。一、man指令Linux的命令有很多参数,可以通过查看联机手册获取帮助。man命令用于查看Linux系统中的手册页,提供了关于各种命令和程序的详细信息。基本语法如下:man[选项]命令......