数据集、源文件可以在Github项目中获得
链接: https://github.com/Raymond-Yang-2001/AndrewNg-Machine-Learing-Homework

1. 主成分分析

1.1 数据降维动机

在某些情况下，机器学习使用的数据可能十分巨大，这种情况下，对数据进行压缩就显得十分必要。假设某一数据集包括两个特征维度：圆形的半径和圆形的面积。根据简单的数学知识，这两个特征实际上是互相关联的，从描述圆特征的角度来说，半径和面积并没有什么不同；在机器学习中，完全可以只是用一个特征进行学习。

这就是数据降维，使得数据中某些“冗余”的特征维度被压缩，从而得到更低位的数据表示——当然，这种压缩是建立在尽可能减小数据信息的损失上的。

1.2 PCA降维目标问题分析

假设我们有如下数据：

如果要对数据进行维度压缩，从二维变成一位，显然水平方向是更适合作为压缩后的维度的。从数学角度来看，是因为水平方向的方差要大于垂直方向，也就意味着，水平方向的特征相对于垂直方向来说，具有更大的特征空间，有更丰富的特征的值的选则，也就是对样本更有区分度。

所以，PCA的优化目标可以总结为如下两个：

1. 降维后同一维度的方差最大
2. 不同维度间的相关度为0

如果不同维度的特征间存在相关度，则证明这两个特征间存在某种内在的关联，从一个特征能推理出另一个特征，显然这没有达到降维的最好目标。

2. PCA数学原理分析

为了同时表示同一维度间的方差和不同维度间的相关性，我们引入了协方差矩阵。设数据维度为$D$，则协方差矩阵为：
$$C={\left[\begin{array}{cccc}cov(x_1,x_1)& cov(x_1,x_2)& \cdots & cov(x_1,x_D)\\ ⋮& ⋮& \ddots & ⋮\\ cov(x_D,x_1)& cov(x_D,x_2)& \cdots & cov(x_D,x_D)\end{array}\right]}_{D\times D}$$
对角线上的元素表明了同一维度的方差，不同维度间的方差表明了不同维度间的相关性。优化后的理想矩阵为：
$$C^{\prime }={\left[\begin{array}{cccc}{\delta }_{11}& 0& \cdots & 0\\ 0& {\delta }_{22}& \cdots & 0\\ 0& 0& \ddots & ⋮\\ 0& 0& \cdots & {\delta }_{RR}\end{array}\right]}_{R\times R}$$
并且，$X$为原数据矩阵，形状为$(N,D)$，$P$为降维变换矩阵。

考虑协方差矩阵的计算：
$$C^{\prime }=\frac{1}{N}(XP{)}^{\top }(XP)=\frac{1}{N}{P}^{\top }{X}^{\top }XP=P^{\top }\left(\frac{1}{N}{X}^{\top }X\right)P=P^{\top }CP$$
要尽量减少不同维度间的相关度，设$P^{\top }P=I$，即$P$的列是正交的。

优化问题如下：
$$\{\begin{array}{rl}& \max \mathrm{tr}({\mathrm{P}}^{\top }\mathrm{CP})\\ & P^{\top }P=I\end{array}$$

采用拉格朗日算子，得到无条件的极值问题：
$$f(P)=\mathrm{tr}({\mathrm{P}}^{\top }\mathrm{CP})+\lambda (\mathrm{I}-{\mathrm{P}}^{\top }\mathrm{P})$$
对$f(P)$求导：
$$\frac{\partial tr(P^{\top }CP)}{\partial P}=\frac{\partial tr(P{P}^{\top }C)}{\partial P}=(P^{\top }C{)}^{\top }=C^{\top }P$$
$$\frac{\partial \lambda (I-P^{\top }P)}{\partial P}=-\lambda P$$
这里用到了迹和矩阵导数的求导的一些知识，不理解得到话可以自行查阅一些线性代数的知识。
$$\frac{\partial f(P)}{\partial P}=C^{\top }P+\lambda P$$
由于$C$是对阵矩阵，所以：
$$\frac{\partial f(P)}{\partial P}=CP-\lambda P$$
$$\frac{\partial f(P)}{\partial P}=0\Rightarrow CP=\lambda P$$
到这里就可以看出，$P$实际上是$C$的特征向量（矩阵）。事实上，对称矩阵特征值不等的特征向量是正交的。

另外，有：
$$P^{\top }CP=P^{\top }\lambda P$$
由于特征值矩阵$\lambda$是一个对角矩阵（除了对角线均为0），所以$C=P^{\top }\lambda P=\lambda P^{\top }P=\lambda$

可以发现，$\lambda$实际上就代表了各个维度的方差，也就是说特征值=方差。选择大方差的维度保留就是保留特征值较大的维度。

2.1 求协方差矩阵的碎碎念

先前我们提到，求对一个$(n\_sample=N,n\_feature=D)$的矩阵$X$求协方差矩阵，公式如下：
$$C=\frac{1}{N}{X}^{\top }X$$
接下来将进行理论推导，为什么协方差要这么求：

从概率论的知识可知，两个随机变量$X,Y$的协方差如下：
$$cov(X,Y)=\mathbb{E}[(X-{\mu }_X)(Y-{\mu }_Y)]$$
其中，${\mu }_X,{\mu }_Y$分别是两个随机变量的均值。

我们在对数据进行PCA之前，先要对数据进行规范化，使得其均值为0，通常采用方法是减去均值，即$X_{new}=X-{\mu }_X$，这样在计算协方差的时候，我们可以得到：
$$cov(X,Y)=\mathbb{E}(XY)$$
这里，$X,Y$都是数据的某一维度。转换成矩阵运算，我们可以得到：
$$Cov=\frac{1}{N}{X}^{\top }X$$

当然，我们如果考虑对样本进行无偏估计，在求期望的时候，均值会改写为$\mathbb{E}(\mathbb{X})=\frac{1}{N-1}{X}_i$，（不熟悉的话可以回顾一下概率论的相关知识），这样我们得到的协方差矩阵是：
$$Cov=\frac{1}{N-1}{X}^{\top }X$$
在numpy库中，采用的计算方法就是这种。

2.2 PCA实现方法

由上文可知，实现PCA的步骤如下：

1. 对输入求协方差矩阵$C$
2. 求矩阵$C$的特征值和特征向量，按照特征值降序排序
3. 保留k个维度，选取前k个特征向量组成矩阵$P$
4. 新的数据矩阵为$C^{\prime }=XP$

3. Python实现

import numpy as np


class PCA:
    """
    ----------------------------------------------------------------------------
    Attributes:
    components_: ndarray with shape of (n_features, n_components)
        Principal axes in feature space, representing the directions of maximum
        variance in the data.

    explained_variance_ : ndarray of shape (n_components,)
        The amount of variance explained by each of the selected components.

    explained_variance_ratio_ : ndarray of shape (n_components,)
        Percentage of variance explained by each of the selected components.
    """
    def __init__(self, n_components):
        self.n_components = n_components
        self.explained_variance_ = None
        self.explained_variance_ratio_ = None
        self.components_ = None

        self.__mean = None

    def fit(self, x):
        """
        :param x: (n,d)
        :return:
        """
        self.__mean = x.mean(axis=0)
        x_norm = x - self.__mean
        x_cov = (x_norm.T @ x_norm) / (x.shape[0] - 1)
        vectors, variance, _ = np.linalg.svd(x_cov)
        # (n_feature, n_components)
        self.components_ = vectors[:, :self.n_components]
        if len(self.components_.shape) == 1:
            self.components_ = np.expand_dims(vectors[:, :self.n_components], axis=1)
        self.explained_variance_ = variance[:self.n_components]
        self.explained_variance_ratio_ = self.explained_variance_ / variance.sum()

    def transform(self, x):
        """
        :param x: (n, n_feature)
        :return:
        """
        if self.__mean is not None:
            x = x - self.__mean
        x_transformed = x @ self.components_
        return x_transformed

数据集可视化

进行PCA降维，可视化数据

from PCA import PCA
m_pca = PCA(1)
m_pca.fit(x)
x_reduced = m_pca.transform(x)
print("The principal axes of components is: {}\n"
      "The variance of each components is: {}\n"
      "The variance ratio of selected components is: {}"
      .format(m_pca.components_.tolist(), m_pca.explained_variance_.tolist(), m_pca.explained_variance_ratio_.tolist()))

The principal axes of components is: [[-0.7690815341368202], [-0.6391506816469459]]
The variance of each components is: [2.1098781795840327]
The variance ratio of selected components is: [0.8706238489732337]

3.1 进行人脸数据压缩

展示人脸数据

faces = loadmat('data/ex7faces.mat')
face_x = faces['X']
face_x.shape

(5000, 1024)

可以看到共5000张人脸数据，维度是1024维度。

使用PCA降维到100维度

from PCA import PCA
face_pca = PCA(n_components=100)
face_pca.fit(face_x)
face_r = face_pca.transform(face_x)

可以看到保留了大部分的人脸特征。