16.1.2 定义和导出

【补充解释】 x = ( x 1 , x 2 , ⋯   , x m ) T x=(x_1,x_2,\cdots,x_m)^T x=(x1​,x2​,⋯,xm​)T是 m m m维随机变量，其中 x x x的每个特征 x i x_i xi​都是一个随机变量。均值向量 μ \mu μ即 m m m维随机变量的每个特征的均值，其中 μ i \mu_i μi​是随机变量 x i x_i xi​的均值。

式16.3和式16.4的推导

由式(16.2)和 Σ = c o v ( x , x ) = E [ ( x − μ ) ( x − μ ) T ] \Sigma = cov(x,x) = E[(x-\mu)(x-\mu)^T] Σ=cov(x,x)=E[(x−μ)(x−μ)T]得
KaTeX parse error: No such environment: align at position 8: \begin{̲a̲l̲i̲g̲n̲}̲ var(y_i) & = E…
式(16.3)得证。同理可证式(16.4)。

总体主成分条件(1)为什么能推出正交变换？

根据定理16.1可知， α 1 , α 2 , ⋯   , α m \alpha_1,\alpha_2,\cdots,\alpha_m α1​,α2​,⋯,αm​分别为 x x x的协方差矩阵 Σ \Sigma Σ的特征值对应的单位特征向量。

因为协方差矩阵是实对称矩阵，所以特征值不等的特征向量正交。而若有特征值相等（有重根），对应的特征向量组成 m m m维空间 R m R^m Rm的一个子空间，子空间的维数等于重根数。在子空间任取一个正交坐标系，这个坐标系的单位向量可作为特征向量。

由此推出：

当 n n n个特征值互不相等时， α 1 , α 2 , ⋯   , α m \alpha_1,\alpha_2,\cdots,\alpha_m α1​,α2​,⋯,αm​唯一且两两正交， α \alpha α为正交矩阵，线性变换为正交变换。

当特征值由重根时， α 1 , α 2 , ⋯   , α m \alpha_1,\alpha_2,\cdots,\alpha_m α1​,α2​,⋯,αm​不唯一，但一定存在两两正交的情况，使 α \alpha α为正交矩阵，使线性变换为正交变换。

16.1.3 主要性质

【补充解释】式(16.19)的推导用到的矩阵的迹的性质

n n n阶方阵 A A A的迹，即方阵 A A A的主对角线上各个元素的总和。矩阵 A A A的迹用 t r ( A ) tr(A) tr(A)表示，矩阵的迹具有如下性质：

  • 矩阵的迹等于矩阵特征值的和 -> 相似矩阵的迹相等
  • t r ( A B ) = t r ( B A ) tr(AB)=tr(BA) tr(AB)=tr(BA)

【补充解释】在式(16.22)的推导中，用到如下引理： ∣ ρ X Y ∣ = 1 |\rho_{XY}|=1 ∣ρXY​∣=1时，表示变量 X X X和 Y Y Y为完全线性相关。

最优化问题(16.7)的拉格朗日函数求导

KaTeX parse error: No such environment: align at position 8: \begin{̲a̲l̲i̲g̲n̲}̲ \frac{\partial…

式(16.13)和式(16.14)

式(16.13)即
$$
\begin{bmatrix}
\sigma_{11} & \sigma_{12} & \cdots & \sigma_{1m} \
\sigma_{21} & \sigma_{22} & \cdots & \sigma_{2m} \
\vdots & \vdots & \ddots & \vdots \
\sigma_{m1} & \sigma_{m2} & \cdots & \sigma_{mm} \
\end{bmatrix}

\begin{bmatrix}
\alpha_{1k} \
\alpha_{2k} \
\vdots \
\alpha_{mk} \
\end{bmatrix}

=

\begin{bmatrix}
\sum_{i=1}^m \sigma_{1i} \alpha_{ik} \
\sum_{i=1}^m \sigma_{2i} \alpha_{ik} \
\vdots \
\sum_{i=1}^m \sigma_{mi} \alpha_{ik} \
\end{bmatrix}

=

\begin{bmatrix}
\lambda_k \alpha_{1k} \
\lambda_k \alpha_{2k} \
\vdots \
\lambda_k \alpha_{mk} \
\end{bmatrix}
式 ( 16.14 ) 即 式(16.14)即 式(16.14)即
\begin{bmatrix}
\sigma_{11} & \sigma_{12} & \cdots & \sigma_{1m} \
\sigma_{21} & \sigma_{22} & \cdots & \sigma_{2m} \
\vdots & \vdots & \ddots & \vdots \
\sigma_{m1} & \sigma_{m2} & \cdots & \sigma_{mm} \
\end{bmatrix}

\begin{bmatrix}
\alpha_{11} & \alpha_{12} & \cdots & \alpha_{1m} \
\alpha_{21} & \alpha_{12} & \cdots & \alpha_{2m} \
\vdots & \vdots & \ddots & \vdots \
\alpha_{m1} & \alpha_{m2} & \cdots & \alpha_{mm} \
\end{bmatrix}

=

\begin{bmatrix}
\alpha_{11} & \alpha_{12} & \cdots & \alpha_{1m} \
\alpha_{21} & \alpha_{12} & \cdots & \alpha_{2m} \
\vdots & \vdots & \ddots & \vdots \
\alpha_{m1} & \alpha_{m2} & \cdots & \alpha_{mm} \
\end{bmatrix}

\begin{bmatrix}
\lambda_1 & 0 & 0 & 0 \
0 & \lambda_2 & \cdots & 0 \
\vdots & \vdots & \ddots & \vdots \
0 & 0 & \cdots & \lambda_m \
\end{bmatrix}

=

\begin{bmatrix}
\lambda_k \alpha_{11} & \lambda_k \alpha_{12} & \cdots & \lambda_k \alpha_{1m} \
\lambda_k \alpha_{21} & \lambda_k \alpha_{22} & \cdots & \lambda_k \alpha_{2m} \
\vdots & \vdots & \ddots & \vdots \
\lambda_k \alpha_{m1} & \lambda_k \alpha_{m2} & \cdots & \lambda_k \alpha_{mm} \
\end{bmatrix}
$$

16.1.4 主成分的个数

定理16.2的证明思路

首先，因为 A A A是正交矩阵，所以可以通过线性变换 C C C，将 A A A映射到 B B B。将 C C C设为变量，则 B B B可通过 C C C确定，同时 Σ y \Sigma_y Σy​也可以通过 C C C确定。此时，求解使 t r ( Σ y ) tr(\Sigma_y) tr(Σy​)取得最大值的线性变换 C C C即可。只要证明 C C C是对角矩阵的前 q q q列，则可以证明 B = A q B=A_q B=Aq​时 t r ( Σ y ) tr(\Sigma_y) tr(Σy​)取得最大值。

16.1.5 规范化变量的总体主成分

【补充说明】规范化随机变量的总体主成分的5条性质，一一对应于前文所述总体主成分的性质。

16.2.1 样本主成分的定义和性质

【勘误补充】式(16.40)中的 x ˉ \bar{x} xˉ应该加粗以表示向量。

【补充解释】式(16.49)中的 X = [ x 1 x 2 ⋯ x n ] X = \begin{bmatrix} x_1 & x_2 & \cdots & x_n\end{bmatrix} X=[x1​​x2​​⋯​xn​​]。之所以协方差矩阵可以写成 S = R = 1 n − 1 X X T S = R = \frac{1}{n-1} X X^T S=R=n−11​XXT，是因为规范化后的所有随机变量的均值均为0。

16.2.2 相关矩阵的特征值分解算法

【勘误补充】在例16.1中，表16.1原本应包含更多位小数，但表中被略去了一部分。这导致仅用两位小数计算出的特征值和特征向量与书中结果不同。仅使用表中的两位小数计算的特征值和特征向量如下：

  • 第1个特征值：2.17016506；对应特征向量：[-0.45990769,-0.4763124,-0.52874972,-0.53106981]；
  • 第2个特征值：0.87100546；对应特征向量：[-0.56790937,-0.49090704,0.47557056，0.45860862]；
  • 第3个特征值：0.56617908；对应特征向量：[-0.6665586,0.71535364,0.11286206,-0.17672284]；
  • 第4个特征值：0.3926504；对应特征向量：[-0.14718523,0.14284941,-0.69391536,0.69022607]。

相关矩阵求解主成分及其因子载荷量和贡献率（原生Python实现）

def pca_by_feature(R, need_accumulative_contribution_rate=0.75):
    """协方差矩阵/相关矩阵求解主成分及其因子载荷量和贡献率（打印到控制台）

    :param R: 协方差矩阵/相关矩阵
    :param need_accumulative_contribution_rate: 需要达到的累计方差贡献率
    :return: None
    """
    n_features = len(R)

    # 求解相关矩阵的特征值和特征向量
    features_value, features_vector = np.linalg.eig(R)

    # 依据特征值大小排序特征值和特征向量
    z = [(features_value[i], features_vector[:, i]) for i in range(n_features)]
    z.sort(key=lambda x: x[0], reverse=True)
    features_value = [z[i][0] for i in range(n_features)]
    features_vector = np.hstack([z[i][1][:, np.newaxis] for i in range(n_features)])

    # 计算所需的主成分数量
    total_features_value = sum(features_value)  # 特征值总和
    need_accumulative_contribution_rate *= total_features_value
    n_principal_component = 0  # 所需的主成分数量
    accumulative_contribution_rate = 0
    while accumulative_contribution_rate < need_accumulative_contribution_rate:
        accumulative_contribution_rate += features_value[n_principal_component]
        n_principal_component += 1

    # 输出单位特征向量和主成分的方差贡献率
    print("【单位特征向量和主成分的方差贡献率】")
    for i in range(n_principal_component):
        print("主成分:", i,
              "方差贡献率:", features_value[i] / total_features_value,
              "特征向量:", features_vector[:, i])

    # 计算各个主成分的因子载荷量：factor_loadings[i][j]表示第i个主成分对第j个变量的相关关系，即因子载荷量
    factor_loadings = np.vstack(
        [[np.sqrt(features_value[i]) * features_vector[j][i] / np.sqrt(R[j][j]) for j in range(n_features)]
         for i in range(n_principal_component)]
    )

    # 输出主成分的因子载荷量和贡献率
    print("【主成分的因子载荷量和贡献率】")
    for i in range(n_principal_component):
        print("主成分:", i, "因子载荷量:", factor_loadings[i])
    print("所有主成分对变量的贡献率:", [np.sum(factor_loadings[:, j] ** 2) for j in range(n_features)])

【测试】例16.1

if __name__ == "__main__":
    X = np.array([[1, 0.44, 0.29, 0.33],
                  [0.44, 1, 0.35, 0.32],
                  [0.29, 0.35, 1, 0.60],
                  [0.33, 0.32, 0.60, 1]])

    pca_by_feature(X)

    # 【单位特征向量和主成分的方差贡献率】
    # 主成分: 0 方差贡献率: 0.542541266192316 特征向量: [-0.45990769 -0.4763124  -0.52874972 -0.53106981]
    # 主成分: 1 方差贡献率: 0.21775136378517398 特征向量: [-0.56790937 -0.49090704  0.47557056  0.45860862]
    # 【主成分的因子载荷量和贡献率】
    # 主成分: 0 因子载荷量: [-0.6775121  -0.70167866 -0.77892661 -0.78234444]
    # 主成分: 1 因子载荷量: [-0.5300166  -0.45815211  0.44383894  0.42800876]
    # 所有主成分对变量的贡献率: [0.7399402466295502, 0.7022563027546291, 0.8037196579404534, 0.7952543125853273]

16.2.3 数据矩阵的奇异值分解算法

数据矩阵奇异值分解进行的主成分分析算法（numpy实现）

def pca_by_svd(X, k):
    """数据矩阵奇异值分解进行的主成分分析算法

    :param X: 样本矩阵X
    :param k: 主成分个数k
    :return:
    """
    n_samples = X.shape[1]

    # 构造新的n×m矩阵
    T = X.T / np.sqrt(n_samples - 1)

    # 对矩阵T进行截断奇异值分解
    U, S, V = np.linalg.svd(T)
    V = V[:, :k]

    # 求k×n的样本主成分矩阵
    return np.dot(V.T, X)

【测试】习题16.1

if __name__ == "__main__":
    X = np.array([[2, 3, 3, 4, 5, 7],
                  [2, 4, 5, 5, 6, 8]])
    X = X.astype("float64")

    # 规范化变量
    avg = np.average(X, axis=1)
    var = np.var(X, axis=1)
    for i in range(X.shape[0]):
        X[i] = (X[i, :] - avg[i]) / np.sqrt(var[i])

    print(pca_by_svd(X, 2))

    # [[-2.02792041 -0.82031104 -0.4330127   0.          0.82031104  2.46093311]
    #  [ 0.2958696  -0.04571437 -0.4330127   0.          0.04571437  0.1371431 ]]