主成分分析（PCA）的数学原理与应用

PCA（主成分分析）主要用于数据降维。

问题：样本 $x_i \in \mathbb R^{m}$ 以列向量表示，现将 $n$ 个样本构成的样本集 $X\in \mathbb R^{m\times n}$ 降维到 $\hat{X} \in \mathbb R^{k\times n}$（$k \leq m$），以达到数据压缩的目的。

不难发现，上述问题的核心是：将 $X$ 中的每个样本点由原来的 $m$ 维坐标系投影到新的 $k$ 维坐标系中。

考虑最简单的二维数据：

import numpy as np
import matplotlib.pyplot as plt

# 绘制数据点
k = 0.3   # 斜率
b = 0.1   # 截距
x1 = np.random.random(30)
x2 = k * x1 + b + np.random.randint(-5, 5, 30) / 100

ax = plt.gca()
ax.set_aspect(1)
plt.scatter(x1, x2, marker='x')

# 中心点
cx = np.mean(x1)
cy = np.mean(x2)

# 绘制第一个维度u1
u1_x = 0.9
u1_y = k * u1_x + b
plt.annotate('', xytext=(cx, cy), xy=(u1_x, u1_y), arrowprops=dict(arrowstyle="->"))
plt.text(u1_x - 0.1, u1_y + 0.01, '$u_1$', fontsize=16)

# 绘制第二个维度u2
u2_x = cx - 0.02
u2_y = (-1 / k) * u2_x + (cy + cx / k)
plt.annotate('', xytext=(cx, cy), xy=(u2_x, u2_y), arrowprops=dict(arrowstyle="->"))
plt.text(u2_x - 0.1, u2_y - 0.01, '$u_2$', fontsize=16)

plt.show()

人能够很直观地看出上述二维数据的2个维度的大致方向，其中最主要的那个方向是 $u_1$ ，因为数据在该方向上的分布最为分散（投影绝对值或方差之和最大）。也可以换个角度看，如果要用1维数据来最大限度地“拟合”二维数据，显然这些二维数据投影在 $u_1$ 这个方向上要比投影在 $u_2$ 上好。因此，$u_1$ 是该数据集的【第1主要成分】，而 $u_2$ 则是【第2主要成分】，这个思想同样适用于更高维的数据。

通过上面的分析可知，当所有的数据点在 $u_1$ 方向上的投影绝对值（或方差）之和最大时，此时的 $u_1$ 就是第1主成分。找 $u_1$ 的问题则可以表示为：

$$\large{\arg\max_{u_1}\left(\sum_{i=1}^n \lvert proj(x_i, u_1) \rvert\right)}\qquad\qquad x_i, u_1\in\mathbb R^{m}$$

其中，$proj(x_i, u_1)$ 为 $x_i$ 在 $u_1$ 上的投影长度。容易求得：

$$\large{proj(x_i, u_1) = \frac{x_i \cdot u_1}{ \lVert{u_1} \rVert_2} = \frac{x_i^T u_1}{ \lVert{u_1} \rVert_2}}$$

因为我们只考虑 $u_1$ 的方向，故设 $u_1$ 具有单位长度，即 $u_1$ 的 $L_2$ 范数为1，则：

$$\begin{align} & \large{\arg\max_{u_1}\left(\sum_{i=1}^n \lvert proj(x_i, u_1)\rvert\right)} \large{=\:} \large{\arg\max_{u_1}\left(\sum_{i=1}^n \lvert x_i^T u_1 \rvert \right)} \qquad s.t. {\lVert u_1 \rVert}_2 = 1\\ \large{=\:}& \large{\arg\max_{u_1}\left(\sum_{i=1}^n \left( {x_i^T u_1} \right) ^2 \right)} \large{=\:} \large{\arg\max_{u_1}\left(\sum_{i=1}^n \left( {x_i^T u_1} \right)^T {x_i^T u_1} \right)} \\ \large{=\:}& \large{\arg\max_{u_1}\left(\sum_{i=1}^n {u_1^T x_i x_i^T u_1} \right)} \large{=\:} \large{\arg\max_{u_1}\left(u_1^T \left( \sum_{i=1}^n {x_i x_i^T} \right) u_1 \right)} \end{align}$$

因为，$\large{X=\left[ \matrix{x_1\quad x_2 \cdots x_n} \right], \qquad X^T=\left[ \matrix{x_1^T \cr x_2^T \cr \vdots \cr x_n^T} \right]}$，所以：

$$\begin{equation} \large{\arg\max_{u_1}\left(u_1^T \left( \sum_{i=1}^m {x_i x_i^T} \right) u_1 \right) = \arg\max_{u_1}\left(u_1^T XX^T u_1 \right)} \end{equation}$$

上式即为我们要最大化的目标函数。若能证明矩阵 $XX^T \in \mathbb R^{m \times m}$ 是（半）正定的，则说明 $u_1^TXX^Tu_1$ 是（半）正定二次型，也就存在最大值。设 $\lambda$ 是对称矩阵 $XX^T$ 的任一特征值，其对应的特征向量是 $v \in \mathbb R^{m}$，则有：

$$\begin{align} & \large{XX^Tv = \lambda v} \\ \large{\Rightarrow} \quad & \large{(XX^Tv)^T v = (\lambda v)^T v} \\ \large{\Rightarrow} \quad & \large{v^TXX^T v = v^T \lambda v} \\ \large{\Rightarrow} \quad & \large{(X^Tv)^TX^Tv = \lambda v^Tv} \\ \large{\Rightarrow} \quad & \large{{\lVert X^Tv \rVert}_2^2 = \lambda {\lVert v \rVert}_2^2} \\ \large{\Rightarrow} \quad & \large{\lambda \geq 0} \end{align}$$

证毕。

下面使用拉格朗日乘子法求 $u_1^TXX^Tu_1$ 在约束 $u_1^Tu_1 = 1$ 下的极值，令：

$$\large{f(u_1) = u_1^TXX^Tu_1 - \alpha(1 - u_1^Tu_1)}$$

对 $u_1$ 求导，并令其为 $0$：

$$\begin{align} & \large{\frac{\partial f}{\partial u_1} = 2XX^Tu_1 - 2\alpha u_1 = 0} \\ \large{\Rightarrow} \quad & \large{XX^Tu_1 = \alpha u_1} \end{align}$$

可见，当 $u_1$ 取矩阵 $XX^T$ 的某个特征值 $\alpha$ 对应的特征向量时，前述目标函数取得最大值。注意，若选取的特征值 $\alpha$ 不同，其对应的特征向量 $u_1$ 也不同，则前述目标函数的最大值也不同。将上式代入目标函数，可得：

$$\large{\arg\max_{u_1}\left(u_1^T XX^T u_1 \right) = \alpha u_1^Tu_1 = \alpha}$$

可见，当 $u_1$ 取最大特征值对应的特征向量时，前述目标函数将取得所有最大值中的最大值，此时的 $u_1$ 即为【第1主成分】。若取前 $k$ 大特征值对应的特征向量 $u_1, u_2, \cdots, u_k$，则样本 $x_i$ 在这 $k$ 个主成分（维度）上的投影则为：

$$\large{\hat{x_i} = \left[ \matrix{x_i^Tu_1\quad x_i^Tu_2 \cdots x_i^Tu_k} \right]} \qquad \qquad k \leq n$$

主成分矩阵：

$$\large{U = \left[ \matrix{u_1\quad u_2 \cdots u_k} \right] \in \mathbb R^{m \times k}} \qquad \qquad k \leq n$$

降维后的数据集：

$$\large{\hat{X} = U^TX \in \mathbb R^{k \times n}}$$

将降维后的 $k$ 维数据集还原到原始的 $m$ 维坐标系：

$$\large{X_{restore} = U\hat{X} = UU^TX \in \mathbb R^{m \times n}}$$

下面通过一个图像压缩的例子来验证PCA。

import math
import numpy as np
from PIL import Image
import matplotlib.pyplot as plt

# 读取图片
X = np.array(Image.open(".\_data\son.jpg").convert('L'), 'f') # M*N
M = X.shape[0]
N = X.shape[1]

XXT = X.dot(X.T) # shape: M*M

e, V = np.linalg.eig(XXT) #特征值、特征向量 
U = V[np.argsort(-e)] # M*M，按特征值从大到小，对特征向量排序

plt.subplots(figsize=(20,10))
r, c, i = (5, 6, 1)

for k in range(M, 0, -math.ceil(M/(r*c))):
    Uk = U[:, 0:k] # M*k
    X_ = Uk.dot(Uk.T).dot(X) # X_restore: M*N
    ratio = 1.0*k*(M+N)/(M*N)
    plt.subplot(r, c, i, title='↓k=%d, ratio=%.1f' % (k, ratio)).axis('off')
    plt.imshow(Image.fromarray(X_))
    i+=1

plt.axis('off')
plt.show()