5. 奇异值分解#
5.1. 概述#
奇异值分解（SVD）是一个强大的数学工具，在数据分析和机器学习中有着广泛的应用。它不仅可以用于最小二乘投影，还是许多统计和机器学习方法的基石。
在这一讲中，我们将首先介绍SVD的基本概念，然后探讨它与以下几个重要领域的关系：
线性代数中的四个基本子空间
最小二乘回归（包括欠定和超定情况）
主成分分析（PCA）
在后续的动态模式分解讲座中，我们会看到如何利用SVD来高效地计算向量自回归（VAR）模型的简化形式。
5.2. 基本设定#
假设我们有一个\(m \times n\)的矩阵\(X\)，其秩为\(p\)。
显然，\(p\)不会超过\(m\)和\(n\)的较小值，即\(p \leq \min(m,n)\)。
在实际应用中，\(X\)通常代表一个数据矩阵：
每一列代表一个观测单位（可以是某个时间点的数据，或是某个个体的信息）
每一行代表一个变量（用来描述观测单位的某个特征或属性）
我们将关注两种情况：
矮胖情况，即\(m << n\)，表示列数（个体）远多于行数（属性）。
高瘦情况，即\(m >> n\)，表示行数（属性）远多于列数（个体）。
我们将在这两种情况下对\(X\)进行奇异值分解。
在 \(m << n\) 的情况下，即个体数量 \(n\) 远大于属性数量 \(m\) 时，我们可以通过对观测值函数取平均来计算联合分布的样本矩。
在这种 \(m << n\) 的情况下，我们将使用奇异值分解来进行主成分分析(PCA)以寻找模式。
在 \(m >> n\) 的情况下，即属性数量 \(m\) 远大于个体数量 \(n\)，且在时间序列环境中 \(n\) 等于数据集 \(X\) 中所覆盖的时间段数量时，我们将采用不同的方法。
在这种情况下，我们仍然使用奇异值分解，但目的是构建动态模态分解(DMD)。
5.3. 奇异值分解#
一个秩为 \(p \leq \min(m,n)\) 的 \(m \times n\) 矩阵 \(X\) 的奇异值分解为：
其中：
且
\(U\) 是 \(X\) 的 \(m \times m\) 正交矩阵，由左奇异向量组成
\(U\) 的列是 \(X X^\top \) 的特征向量
\(V\) 是 \(X\) 的 \(n \times n\) 正交矩阵，由右奇异向量组成
\(V\) 的列是 \(X^\top X\) 的特征向量
\(\Sigma\) 是一个 \(m \times n\) 矩阵，其主对角线上的前 \(p\) 个位置是正数 \(\sigma_1, \sigma_2, \ldots, \sigma_p\)，称为奇异值；\(\Sigma\) 的其余元素都为零
这 \(p\) 个奇异值是 \(m \times m\) 矩阵 \(X X^\top \) 以及 \(n \times n\) 矩阵 \(X^\top X\) 的特征值的正平方根
我们约定，当 \(U\) 是复值矩阵时，\(U^\top \) 表示 \(U\) 的共轭转置或厄米特转置，即 \(U_{ij}^\top \) 是 \(U_{ji}\) 的复共轭。
类似地，当 \(V\) 是复值矩阵时，\(V^\top\) 表示 \(V\) 的共轭转置或厄米特转置
矩阵 \(U,\Sigma,V\) 通过以下方式对向量进行线性变换：
用酉矩阵 \(U\) 和 \(V\) 乘以向量会使其旋转，但保持向量之间的角度和向量的长度不变。
用对角矩阵 \(\Sigma\) 乘以向量会保持向量之间的角度不变，但会重新缩放向量。
因此，表示式 (5.1) 表明，用 \(m \times n\) 矩阵 \(X\) 乘以 \(n \times 1\) 向量 \(y\) 相当于按顺序执行以下三个乘法运算：
通过计算 \(V^\top y\) 来旋转 \(y\)
通过乘以 \(\Sigma\) 来重新缩放 \(V^\top y\)
通过乘以 \(U\) 来旋转 \(\Sigma V^\top y\)
\(m \times n\) 矩阵 \(X\) 的这种结构为构建系统开启了大门
数据编码器和解码器。
因此，
\(V^\top y\) 是一个编码器
\(\Sigma\) 是一个应用于编码数据的运算符
\(U\) 是一个解码器，用于处理将运算符 \(\Sigma\) 应用于编码数据后的输出
我们将在本讲稍后研究动态模态分解时应用这些概念。
未来路线
我们上面描述的是所谓的完全 SVD。
在完全 SVD中，\(U\)、\(\Sigma\) 和 \(V\) 的形状分别为 \(\left(m, m\right)\)、\(\left(m, n\right)\)、\(\left(n, n\right)\)。
稍后我们还将描述经济型或简化 SVD。
在研究简化 SVD之前，我们将进一步讨论完全 SVD的性质。
5.4. 四个基本子空间#
让 \({\mathcal C}\) 表示列空间，\({\mathcal N}\) 表示零空间，\({\mathcal R}\) 表示行空间。
让我们首先回顾一下秩为 \(p\) 的 \(m \times n\) 矩阵 \(X\) 的四个基本子空间。
列空间\(X\)，记作\({\mathcal C}(X)\)，是\(X\)的列向量的张成空间，即所有可以写成\(X\)的列向量的线性组合的向量\(y\)。其维数为\(p\)。
零空间\(X\)，记作\({\mathcal N}(X)\)，包含所有满足\(Xy=0\)的向量\(y\)。其维数为\(n-p\)。
行空间\(X\)，记作\({\mathcal R}(X)\)，是\(X^\top\)的列空间。它包含所有可以写成\(X\)的行向量的线性组合的向量\(z\)。其维数为\(p\)。
左零空间\(X\)，记作\({\mathcal N}(X^\top)\)，包含所有满足\(X^\top z=0\)的向量\(z\)。其维数为\(m-p\)。
对于矩阵\(X\)的完全奇异值分解，左奇异向量矩阵\(U\)和右奇异向量矩阵\(V\)包含了所有四个子空间的正交基。
它们形成两对正交子空间，我们现在来描述。
令\(u_i, i = 1, \ldots, m\)为\(U\)的\(m\)个列向量，令
设 \(v_i, i = 1, \ldots, n\) 为 \(V\) 的 \(n\) 个列向量。
让我们将 X 的完整奇异值分解写作
其中 \(\Sigma_p\) 是一个 \(p \times p\) 对角矩阵，对角线上是 \(p\) 个奇异值，且
表示式 (5.2) 意味着
或
或
方程 (5.4) 说明了变换 \(X\) 如何将一对正交单位向量 \(v_i, v_j\)（其中 \(i\) 和 \(j\) 都小于或等于 \(X\) 的秩 \(p\)）映射到一对正交单位向量 \(u_i, u_j\)。
方程 (5.3) 表明
对表示式 (5.2) 两边取转置得到
或
或
注意方程 (5.6) 表明变换 \(X^\top\) 将一对不同的正交单位向量 \(u_i, u_j\)（其中 \(i\) 和 \(j\) 都小于或等于 \(X\) 的秩 \(p\)）映射到一对不同的正交单位向量 \(v_i, v_j\)。
方程 (5.5) 表明：
因此，方程组 (5.3) 和 (5.5) 共同描述了 \(X\) 的四个基本子空间，如下所示：
由于 \(U\) 和 \(V\) 都是正交矩阵，集合 (5.7) 表明
\(U_L\) 是 \(X\) 列空间的标准正交基
\(U_R\) 是 \(X^\top\) 零空间的标准正交基
\(V_L\) 是 \(X\) 行空间的标准正交基
\(V_R\) 是 \(X\) 零空间的标准正交基
我们通过执行(5.2)右侧要求的乘法并读取结果，已经验证了(5.7)中的四个声明。
(5.7)中的声明以及\(U\)和\(V\)都是酉矩阵（即正交矩阵）这一事实意味着：
\(X\)的列空间与\(X^\top\)的零空间正交
\(X\)的零空间与\(X\)的行空间正交
这些性质有时用以下两对正交补空间来描述：
\({\mathcal C}(X)\)是\({\mathcal N}(X^\top)\)的正交补
\({\mathcal R}(X)\)是\({\mathcal N}(X)\)的正交补
让我们看一个例子。
import numpy as np
import numpy.linalg as LA
import matplotlib.pyplot as plt
import matplotlib as mpl
FONTPATH = "fonts/SourceHanSerifSC-SemiBold.otf"
mpl.font_manager.fontManager.addfont(FONTPATH)
plt.rcParams['font.family'] = ['Source Han Serif SC']
导入这些模块后，让我们来看示例。
np.set_printoptions(precision=2)
# 定义矩阵
A = np.array([[1, 2, 3, 4, 5],
[2, 3, 4, 5, 6],
[3, 4, 5, 6, 7],
[4, 5, 6, 7, 8],
[5, 6, 7, 8, 9]])
# 计算矩阵的奇异值分解
U, S, V = np.linalg.svd(A,full_matrices=True)
# 计算矩阵的秩
rank = np.linalg.matrix_rank(A)
# 打印矩阵的秩
print("矩阵的秩:\n", rank)
print("S: \n", S)
# 计算四个基本子空间
row_space = U[:, :rank]
col_space = V[:, :rank]
null_space = V[:, rank:]
left_null_space = U[:, rank:]
print("U:\n", U)
print("列空间:\n", col_space)
print("左零空间:\n", left_null_space)
print("V.T:\n", V.T)
print("行空间:\n", row_space.T)
print("右零空间:\n", null_space.T)
矩阵的秩:
2
S:
[2.69e+01 1.86e+00 8.62e-16 5.26e-16 2.77e-17]
U:
[[-0.27 -0.73 -0.53 -0.34 0.03]
[-0.35 -0.42 0.24 0.8 -0.07]
[-0.43 -0.11 0.67 -0.41 0.43]
[-0.51 0.19 0.09 -0.22 -0.8 ]
[-0.59 0.5 -0.46 0.17 0.4 ]]
列空间:
[[-0.27 -0.35]
[ 0.73 0.42]
[-0.03 0.16]
[-0.51 0.82]
[ 0.37 0.06]]
左零空间:
[[-0.53 -0.34 0.03]
[ 0.24 0.8 -0.07]
[ 0.67 -0.41 0.43]
[ 0.09 -0.22 -0.8 ]
[-0.46 0.17 0.4 ]]
V.T:
[[-0.27 0.73 -0.03 -0.51 0.37]
[-0.35 0.42 0.16 0.82 0.06]
[-0.43 0.11 0.25 -0.23 -0.83]
[-0.51 -0.19 -0.84 0.04 -0.02]
[-0.59 -0.5 0.46 -0.12 0.41]]
行空间:
[[-0.27 -0.35 -0.43 -0.51 -0.59]
[-0.73 -0.42 -0.11 0.19 0.5 ]]
右零空间:
[[-0.43 0.11 0.25 -0.23 -0.83]
[-0.51 -0.19 -0.84 0.04 -0.02]
[-0.59 -0.5 0.46 -0.12 0.41]]
5.5. Eckart-Young定理#
假设我们要构造一个\(m \times n\)矩阵\(X\)的最佳秩\(r\)近似。
这里的最佳，指的是在所有秩为\(r < p\)的矩阵中，找到一个矩阵\(X_r\)使得以下范数最小：
其中\(|| \cdot ||\)表示矩阵\(X\)的范数，且\(X_r\)属于所有维度为\(m \times n\)的秩\(r\)矩阵空间。
一个\(m \times n\)矩阵\(X\)的三种常用矩阵范数可以用\(X\)的奇异值表示：
谱范数或\(l^2\)范数 \(|| X ||_2 = \max_{||y|| \neq 0} \frac{||X y ||}{||y||} = \sigma_1\)
Frobenius范数 \(||X ||_F = \sqrt{\sigma_1^2 + \cdots + \sigma_p^2}\)
核范数 \( || X ||_N = \sigma_1 + \cdots + \sigma_p \)
Eckart-Young定理指出，对于这三种范数，最佳的秩\(r\)矩阵是相同的，等于：
这是一个非常强大的定理，它表明我们可以将一个非满秩的 \(m \times n\) 矩阵 \(X\) 通过SVD分解，用一个满秩的 \(p \times p\) 矩阵来最佳近似。
此外，如果这些 \(p\) 个奇异值中有些携带的信息比其他的更多，而我们想用最少的数据获得最多的信息，我们可以取按大小排序的 \(r\) 个主要奇异值。
在介绍主成分分析时，我们会对此进行更详细的讨论。
你可以在这里阅读关于Eckart-Young定理及其应用的内容。
在讨论主成分分析(PCA)和动态模态分解(DMD)时，我们将会用到这个定理。
5.6. 完全SVD和简化SVD#
到目前为止，我们描述的是完全SVD的性质，其中 \(U\)、\(\Sigma\) 和 \(V\) 的形状分别为 \(\left(m, m\right)\)、\(\left(m, n\right)\)、\(\left(n, n\right)\)。
有一种替代性的矩阵分解记法，称为经济型或简化型 SVD，其中 \(U, \Sigma\) 和 \(V\) 的形状与完全SVD中的不同。
注意，因为我们假设 \(X\) 的秩为 \(p\)，所以只有 \(p\) 个非零奇异值，其中 \(p=\textrm{rank}(X)\leq\min\left(m, n\right)\)。
简化型 SVD利用这一事实，将 \(U\)、\(\Sigma\) 和 \(V\) 表示为形状分别为 \(\left(m, p\right)\)、\(\left(p, p\right)\)、\(\left(n, p\right)\) 的矩阵。
你可以在这里了解简化型和完全型SVD https://numpy.org/doc/stable/reference/generated/numpy.linalg.svd.html
对于完全型SVD，
但这些性质并非都适用于简化型 SVD。
哪些性质成立取决于我们处理的是高瘦型矩阵还是矮胖型矩阵。
在高瘦型情况下，即 \(m > > n\)，对于简化型 SVD
在短胖情况下(即 \(m < < n\)),对于简化SVD
当我们研究动态模态分解时,我们需要记住这些性质,因为我们会使用简化SVD来计算一些DMD表示。
让我们做一个练习来比较完全和简化SVD。
回顾一下,
在完全SVD中
\(U\) 是 \(m \times m\) 维
\(\Sigma\) 是 \(m \times n\) 维
\(V\) 是 \(n \times n\) 维
在简化SVD中
\(U\) 是 \(m \times p\) 维
\(\Sigma\) 是 \(p \times p\) 维
\(V\) 是 \(n \times p\) 维
首先,让我们研究一个 \(m = 5 > n = 2\) 的情况。
(这是我们在研究动态模态分解时会遇到的高瘦情况的一个小例子。)
import numpy as np
X = np.random.rand(5,2)
U, S, V = np.linalg.svd(X,full_matrices=True) # 完整SVD
Uhat, Shat, Vhat = np.linalg.svd(X,full_matrices=False) # 简化SVD
print('U, S, V =')
U, S, V
U, S, V =
(array([[-0.22, 0.67, 0.26, -0.6 , 0.27],
[-0.63, -0.1 , -0.46, -0.32, -0.53],
[-0.24, -0.45, 0.81, -0.18, -0.22],
[-0.65, 0.26, 0.15, 0.68, 0.17],
[-0.29, -0.51, -0.21, -0.22, 0.75]]),
array([1.85, 0.38]),
array([[-0.63, -0.78],
[-0.78, 0.63]]))
print('Uhat, Shat, Vhat = ')
Uhat, Shat, Vhat
Uhat, Shat, Vhat =
(array([[-0.22, 0.67],
[-0.63, -0.1 ],
[-0.24, -0.45],
[-0.65, 0.26],
[-0.29, -0.51]]),
array([1.85, 0.38]),
array([[-0.63, -0.78],
[-0.78, 0.63]]))
rr = np.linalg.matrix_rank(X)
print(f'X的秩 = {rr}')
X的秩 = 2
性质：
当\(U\)通过完全SVD构造时，\(U^\top U = I_{m\times m}\) 且 \(U U^\top = I_{m \times m}\)
当\(\hat U\)通过简化SVD构造时，虽然\(\hat U^\top \hat U = I_{p\times p}\)，但\(\hat U \hat U^\top \neq I_{m \times m}\)
我们通过以下代码单元来说明这些性质。
UTU = U.T@U
UUT = U@U.T
print('UUT, UTU = ')
UUT, UTU
UUT, UTU =
(array([[ 1.00e+00, -1.39e-16, -1.67e-16, -2.08e-16, -2.50e-16],
[-1.39e-16, 1.00e+00, -5.55e-17, -6.94e-17, 0.00e+00],
[-1.67e-16, -5.55e-17, 1.00e+00, -9.02e-17, 0.00e+00],
[-2.08e-16, -6.94e-17, -9.02e-17, 1.00e+00, -8.33e-17],
[-2.50e-16, 0.00e+00, 0.00e+00, -8.33e-17, 1.00e+00]]),
array([[ 1.00e+00, 1.94e-16, 1.73e-16, -1.11e-16, 8.33e-17],
[ 1.94e-16, 1.00e+00, 0.00e+00, 8.33e-17, -5.55e-17],
[ 1.73e-16, 0.00e+00, 1.00e+00, 6.25e-17, -2.78e-17],
[-1.11e-16, 8.33e-17, 6.25e-17, 1.00e+00, 5.55e-17],
[ 8.33e-17, -5.55e-17, -2.78e-17, 5.55e-17, 1.00e+00]]))
UhatUhatT = Uhat@Uhat.T
UhatTUhat = Uhat.T@Uhat
print('UhatUhatT, UhatTUhat= ')
UhatUhatT, UhatTUhat
UhatUhatT, UhatTUhat=
(array([[ 0.5 , 0.07, -0.25, 0.32, -0.28],
[ 0.07, 0.4 , 0.19, 0.38, 0.23],
[-0.25, 0.19, 0.26, 0.03, 0.3 ],
[ 0.32, 0.38, 0.03, 0.49, 0.05],
[-0.28, 0.23, 0.3 , 0.05, 0.35]]),
array([[1.00e+00, 1.94e-16],
[1.94e-16, 1.00e+00]]))
注释：
上述代码展示了两种不同的SVD计算方式:
full_matrices=True计算完整的SVD分解full_matrices=False计算简化的SVD分解,只保留非零奇异值对应的部分
完整和简化的奇异值分解都能准确地分解一个 \(m \times n\) 矩阵 \(X\)
当我们在后面学习动态模态分解时，记住在这种高瘦矩阵情况下完整和简化奇异值分解的上述性质将很重要。
现在让我们来看一个矮胖矩阵的情况。
为了说明这种情况，我们将设置 \(m = 2 < 5 = n\)，并计算完整和简化的奇异值分解。
import numpy as np
X = np.random.rand(2,5)
U, S, V = np.linalg.svd(X,full_matrices=True) # 完整SVD
Uhat, Shat, Vhat = np.linalg.svd(X,full_matrices=False) # 简化SVD
print('U, S, V = ')
U, S, V
U, S, V =
(array([[ 0.53, -0.85],
[ 0.85, 0.53]]),
array([1.44, 0.66]),
array([[ 0.61, 0.31, 0.09, 0.36, 0.63],
[-0.11, 0.31, 0.07, -0.84, 0.43],
[-0.05, -0.11, 0.99, 0.02, -0.06],
[-0.7 , 0.58, 0.03, 0.39, 0.17],
[-0.35, -0.68, -0.06, 0.11, 0.63]]))
print('Uhat, Shat, Vhat = ')
Uhat, Shat, Vhat
Uhat, Shat, Vhat =
(array([[ 0.53, -0.85],
[ 0.85, 0.53]]),
array([1.44, 0.66]),
array([[ 0.61, 0.31, 0.09, 0.36, 0.63],
[-0.11, 0.31, 0.07, -0.84, 0.43]]))
让我们验证我们的简化SVD是否准确表示\(X\)
SShat=np.diag(Shat)
np.allclose(X, Uhat@SShat@Vhat)
True
5.7. 极分解#
矩阵 \(X\) 的简化奇异值分解与其极分解相关
其中
这里
\(S\) 是一个 \(m \times m\) 对称矩阵
\(Q\) 是一个 \(m \times n\) 正交矩阵
在我们的简化SVD中
\(U\) 是一个 \(m \times p\) 正交矩阵
\(\Sigma\) 是一个 \(p \times p\) 对角矩阵
\(V\) 是一个 \(n \times p\) 正交矩阵
5.8. 应用：主成分分析(PCA)#
让我们从 \(n >> m\) 的情况开始，即个体数量 \(n\) 远大于属性数量 \(m\) 的情况。
在 \(n >> m\) 的情况下，矩阵 \(X\) 是矮胖型的，这与后面要讨论的 \(m >> n\) 情况下的高瘦型相对。
我们将 \(X\) 视为一个 \(m \times n\) 的数据矩阵：
其中对于 \(j = 1, \ldots, n\)，列向量 \(X_j = \begin{bmatrix}x_{1j}\\x_{2j}\\\vdots\\x_{mj}\end{bmatrix}\) 是变量 \(\begin{bmatrix}X_1\\X_2\\\vdots\\X_m\end{bmatrix}\) 的观测值向量。
在时间序列分析中，列索引 \(j\) 代表不同的时间点，而行索引代表不同的随机变量。
在横截面分析中，列索引 \(j\) 代表不同的个体，而行索引代表它们的不同属性。
SVD 是一种将矩阵分解为有用组件的方法，类似于极分解、特征分解等。而 PCA 则是一种基于 SVD 的数据分析方法，它通过一系列统计步骤来捕捉数据中最重要的模式，帮助我们更好地理解和可视化数据。
让我们来看看 PCA 的具体步骤：
第1步：数据标准化
由于数据可能包含不同单位和尺度的变量，我们首先需要标准化数据：
首先计算 \(X\) 的每一行的平均值。
然后用这些平均值创建一个平均值矩阵：
从原始矩阵中减去平均值矩阵以创建一个均值中心化矩阵：
第2步：计算协方差矩阵
为了研究变量之间的关系，我们计算中心化数据的协方差矩阵：
第3步：分解协方差矩阵
由于矩阵\(C\)是正定的，我们可以对其进行特征值分解，找出其特征值，并按降序重新排列特征值和特征向量矩阵。
\(C\)的特征值分解可以通过分解\(B\)来得到。由于\(B\)不是方阵，我们对\(B\)进行SVD分解：
然后我们可以重新排列矩阵\(U\)和\(\Sigma\)中的列，使奇异值按降序排列。
第4步：选择主成分
我们现在可以根据想要保留的方差量来决定选择多少个奇异值（例如，保留95%的总方差）。
我们可以通过计算前\(r\)个因子包含的方差除以总方差来获得百分比：
第5步：创建主成分得分矩阵：
这个矩阵的每一列代表一个主成分，其中包含了原始数据在新坐标系下的投影。
5.9. PCA与SVD的关系#
让我们来探讨SVD与PCA之间的关系。假设我们有一个数据矩阵\(X\)，其中所有变量的样本均值都为零。
\(X\)的SVD分解可以写成:
其中矩阵\(U\)和\(V\)分别由列向量组成:
在方程(5.9)中，每个\(m \times n\)矩阵\(U_{j}V_{j}^\top \)显然是秩1的。
因此，我们有
在时间序列分析的背景下，这个分解有着重要的含义:
\( \textrm{对于每个} \ k=1, \ldots, n \)，对象 \(\lbrace V_{kj} \rbrace_{j=1}^n\) 是第\(k\)个主成分的时间序列
\(U_j = \begin{bmatrix}U_{1k}\\U_{2k}\\\ldots\\U_{mk}\end{bmatrix} \ k=1, \ldots, m\) 是变量\(X_i\)在第\(k\)个主成分上的载荷向量，其中\(i=1, \ldots, m\)
对于每个\(k=1, \ldots, p\)，\(\sigma_k\)是第\(k\)个主成分的强度，这里的强度指的是对\(X\)的整体协方差的贡献。
5.10. 基于特征值分解的PCA#
现在我们使用样本协方差矩阵的特征分解来进行PCA。
设\(X_{m \times n}\)为我们的\(m \times n\)数据矩阵。
假设所有变量的样本均值都为零。
我们可以通过减去样本均值的预处理来确保这一点。
定义样本协方差矩阵\(\Omega\)为
其特征值分解为:
这里\(P\)是\(m \times m\)特征向量矩阵,而\(\Lambda\)是特征值对角矩阵。
我们可以将\(X\)表示为:
其中
且
我们可以验证
因此数据矩阵\(X\)可以写成:
为了将前面的表示与我们之前通过SVD获得的PCA相协调，我们首先注意到\(\epsilon_j^2=\lambda_j\equiv\sigma^2_j\)。
现定义\(\tilde{\epsilon_j} = \frac{\epsilon_j}{\sqrt{\lambda_j}}\)， 这意味着\(\tilde{\epsilon}_j\tilde{\epsilon}_j^\top =1\)。
因此
这与SVD分解
是等价的,只要我们令:
\(U_j=P_j\) (变量在第j个主成分上的载荷向量)
\({V_k}^{T}=\tilde{\epsilon_k}\) (第k个主成分)
由于计算数据矩阵\(X\)的\(P\)和\(U\)有不同的算法，根据所使用的算法，我们可能会得到符号差异或特征向量顺序的不同。
我们可以通过以下方式解决关于\(U\)和\(P\)的这些歧义：
将特征值和奇异值按降序排列
在\(P\)和\(U\)中强制使对角线为正，并相应地调整\(V^\top\)中的符号
5.11. 联系与总结#
为了更好地理解前面的内容，让我们把这些公式联系起来，看看它们之间的关系。
首先，对于一个\(m \times n\)矩阵的奇异值分解（SVD）：
我们可以计算：
将这个表达式(5.12)与前面的方程(5.11)对比，我们可以发现SVD分解中的\(U\)矩阵实际上就是\(XX^\top\)的特征向量矩阵\(P\)，而\(\Sigma \Sigma^\top\)就对应着特征值矩阵\(\Lambda\)。
类似地，我们再来看：
这个结果告诉我们，SVD中的\(V\)矩阵其实就是\(X^\top X\)的特征向量矩阵。
把这些发现整理一下，我们可以得到样本协方差矩阵的特征分解：
这里的\(P\)是一个正交矩阵。而从\(X\)的SVD分解，我们也知道：
其中\(U\)同样是一个正交矩阵。
这两个表达式告诉我们\(P = U\)，因此\(X\)可以表示为：
进一步推导可得：
注意上述推导意味着
这样所有部分都完美契合。
下面我们定义一个DecomAnalysis类，用于对给定的数据矩阵X进行PCA和SVD分析。
class DecomAnalysis:
"""
用于进行PCA和SVD分析的类。
X: 数据矩阵
r_component: 最佳近似所选择的秩
"""
def __init__(self, X, r_component=None):
self.X = X
self.Ω = (X @ X.T)
self.m, self.n = X.shape
self.r = LA.matrix_rank(X)
if r_component:
self.r_component = r_component
else:
self.r_component = self.m
def pca(self):
𝜆, P = LA.eigh(self.Ω) # P的列是特征向量
ind = sorted(range(𝜆.size), key=lambda x: 𝜆[x], reverse=True)
# 按特征值排序
self.𝜆 = 𝜆[ind]
P = P[:, ind]
self.P = P @ diag_sign(P)
self.Λ = np.diag(self.𝜆)
self.explained_ratio_pca = np.cumsum(self.𝜆) / self.𝜆.sum()
# 计算N乘T的主成分矩阵
self.𝜖 = self.P.T @ self.X
P = self.P[:, :self.r_component]
𝜖 = self.𝜖[:self.r_component, :]
# 转换数据
self.X_pca = P @ 𝜖
def svd(self):
U, 𝜎, VT = LA.svd(self.X)
ind = sorted(range(𝜎.size), key=lambda x: 𝜎[x], reverse=True)
# 按特征值排序
d = min(self.m, self.n)
self.𝜎 = 𝜎[ind]
U = U[:, ind]
D = diag_sign(U)
self.U = U @ D
VT[:d, :] = D @ VT[ind, :]
self.VT = VT
self.Σ = np.zeros((self.m, self.n))
self.Σ[:d, :d] = np.diag(self.𝜎)
𝜎_sq = self.𝜎 ** 2
self.explained_ratio_svd = np.cumsum(𝜎_sq) / 𝜎_sq.sum()
# 按使用的成分数量切分矩阵
U = self.U[:, :self.r_component]
Σ = self.Σ[:self.r_component, :self.r_component]
VT = self.VT[:self.r_component, :]
# 转换数据
self.X_svd = U @ Σ @ VT
def fit(self, r_component):
# pca
P = self.P[:, :r_component]
𝜖 = self.𝜖[:r_component, :]
# 转换数据
self.X_pca = P @ 𝜖
# svd
U = self.U[:, :r_component]
Σ = self.Σ[:r_component, :r_component]
VT = self.VT[:r_component, :]
# 转换数据
self.X_svd = U @ Σ @ VT
def diag_sign(A):
"计算矩阵A对角线元素的符号"
D = np.diag(np.sign(np.diag(A)))
return D
我们还定义一个函数来打印信息，以便比较不同算法得到的分解结果。
def compare_pca_svd(da):
"""
比较PCA和SVD的结果。
"""
da.pca()
da.svd()
print('特征值和奇异值\n')
print(f'λ = {da.λ}\n')
print(f'σ^2 = {da.σ**2}\n')
print('\n')
# 载荷矩阵
fig, axs = plt.subplots(1, 2, figsize=(14, 5))
plt.suptitle('载荷')
axs[0].plot(da.P.T)
axs[0].set_title('P')
axs[0].set_xlabel('m')
axs[1].plot(da.U.T)
axs[1].set_title('U')
axs[1].set_xlabel('m')
plt.show()
# 主成分
fig, axs = plt.subplots(1, 2, figsize=(14, 5))
plt.suptitle('主成分')
axs[0].plot(da.ε.T)
axs[0].set_title('ε')
axs[0].set_xlabel('n')
axs[1].plot(da.VT[:da.r, :].T * np.sqrt(da.λ))
axs[1].set_title(r'$V^\top *\sqrt{\lambda}$')
axs[1].set_xlabel('n')
plt.show()
5.12. 练习#
练习 5.1
在普通最小二乘法(OLS)中，我们学会计算 \( \hat{\beta} = (X^\top X)^{-1} X^\top y \)，但在某些情况下，比如当我们遇到共线性或欠定系统时：即短而宽的矩阵。
在这些情况下，\( (X^\top X) \)矩阵不可逆（其行列式为零）或病态（其行列式非常接近零）。
我们可以改用所谓的伪逆，即创建一个满秩的逆矩阵近似，以此来计算 \( \hat{\beta} \)。
根据Eckart-Young定理，构建伪逆矩阵 \( X^{+} \) 并用它来计算 \( \hat{\beta} \)。
解答 练习 5.1
我们可以使用SVD来计算伪逆：
对X求逆，我们得到：
其中：
最后：
关于PCA应用于分析智力测试结构的示例，请参见本讲座 多元正态分布。
查看该讲座中描述和说明经典因子分析模型的部分。
如前所述，在后续关于 动态模态分解 的讲座中，我们将描述SVD如何提供快速计算一阶向量自回归(VARs)的降阶近似的方法。