主成分分析的笔记与思考

主成分分析（Principle Component Analysis）

主成分分析（Principle Component Analysis），旨在找到高维数据中的主成分，并利用这些主成分来表示原始高维数据，从而达到降维的目的。

例如，一个三维空间中有$n$个样本数据，每个样本$x_i$的特征可以表示为$(x,y,z)$，因此$X\in {\mathbb{R}}^{n\times 3}$,但实际上在三维空间中它只出现在一个平面里，因此，如果我们可以通过坐标系旋转使数据所在的平面与三维空间的一个平面$x^{\prime },y^{\prime }$或$y^{\prime },z^{\prime }$或$x^{\prime },z^{\prime }$平面重合，就可以利用这的任意两个维度表示所有的原始数据，并且没有任何信息的损失，达到三维降到二维的效果。$x^{\prime },y^{\prime }$等就是我们要求的主成分。

如图，在二维空间里的中心化后的一组数据，我们可以大概看出主成分所在的轴（斜着的红线），在这个轴上，数据分布的更加分散，这也意味着数据在这个方向上的方差更大。

通常我们认为，数据的方差越大越好，噪声的方差越小越好，因此主成分分析的目标就是 最大化投影方差，使得数据点在主轴上投影的方差最大。

对于一组给定的数据$v_1,...,v_n,v_i\in {\mathbb{R}}^m$，中心化后$x_1,x_2,...,x_n=v_1-\mu ,v_2-\mu ,...,v_n-\mu$，其中$\mu =\frac{1}{n}{\sum }_{i=1}^n{v}_i$。

由于向量内积的几何意义可以表示为第一个向量投影到第二个向量的长度，因此向量$x_i$在$w$（单位方向向量）上的投影坐标可以表示为$⟨x_i,w⟩=x_i^{T}w$。所以目标是找到一个投影方向$w$，使得$x_1,x_2,...,x_n$在$w$上的投影长度（方差）尽可能大。

由于中心化后，投影的均值为$0$，（因为中心化后${\mu }^{\prime }=\frac{1}{n}{\sum }_{i=1}^n{x}_i^{T}w=(\frac{1}{n}{\sum }_{i=1}^n{x}_i^T)w=(\frac{1}{n}{\sum }_{i=1}^n(v_i-\mu {)}^T)w=0$），因此投影后的方差就可以写为

$$\begin{gathered} D(x)=\frac{1}{n}\sum _{i=1}^n(x_i^{T}w{)}^2=\frac{1}{n}\sum _{i=1}^n(x_i^{T}w{)}^T(x_i^{T}w) \\ =\frac{1}{n}\sum _{i=1}^n{w}^T{x}_i{x}_i^{T}w \\ =w^T(\frac{1}{n}\sum _{i=1}^n{x}_i{x}_i^T)w \end{gathered}$$

其中$\frac{1}{n}{\sum }_{i=1}^n{x}_i{x}_i^T$如果熟悉概率论的话，不难看出是样本的协方差矩阵；我们将其写作$\Sigma$。$w$是单位向量，所以$w^{T}w=1$。因此我们要求最大化投影方差，可以表示为

$$\begin{gathered} \max w^T\Sigma w \\ \mathrm{s}.\mathrm{t}.w^{T}w=1 \end{gathered}$$

求解该优化问题，我们使用拉格朗日乘子法，
$$\begin{gathered} F(w)=w^T\Sigma w-\lambda (w^{T}w-1) \\ \frac{\partial F}{\partial w}=0 \\ \Sigma w-\lambda w=0 \end{gathered}$$
得$\Sigma w=\lambda w$，此时$D(x)=w^T\Sigma w=w^T\lambda w=\lambda$。根据线性代数的知识，投影后的方差就是协方差矩阵的特征值，要求的最大化方差也就是协方差矩阵的最大特征值，最佳投影方向就是最大特征值对应的特征向量。次佳投影方向位于最佳投影方向的正交空间，也就是次大特征值对应的特征向量（特征向量彼此正交），以此类推。因此，我们可以得到主成分分析的求解步骤如下：

---

Algorithm: Principle Component Analysis

• 对样本数据进行中心化处理
• 求样本的协方差矩阵
• 对协方差矩阵进行特征值分解，将特征值从大到小排列
• 取前d个最大的特征值对应的特征向量$w_1,w_2,...w_d,W=[w_1,...,w_d]\in {\mathbb{R}}^{m\times d}$，设$z_i$是$x_i$降维后的表征，$z_i=\left[\begin{array}{c}{w}_1^T{x}_i\\ w_2^T{x}_i\\ ⋮\\ w_d^T{x}_i\end{array}\right]=W^T{x}_i$。

---

python的代码可以参考如下：

def pca(xs: np.ndarray, n_pc: int = 2) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
    """
    Implement PCA via eigen-decomposition
    :param xs: a data matrix with (N, D), N is the number of samples, D is the dimension of sample space
    :param n_pc: the number of principal components we would like to output
    :return:
        the matrix containing top-k principal components, with size (D, n_pc)
        the vector indicating the top-k eigenvalues, with size (n_pc,)
        the data recovered from the projections along the principal components, with size (N, D)
        the zero-mean data with size (N, D)
    """
    xs_zero_mean = xs - np.mean(xs, axis = 0)
    X = xs.T @ xs / (xs.shape[0] - 1) # calculate the covariance matrix
    eig, w = np.linalg.eig(X) # Eigen-decomposition
    idx = eig.argsort()[: : -1] # to sort the eig, we get the index
    eig = eig[idx] 
    w = w[:, idx]
    W = w[: , 0 : n_pc] # the matrix containing top-k principal components, with size (D, n_pc)
    EIG = eig[0 : n_pc] # the vector indicating the top-k eigenvalues
    xs_r = xs_zero_mean @ W @ W.T + np.mean(xs, axis = 0)
    return W, EIG, xs_r, xs_zero_mean