# 主成分分析(Principal Component Analysis, PCA)

## 符号定义

|符号|含义|
|:-:|:-:|
|$X$|样本集合|
|$\pmb{x}$|数据点|
|$d$|数据点维度|
|$k$|降维后维度|
|$\hat{\pmb{x}}$|降维后数据点|
|$\pmb{y}$|降维后在k维子空间的坐标|

## 维度、内在维度与降维

对于一个多维的数据点$\pmb{x}\in \mathcal{R}^d$，若需要$d$个值来确切定义这一数据点，则其内在维度(inherent dimensionality)与其自然维度相同。

但是在实际数据中，往往不同维度的数据间有一定的关系，即对于一个多维的数据点 $ \pmb{x}\in \mathcal{R}^d $ ，往往仅需要 $ k<d $ 个值来定义这一数据点，此时内在维度要小于自然维度。

上述分析显然是存在问题的，因为实际数据中有噪声。噪声会导致一个多维的数据点 $ \pmb{x}\in \mathcal{R}^d $ 必须使用$d$个值来定义。但是噪声往往是无意义的，甚至在后续的分析中是需要被排除的，因此可以合理的认为噪声不影响内在维度的确定。

降维就是在寻找这么一个变换，将原本为$d$个维度的数据转换为$k$维。

降维显然有诸多好处：

* **降低数据对计算资源以及存储资源的需求**：这是显而易见的
* **凸显特征**：降维后能更好的体现原始数据在某一特定观点下的关键特征。无论是用于可视化（降至3维及以下），还是用于后续的分析均是有帮助的

## 降到零维子空间

若降维后为零维子空间，那么所有的样本均会降维到同一个点，设该点为$\pmb{m}$

对于样本集合$X=\{ \pmb{x_1}, \pmb{x_2}, \cdots, \pmb{x_n} \}$

现在即需要寻找点$\pmb{m}$，使得该点在一定的条件下是样本集合$X$降维到零维子空间后的最佳点

现在定义条件为：点$\pmb{m}$与样本集合中所有的数据点平均距离最小

有：

$$
\begin{equation}
    \pmb{m} = \argmin\limits_{m}\frac{1}{N}\sum_{i=1}^N||\pmb{x_i}-\pmb{m}||^2
\end{equation}
$$

寻找上述待优化问题的极小值点

$$
\begin{equation}
    \frac{\partial{J}}{\partial{\pmb{m}}} = \frac{2}{N}\sum_{i=1}^N(\pmb{m}-\pmb{x_i}) = 0
\end{equation}
$$

可以得到

$$
\begin{equation}
    \pmb{m} = \frac{1}{N}\sum_{i=1}^N\pmb{x_i}
\end{equation}
$$

即样本集合$X$在平均距离最小的条件下得到的最佳零维子空间降维结果为所有数据点的均值

## 降维到一维子空间

对于样本集合$X$中的任意数据点，其降维到一维子空间后的表示为：

$$
\begin{equation}
    \hat{\pmb{x_i}} = \pmb{x_{0}} + a_i\pmb{w}, i = 1, 2, \cdots, n
\end{equation}
$$

其中$\pmb{x_{0}}$和$\pmb{w}$由降维后的一维子空间决定，$a_i$由待降维的数据点决定，显然上述结果为一条直线。其中

注意到，最佳的零维子空间表示为所有数据点的均值，因此可以合理的将$\pmb{x_{0}}$设置为$\bar{\pmb{x}}$

即：
$$
\begin{equation}
    \hat{\pmb{x_i}} = \bar{\pmb{x}} + a_i\pmb{w}, i = 1, 2, \cdots, n
\end{equation}
$$

同样的，定义条件为：降维前后的平均距离最小


定义损失函数：
$$
\begin{equation}
    \begin{split}
    \mathcal{L} 
    &= \frac{1}{N}\sum_{i=1}^N||\pmb{x_i} - (\bar{\pmb{x}} + a_i\pmb{w})||^2\\
    &= \frac{1}{N}\sum_{i=1}^N||a_i\pmb{w}-(\pmb{x_i}-\bar{\pmb{x}})||^2\\
    &= \frac{1}{N}\sum_{i=1}^N(a_i^2||\pmb{w}||^2+||\pmb{x_i}-\bar{\pmb{x}}||^2 - 2a_i\pmb{w}^T(\pmb{x_i}-\bar{\pmb{x}}))
    \end{split}
\end{equation}
$$

上式分别对$a_i$和$\pmb{w}$求偏导有

$$
\begin{equation}
    \frac{\partial{J}}{\partial{a_i}} = \frac{2}{N}(a_i||\pmb{w}||^2-\pmb{w}^T(\pmb{x_i}-\bar{\pmb{x}})) = 0
\end{equation}
$$

$$
\begin{equation}
    \frac{\partial{J}}{\partial{\pmb{w}}} = \frac{2}{N}\sum_{i=1}^N(a_i^2\pmb{w}-a_i(\pmb{x_i}-\bar{\pmb{x}})) = 0
\end{equation}
$$

由损失函数对$a_i$的偏导可以得到

$$
\begin{equation}
    a_i = \frac{\pmb{w}^T(\pmb{x_i}-\bar{\pmb{x}})}{||\pmb{w}||^2}
\end{equation}
$$

若$\pmb{w}$为单位向量，上式显然是向量$\pmb{x_i}-\bar{\pmb{x}}$在$\pmb{w}$上投影的带符号长度

注意到

$$
\begin{equation}
    a_i\pmb{w} = \frac{\pmb{w}^T(\pmb{x_i}-\bar{\pmb{x}})}{||\pmb{w}||^2} \pmb{w} = \frac{(c\pmb{w})^T(\pmb{x_i}-\bar{\pmb{x}})}{||c\pmb{w}||^2} c\pmb{w},c \neq 0 \ and \ c\in \mathcal{R} 
\end{equation}
$$

即$\pmb{w}$乘以任意非零常数均不会改变降维后的结果

不妨设$\pmb{w} = \frac{\pmb{w}}{||\pmb{w}||_2}$，此时$\pmb{w}$为单位向量

因此有

$$
\begin{equation}
    a_i = \pmb{w}^T(\pmb{x_i}-\bar{\pmb{x}}) = (\pmb{x_i}-\bar{\pmb{x}})^T\pmb{w}
\end{equation}
$$

将式-11带入到损失函数有

$$
\begin{equation}
    \mathcal{L} = \frac{1}{N}\sum_{i=1}^N(||\pmb{x_i}-\bar{\pmb{x}}||^2 - a_i^2)
\end{equation}
$$

显然，当$\sum_{i=1}^Na_i^2$取得最大值时损失函数有最小值，即$\min\mathcal{L}$和$\max{\sum_{i=1}^Na_i^2}$等价

又由式-8有

$$
\begin{equation}
    \begin{split}
        \frac{\sum_{i=1}^Na_i^2}{N}\pmb{w} 
        &= \frac{1}{N}\sum_{i=1}^Na_i(\pmb{x_i}-\bar{\pmb{x}}) \\
        &= \frac{\sum_{i=1}^N(\pmb{x_i}-\bar{\pmb{x}})(\pmb{x_i}-\bar{\pmb{x}})^T\pmb{w}}{N} \\
        &= Cov(X)\pmb{w}
    \end{split}
\end{equation}
$$

其中$Cov(X)$为协方差矩阵

式-13成立表明$\pmb{w}$是协方差矩阵$Cov(X)$的特征向量，并且$\frac{\sum_{i=1}^Na_i^2}{N}$为对应的特征值

又结合式-12可知，需要最大化$\sum_{i=1}^Na_i^2$

因此，当降维到一维子空间时，$\pmb{w}$为协方差矩阵$Cov(X)$最大特征值对应的特征向量

考虑到协方差矩阵为实对称矩阵，设协方差矩阵$Cov(X)$的特征向量为$\pmb{\xi_1}, \pmb{\xi_2}, \cdots, \pmb{\xi_d}$，其对应的特征值依次增大

因此降维到一维子空间的表示为

$$
\begin{equation}
    \hat{\pmb{x_i}} = \bar{\pmb{x}} + (\pmb{\xi_1}^T(\pmb{x_i}-\bar{\pmb{x}}))\pmb{\xi_1}
\end{equation}
$$

## k维子空间

同上分析，可以推测k维子空间实际上就是协方差矩阵$Cov(X)$的前k个最大特征值对应的特征向量构成的子空间

对应的k维子空间表示为

$$
\begin{equation}
    \hat{\pmb{x_i}} = \bar{\pmb{x}} + (\pmb{\xi_1}^T(\pmb{x_i}-\bar{\pmb{x}}))\pmb{\xi_1} + (\pmb{\xi_2}^T(\pmb{x_i}-\bar{\pmb{x}}))\pmb{\xi_2} + \cdots + (\pmb{\xi_k}^T(\pmb{x_i}-\bar{\pmb{x}}))\pmb{\xi_k}
\end{equation}
$$

可以得到在k维子空间的坐标为

$$
\begin{equation}
    \pmb{y} = (\pmb{\xi_1}^T(\pmb{x_i}-\bar{\pmb{x}}), \pmb{\xi_2}^T(\pmb{x_i}-\bar{\pmb{x}}), \cdots, \pmb{\xi_k}^T(\pmb{x_i}-\bar{\pmb{x}}))^T
\end{equation}
$$

## 算法流程
* 定义数据$X=\{\pmb{x_1}, \pmb{x_2}, \cdots, \pmb{x_N}\}$以及一个新的维度$k, k<d$
* 计算均值向量$\bar{\pmb{x}}=\frac{1}{N}\sum_{i=1}^N\pmb{x_i}$
* 计算协方差矩阵
$$
\begin{equation}
    Cov(\pmb{x}) = \frac{1}{N}\sum_{i=1}^N(\pmb{x_i}-\bar{\pmb{x}})(\pmb{x_i}-\bar{\pmb{x}})^T
\end{equation}
$$
* 对协方差矩阵进行谱分解，得到d个特征向量及特征值
* 则任意$\pmb{x}\in\mathcal{R}^d$其近似为
$$
\begin{equation}
    \hat{\pmb{x_i}} = \bar{\pmb{x}} + (\pmb{\xi_1}^T(\pmb{x_i}-\bar{\pmb{x}}))\pmb{\xi_1} + (\pmb{\xi_2}^T(\pmb{x_i}-\bar{\pmb{x}}))\pmb{\xi_2} + \cdots + (\pmb{\xi_k}^T(\pmb{x_i}-\bar{\pmb{x}}))\pmb{\xi_k}
\end{equation}
$$
其对应的坐标为
$$
\begin{equation}
    \pmb{y} = (\pmb{\xi_1}^T(\pmb{x_i}-\bar{\pmb{x}}), \pmb{\xi_2}^T(\pmb{x_i}-\bar{\pmb{x}}), \cdots, \pmb{\xi_k}^T(\pmb{x_i}-\bar{\pmb{x}}))^T
\end{equation}
$$



```python
import numpy as np
from sklearn.decomposition import PCA
from sklearn.datasets import make_blobs
import matplotlib.pyplot as plt
```


```python
# 数据数量
sample_num = 1000
# 方差
std_ = 0.1

x_list, _ = make_blobs(n_samples=sample_num, centers=1, 
                            cluster_std=std_, random_state=1024, center_box=(-2, 2))

transformation = [[0.60834549, -0.63667341], [-0.40887718, 0.85253229]]
x_list = np.dot(x_list, transformation)
```


```python
class MyPCA(object):

    def __init__(self, input_dim, reduction_dim):

        self.input_dim = input_dim
        self.reduction_dim = reduction_dim

        self.pca_matrix = None
        self.avg_vector = None

        assert self.reduction_dim <= self.input_dim
    
    def fit(self, data):
        
        input_data = np.array(data)
        assert input_data.shape[1] == self.input_dim

        # cal mean vector
        self.avg_vector = np.mean(input_data, 0)
        # cal cov matrix
        conv_matrix = np.cov(input_data.T)
        # cal top k
        _, topk_e_vecs_mat = self._topk_eigen(conv_matrix, self.reduction_dim)
        self.pca_matrix = topk_e_vecs_mat.T
    
    def transform(self, x):
        input_vec = np.array(x)
        input_vec = np.reshape(x, (self.input_dim))
        return np.matmul(self.pca_matrix, input_vec-self.avg_vector)
    
    def fit_transform(self, data):
        
        # fit
        self.fit(data)
        
        reduction_vecs = list()

        for input_vec in data:
            reduction_vecs.append(self.transform(input_vec))

        return np.array(reduction_vecs)
    
    def _topk_eigen(self, mat, k):
        e_vals, e_vecs = np.linalg.eig(mat)
        sorted_ind = np.argsort(e_vals)

        topk_e_vals = e_vals[sorted_ind[:-1-k:-1]]
        topk_e_vecs = e_vecs[:, sorted_ind[:-1-k:-1]]
        return topk_e_vals, topk_e_vecs

```


```python
# pca

# sklearn pca
sklearn_pca = PCA(n_components=1)
sklearn_pred_x = sklearn_pca.fit_transform(x_list) * sklearn_pca.components_[0] + sklearn_pca.mean_

# my pca
my_pca = MyPCA(2, 1)
my_pred_x = my_pca.fit_transform(x_list) * my_pca.pca_matrix[0] + my_pca.avg_vector

print("sklearn transfer matrix: \n{}".format(sklearn_pca.components_))
print("My transfer matrix: \n{}".format(sklearn_pca.components_))

print("average error:{}".format(np.mean(np.abs(my_pred_x - sklearn_pred_x))))

plt.figure(figsize=(20, 10))
plt.subplot(1, 2, 1)
plt.title("Projection of training data using Sklearn PCA", fontsize=20)
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.scatter(x_list[:, 0], x_list[:, 1], c="#2b5f8e", label="training data")
plt.scatter(sklearn_pred_x[:, 0], sklearn_pred_x[:, 1], c="#cc3536", label="PCA Result")
plt.legend()
plt.subplot(1, 2, 2)
plt.title("Projection of training data using My PCA", fontsize=20)
plt.xticks(fontsize=20)
plt.yticks(fontsize=20)
plt.scatter(x_list[:, 0], x_list[:, 1], c="#2b5f8e", label="training data")
plt.scatter(my_pred_x[:, 0], my_pred_x[:, 1], c="#cc3536", label="PCA Result")
plt.legend()
plt.show()
```

e03c55f9cc75e305dd31830e55d27d1f52a45030
ipynb
Jupyter Notebook
01_PCA/PCA.ipynb
koolo233/dimensionality_reduction_python
452a927772c546f68d6a63e96cdb017b23e4077c
MIT
01_PCA/PCA.ipynb
koolo233/dimensionality_reduction_python
452a927772c546f68d6a63e96cdb017b23e4077c
MIT
01_PCA/PCA.ipynb
koolo233/dimensionality_reduction_python
452a927772c546f68d6a63e96cdb017b23e4077c
MIT
Qwen/Qwen-72B
1. YES
2. YES