第五章 降维

1 常见降维算法

• 线性降维方法：各个属性（特征）独立无关
  • 奇异值分解SVD
  • 主成分分析PCA
  • 因子分析FA
  • 线性判别分析LDA
• 非线性降维方法：各个属性（特征）存在较强的相关性
  • 多维尺度变换MDS
  • 等距特征映射IsoMap
  • t分布随机邻域嵌入t-SNE
  • 统一流行逼近与投影UMAP

---

2 奇异值分解SVD

from numpy import linalg.svd

一种将高维矩阵表示为三个低维矩阵乘积的矩阵分解方法

算法描述

$${\mathbf{A}}_{(N\times p)}={\mathbf{U}}_{(N\times N)}{\mathbf{\Sigma }}_{(N\times p)}{\mathbf{V}}_{(p\times p)}^T$$

其中${\mathbf{U}}_{(N\times N)}$为左奇异矩阵，${\mathbf{V}}_{(p\times p)}$为右奇异矩阵，且$\mathbf{U}、\mathbf{V}$均为正交矩阵；${\mathbf{\Sigma }}_{(N\times p)}$为奇异值矩阵，其主对角线上的值称为奇异值，其他位置元素都为0

算法步骤

1. 求解矩阵${\mathbf{A}}^T\mathbf{A}$的特征值${\lambda }_1,{\lambda }_2,\cdot \cdot \cdot ,{\lambda }_p({\lambda }_1\ge {\lambda }_2\ge \cdot \cdot \cdot \ge {\lambda }_p)$，以及对应的正交化特征向量矩阵$({\mathbf{v}}_1,{\mathbf{v}}_2,...,{\mathbf{v}}_{\mathbf{p}})$，得到右奇异矩阵$\mathbf{V}=({\mathbf{v}}_1,{\mathbf{v}}_2,...,{\mathbf{v}}_{\mathbf{p}})$
2. 求解矩阵${\mathbf{AA}}^T$的特征值${\omega }_1,{\omega }_2,...,{\omega }_N({\omega }_1\ge {\omega }_2\ge \cdot \cdot \cdot \ge {\omega }_N)$，以及对应的正交化特征向量矩阵$({\mathbf{u}}_1,{\mathbf{u}}_2,...,{\mathbf{u}}_{\mathbf{N}})$，得到右奇异矩阵$\mathbf{U}=({\mathbf{u}}_1,{\mathbf{u}}_2,...,{\mathbf{u}}_{\mathbf{N}})$
3. 求出奇异值矩阵${\mathbf{\Sigma }}_{(N\times p)}={\mathbf{\Sigma }}_{(N\times p)}=\left(\begin{array}{cc}{\mathbf{\Sigma }}_1& 0\\ 0& 0\end{array}\right)$，${\mathbf{\Sigma }}_1=diag(\sqrt{{\lambda }_1},\sqrt{{\lambda }_2},\cdot \cdot \cdot ,\sqrt{{\lambda }_p})$，${\lambda }_i$为矩阵${\mathbf{A}}^{\mathbf{T}}\mathbf{A}$的特征值
4. 得到奇异值分解结果$\mathbf{A}=\mathbf{U}\mathbf{\Sigma }{\mathbf{V}}^T$

截断奇异值分解Truncated SVD

SVD算法在分解过程中，仅保留最大的d个奇异值对应的部分，其中d远远小于p（原始数据维数）
$${\mathbf{A}}_{(N\times p)}={\mathbf{U}}_{(N\times N)}{\mathbf{\Sigma }}_{(N\times p)}{\mathbf{V}}_{(p\times p)}^T\approx {\mathbf{U}}_{(N\times d)}{\mathbf{\Sigma }}_{(d\times d)}{\mathbf{V}}_{(d\times p)}^T$$

算法不足

使用SVD算法对数据矩阵分解后得到的结果（奇异值和奇异矩阵），无法直接对应于原始数据的特征，导致降维后的数据可读性与可解释性较低

---

3 主成分分析PCA

from sklearn.decomposition import PCA

采用线性变换将一组统计变量（指标）变换成为新的一组综合变量，转换后的这组变量称为主成分，且主成分之间互不相关，其着重解释各变量的总方差

为消除由于量纲不同可能带来的影响，一般在主成分分析中会首先对原始数据进行标准化处理

标准化处理：from sklearn.preprocessing import StandardScaler

算法描述

$$\mathbf{F}={\mathbf{A}}^T\mathbf{X}$$

其中$\mathbf{F}=(F_1,F_2,\cdot \cdot \cdot ,F_p{)}^T,\mathbf{A}=({\mathbf{A}}_1,{\mathbf{A}}_2,\cdot \cdot \cdot ,{\mathbf{A}}_{\mathbf{P}})=\left(\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1p}\\ a_{21}& a_{22}& \cdots & a_{2p}\\ ⋮& ⋮& \ddots & ⋮\\ a_{p1}& a_{p2}& \cdots & a_{pp}\end{array}\right)$

算法步骤

1. 求解标准化后的数据矩阵$\mathbf{X}$的协方差矩阵$\mathbf{\Sigma }$的特征值，并进行降序排列${\lambda }_1,{\lambda }_2,\cdot \cdot \cdot ,{\lambda }_p({\lambda }_1\ge {\lambda }_2\ge \cdot \cdot \cdot \ge {\lambda }_p)$
2. 求解各特征值对应的单位化特征向量矩阵$\mathbf{A}=({\mathbf{A}}_1,{\mathbf{A}}_2,\cdot \cdot \cdot ,{\mathbf{A}}_{\mathbf{p}})$
3. 确定主成分个数$m$。前$m$个主成分的累计方差贡献率为${\sum }_{i=1}^m{\lambda }_i/{\sum }_{j=1}^p{\lambda }_j$，一般可要求累计方差贡献率大于$80\%$
4. 得到降维后的主成分$\mathbf{F}={\mathbf{A}}^T\mathbf{X}$

算法性质

• 每个主成分的系数平方和为1
• 主成分之间相互独立
• 主成分的方差依次递减

算法评价

PCA算法适用于变量之间有较强的相关性的数据，如果原始数据相关性较弱时，方差小的成分可能含有影响样本差异的重要信息，降维丢弃可能对后续数据处理有影响

---

4 因子分析FA

from factor_analyzer import FactorAnalyzer

通过使用少数公共因子与特殊因子来描述原有的变量，并着重分析各变量之间的协方差，以此实现数据的降维

主成分分析 VS 因子分析

主成分分析：想把现有的变量变为少数几个新的变量（新的变量几乎带有原来所有变量的信息）

因子分析：需要寻找潜在的因子，并对这些因子进行解释

与主成分分析相比，因子分析中可以对因子进行旋转，在解释和可读性方面更加有优势

算法描述

$$\mathbf{X}-\mathbf{\mu }=\mathbf{AF}+\mathbf{\epsilon }$$

其中$\mathbf{F}=(f_1,f_2,\cdot \cdot \cdot ,f_m{)}^T,\mathbf{\epsilon }=({\epsilon }_1,{\epsilon }_2,\cdot \cdot \cdot ,{\epsilon }_p{)}^T$，因子载荷矩阵$\mathbf{A}=({\mathbf{A}}_1,{\mathbf{A}}_2,\cdot \cdot \cdot ,{\mathbf{A}}_{\mathbf{P}})=\left(\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1m}\\ a_{21}& a_{22}& \cdots & a_{2m}\\ ⋮& ⋮& \ddots & ⋮\\ a_{p1}& a_{p2}& \cdots & a_{pm}\end{array}\right)$

记载荷矩阵$\mathbf{A}$中第$i$行元素的平方和$h_i^2$为第$i$共同度：$h_i^2$越大，则该变量能被公共因子解释的程度越高

记载荷矩阵$\mathbf{A}$中第$j$列元素的平方和$g_j^2$：$g_j^2$越大，则该公共因子对于解释原始变量的作用就越大

算法步骤

1. 进行适用性检验。对数据矩阵$\mathbf{D}$进行标准化处理后，求得协方差矩阵$\mathbf{\Sigma }$，利用KMO检验和Bartlett球性检验方法判断是否可以对变量进行因子分析

   适用性检验：

   检验特征之间是否具有相关性，如果相关性太低则不适合使用因子分析

   • KMO检验：$KMO-value>0.6$

     from factor_analyzer.factor_analyzer import calculate_kmo

   • Bartlett球性检验：$p-value<0.05$

     from factor_analyzer.factor_analyzer import calculate_bartlett_sphericity

2. 确定公共因子个数。根据原有变量求解协方差矩阵并求解特征值及特征向量，可根据因子特征值大于1或累计特征值所占百分比大于$80\%$，确定公共因子个数$m$，求得对应的特征向量

3. 求解初始因子载荷矩阵。根据公式$\mathbf{A}=\mathbf{U}{\mathbf{\Lambda }}^{1/2}=(\sqrt{{\lambda }_1}{\mathbf{u}}_1,\sqrt{{\lambda }_2}{\mathbf{u}}_2,\cdots ,\sqrt{{\lambda }_m}{\mathbf{u}}_{\mathbf{m}})$，求得初始因子载荷矩阵

4. 进行因子旋转${\mathbf{A}}^{\ast }=\mathbf{AQ}$。采用因子旋转方法将因子载荷矩阵正交旋转，使得因子变量更具有可解读性

   因子旋转：

   基于公共因子模型的选择不定性，可以对载荷矩阵进行旋转。因子旋转是指用一个正交矩阵$\mathbf{Q}$对载荷矩阵$\mathbf{A}$旋转，记旋转因子载荷矩阵$A^{\ast }=AQ$，此时正交因子模型可写为：
   $$\mathbf{X}-\mathbf{\mu }={\mathbf{A}}^{\ast }({\mathbf{Q}}^T\mathbf{F})+\mathbf{\epsilon }$$

   因子旋转方法：

   • 四次方最大法：因子负荷平方和最大化
   • 方差最大法：方差最大化
   • 等量最大法：前两者权衡

5. 计算因子得分。采用最小二乘法计算因子得分，并对因子进行命名

算法性质

• 公共因子和特殊因子之间不相关
• 各公共因子之间不相关，且方差均为1
• 各特殊因子之间不相关且方差不相等

算法评价

通过识别数据中的共性和模式，可以将大量的变量归纳为较少数量的因子

---

5 线性判别分析LDA

from sklearn.discriminant_analysis import LinearDiscriminantAnalysis

旨在找到原始数据特征变量的一个线性组合，然后将数据投影到一个新的空间。这使得不同类别数据投影点的中心尽可能的远，同类别数据的投影点的中心尽可能的近，即类间方差最大化，类内方差最小化

算法步骤

1. 构造最大化目标函数$J(\mathbf{w})=\prod _{i=1}^d\frac{{\mathbf{w}}_{\mathbf{i}}^{\mathbf{T}}{\mathbf{S}}_{\mathbf{b}}{\mathbf{w}}_{\mathbf{i}}}{{\mathbf{w}}_{\mathbf{i}}^{\mathbf{T}}{\mathbf{S}}_{\mathbf{w}}{\mathbf{w}}_{\mathbf{i}}}$

   • 类间散度矩阵${\mathbf{S}}_{\mathbf{b}}={\sum }_{j=1}^k{N}_j({\mathbf{\mu }}_{\mathbf{j}}-\mathbf{\mu })({\mathbf{\mu }}_{\mathbf{j}}-\mathbf{\mu }{)}^T$

   • 类内散度矩阵${\mathbf{S}}_{\mathbf{w}}={\sum }_{j=1}^k{\sum }_{x\in {\mathbf{X}}_{\mathbf{j}}}(\mathbf{x}-{\mathbf{\mu }}_{\mathbf{j}})(\mathbf{x}-{\mathbf{\mu }}_{\mathbf{j}}{)}^T$

2. 计算类内散度矩阵${\mathbf{S}}_{\mathbf{w}}$与类间散度矩阵${\mathbf{S}}_{\mathbf{b}}$，求出矩阵${{\mathbf{S}}_{\mathbf{w}}}^{-1}{\mathbf{S}}_{\mathbf{b}}$

   矩阵不可逆时：

   • 对矩阵${\mathbf{S}}_{\mathbf{w}}$进行SVD分解。得${\mathbf{S}}_{\mathbf{w}}=\mathbf{U}\mathbf{\Sigma }{\mathbf{V}}^T$，进而${{\mathbf{S}}_{\mathbf{w}}}^{-1}=\mathbf{V}{\mathbf{\Sigma }}^{-1}{\mathbf{U}}^T$
   • 先用PCA对数据进行降维，使得降维后的矩阵${{\mathbf{S}}_{\mathbf{w}}}^{\ast }$可逆。$({{\mathbf{S}}_{\mathbf{w}}}^{\ast }{)}^{-1}={\sum }_{j=1}^k{N}_j\mathbf{w}({\mathbf{\mu }}_{\mathbf{j}}-\mathbf{\mu }{)}^{-1}$

3. 通过对矩阵${{\mathbf{S}}_{\mathbf{w}}}^{-1}{\mathbf{S}}_{\mathbf{b}}$进行特征值求解得到投影方向$\mathbf{w}$

4. 计算投影后的数据点${\mathbf{z}}_{\mathbf{i}}={\mathbf{w}}^T{\mathbf{x}}_{\mathbf{i}}$，将原始数据投影到直线$y$上

算法评价

LDA算法寻找满足最大化类间距离和最小化类内距离的最佳投影方向，在低维空间中能够区分不同类别的数据

但LDA算法的计算复杂度较高，而且当样本数量较少或者类别之间的方差差异很大时，容易出现过拟合现象

LDA算法能够找到一个投影方向，使得不同类别的数据点在投影后能够有很好的可分性，利用类别信息来最大程度保留数据的判别性能、

---

6 多维尺度变换MDS

from sklearn.manifold import MDS

在尽可能保持数据点之间的相似性不变的前提下，将高维坐标的数据点投影到低维的坐标系中的一种降维方法

算法步骤

1. 计算数据集的距离矩阵Dd(N×N)={dij}=∥xi−xj∥\boldsymbol{D}_{d(N\times N)}=\{d_{ij}\}=\parallel{x_i-x_j}\parallelDd(N×N)​={dij​}=∥xi​−xj​∥，求出$d_{ij}^2,d_{\cdot j}^2,d_{i\cdot }^2,d_{\cdot \cdot }^2$
2. 根据公式$b_{ij}=-\frac{1}{2}(d_{ij}^2-d_{\cdot j}^2-d_{i\cdot }^2+d_{\cdot \cdot }^2)$，得到内积矩阵$\mathbf{B}$
3. 对矩阵$\mathbf{B}$进行特征值分解，得$\mathbf{B}=\mathbf{V}\mathbf{\Lambda }{\mathbf{V}}^{\mathbf{T}}$，选择前$k$个最大的特征值${\mathbf{\Lambda }}_{\mathbf{k}}$及对应的特征向量${\mathbf{V}}_{\mathbf{k}}$
4. 计算降维后的矩阵${\mathbf{Z}}_{(N\times k)}={\mathbf{V}}_{\mathbf{k}}{\mathbf{\Lambda }}_{\mathbf{k}}^{1/2}$

算法评价

MDS算法降维后的数据布局反映原始数据的相似性和差异性，但使用的前提是数据点之间的距离可以用线性函数表示，这里采用欧式距离作为数据点之间的距离度量

---

7 等距特征映射IsoMap

from sklearn.manifold import Isomap

一种基于全局特征保持的流形学习算法

MDS VS IsoMap

• MDS算法：两点之间的距离为欧式距离（直线距离）
• IsoMap算法：两点之间的距离为测底线距离

算法步骤

1. 根据输入的高维数据集，计算样本之间的距离矩阵

2. 根据指定近邻点的个数或者设定的距离与之构建近邻图

   近邻图：

   • $k$近邻图：指定近邻点的个数$k$
   • $\epsilon$近邻图：制定距离阈值$\epsilon$

3. 通过最短路径算法计算邻接图中数据点之间的最短距离，得到数据点的实际侧低距离矩阵

   常用最短路径算法：

   • Dijkstra算法
   • Floyd算法
   • Bellman-Ford算法
   • SPFA算法

4. 通过MDS算法对侧低距离矩阵进行降维

算法评价

IsoMap算法能够在降维的同时尽量保持原始数据之间的相对位置关系，在保留原始数据的总体结构和流形结构特性的同时进行降维，在较低维数下对数据进行可视化展示。

IsoMap算法计算复杂度和空间复杂度都较高，不适用于数量级较大的数据集

---

8 t分布随机邻域嵌入t-SNE

from sklearn.manifold import TSNE

关注数据的局部结构流形学习方法

用条件概率度量数据点之间的相似性，并最小化高维空间和低维空间中数据点相似度的差异来实现数据的非线性降维

算法描述

高维空间：

采用高斯分布构建条件分布函数$P_i$，$P_i$表示给定数据点$x_i$时，其他所有数据点的条件概率分布。以$x_i$为中心构建方差为${\sigma }_i$的高斯分布，$x_i$以条件概率$p_{j|i}$选择$x_j$作为它的邻近点，定义$p_{j|i}$为
$$p_{j|i}=\frac{\exp \left(-\parallel x_i-x_j{\parallel }^2/2{\sigma }_i^2\right)}{{\sum }_{k\ne i}\exp \left(-\parallel x_i-x_k{\parallel }^2/2{\sigma }_i^2\right)}$$
SNE算法只关注不同数据点之间的相似度，所以认为$p_{i|i}=0$

低维空间：

采用自由度为1的t分布构建条件概率分布函数$Q_i$，同理定义$q_{j|i}$为
$$q_{j|i}=\frac{(1+\parallel y_i-y_j{\parallel }^2{)}^{-1}}{{\sum }_{k\ne l}(1+\parallel y_k+y_l{\parallel }^2{)}^{-1}}$$
目标函数：

采用KL散度度量概率分布的相似性，运用梯度下降算法求解目标函数C
$$C=\sum _{i}KL(P_i\parallel Q_i)=\sum _i\sum _j{p}_{j|i}\log \frac{p_{j|i}}{q_{j|i}}$$

改进随机邻域嵌入SNE算法

SNE算法的局限性	t-SNE算法的改进
条件不对称问题：$p_{j	i}\ne p_{i
**数据”拥挤“问题：**SNE算法使用高斯函数构建条件概率分布，但由于高斯分布在远离均值处呈现指数式下降趋势，导致数据”拥挤“问题的出现	t-SNE算法在低维空间中使用t分布构建条件概率函数

算法步骤

1. 对于给定的高维数据集，构建条件概率$p_{j|i}=\frac{\exp \left(-\parallel x_i-x_j{\parallel }^2/2{\sigma }_i^2\right)}{{\sum }_{k\ne i}\exp \left(-\parallel x_i-x_k{\parallel }^2/2{\sigma }_i^2\right)}$
2. 在低维空间中，引入t分布函数构建条件概率$q_{j|i}=\frac{(1+\parallel y_i-y_j{\parallel }^2{)}^{-1}}{{\sum }_{k\ne l}(1+\parallel y_k+y_l{\parallel }^2{)}^{-1}}$
3. 计算联合概率$p_{ij}=\frac{p_{i|j}+p_{j|i}}{2n}$，令$p_{i|j},p_{j|i}:=p_{ij}$
4. 使用梯度下降算法求解$\frac{\partial C}{\partial y_i}=4{\sum }_j(p_{j|i}-q_{j|i})(y_i-y_j)(1+\parallel y_i-y_j{\parallel }^2{)}^{-1}$，获得数据点在低维空间中的分布$q_{j|i}$，以最小化目标函数$C$

算法评价

t-SNE算法在处理数据时对异常值较为敏感，且其计算复杂度较高，因此不适用于处理高维稀疏数据和大规模数据

由于其局部优化性质，t-SNE算法在降维过程中通常无法保留全局数据结构，这使其保留全局数据结构，这使其在需要完整保留数据集的整体分布和关系的情况下可能表现不佳

当高维空间中两个数据点距离较近，映射到低维空间后距离较远时，将得到一个很高的惩罚值；但如果高维空间中两个数据点距离较远，映射到低维空间距离较近时，惩罚值反而很低

---

9 统一流行逼近与投影UMAP

from umap import UMAP

基本思想是先构建高维空间的局部结构，然后通过梯度下降算法最优化目标函数，以保持低维空间中数据点的局部关系和全局结构

算法描述

高维空间：

计算点$x_i$与最近$k$个点之间的距离为$d(x_i,x_{ij})$，其中$j=1,2,\cdots ,k$，记$x_i$的最近邻值${\rho }_i$为
ρi=min⁡{d(xi,xij)∣1≤j≤k,d(xi,xij)>0} \rho_i=\min\{d(x_i,x_{ij})|1 \le j \le k,d(x_i,x_{ij})>0\} ρi​=min{d(xi​,xij​)∣1≤j≤k,d(xi​,xij​)>0}
将数据点$x_i$选择数据点$x_{ij}$作为其最近邻的概率记为$p_{j|i}$
$$p_{j|i}=e^{-\frac{d(x_i,x_{ij})-{\rho }_i}{{\sigma }_i}}$$
其中${\sigma }_i$为归一化因子，可由如下公式计算得出
$$\sum _{j=1}^k{e}^{-\frac{d(x_i,x_{ij})-{\rho }_i}{{\sigma }_i}}={\log }_{2}k$$
用联合概率代替条件概率
$$\begin{gathered} p_{ij}=(p_{j|i}+p_{i|j})-p_{j|i}{p}_{i|j} \\ {p}_{j|i},p_{i|j}:=p_{ij} \end{gathered}$$
低维空间：

引入$a,b$两个参数计算概率分布函数，通常取参数$a\approx 1.929,b\approx 0.7915$，$y_i$和$y_j$之间的条件概率分布表示为
$$q_{ij}={\left(1+a(\parallel y_i-y_j{\parallel }^2{)}^b\right)}^{-1}$$
目标函数：
$$C=\sum _{i\ne j}{p}_{ij}\log (\frac{p_{ij}}{q_{ij}})+(1-p_{ij})\log (\frac{1-p_{ij}}{1-q_{ij}})$$
通过梯度下降等优化方法最小化目标函数$C$

算法评价

适用于高维数据的降维，尤其是当数据具有非线性结构和复杂的局部关系时

对于邻域大小、最小距离和其他超参数的选择非常敏感，这些选择会影响最终的降维结果