# 用Python验证特征值重数定理：从矩阵构建到Sympy实战 你是否曾在学习线性代数时，对特征值的“代数重数”和“几何重数”这两个概念感到困惑？教科书上的证明虽然严谨，但总感觉隔着一层抽象的薄纱，难以触及本质。作为一名数据科学或机器学习从业者，我们习惯于与数据和代码打交道，通过“运行一下”来加深理解。今天，我们就换个思路，不满足于纸面推导，而是直接动手，用Python和Sympy库，在Jupyter Notebook的交互环境中，亲手构建矩阵、计算特征值、求解特征空间，从而直观地验证那个重要的定理：**特征值的几何重数不超过其代数重数**。这个过程不仅能帮你彻底吃透理论，更能让你掌握一套将抽象数学定理转化为可验证、可探索代码的实战方法，这在处理实际模型中的奇异值分解、主成分分析稳定性等问题时，将是一把利器。 ## 1. 重数定理：从抽象定义到可计算概念 在深入代码之前，我们必须清晰地界定战场。特征值问题，简单说，就是寻找一个标量λ和一个非零向量v，使得对于给定矩阵A，满足方程 **A v = λ v**。这个λ就是特征值，v就是对应的特征向量。 * **代数重数**：听起来很“代数”。它指的是特征多项式 `det(A - λI) = 0` 中，根 λ 的重数。比如，如果 `(λ-3)^2 * (λ+1) = 0`，那么特征值3的代数重数就是2，特征值-1的代数重数就是1。它描述的是这个特征值在多项式意义上的“出现次数”。 * **几何重数**：听起来很“几何”。它指的是对应某个特征值λ的所有特征向量所张成的空间（即特征空间）的维数。换句话说，就是矩阵 `(A - λI)` 的零空间的维数，即 `nullity(A - λI)`。它描述的是有多少个线性无关的特征向量“指向”这个特征值。 > 注意：一个常见的误解是认为代数重数直接等于线性无关特征向量的个数。实际上，几何重数才是这个个数，而代数重数只是一个多项式意义上的计数，它总是大于等于几何重数。 为什么这个不等式 `几何重数 ≤ 代数重数` 如此重要？在机器学习中，例如PCA（主成分分析），我们依赖协方差矩阵的特征分解。如果一个特征值的几何重数小于代数重数（即矩阵是“亏损”的），意味着无法找到足够多的线性无关特征向量来构成一组完备基，这可能会在数值计算中引发稳定性问题，或者在试图对角化矩阵时遭遇失败。理解并能够验证这一定理，是诊断此类潜在问题的第一步。 为了后续用Sympy验证，我们先明确计算路径： 1. **计算代数重数**：求解特征多项式，因式分解，统计特定特征值的指数。 2. **计算几何重数**：对特定特征值λ，计算矩阵 `(A - λI)` 的秩，然后利用公式：`几何重数 = n - rank(A - λI)`，其中n是矩阵的阶数。 ## 2. 构建Sympy实验环境与示例矩阵 理论清晰后，我们进入实战环节。我们将使用Sympy，这是一个用于符号数学的Python库。与NumPy不同，Sympy能进行精确的符号计算，避免浮点数误差对重数判断的干扰，非常适合教学和理论验证。 首先，确保你的环境已安装Sympy。在Jupyter Notebook中，我们开始初始化并构建一个精心设计的示例矩阵。这个矩阵需要能同时展示“几何重数等于代数重数”和“几何重数小于代数重数”两种情况。 ```python import sympy as sp # 初始化符号计算环境 sp.init_printing(use_unicode=True) # 构建一个4x4的示例矩阵A # 我们设计它有一个代数重数为2，几何重数为1的特征值（例如2） # 以及另一个代数重数为2，几何重数也为2的特征值（例如3） # 这样可以在一个矩阵中对比两种情形 A = sp.Matrix([ [2, 1, 0, 0], [0, 2, 0, 0], [0, 0, 3, 0], [0, 0, 0, 3] ]) print("示例矩阵 A:") A ``` 运行上述代码，你会看到一个分块上三角矩阵。左上角的2x2块是一个若尔当块（Jordan Block），其特征值2的代数重数为2，但几何重数仅为1。右下角的2x2块是一个对角矩阵，其特征值3的代数重数和几何重数均为2。这个矩阵是我们验证定理的完美沙盒。 > 提示：在Jupyter中，`sp.init_printing()` 能让矩阵和公式以更美观的LaTeX格式渲染，极大地提升可读性。 接下来，我们计算它的特征多项式，这是求代数重数的起点。 ```python # 定义符号变量 lambda λ = sp.symbols('λ') # 计算特征多项式：det(A - λ * I) I = sp.eye(4) # 4阶单位矩阵 char_poly = (A - λ * I).det() char_poly_simplified = sp.factor(char_poly) print("特征多项式 det(A - λI) = ") char_poly_simplified ``` Sympy会输出因式分解后的结果：`(λ - 3)**2 * (λ - 2)**2`。这清晰地告诉我们： * 特征值 `λ1 = 2` 的代数重数为 **2**。 * 特征值 `λ2 = 3` 的代数重数为 **2**。 ## 3. 手动计算与验证：几何重数的求解 现在，我们进入核心环节：分别计算两个特征值的几何重数。回忆一下，几何重数等于 `n - rank(A - λI)`。 ### 3.1 计算特征值 λ=2 的几何重数 ```python # 对于特征值 λ = 2 λ_val = 2 M = A - λ_val * I # 计算 A - 2I print(f"矩阵 A - {λ_val}I:") display(M) # 计算矩阵的秩 rank_M = M.rank() print(f"秩 rank(A - {λ_val}I) = {rank_M}") # 计算几何重数 geo_mult = 4 - rank_M # 矩阵是4x4的 print(f"特征值 {λ_val} 的几何重数 = 4 - {rank_M} = {geo_mult}") # 我们也可以直接求出特征空间（零空间）的基，其向量个数即几何重数 nullspace = M.nullspace() print(f"\n特征空间（零空间）的基向量（即特征向量）：") for i, vec in enumerate(nullspace): print(f"v{i+1}:", vec.T) # .T 是为了转置显示为行向量，便于阅读 ``` 观察输出。`A - 2I` 的秩为3，因此几何重数为 `4-3=1`。同时，`nullspace()` 方法只返回了一个基向量，例如 `[1, 0, 0, 0]`。这验证了我们的设计：对于特征值2，**代数重数(2) > 几何重数(1)**。 ### 3.2 计算特征值 λ=3 的几何重数 ```python # 对于特征值 λ = 3 λ_val = 3 M = A - λ_val * I # 计算 A - 3I print(f"矩阵 A - {λ_val}I:") display(M) rank_M = M.rank() print(f"秩 rank(A - {λ_val}I) = {rank_M}") geo_mult = 4 - rank_M print(f"特征值 {λ_val} 的几何重数 = 4 - {rank_M} = {geo_mult}") nullspace = M.nullspace() print(f"\n特征空间（零空间）的基向量：") for i, vec in enumerate(nullspace): print(f"v{i+1}:", vec.T) ``` 这次，`A - 3I` 的秩为2，因此几何重数为 `4-2=2`。`nullspace()` 返回了两个线性无关的向量，例如 `[0, 0, 1, 0]` 和 `[0, 0, 0, 1]`。这展示了另一种情况：对于特征值3，**代数重数(2) = 几何重数(2)**。 我们将两次计算的结果汇总对比： | 特征值 (λ) | 代数重数 (来自特征多项式) | 矩阵 (A - λI) 的秩 | 几何重数 (n - rank) | 不等式关系 | | :--- | :---: | :---: | :---: | :--- | | 2 | 2 | 3 | 1 | 代数重数 > 几何重数 | | 3 | 2 | 2 | 2 | 代数重数 = 几何重数 | 这个表格直观地证实了定理：在所有情况下，几何重数都未超过代数重数。 ## 4. 自动化校验与常见数值陷阱剖析 手动验证一两个矩阵固然有启发性，但作为一个严谨的实践者，我们更希望有一个自动化工具来批量验证或检查随机矩阵。同时，我们必须警惕从理论世界踏入数值计算世界时可能遇到的“陷阱”。 ### 4.1 编写通用验证函数 我们可以编写一个函数，输入一个Sympy矩阵，自动计算其每个不同特征值的代数和几何重数，并进行比较。 ```python def verify_eigenvalue_multiplicity_theorem(matrix): """ 验证给定方阵的特征值重数定理。 返回一个字典，包含每个特征值的代数和几何重数，以及验证结果。 """ n = matrix.shape[0] if n != matrix.shape[1]: raise ValueError("输入必须是方阵。") # 计算特征多项式并因式分解 λ = sp.symbols('λ') char_poly = sp.factor((matrix - λ * sp.eye(n)).det()) print(f"特征多项式: {char_poly}") # 从因式分解中提取特征值及其代数重数 # 注意：sympy返回的因子形式可能是 (λ - a)**m from sympy import Pow algebraic_mults = {} if isinstance(char_poly, sp.Mul): factors = char_poly.args else: factors = (char_poly,) for factor in factors: if isinstance(factor, Pow): base, exp = factor.args if base.has(λ): # 形式为 (λ - a)**m eigenvalue = sp.solve(base, λ)[0] algebraic_mults[eigenvalue] = exp elif factor.has(λ) and factor.is_Poly: # 处理不可约的一次因子 (λ - a) roots = sp.solve(factor, λ) for root in roots: algebraic_mults[root] = algebraic_mults.get(root, 0) + 1 results = {} for eig_val, alg_mult in algebraic_mults.items(): # 计算几何重数 M = matrix - eig_val * sp.eye(n) geo_mult = n - M.rank() theorem_holds = (geo_mult <= alg_mult) results[eigenvalue] = { '代数重数': alg_mult, '几何重数': geo_mult, '定理成立？': theorem_holds } print(f"\n特征值: {eig_val}") print(f" 代数重数: {alg_mult}") print(f" 几何重数: {geo_mult}") print(f" 验证 (几何≤代数): {theorem_holds}") all_hold = all(info['定理成立？'] for info in results.values()) print(f"\n所有特征值均满足定理: {all_hold}") return results, all_hold # 用我们之前的矩阵A测试 print("=== 验证矩阵A ===") results_A, hold_A = verify_eigenvalue_multiplicity_theorem(A) ``` 这个函数提供了自动化的验证流程。对于更复杂的矩阵，它能高效地完成检查。 ### 4.2 数值计算中的陷阱与应对策略 当我们离开Sympy的精确符号计算，转向像NumPy/SciPy这样的数值计算库处理真实世界数据时，问题会变得微妙。 * **陷阱一：特征值的重复与近似**。数值算法（如`numpy.linalg.eig`）计算出的特征值可能存在微小误差。理论上重数为2的特征值，在数值上可能显示为两个非常接近但不完全相等的数（如 3.0000001 和 2.9999999）。此时，直接判断它们的重数会失败。 * **陷阱二：秩的数值判定**。计算 `几何重数 = n - rank(A - λI)` 时，数值库计算的秩是基于矩阵的奇异值分解（SVD）和一个阈值（如 `tol=1e-12`）来判断的。对于病态矩阵或特征值有微小误差时，`A - λI` 可能不是真正的奇异矩阵，导致算出的秩错误，从而影响几何重数。 **应对策略：** 1. **聚类近似特征值**：对于计算出的特征值列表，先进行聚类。将距离小于某个容差（例如 `tol=1e-10`）的特征值视为同一个。 ```python import numpy as np from scipy.linalg import eig # 假设 eigvals 是从 numpy.linalg.eig 得到的特征值数组 eigvals = np.array([3.0, 3.000000001, 2.0, 1.999999999]) tol = 1e-9 clustered = [] for val in eigvals: if not clustered: clustered.append([val]) else: # 检查是否与已有聚类中心足够接近 assigned = False for cluster in clustered: if abs(val - np.mean(cluster)) < tol: cluster.append(val) assigned = True break if not assigned: clustered.append([val]) # 每个 cluster 代表一个数值近似下的特征值 print("聚类后的特征值组:", clustered) ``` 2. **使用奇异值分解（SVD）判断几何重数**：对于疑似特征值λ，计算矩阵 `B = A - λI` 的奇异值。几何重数（零空间维数）等于奇异值中小于某个阈值的个数。这比直接求秩更稳健。 ```python import numpy as np from scipy.linalg import svd def numerical_geometric_multiplicity(A, eigval_approx, tol=1e-12): B = A - eigval_approx * np.eye(A.shape[0]) s = svd(B, compute_uv=False) # 只计算奇异值 # 统计小于阈值的奇异值个数，即零空间维数 geo_mult = np.sum(s < tol) return geo_mult ``` 3. **理解工具限制**：对于高度病态或接近亏损的矩阵，数值方法可能无法给出可靠的重数信息。此时，可能需要结合问题背景知识，或使用更高精度的计算库（如`mpmath`）。 在Jupyter Notebook中探索这些陷阱非常有益。你可以尝试用NumPy生成一个接近我们示例矩阵A的数值矩阵，然后观察`eig`函数返回的特征值和特征向量，并与Sympy的精确结果对比，亲身感受数值误差带来的挑战。这种经验能让你在未来处理实际数据，例如近乎共线的协方差矩阵时，保持必要的警惕。 通过这一整套从理论理解、符号计算验证到数值陷阱剖析的流程，我们不仅牢牢掌握了特征值重数定理，更获得了一套可迁移的方法论：如何用计算工具去探索和验证数学理论。下次当你怀疑一个优化问题的Hessian矩阵是否可对角化，或者分析动力系统稳定性时，你会知道，除了纸笔，你还有Python这个强大的实验伙伴。记住，在数据科学里，能跑通的代码，往往是最有说服力的证明。