机器学习之叁：PCA 与 LDA

Posted on 繁/简： set

准备写这篇笔记的时候正好刷到了 MeAqua 的突击联动录播。记忆的锚点一瞬展开；我恍惚看向窗外，对面是五年前那个狭小的出租屋，夏日午后的阳光斜斜照在书桌上。摆头的风扇，被手机压住的试卷的边角不时被吹起。耳边传来的还是那么令人安心的吵闹嬉笑声，但那个抱着屏幕傻笑的少年已经不见踪影。


  This article is a self-administered course note.

  It will NOT cover any exam or assignment related content.


Feature Extraction

区分两个概念：Feature Extraction 与 Feature Selection。

两者都属于 Dimensionality Reduction 的范畴，但前者将自动把数据投影到一个全新的特征空间，后者通常需要特征工程的参与 (人为的设计特征的选取方案)，保留大量原有的特征。

接下来介绍的 PCA 与 LDA，都属于 Feature Extraction 的算法。


Principal Component Analysis

(研究了一天。线代不好地动山摇) Pre-requisite:

  • 协方差矩阵 (Covariance Matrix)
  • 简单的矩阵求导 (Derivative of a Matrix)
    • \(\frac{d}{dx} Ax=A\)
    • \(\frac{d}{dx} x^T Ax=2Ax\)
  • 拉格朗日乘子法 (Lagrange Multiplier)

PCA 的目的是在 \(d\) 维特征空间中找到 \(k\) (\(k<<d\)) 个单位正交基 \(\{u_1,u_2,...,u_k\}\)，使得投影到这 \(k\) 维空间的数据点的方差最大。这能够实现由 \(d\) 到 \(k\)​ 降维的同时，尽可能的保留较多的信息 (方差越大，熵越大)。

Shape 参考表：

  • 列向量 \(x_n:d\times 1\)。​
  • 列向量 \(u_j:d\times 1\)。​
  • 协方差矩阵 \(S:d\times d\)，其中 \(S_{i,j}\) 是第 \(i\) 维与第 \(j\)​ 维间的协方差。
  • 主成分矩阵/转换矩阵/降维矩阵/投影矩阵 \(\mathbf{U}=\{u_1,u_2,...,u_k\}:d\times k\)。

由于 \(\{u_1,u_2,...,u_k\}\) 是单位正交基，第 \(n\) 个数据 \(x_n\) 在第 \(j\) 个基上的投影是 \(u_j^{T}x_n\)。那么数据集在第 \(j\) 维的方差是： \[ \begin{aligned} \frac{1}{N} \sum_{n=1}^{N} (u_j^{T}x_n-u_j^{T}\overline{x})^2 &=\frac{1}{N}\sum_{n=1}^N[u_j^T (x_n-\overline{x})]^2\\ &=\frac{1}{N}\sum_{n=1}^N u_j^T (x_n-\overline{x})(x_n-\overline{x})^T u_j\\ &=u_j^{T}[\frac{1}{N}\sum_{n=1}^N(x_n-\overline{x})(x_n-\overline{x})^T]u_j \\ &= u_j^T Su_j \end{aligned} \] 在代码实现中，常常先对 \(\mathbf{X}\) 应用 standardization 得到 \(\mathbf{X}'\)，这里就可以直接使用 \(\frac{1}{N}\sum x_n' x_n'^{T}\) 构建协方差矩阵。标准化的另一个好处是消除 PCA 对均值的敏感性。

建立最优化模型：求一系列 \(\mathbf{U}=\{u_1,u_2,...,u_k\}\)​，使得： \[ \max f(\{u_1,u_2,...,u_k\})=\sum_{j=1}^k u_j^T Su_j \ \ s.t. \ \ u_j^Tu_j=1 \] 对每一维应用拉格朗日乘子，在第 \(j\) 维有： \[ \begin{aligned} \frac{df}{du_j}+\lambda_j(1-u_j^Tu_j)&=0 \\ 2Su_j-2\lambda_j u_j&=0 \\ Su_j&=\lambda_ju_j \end{aligned} \] 等等，这不是我们熟悉的 Eigen Decomposition 吗？这么看来，\(u_j\) 是协方差矩阵 \(S\) 的一个特征向量，\(\lambda_j\) 是它对应的特征值。把所有 \(k\) 维都考虑进来，写成矩阵形式： \[ \begin{bmatrix} | & | &{...}&| \\ Su_1 & Su_2 &{...}& Su_k \\ |&|&{...}&|\end{bmatrix}= \begin{bmatrix} | & | &{...}&| \\ \lambda_1 u_1 & \lambda_2u_2&{...}& \lambda_k u_k \\ | & | &{...}& |\end{bmatrix} \] 左右两个矩阵的 shape 都是 \(d\times k\)。同时左乘一个 \(\mathbf{U}^T\)，再求迹，有： \[ \begin{aligned} \text{tr}\begin{pmatrix} \begin{bmatrix}-&u_1&-\\ -&u_2&-\\ -&:&-\\ -&u_k&-\end{bmatrix}\begin{bmatrix} | & | &{...}&| \\ Su_1 & Su_2 &{...}& Su_k \\ |&|&{...}&|\end{bmatrix} \end{pmatrix} &= \text{tr}\begin{pmatrix} \begin{bmatrix}-&u_1&-\\ -&u_2&-\\ -&:&-\\ -&u_k&-\end{bmatrix} \begin{bmatrix} | & | &{...}&| \\ \lambda_1 u_1 & \lambda_2u_2&{...}& \lambda_k u_k \\ | & | &{...}& |\end{bmatrix}\end{pmatrix} \\ \sum_{j=1}^k u_j^TSu_j&=\sum_{j=1}^k \lambda_j u_j^Tu_j \ \ s.t. \ \ u_j^Tu_j=1 \\ \sum_{j=1}^k u_j^TSu_j&=\lambda_1+\lambda_2+{...}+\lambda_k \end{aligned} \] 对比一下上面我们建立的最优化模型：等式左边不就是我们要求最大化的项吗？

转化过后，这个问题就很好解决了：特征分解协方差矩阵 \(S\) 后，再做单位正交化；从大到小排序所有特征值，取前 \(k\) 个，就获得了 \(\{\lambda_1,\lambda_2,...,\lambda_k\}\)。这 \(k\) 个特征值对应的单位正交特征向量，就是要求的 \(\mathbf{U}=\{u_1,u_2,...,u_k\}\)​。

用计算好的主成分矩阵乘上原 \(d\) 维数据集 \(\mathbf{U}^T \mathbf{X}\)，就得到了降到 \(k\)​​ 维且方差最大的新数据集。

代码实现：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
import pandas as pd
from sklearn.preprocessing import StandardScaler

X, y = df_wine.iloc[:, 1:].values, df_wine.iloc[:, 0].values
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, stratify=y)

# step 1: standardization
scaler = StandardScaler()
X_train_std = scaler.fit_transform(X_train)
X_test_std = scaler.fit_transform(X_test)

# step 2: construct the covariance matrix
cov_mat = np.cov(X_train_std.T)

# step 3: eigen decomposition
eigen_Vals, eigen_vecs = np.linalg.eig(cov_mat)

# step 4: sort eigenvalues
eigen_pairs = [(np.abs(eigen_vals[i]), eigen_vecs[:, i]) for i in range(len(eigen_vals))]
eigen_pairs.sort(key=lambda k: k[0], reverse=True)

# step 5: construct principal component matrix, k=2
w = np.hstack((eigen_pairs[0][1][:, np.newaxis],
               eigen_pairs[1][1][:, np.newaxis]))

# step 6: apply projection
X_train_pca = X_train_std.dot(w)
X_test_pca = X_test_std.dot(w)

另外一种就是直接调包啦。

1
2
3
4
5
from sklearn.decomposition import PCA
pca = PCA(n_components=2)
X_train_pca = pca.fit_transform(X_train_std)
X_test_pca = pca.fit_transform(X_test_pca)
# pca.explained_variance_ratio_: check explained variance for each principal component


Linear Discriminant Analysis

PCA 是 unsupervised 的降维算法，LDA 则是 supervised 的，需要得到数据的 labels。

不同于 PCA 最大化方差，LDA 希望原始数据在投影到低维后能够使得同类的数据尽可能聚集，不同类的数据尽可能分散，即最大化 Class Separability。

LDA 假设数据 normally distributed，且每一类数据都有相同的协方差矩阵；由此可见，如果不同类的数据之间方差差异较大，PCA 效果好；均值差异较大，LDA 效果好。

LDA 的数学推导和 PCA 很相似，还是构建最优化模型 - 拉格朗日乘子 - 特征分解这一系列过程。

LDA 的重心在于均值，因此它使用散度矩阵 (scatter matrix) 而非协变量矩阵来描述维度间的关系。

Shape 参考表：

  • 列向量 \(x:1\times d\)​。
  • 列向量 - 全局均值向量 \(\overline{m}:1\times d\)。
  • 列向量 - 第 \(i\) 类数据的均值向量 \(\overline{m_i}:1\times d\)。
  • 第 \(i\) 类数据的散度矩阵 \(\mathbf{S}_i:d\times d\)​。
  • 类内散度矩阵 (within-class SM) \(\mathbf{S}_w\) 与类间散度矩阵 (between-class SM) \(\mathbf{S}_b:d\times d\)。

类内与类间散度矩阵的求法。 \[ \begin{aligned} \overline{m_i}&=\frac{1}{n_i}\sum_{x\in D_i}x, \ \overline{m}=\frac{1}{N}\sum x \\ \mathbf{S}_i&=\sum_{x\in D_i}(x-\overline{m_i})(x-\overline{m_i})^T \\ \mathbf{S}_W&=\sum_{i=1}^c \mathbf{S}_i \\ \mathbf{S}_B&=\sum_{i=1}^c (\overline{m_i}-\overline{m})(\overline{m_i}-\overline{m})^T \end{aligned} \] 定义 Class Separability 为 \(\mathbf{S}_W^{-1}\mathbf{S}_B\)：最大化类间散布 (分子)，最小化类内散布 (分母)。构建最优化模型：求一系列 \(\mathbf{U}=\{u_1,u_2,...,u_k\}\)，使得： \[ \max f(\{u_1,u_2,...,u_k\})=\sum_{j=1}^k u_j^T \mathbf{S}_W^{-1}\mathbf{S}_Bu_j \ \ s.t. \ \ u_j^Tu_j=1 \] 接下来的推导就和 PCA 一样了，把协方差矩阵 \(S\) 替换为 Class Separability 矩阵 \(\mathbf{S}_W^{-1}\mathbf{S}_B\) 即可：特征分解后，取前 \(k\)​ 大特征值对应的单位正交特征向量构建转换矩阵。

1
2
3
4
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA
lda = LDA(n_components=2)
X_train_lda = lda.fit_transform(X_train_std, y_train)
X_test_lda = lda.transform(X_test_std)


Reference

  This article is a self-administered course note.

  References in the article are from corresponding course materials if not specified.

Course info: COMP3314. Lecturer: Prof. Yizhou Yu.

Course textbook: Mathematics for Machine Learning. M. Deisenroth, A. Faisal & C.Ong.

参考：

-----------------------------------そして、次の曲が始まるのです。-----------------------------------