22.2. 特征分解¶ 在 SageMaker Studio Lab 中打开 Notebook
特征值（Eigenvalues）通常是我们在学习线性代数时会遇到的最有用的概念之一。然而，对于初学者来说，很容易忽视它们的重要性。下面，我们将介绍特征分解，并试图传达它为何如此重要的感觉。
假设我们有一个矩阵 \(A\)，其元素如下：
如果我们将 \(A\) 应用于任何向量 \(\mathbf{v} = [x, y]^\top\)，我们会得到一个向量 \(\mathbf{A}\mathbf{v} = [2x, -y]^\top\)。这有一个直观的解释：在 \(x\) 方向上将向量拉伸两倍宽，然后在 \(y\) 方向上翻转它。
然而，对于某些向量，有些东西保持不变。也就是说，\([1, 0]^\top\) 被映射到 \([2, 0]^\top\)，而 \([0, 1]^\top\) 被映射到 \([0, -1]^\top\)。这些向量仍然在同一条线上，唯一的改变是矩阵分别将它们拉伸了 \(2\) 倍和 \(-1\) 倍。我们称这样的向量为特征向量（eigenvectors），它们被拉伸的因子称为特征值（eigenvalues）。
一般地，如果我们能找到一个数 \(\lambda\) 和一个向量 \(\mathbf{v}\) 使得：
我们称 \(\mathbf{v}\) 是 \(A\) 的一个特征向量，\(\lambda\) 是一个特征值。
22.2.1. 寻找特征值¶
让我们来弄清楚如何找到它们。通过从两边减去 \(\lambda \mathbf{v}\)，然后提取出向量，我们看到上面的式子等价于：
为了使 (22.2.3) 成立，我们看到 \((\mathbf{A} - \lambda \mathbf{I})\) 必须将某个方向压缩到零，因此它是不可逆的，所以它的行列式为零。因此，我们可以通过找到使得 \(\det(\mathbf{A}-\lambda \mathbf{I}) = 0\) 的 \(\lambda\) 来找到特征值。一旦我们找到特征值，我们就可以解 \(\mathbf{A}\mathbf{v} = \lambda \mathbf{v}\) 来找到相关的特征向量。
22.2.1.1. 一个例子¶
让我们用一个更具挑战性的矩阵来看这个问题：
如果我们考虑 \(\det(\mathbf{A}-\lambda \mathbf{I}) = 0\)，我们看到这等价于多项式方程 \(0 = (2-\lambda)(3-\lambda)-2 = (4-\lambda)(1-\lambda)\)。因此，两个特征值是 \(4\) 和 \(1\)。为了找到相关的向量，我们然后需要解：
我们可以分别用向量 \([1, -1]^\top\) 和 \([1, 2]^\top\) 来解这个问题。
我们可以使用内置的 numpy.linalg.eig 例程在代码中检查这一点。
%matplotlib inline
import torch
from IPython import display
from d2l import torch as d2l
torch.linalg.eig(torch.tensor([[2, 1], [2, 3]], dtype=torch.float64))
torch.return_types.linalg_eig(
eigenvalues=tensor([1.+0.j, 4.+0.j], dtype=torch.complex128),
eigenvectors=tensor([[-0.7071+0.j, -0.4472+0.j],
[ 0.7071+0.j, -0.8944+0.j]], dtype=torch.complex128))
%matplotlib inline
import numpy as np
from IPython import display
from d2l import mxnet as d2l
np.linalg.eig(np.array([[2, 1], [2, 3]]))
(array([1., 4.]),
array([[-0.70710678, -0.4472136 ],
[ 0.70710678, -0.89442719]]))
%matplotlib inline
import tensorflow as tf
from IPython import display
from d2l import tensorflow as d2l
tf.linalg.eig(tf.constant([[2, 1], [2, 3]], dtype=tf.float64))
(<tf.Tensor: shape=(2,), dtype=complex128, numpy=array([1.+0.j, 4.+0.j])>,
<tf.Tensor: shape=(2, 2), dtype=complex128, numpy=
array([[-0.70710678+0.j, -0.4472136 +0.j],
[ 0.70710678+0.j, -0.89442719+0.j]])>)
注意，numpy 将特征向量归一化为长度为一，而我们取的特征向量是任意长度的。此外，符号的选择是任意的。然而，计算出的向量与我们手动找到的具有相同特征值的向量是平行的。
22.2.2. 分解矩阵¶
让我们将前面的例子再推进一步。令
是列为矩阵 \(\mathbf{A}\) 的特征向量的矩阵。令
是对角线上有相关特征值的矩阵。那么，特征值和特征向量的定义告诉我们
矩阵 \(W\) 是可逆的，所以我们可以在两边右乘 \(W^{-1}\)，我们看到我们可以写成
在下一节中，我们将看到这样做的一些好的结果，但现在我们只需要知道，只要我们能找到一整套线性无关的特征向量（使得 \(W\) 可逆），这样的分解就存在。
22.2.3. 特征分解的运算¶
特征分解 (22.2.9) 的一个好处是，我们可以将我们通常遇到的许多运算用特征分解清晰地表示出来。作为第一个例子，考虑
这告诉我们，对于矩阵的任何正整数次幂，特征分解只需将特征值提升到相同的次幂即可。对于负次幂也可以证明同样的情况，所以如果我们想求矩阵的逆，我们只需要考虑
换句话说，只需将每个特征值取倒数。只要每个特征值都非零，这就成立，所以我们看到可逆等同于没有零特征值。
事实上，进一步的工作可以表明，如果 \(\lambda_1, \ldots, \lambda_n\) 是一个矩阵的特征值，那么该矩阵的行列式是
也就是所有特征值的乘积。这在直观上是有道理的，因为无论 \(\mathbf{W}\) 做了什么拉伸，\(W^{-1}\) 都会把它撤销，所以最终发生的唯一拉伸是通过乘以对角矩阵 \(\boldsymbol{\Sigma}\)，它将体积拉伸了对角元素的乘积。
最后，回想一下秩是矩阵线性无关列的最大数量。通过仔细研究特征分解，我们可以看到秩与 \(\mathbf{A}\) 的非零特征值的数量相同。
例子可以继续，但希望要点是清楚的：特征分解可以简化许多线性代数计算，并且是许多数值算法和我们在线性代数中所做的许多分析的基础操作。
22.2.4. 对称矩阵的特征分解¶
并不总是可能找到足够的线性无关的特征向量来使上述过程工作。例如矩阵
只有一个特征向量，即 \((1, 0)^\top\)。为了处理这样的矩阵，我们需要比我们能涵盖的更高级的技术（例如若尔当标准型或奇异值分解）。我们通常需要将注意力限制在那些我们可以保证存在一整套特征向量的矩阵上。
最常遇到的族是对称矩阵，即那些 \(\mathbf{A} = \mathbf{A}^\top\) 的矩阵。在这种情况下，我们可以取 \(W\) 为一个正交矩阵——一个其列都是长度为1且相互垂直的向量的矩阵，其中 \(\mathbf{W}^\top = \mathbf{W}^{-1}\)——并且所有的特征值都将是实数。因此，在这种特殊情况下，我们可以将 (22.2.9) 写成
22.2.5. 盖尔圆定理¶
特征值通常很难直观地推理。如果给出一个任意矩阵，在不计算它们的情况下，关于特征值是什么几乎没有什么可说的。然而，有一个定理可以使其在最大值在对角线上的情况下很容易地进行很好的近似。
令 \(\mathbf{A} = (a_{ij})\) 为任意方阵（\(n\times n\)）。我们将定义 \(r_i = \sum_{j \neq i} |a_{ij}|\)。令 \(\mathcal{D}_i\) 表示复平面中以 \(a_{ii}\) 为中心，半径为 \(r_i\) 的圆盘。那么，\(\mathbf{A}\) 的每个特征值都包含在某个 \(\mathcal{D}_i\) 中。
这可能有点难以理解，所以让我们看一个例子。考虑矩阵
我们有 \(r_1 = 0.3\)，\(r_2 = 0.6\)，\(r_3 = 0.8\) 和 \(r_4 = 0.9\)。该矩阵是对称的，所以所有特征值都是实数。这意味着我们所有的特征值都将在以下范围之一：
进行数值计算表明，特征值大约是 \(0.99\), \(2.97\), \(4.95\), \(9.08\)，都舒适地在提供的范围内。
A = torch.tensor([[1.0, 0.1, 0.1, 0.1],
[0.1, 3.0, 0.2, 0.3],
[0.1, 0.2, 5.0, 0.5],
[0.1, 0.3, 0.5, 9.0]])
v, _ = torch.linalg.eig(A)
v
tensor([0.9923+0.j, 9.0803+0.j, 4.9539+0.j, 2.9734+0.j])
A = np.array([[1.0, 0.1, 0.1, 0.1],
[0.1, 3.0, 0.2, 0.3],
[0.1, 0.2, 5.0, 0.5],
[0.1, 0.3, 0.5, 9.0]])
v, _ = np.linalg.eig(A)
v
array([9.08033648, 0.99228545, 4.95394089, 2.97343718])
A = tf.constant([[1.0, 0.1, 0.1, 0.1],
[0.1, 3.0, 0.2, 0.3],
[0.1, 0.2, 5.0, 0.5],
[0.1, 0.3, 0.5, 9.0]])
v, _ = tf.linalg.eigh(A)
v
<tf.Tensor: shape=(4,), dtype=float32, numpy=array([0.99228525, 2.9734395 , 4.953943 , 9.080336 ], dtype=float32)>
通过这种方式，特征值可以被近似，并且在对角线元素明显大于所有其他元素的情况下，近似值将相当准确。
这是一件小事，但对于像特征分解这样复杂而微妙的话题，我们能获得任何直观的理解都是好的。
22.2.6. 一个有用的应用：迭代映射的增长¶
现在我们原则上理解了特征向量是什么，让我们看看如何用它们来深入理解神经网络行为的核心问题：正确的权重初始化。
22.2.6.1. 特征向量作为长期行为¶
对深度神经网络初始化的完整数学研究超出了本文的范围，但我们可以在这里看到一个玩具版本，以了解特征值如何帮助我们理解这些模型的工作方式。我们知道，神经网络通过将线性变换层与非线性操作交错来运作。为简单起见，这里我们假设没有非线性，并且变换是单个重复的矩阵操作 \(A\)，因此我们模型的输出是
当这些模型被初始化时，\(A\) 被取为一个具有高斯分布元素的随机矩阵，所以让我们来创建一个。具体来说，我们从一个均值为零、方差为一的高斯分布的 \(5 \times 5\) 矩阵开始。
torch.manual_seed(42)
k = 5
A = torch.randn(k, k, dtype=torch.float64)
A
tensor([[ 0.2996, 0.2424, 0.2832, -0.2329, 0.6712],
[ 0.7818, -1.7903, -1.7484, 0.1735, -0.1182],
[-1.7446, -0.4695, 0.4573, 0.5177, -0.2771],
[-0.6641, 0.6551, 0.2616, -1.5265, -0.3311],
[-0.6378, 0.1072, 0.7096, 0.3009, -0.2869]], dtype=torch.float64)
np.random.seed(8675309)
k = 5
A = np.random.randn(k, k)
A
array([[ 0.58902366, 0.73311856, -1.1621888 , -0.55681601, -0.77248843],
[-0.16822143, -0.41650391, -1.37843129, 0.74925588, 0.17888446],
[ 0.69401121, -1.9780535 , -0.83381434, 0.56437344, 0.31201299],
[-0.87334496, 0.15601291, -0.38710108, -0.23920821, 0.88850104],
[ 1.29385371, -0.76774106, 0.20131613, 0.91800842, 0.38974115]])
k = 5
A = tf.random.normal((k, k), dtype=tf.float64)
A
<tf.Tensor: shape=(5, 5), dtype=float64, numpy=
array([[-0.97812483, 1.20445062, 1.11533381, -0.18427915, -1.77590971],
[-1.15848523, 1.17169049, -1.19847285, 0.55904217, -1.11710177],
[ 0.55927053, 0.66065268, -1.0880367 , -0.24903121, 0.58888241],
[ 0.4700517 , 1.03649145, -0.55023977, -0.33531388, -0.80425649],
[ 0.83811055, 0.16131817, 1.07698088, 1.71776595, -0.6111541 ]])>
22.2.6.2. 在随机数据上的行为¶
为了简化我们的玩具模型，我们将假设我们输入的数据向量 \(\mathbf{v}_{in}\) 是一个随机的五维高斯向量。让我们想想我们希望发生什么。作为背景，让我们考虑一个通用的机器学习问题，我们试图将输入数据，如一张图片，转换成一个预测，比如这张图片是猫的概率。如果重复应用 \(\mathbf{A}\) 将一个随机向量拉伸得很长，那么输入的微小变化将被放大成输出的巨大变化——输入图像的微小修改将导致截然不同的预测。这似乎不对！
另一方面，如果 \(\mathbf{A}\) 将随机向量收缩得更短，那么经过许多层之后，向量基本上会收缩为无，输出将不依赖于输入。这也显然是不对的！
我们需要在增长和衰减之间走一条窄路，以确保我们的输出随输入而变化，但变化不大！
让我们看看当我们重复地将我们的矩阵 \(\mathbf{A}\) 与一个随机输入向量相乘，并跟踪其范数时会发生什么。
# Calculate the sequence of norms after repeatedly applying `A`
v_in = torch.randn(k, 1, dtype=torch.float64)
norm_list = [torch.norm(v_in).item()]
for i in range(1, 100):
v_in = A @ v_in
norm_list.append(torch.norm(v_in).item())
d2l.plot(torch.arange(0, 100), norm_list, 'Iteration', 'Value')
# Calculate the sequence of norms after repeatedly applying `A`
v_in = np.random.randn(k, 1)
norm_list = [np.linalg.norm(v_in)]
for i in range(1, 100):
v_in = A.dot(v_in)
norm_list.append(np.linalg.norm(v_in))
d2l.plot(np.arange(0, 100), norm_list, 'Iteration', 'Value')
# Calculate the sequence of norms after repeatedly applying `A`
v_in = tf.random.normal((k, 1), dtype=tf.float64)
norm_list = [tf.norm(v_in).numpy()]
for i in range(1, 100):
v_in = tf.matmul(A, v_in)
norm_list.append(tf.norm(v_in).numpy())
d2l.plot(tf.range(0, 100), norm_list, 'Iteration', 'Value')
范数正在失控地增长！事实上，如果我们取商的列表，我们会看到一个模式。
# Compute the scaling factor of the norms
norm_ratio_list = []
for i in range(1, 100):
norm_ratio_list.append(norm_list[i]/norm_list[i - 1])
d2l.plot(torch.arange(1, 100), norm_ratio_list, 'Iteration', 'Ratio')
# Compute the scaling factor of the norms
norm_ratio_list = []
for i in range(1, 100):
norm_ratio_list.append(norm_list[i]/norm_list[i - 1])
d2l.plot(np.arange(1, 100), norm_ratio_list, 'Iteration', 'Ratio')
# Compute the scaling factor of the norms
norm_ratio_list = []
for i in range(1, 100):
norm_ratio_list.append(norm_list[i]/norm_list[i - 1])
d2l.plot(tf.range(1, 100), norm_ratio_list, 'Iteration', 'Ratio')
如果我们看一下上述计算的最后一部分，我们会看到随机向量被拉伸了 1.974459321485[...] 倍，其中结尾部分有一点变化，但拉伸因子是稳定的。
22.2.6.3. 回归到特征向量¶
我们已经看到特征向量和特征值对应于某物被拉伸的量，但那是针对特定的向量和特定的拉伸。让我们来看看它们对于 \(\mathbf{A}\) 是什么。这里有一点需要注意：事实证明，要看到它们全部，我们需要使用复数。你可以把这些看作是拉伸和旋转。通过取复数的模（实部和虚部平方和的平方根），我们可以测量那个拉伸因子。让我们也对它们进行排序。
# Compute the eigenvalues
eigs = torch.linalg.eig(A).eigenvalues.tolist()
norm_eigs = [torch.abs(torch.tensor(x)) for x in eigs]
norm_eigs.sort()
print(f'norms of eigenvalues: {norm_eigs}')
norms of eigenvalues: [tensor(0.3490), tensor(1.1296), tensor(1.1296), tensor(1.1828), tensor(2.4532)]
# Compute the eigenvalues
eigs = np.linalg.eigvals(A).tolist()
norm_eigs = [np.absolute(x) for x in eigs]
norm_eigs.sort()
print(f'norms of eigenvalues: {norm_eigs}')
norms of eigenvalues: [0.8786205280381857, 1.2757952665062624, 1.4983381517710659, 1.4983381517710659, 1.974459321485074]
# Compute the eigenvalues
eigs = tf.linalg.eigh(A)[0].numpy().tolist()
norm_eigs = [tf.abs(tf.constant(x, dtype=tf.float64)) for x in eigs]
norm_eigs.sort()
print(f'norms of eigenvalues: {norm_eigs}')
norms of eigenvalues: [<tf.Tensor: shape=(), dtype=float64, numpy=0.30955595104947614>, <tf.Tensor: shape=(), dtype=float64, numpy=1.3937271988640285>, <tf.Tensor: shape=(), dtype=float64, numpy=2.0096259821132922>, <tf.Tensor: shape=(), dtype=float64, numpy=2.2404280282645432>, <tf.Tensor: shape=(), dtype=float64, numpy=3.1559123193686784>]
22.2.6.4. 一个观察¶
我们在这里看到了一些意想不到的事情：我们之前为我们的矩阵 \(\mathbf{A}\) 应用于随机向量的长期拉伸确定的那个数，恰好（精确到十三位小数！）是 \(\mathbf{A}\) 的最大特征值。这显然不是巧合！
但是，如果我们现在从几何角度思考正在发生什么，这就开始说得通了。考虑一个随机向量。这个随机向量指向每个方向一点点，所以特别是，它至少有一点点指向与 \(\mathbf{A}\) 的最大特征值相关的特征向量相同的方向。这非常重要，以至于它被称为主特征值和主特征向量。在应用 \(\mathbf{A}\) 之后，我们的随机向量在每个可能的方向上都被拉伸，这与每个可能的特征向量相关，但它在与这个主特征向量相关的方向上被拉伸得最多。这意味着在应用 \(A\) 之后，我们的随机向量更长，并且指向一个更接近于与主特征向量对齐的方向。在多次应用矩阵之后，与主特征向量的对齐变得越来越近，直到实际上，我们的随机向量已经转换为主特征向量！事实上，这个算法是所谓的幂迭代法的基础，用于寻找矩阵的最大特征值和特征向量。详情请见，例如，（Golub and Van Loan, 1996）。
22.2.6.5. 修正归一化¶
现在，根据上面的讨论，我们得出结论，我们不希望随机向量被拉伸或挤压，我们希望随机向量在整个过程中保持大约相同的大小。为此，我们现在用这个主特征值重新缩放我们的矩阵，使得最大的特征值现在变成了一。让我们看看在这种情况下会发生什么。
# Rescale the matrix `A`
A /= norm_eigs[-1]
# Do the same experiment again
v_in = torch.randn(k, 1, dtype=torch.float64)
norm_list = [torch.norm(v_in).item()]
for i in range(1, 100):
v_in = A @ v_in
norm_list.append(torch.norm(v_in).item())
d2l.plot(torch.arange(0, 100), norm_list, 'Iteration', 'Value')
# Rescale the matrix `A`
A /= norm_eigs[-1]
# Do the same experiment again
v_in = np.random.randn(k, 1)
norm_list = [np.linalg.norm(v_in)]
for i in range(1, 100):
v_in = A.dot(v_in)
norm_list.append(np.linalg.norm(v_in))
d2l.plot(np.arange(0, 100), norm_list, 'Iteration', 'Value')
# Rescale the matrix `A`
A /= norm_eigs[-1]
# Do the same experiment again
v_in = tf.random.normal((k, 1), dtype=tf.float64)
norm_list = [tf.norm(v_in).numpy()]
for i in range(1, 100):
v_in = tf.matmul(A, v_in)
norm_list.append(tf.norm(v_in).numpy())
d2l.plot(tf.range(0, 100), norm_list, 'Iteration', 'Value')
我们也可以像之前一样绘制连续范数之间的比率，并看到它确实稳定下来了。
# Also plot the ratio
norm_ratio_list = []
for i in range(1, 100):
norm_ratio_list.append(norm_list[i]/norm_list[i-1])
d2l.plot(torch.arange(1, 100), norm_ratio_list, 'Iteration', 'Ratio')
# Also plot the ratio
norm_ratio_list = []
for i in range(1, 100):
norm_ratio_list.append(norm_list[i]/norm_list[i-1])
d2l.plot(np.arange(1, 100), norm_ratio_list, 'Iteration', 'Ratio')
# Also plot the ratio
norm_ratio_list = []
for i in range(1, 100):
norm_ratio_list.append(norm_list[i]/norm_list[i-1])
d2l.plot(tf.range(1, 100), norm_ratio_list, 'Iteration', 'Ratio')
22.2.7. 讨论¶
我们现在看到了我们所期望的！在通过主特征值对矩阵进行归一化之后，我们看到随机数据没有像以前那样爆炸，而是最终平衡到一个特定的值。如果能够从第一性原理出发做这些事情会很好，事实证明，如果我们深入研究其数学原理，我们可以看到，一个具有独立均值为零、方差为一的高斯分布元素的大型随机矩阵的最大特征值平均约为 \(\sqrt{n}\)，在我们的例子中是 \(\sqrt{5} \approx 2.2\)，这是由于一个被称为圆律的迷人事实 （Ginibre, 1965）。随机矩阵的特征值（以及一个相关的对象称为奇异值）与神经网络的正确初始化之间的关系已被证明有深刻的联系，正如在 Pennington et al. (2017) 和后续工作中讨论的那样。
22.2.8. 总结¶
特征向量是被矩阵拉伸而不改变方向的向量。
特征值是特征向量被矩阵作用后拉伸的量。
矩阵的特征分解可以使许多运算简化为对特征值的运算。
盖尔圆定理可以提供矩阵特征值的近似值。
迭代矩阵幂的行为主要取决于最大特征值的大小。这种理解在神经网络初始化的理论中有许多应用。
22.2.9. 练习¶
下面矩阵的特征值和特征向量是什么？
(22.2.21)¶\[\begin{split}\mathbf{A} = \begin{bmatrix} 2 & 1 \\ 1 & 2 \end{bmatrix}?\end{split}\]下面矩阵的特征值和特征向量是什么？与前一个例子相比，这个例子有什么奇怪之处？
(22.2.22)¶\[\begin{split}\mathbf{A} = \begin{bmatrix} 2 & 1 \\ 0 & 2 \end{bmatrix}.\end{split}\]不计算特征值，下面矩阵的最小特征值是否可能小于 \(0.5\)？注意：这个问题可以在脑中完成。
(22.2.23)¶\[\begin{split}\mathbf{A} = \begin{bmatrix} 3.0 & 0.1 & 0.3 & 1.0 \\ 0.1 & 1.0 & 0.1 & 0.2 \\ 0.3 & 0.1 & 5.0 & 0.0 \\ 1.0 & 0.2 & 0.0 & 1.8 \end{bmatrix}.\end{split}\]