矩阵求特征值的方法

求解矩阵的特征值 (Eigenvalue) 是线性代数中的一个核心问题，在物理学、工程学、计算机科学（如图论、机器学习）、经济学等众多领域都有广泛应用。特征值揭示了矩阵所代表的线性变换在特定方向上的“伸缩”特性。寻找特征值的方法多种多样，从理论定义到高效的数值迭代算法，各有侧重和适用范围。

最基础也是理论上最直接的方法是特征方程法 (Characteristic Equation Method)。根据定义，如果 λ 是方阵 A 的一个特征值，那么必然存在一个非零向量 x（称为特征向量 Eigenvector），使得 Ax = λx 成立。这个等式可以改写为 (A – λI)x = 0，其中 I 是单位矩阵。由于 x 是非零向量，这意味着矩阵 (A – λI) 是奇异的（不可逆），其行列式 (Determinant) 必须为零。

因此，求解特征值的核心就是解这个方程：

det(A – λI) = 0

这个方程称为矩阵 A 的特征方程。左边 det(A – λI) 是一个关于 λ 的多项式，称为特征多项式 (Characteristic Polynomial)。对于 n×n 的矩阵 A，这是一个关于 λ 的 n 次多项式。根据代数基本定理 (Fundamental Theorem of Algebra)，这个 n 次多项式在复数域内恰好有 n 个根（计入重数），这些根就是矩阵 A 的全部特征值。

虽然特征方程法在理论上清晰明了，对于低阶（如 2×2 或 3×3）矩阵，手算求解特征值是可行的。但对于高阶矩阵，这种方法在实践中面临巨大挑战：

1. 计算量巨大：直接计算 n 阶矩阵的行列式本身就需要大量运算（复杂度约为 O(n!) 或使用高斯消元法的 O(n³)），而得到特征多项式后，求解高次多项式的根在数值上也是一个难题，没有通用的精确解析解（对于 n≥5）。

2. 数值不稳定性：多项式的根对其系数非常敏感。在实际计算中，由于舍入误差，计算出的特征多项式系数可能与真实值有微小偏差，但这可能导致计算出的根（特征值）与真实特征值相差甚远，即存在数值不稳定性 (Numerical Instability)。

因此，在实际工程和科学计算中，特别是对于大型矩阵，很少直接使用特征方程法。取而代ل的是各种迭代算法 (Iterative Algorithms)。

幂法 (Power Iteration) 是其中一种相对简单的迭代方法，主要用于计算矩阵按模最大 (Largest Magnitude) 的特征值，即占优特征值 (Dominant Eigenvalue)，以及对应的特征向量。其基本思想是：任取一个非零的初始向量 x₀（通常选择随机向量或全1向量），然后进行迭代计算：

x_k+1 = Ax_k

在每次迭代后，通常需要对向量进行归一化 (Normalization)，以防止其分量过大或过小，例如 x_k+1 = Ax_k / ||Ax_k||。当矩阵 A 满足一定条件（例如，存在唯一的按模最大的特征值，且初始向量 x₀ 在对应特征向量方向上有非零分量），序列 x_k 的方向会逐渐收敛到占优特征值对应的特征向量 v₁ 的方向。此时，可以通过瑞利商 (Rayleigh Quotient) 来估计占优特征值 λ₁：

λ₁ ≈ (x_kᵀ A x_k) / (x_kᵀ x_k) ≈ (x_k+1ᵀ x_k) / (x_kᵀ x_k) (使用归一化向量时简化)

幂法的优点是简单易实现，尤其适用于只需要最大特征值的场景。缺点是它通常只能找到一个特征值，且收敛速度依赖于占优特征值与其他特征值模长之比，如果两者接近，收敛会很慢。

与幂法相对应的是反幂法 (Inverse Iteration)。它用于寻找距离某个位移 (Shift) σ 最近的特征值。其核心思想是对矩阵 (A – σI)⁻¹ 应用幂法。因为如果 λ 是 A 的特征值，那么 (λ – σ)⁻¹ 是 (A – σI)⁻¹ 的特征值。如果 λ 是距离 σ 最近的 A 的特征值，那么 (λ – σ)⁻¹ 就是 (A – σI)⁻¹ 按模最大的特征值。迭代过程为：

解线性方程组 (A – σI) y_k+1 = x_k

x_k+1 = y_k+1 / ||y_k+1||

当迭代收敛时，x_k 收敛到对应特征向量的方向，而特征值 λ 可以通过瑞利商或其他方式估计出来，它将非常接近 σ。特别地，当 σ = 0 时（即不进行位移），反幂法寻找的是矩阵 A 按模最小的特征值，因为此时对应的是 A⁻¹ 的占优特征值。反幂法的优势在于，如果位移 σ 选择得当（非常接近某个真实特征值），其收敛速度 (Convergence Speed) 会非常快（通常是二次或三次收敛）。缺点是每次迭代都需要求解一个线性方程组，这对于大型稠密矩阵来说计算成本较高，但对于稀疏矩阵或可以通过快速方法求解线性系统的矩阵则非常有效。

对于需要求解所有特征值的场景，尤其是对于中等规模的稠密矩阵 (Dense Matrix)，QR 算法 (QR Algorithm) 是目前最为常用和稳健 (Widely Used and Robust) 的方法。QR 算法基于矩阵的 QR 分解 (QR Decomposition)，即任何方阵 A 都可以分解为一个正交矩阵 Q (QᵀQ = I) 和一个上三角矩阵 R 的乘积：A = QR。

基本的 QR 算法迭代过程如下：

1. 初始化 A₀ = A

2. 对 k = 0, 1, 2, …

a. 对 A_k 进行 QR 分解：A_k = Q_k R_k

b. 计算下一个矩阵：A_{k+1} = R_k Q_k

可以证明，通过这种迭代方式生成的矩阵序列 {A_k} 与初始矩阵 A 是相似 (Similar) 的（A_{k+1} = R_k Q_k = Q_k⁻¹ A_k Q_k = Q_kᵀ A_k Q_k），因此它们拥有相同的特征值。在相当广泛的条件下，当 k 趋于无穷时，A_k 会收敛到一个上三角矩阵 (Upper Triangular Matrix) 或准上三角矩阵 (Quasi-Triangular Matrix，也称为实舒尔形式 Real Schur Form)，其对角线上的元素（或 2×2 块的特征值）就是原矩阵 A 的特征值。

为了加速收敛，实践中通常使用带位移的 QR 算法 (QR Algorithm with Shifts)。通过在每步迭代中引入合适的位移 σ_k，对 (A_k – σ_k I) 进行 QR 分解，可以显著加快收敛速度，特别是对于接近收敛的特征值。还有诸如双重位移 (Double Shift) 等技术用于处理复数特征值的情况，同时保持计算过程在实数域内进行。QR 算法是现代数值线性代数软件库（如 LAPACK）中计算一般矩阵特征值的标准核心算法。

对于对称矩阵 (Symmetric Matrix) 或埃尔米特矩阵 (Hermitian Matrix)，其特征值都是实数，并且存在一组正交的特征向量。针对这类矩阵，存在更高效或更简单的算法。雅可比方法 (Jacobi Method) 是一种经典的迭代方法，它通过一系列平面旋转 (Plane Rotations，也称 Givens 旋转) 来逐步对角化 (Diagonalize) 矩阵。每次旋转选择一个非对角元素 a_ij，然后构造一个旋转矩阵 J(i, j, θ)，使得经过相似变换 Jᵀ A J 后，原来 (i, j) 和 (j, i) 位置的元素变为零。不断重复这个过程，迭代地将非对角元素“扫”到对角线上。当所有非对角元素足够小时，对角线上的元素就近似为矩阵的特征值。雅可比方法相对容易理解和实现，并且对于某些类型的对称矩阵表现良好，精度也高。

对于大型稀疏矩阵 (Large Sparse Matrix)，上述的 QR 算法和雅可比方法可能因为填充（fill-in，即迭代过程中产生新的非零元素）或计算量过大而变得不适用。这时，基于Krylov 子空间 (Krylov Subspace) 的迭代方法成为主流，例如Arnoldi 算法（对于一般矩阵）和Lanczos 算法（对于对称/埃尔米特矩阵）。这些方法并不直接计算矩阵的变换，而是构造一个与特定向量相关的 Krylov 子空间（span{v, Av, A²v, …}），然后在这个低维子空间上寻找原矩阵特征值的近似。这些方法特别适合于寻找矩阵的部分特征值（例如，模最大或最小的几个特征值），并且能很好地利用矩阵的稀疏性，因为它们主要依赖于矩阵向量乘积 (Matrix-Vector Product) 运算。

总结来说，矩阵特征值的计算方法呈现多样性：

特征方程法是理论基础，但实践中数值性能差。

幂法和反幂法简单，适用于计算特定（最大/最小/最接近给定值）的单个或少数几个特征值。

QR 算法是计算稠密矩阵所有特征值的工业标准，稳健且高效，尤其结合位移策略后。

雅可比方法适用于对称矩阵的对角化，精度高。

Arnoldi/Lanczos 算法则是在处理大型稀疏矩阵时求解部分特征值的有力武器。

选择哪种方法取决于矩阵的规模 (Size)、结构 (Structure，如对称性、稀疏性)、以及需要求解的特征值范围 (all eigenvalues or specific ones)。理解这些方法的原理、优缺点和适用场景，对于在科学研究和工程实践中有效解决相关问题至关重要。