Hessenberg-QR算法计算非对称矩阵特征值
字数 2707 2025-12-09 11:44:52

Hessenberg-QR算法计算非对称矩阵特征值

题目描述
计算一个实方阵 \(A \in \mathbb{R}^{n \times n}\) 的全部特征值。这是一个非对称矩阵的特征值问题。Hessenberg-QR算法是解决这类问题的经典数值方法，其核心思想是先通过相似变换将 \(A\) 化为上Hessenberg矩阵（次对角线以下为零的矩阵），然后对Hessenberg矩阵应用带位移的QR迭代，高效地收敛到特征值。

解题过程

第一步：理解目标与基本思路
矩阵 \(A\) 的特征值是满足 \(A v = \lambda v\) 的标量 \(\lambda\)。直接对 \(A\) 做QR迭代是低效的，因为每次QR分解的运算量是 \(O(n^3)\)。Hessenberg-QR算法的优化分为两步：

  1. 预处理阶段：将 \(A\) 化为上Hessenberg矩阵 \(H\)（即当 \(i > j+1\) 时 \(h_{ij} = 0\)），这个过程通过相似变换完成，不改变特征值。
  2. 迭代阶段：对 \(H\) 进行带位移的QR迭代，此时每次迭代的运算量降为 \(O(n^2)\)，且可通过位移策略加速收敛。

相似变换保持特征值不变：若 \(A = P H P^{-1}\)，则 \(A\) 与 \(H\) 的特征值相同。这里 \(P\) 是正交矩阵（保证数值稳定性）。

第二步：化矩阵为上Hessenberg形式
给定矩阵 \(A \in \mathbb{R}^{n \times n}\)，我们希望找到正交矩阵 \(P\)（一系列Householder反射器的乘积），使得 \(H = P A P^T\) 是上Hessenberg矩阵。

具体步骤（以 \(n=5\) 为例，逐列消元）：

  1. 对第1列，考虑元素 \(a_{21}, a_{31}, \dots, a_{51}\)，构造Householder反射器 \(P_1\) 使得 \(a_{31}\) 到 \(a_{51}\) 变为零，但保持 \(a_{21}\) 不变（以保留次对角线）。将 \(P_1\) 同时左乘和右乘到 \(A\) 上：

\[ A^{(1)} = P_1 A P_1^T \]

右乘 \(P_1^T\) 是为了保持相似性，但会影响已引入的零元吗？实际上右乘只会影响行，不破坏列中已引入的零元结构。
2. 对第2列，类似地对元素 \(a_{32}, a_{42}, a_{52}\) 构造 \(P_2\)，将 \(a_{42}, a_{52}\) 化为零。注意 \(P_2\) 的前两行和前两列是单位矩阵，以保护第1列已引入的零。然后更新：

\[ A^{(2)} = P_2 A^{(1)} P_2^T \]

  1. 重复直到第 \(n-2\) 列。最终得到上Hessenberg矩阵 \(H = A^{(n-2)}\)。

运算量：约 \(\frac{10}{3} n^3\) 浮点运算，但只需做一次。

第三步：理解QR迭代对上Hessenberg矩阵的加速
标准QR迭代：对矩阵 \(B_k\) 做QR分解 \(B_k = Q_k R_k\)，然后反向乘积得 \(B_{k+1} = R_k Q_k\)，如此重复，\(B_k\) 会收敛到上三角矩阵（特征值在对角线）。
如果 \(B_0\) 是上Hessenberg矩阵，可以证明：

  • 在QR分解中，\(Q_k\) 也是上Hessenberg矩阵（可利用Givens旋转高效计算）。
  • \(B_{k+1}\) 保持上Hessenberg形。
    因此每次迭代只需对次对角线非零的部分操作，运算量从 \(O(n^3)\) 降为 \(O(n^2)\)。

第四步：引入位移策略加速收敛
QR迭代的收敛速度依赖于相邻特征值的比值。对实矩阵，特征值可能是复数，但复特征值以共轭对出现。为了加速，采用双重隐式位移QR迭代（Francis double shift）：

  1. 从当前上Hessenberg矩阵 \(H\) 取出右下角 \(2\times2\) 子矩阵，计算它的两个特征值 \(\mu_1, \mu_2\)（可能是复数）。
  2. 计算位移量：实际使用 \(\mu_1, \mu_2\) 作为两次位移，但为了避免复数运算，在一次迭代中隐式完成。
  3. 构造一个正交变换矩阵 \(Q\)（通过计算 \(H\) 的第一列与单位矩阵的线性组合），使得 \(Q\) 能将 \(H\) 变为新的上Hessenberg矩阵，且这个变换等价于做了两次带位移的QR迭代，但完全在实数运算下完成。

隐式位移的关键：通过一个精心设计的第一列变换，利用隐式Q定理保证最终矩阵与显式做两次位移QR迭代得到的矩阵相同（在舍入误差范围内）。

第五步：迭代直至收敛
每次迭代后，检查次对角线元素 \(h_{i+1,i}\) 是否可视为零（若 \(|h_{i+1,i}| \le \epsilon (|h_{ii}|+|h_{i+1,i+1}|)\)）。如果可视为零，则将矩阵分裂为两个更小的Hessenberg子矩阵，分别处理。

  • 如果 \(1 \times 1\) 块分离，得到一个实特征值。
  • 如果 \(2 \times 2\) 块分离，计算它的两个特征值（可能是实或共轭复对）。
    继续迭代直至所有特征值被提取。

第六步：整体算法流程

  1. 输入实方阵 \(A\)。
  2. 用Householder变换化 \(A\) 为上Hessenberg矩阵 \(H\)。
  3. 重复：
    • 如果 \(H\) 的次对角线元足够小，分割 \(H\)。
    • 对当前Hessenberg块，用双重隐式位移做一次QR迭代，更新 \(H\)。
  4. 直到所有特征值被求出（对角块为 \(1\times1\) 或 \(2\times2\)）。

输出：所有特征值（实数或共轭复数对）。

总结
Hessenberg-QR算法将非对称矩阵特征值计算分解为“预处理+迭代”两阶段，通过上Hessenberg化降低迭代成本，再通过隐式位移QR迭代高效、稳定地逼近特征值。这是EISPACK、LAPACK等库中 dgeev 等例程的核心算法。

Hessenberg-QR算法计算非对称矩阵特征值 题目描述 计算一个实方阵 \( A \in \mathbb{R}^{n \times n} \) 的全部特征值。这是一个非对称矩阵的特征值问题。Hessenberg-QR算法是解决这类问题的经典数值方法，其核心思想是先通过相似变换将 \( A \) 化为上Hessenberg矩阵（次对角线以下为零的矩阵），然后对Hessenberg矩阵应用带位移的QR迭代，高效地收敛到特征值。 解题过程 第一步：理解目标与基本思路 矩阵 \( A \) 的特征值是满足 \( A v = \lambda v \) 的标量 \( \lambda \)。直接对 \( A \) 做QR迭代是低效的，因为每次QR分解的运算量是 \( O(n^3) \)。Hessenberg-QR算法的优化分为两步： 预处理阶段 ：将 \( A \) 化为上Hessenberg矩阵 \( H \)（即当 \( i > j+1 \) 时 \( h_ {ij} = 0 \)），这个过程通过相似变换完成，不改变特征值。 迭代阶段 ：对 \( H \) 进行带位移的QR迭代，此时每次迭代的运算量降为 \( O(n^2) \)，且可通过位移策略加速收敛。 相似变换保持特征值不变 ：若 \( A = P H P^{-1} \)，则 \( A \) 与 \( H \) 的特征值相同。这里 \( P \) 是正交矩阵（保证数值稳定性）。 第二步：化矩阵为上Hessenberg形式 给定矩阵 \( A \in \mathbb{R}^{n \times n} \)，我们希望找到正交矩阵 \( P \)（一系列Householder反射器的乘积），使得 \( H = P A P^T \) 是上Hessenberg矩阵。 具体步骤 （以 \( n=5 \) 为例，逐列消元）： 对第1列，考虑元素 \( a_ {21}, a_ {31}, \dots, a_ {51} \)，构造Householder反射器 \( P_ 1 \) 使得 \( a_ {31} \) 到 \( a_ {51} \) 变为零，但保持 \( a_ {21} \) 不变（以保留次对角线）。将 \( P_ 1 \) 同时左乘和右乘到 \( A \) 上： \[ A^{(1)} = P_ 1 A P_ 1^T \] 右乘 \( P_ 1^T \) 是为了保持相似性，但会影响已引入的零元吗？实际上右乘只会影响行，不破坏列中已引入的零元结构。 对第2列，类似地对元素 \( a_ {32}, a_ {42}, a_ {52} \) 构造 \( P_ 2 \)，将 \( a_ {42}, a_ {52} \) 化为零。注意 \( P_ 2 \) 的前两行和前两列是单位矩阵，以保护第1列已引入的零。然后更新： \[ A^{(2)} = P_ 2 A^{(1)} P_ 2^T \] 重复直到第 \( n-2 \) 列。最终得到上Hessenberg矩阵 \( H = A^{(n-2)} \)。 运算量 ：约 \( \frac{10}{3} n^3 \) 浮点运算，但只需做一次。 第三步：理解QR迭代对上Hessenberg矩阵的加速 标准QR迭代：对矩阵 \( B_ k \) 做QR分解 \( B_ k = Q_ k R_ k \)，然后反向乘积得 \( B_ {k+1} = R_ k Q_ k \)，如此重复，\( B_ k \) 会收敛到上三角矩阵（特征值在对角线）。 如果 \( B_ 0 \) 是上Hessenberg矩阵，可以证明： 在QR分解中，\( Q_ k \) 也是上Hessenberg矩阵（可利用Givens旋转高效计算）。 \( B_ {k+1} \) 保持上Hessenberg形。 因此每次迭代只需对次对角线非零的部分操作，运算量从 \( O(n^3) \) 降为 \( O(n^2) \)。 第四步：引入位移策略加速收敛 QR迭代的收敛速度依赖于相邻特征值的比值。对实矩阵，特征值可能是复数，但复特征值以共轭对出现。为了加速，采用 双重隐式位移QR迭代 （Francis double shift）： 从当前上Hessenberg矩阵 \( H \) 取出右下角 \( 2\times2 \) 子矩阵，计算它的两个特征值 \( \mu_ 1, \mu_ 2 \)（可能是复数）。 计算位移量：实际使用 \( \mu_ 1, \mu_ 2 \) 作为两次位移，但为了避免复数运算，在一次迭代中隐式完成。 构造一个正交变换矩阵 \( Q \)（通过计算 \( H \) 的第一列与单位矩阵的线性组合），使得 \( Q \) 能将 \( H \) 变为新的上Hessenberg矩阵，且这个变换等价于做了两次带位移的QR迭代，但完全在实数运算下完成。 隐式位移的关键 ：通过一个精心设计的第一列变换，利用隐式Q定理保证最终矩阵与显式做两次位移QR迭代得到的矩阵相同（在舍入误差范围内）。 第五步：迭代直至收敛 每次迭代后，检查次对角线元素 \( h_ {i+1,i} \) 是否可视为零（若 \( |h_ {i+1,i}| \le \epsilon (|h_ {ii}|+|h_ {i+1,i+1}|) \)）。如果可视为零，则将矩阵分裂为两个更小的Hessenberg子矩阵，分别处理。 如果 \( 1 \times 1 \) 块分离，得到一个实特征值。 如果 \( 2 \times 2 \) 块分离，计算它的两个特征值（可能是实或共轭复对）。 继续迭代直至所有特征值被提取。 第六步：整体算法流程 输入实方阵 \( A \)。 用Householder变换化 \( A \) 为上Hessenberg矩阵 \( H \)。 重复 ： 如果 \( H \) 的次对角线元足够小，分割 \( H \)。 对当前Hessenberg块，用双重隐式位移做一次QR迭代，更新 \( H \)。 直到所有特征值被求出（对角块为 \( 1\times1 \) 或 \( 2\times2 \)）。 输出 ：所有特征值（实数或共轭复数对）。 总结 Hessenberg-QR算法将非对称矩阵特征值计算分解为“预处理+迭代”两阶段，通过上Hessenberg化降低迭代成本，再通过隐式位移QR迭代高效、稳定地逼近特征值。这是EISPACK、LAPACK等库中 dgeev 等例程的核心算法。