用Python实战验证Hamilton-Cayley定理：当矩阵遇见自己的命运方程式

第一次听说Hamilton-Cayley定理时，我的反应和大多数理工科学生一样——这个宣称"每个矩阵都是自己特征多项式根"的命题，听起来就像在说"每个人最终都会成为自己最讨厌的样子"。直到我用Python把它计算出来，看着屏幕上那个逼近零矩阵的结果，才真正理解这个线性代数中最优雅的"自我指涉"命题。本文将带你用NumPy进行数值验证，用SymPy进行符号推导，最后还会探讨这个定理如何实际帮我们简化矩阵运算。

1. 准备工作：理解定理与搭建环境

Hamilton-Cayley定理的核心内容是：对于任意n阶方阵A，如果它的特征多项式是φ(λ)=det(λI-A)，那么将矩阵A代入这个多项式得到的φ(A)结果必然是零矩阵。换句话说，矩阵满足自己的特征方程。

需要准备的Python环境：

# 安装必要库
pip install numpy sympy matplotlib

关键概念快速回顾：

• 特征多项式：det(λI-A)得到的关于λ的n次多项式
• 矩阵代入多项式：将A⁰视为单位矩阵，A¹就是A本身，A²是矩阵乘法，以此类推
• 零矩阵：所有元素都为0的矩阵

2. 数值验证：用NumPy让定理"显形"

让我们先从一个具体的3×3矩阵开始，用NumPy进行实际计算验证：

import numpy as np

# 定义一个随机矩阵（可替换为你感兴趣的任意矩阵）
A = np.array([[2, -1, 0],
              [-1, 2, -1],
              [0, -1, 2]])

# 计算特征多项式系数
coeffs = np.poly(A)  # 返回从高次到低次的系数

# 验证Hamilton-Cayley定理
result = np.zeros_like(A)  # 准备累加器
n = A.shape[0]
for i, c in enumerate(coeffs):
    power = n - i  # 当前项的次数
    result += c * np.linalg.matrix_power(A, power)

print("验证结果（应接近零矩阵）：\n", result)

典型输出示例：

验证结果（应接近零矩阵）：
 [[ 3.55271368e-15  0.00000000e+00  0.00000000e+00]
 [ 0.00000000e+00  3.55271368e-15  0.00000000e+00]
 [ 0.00000000e+00  0.00000000e+00  3.55271368e-15]]

注意：由于浮点数计算精度限制，结果中出现的极小值（如1e-15量级）可视为计算误差，实际理论值应为严格的零矩阵。

3. 符号推导：用SymPy展示精确证明

为了消除数值计算中的精度问题，我们可以使用SymPy进行符号计算，这相当于在数学上做严格推导：

from sympy import symbols, Matrix, eye, det, simplify

# 定义符号矩阵
a, b, c, d = symbols('a b c d')
A = Matrix([[a, b], [c, d]])

# 计算特征多项式
lambda_ = symbols('lambda')
char_poly = det(lambda_ * eye(2) - A)
expanded_poly = char_poly.expand()

# 提取多项式系数
poly_terms = expanded_poly.as_ordered_terms()
coeff_dict = expanded_poly.as_coefficients_dict()

# 代入矩阵计算
I = eye(2)  # 单位矩阵
result = coeff_dict[lambda_**2] * A**2 + coeff_dict[lambda_**1] * A + coeff_dict[lambda_**0] * I

print("符号验证结果：\n", simplify(result))

这段代码将对任意2×2符号矩阵进行验证，输出结果确实是零矩阵。你可以修改矩阵维度，虽然计算量会增加，但结论依然成立。

4. 实际应用：矩阵高次幂的简化计算

Hamilton-Cayley定理不仅是个漂亮的数学结论，它还能帮我们实际简化矩阵运算。例如计算矩阵的高次幂时：

传统方法的问题：

• 直接连乘：时间复杂度O(n³logk)（对k次幂）
• 对角化方法：需要求特征值和特征向量，可能遇到数值不稳定

利用Hamilton-Cayley的简化方法：

1. 根据特征多项式，最高次项可以表示为低次项的线性组合
2. 这意味着任何A^k都可以表示为I, A, A²,...,A^(n-1)的线性组合
3. 大大降低了计算复杂度

def matrix_power_hc(A, k):
    """利用Hamilton-Cayley定理计算矩阵幂"""
    n = A.shape[0]
    coeffs = np.poly(A)[:-1]  # 去掉最高次项系数1
    
    # 初始化幂矩阵列表：I, A, A²,...,A^(n-1)
    powers = [np.eye(n)]
    for i in range(1, n):
        powers.append(np.dot(powers[-1], A))
    
    # 计算递推系数
    def get_coeff(m):
        if m < n:
            return (m == k) * np.eye(n)
        else:
            return -sum(coeff * get_coeff(m - n + i) 
                       for i, coeff in enumerate(coeffs))
    
    return sum(np.kron(get_coeff(k - i), power) 
               for i, power in enumerate(powers[:min(k+1, n)]))

提示：这个方法特别适合需要重复计算同一矩阵不同幂次的情况，因为可以将I到A^(n-1)预先计算存储。

5. 可视化验证：特征值与矩阵代入的关系

为了更直观理解为什么Hamilton-Cayley定理成立，我们可以可视化特征值与矩阵代入的关系：

import matplotlib.pyplot as plt

# 计算特征值
eigenvalues = np.linalg.eigvals(A)

# 生成特征多项式函数
def char_poly_func(x):
    return np.prod(x - eigenvalues)

# 矩阵代入后的范数变化
matrix_values = []
x_range = np.linspace(min(eigenvalues)-1, max(eigenvalues)+1, 100)
for x in x_range:
    substituted = np.sum([coeff * np.linalg.matrix_power(A, i) 
                         for i, coeff in enumerate(coeffs[::-1])], axis=0)
    matrix_values.append(np.linalg.norm(substituted))

# 绘制对比图
plt.figure(figsize=(10, 5))
plt.plot(x_range, [char_poly_func(x) for x in x_range], label='标量特征多项式')
plt.plot(x_range, matrix_values, label='矩阵代入的范数')
plt.scatter(eigenvalues, [0]*len(eigenvalues), c='r', label='特征值位置')
plt.axhline(0, color='gray', linestyle='--')
plt.legend()
plt.title("Hamilton-Cayley定理可视化验证")
plt.xlabel("λ")
plt.ylabel("值/范数")
plt.show()

这个可视化会显示：

• 标量特征多项式在特征值处过零点
• 矩阵代入的范数在所有位置都接近零（在特征值处可能更明显）
• 直观验证φ(A)≈O的结论

6. 常见问题与调试技巧

在实际验证过程中，可能会遇到以下典型问题：

问题1：数值验证结果不接近零矩阵

• 检查矩阵是否确实可对角化（尝试 np.linalg.eig 看特征向量矩阵是否可逆）
• 尝试减小矩阵元素的大小（大数值会放大浮点误差）
• 考虑使用更高精度的数据类型（如 np.float128 ）

问题2：符号计算速度太慢

• 对于大于3×3的矩阵，符号计算会变得非常耗时
• 可以尝试：

  from sympy import simplify
  simplify(result, ratio=1.7)  # 调整简化强度
  

问题3：特征多项式系数提取错误

• NumPy的 poly 函数返回的系数顺序是从高次到低次
• 确保在矩阵代入时正确对应：

  # 正确顺序：
  # coeffs[0]*A^n + coeffs[1]*A^(n-1) + ... + coeffs[n]*I
  

7. 扩展应用：矩阵函数的计算

Hamilton-Cayley定理的一个重要应用是计算矩阵函数，如矩阵指数e^A：

基本原理：

1. 任何解析函数f(x)都可以表示为特征多项式的商加上余项
2. 余项次数小于n，因此f(A)的计算简化为低次幂的线性组合

def matrix_exp_hc(A, terms=10):
    """利用Hamilton-Cayley定理计算矩阵指数"""
    n = A.shape[0]
    coeffs = np.poly(A)[:-1]  # 特征多项式系数
    
    # 预先计算I到A^(n-1)
    powers = [np.eye(n)]
    for i in range(1, n):
        powers.append(powers[-1] @ A)
    
    # 计算泰勒级数余项系数
    def residual_coeff(k):
        fact = 1
        for i in range(1, k+1):
            fact *= i
        return 1/fact
    
    # 组合结果
    return sum(residual_coeff(i) * powers[i] 
               for i in range(min(terms, n)))

这个方法相比直接泰勒展开，在保持精度的同时大大减少了计算量。类似的方法还可以应用于计算矩阵三角函数、矩阵对数等。