Hessenberg矩阵的隐式位移QR迭代算法
字数 3963 2025-12-16 06:31:36

Hessenberg矩阵的隐式位移QR迭代算法

我来为你讲解Hessenberg矩阵的隐式位移QR迭代算法。这个算法是计算一般矩阵所有特征值的核心方法，结合了Hessenberg化减少计算量、位移技术加速收敛，以及隐式QR技巧保证数值稳定性。

题目描述

给定一个 \(n \times n\) 的实矩阵 \(A\)，我们需要计算它的所有特征值（可能是实数或共轭复数对）。直接对原矩阵进行QR迭代计算量巨大（\(O(n^3)\) 每步）。高效的算法是：

  1. 先将 \(A\) 通过相似变换化为上Hessenberg形式 \(H\)（次对角线以下全为0），这只需 \(O(n^3)\) 一次。
  2. 然后对上Hessenberg矩阵 \(H\) 应用带位移的隐式QR迭代，逐步将其化为实Schur形式（对角块为1×1或2×2的块上三角矩阵），从而求出所有特征值。

循序渐进讲解

第1步：为什么要用Hessenberg矩阵？

对于一个一般的稠密矩阵 \(A\)，完整的QR分解每步需要 \(O(n^3)\) 运算量。如果我们能把 \(A\) 相似变换（不改变特征值）成一个特殊形状的矩阵，就能大幅减少每步QR迭代的计算量。

  • 上Hessenberg矩阵 \(H\) 的形式如下：

\[H = \begin{bmatrix} * & * & * & \dots & * \\ * & * & * & \dots & * \\ 0 & * & * & \dots & * \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & \dots & 0 & * & * \end{bmatrix} \]

即所有元素 \(h_{ij} = 0\) 当 \(i > j+1\)。简单说，只有主对角线、上次对角线和上三角区域有非零元。

  • 关键好处：对上Hessenberg矩阵进行QR迭代，每步的计算量仅为 \(O(n^2)\)，而不是 \(O(n^3)\)。

将任意矩阵约化为上Hessenberg形式，可以使用一系列Householder变换（或Givens旋转），这是一个稳定的 \(O(n^3)\) 预处理步骤。

第2步：基本的QR迭代算法思想

  1. 从 \(A_0 = A\) 开始。
  2. 对 \(k = 0, 1, 2, \dots\)，重复：
    • 计算 \(A_k = Q_k R_k\) （QR分解）
    • 令 \(A_{k+1} = R_k Q_k\)
  3. 在满足一定条件下，\(A_k\) 会收敛到一个上三角矩阵（或实Schur形式），其对角线（或2×2块的对角线）元素就是特征值。

为什么有效？ 因为 \(A_{k+1} = Q_k^T A_k Q_k\)，所以所有 \(A_k\) 都相似，特征值不变。这个过程本质上是幂法的推广。

问题：基本QR迭代收敛可能很慢。

第3步：位移技术加速收敛

为了加速收敛，我们引入位移（shift）。其思想是：如果 \(\mu\) 是矩阵特征值的一个好的近似，那么对 \((A - \mu I)\) 进行QR迭代，收敛会快得多。

带单位移的QR迭代步骤（假设 \(A_k\) 已经是上Hessenberg矩阵）：

  1. 选择一个位移量 \(\mu_k\)（例如，取 \(A_k\) 的右下角元素，即Rayleigh商位移）。
  2. 计算 \(A_k - \mu_k I = Q_k R_k\)。
  3. 计算 \(A_{k+1} = R_k Q_k + \mu_k I\)。
  4. 可以证明 \(A_{k+1} = Q_k^T A_k Q_k\)，仍是上Hessenberg矩阵。

通过巧妙选择位移（如双位移策略处理复共轭特征值对），可以使收敛速度达到二次甚至三次。

第4步：隐式QR技巧——算法的核心

直接按公式 \(A_{k+1} = Q_k^T A_k Q_k\) 计算有两个问题：1) 需要显式计算QR分解；2) 会引入不必要的舍入误差。

隐式QR定理（Implicit Q Theorem） 是关键。它告诉我们：如果一个正交矩阵 \(Q\) 能把一个上Hessenberg矩阵 \(H\) 变换成另一个上Hessenberg矩阵 \(\tilde{H}\)（假设两者都有正上次对角元），并且 \(Q\) 的第一列已知，那么 \(Q\) 本质上唯一确定，\(\tilde{H}\) 也唯一确定（只差对角块的排列）。

隐式QR迭代的步骤（以单位移为例）：

  1. 计算第一列：我们想要一个正交矩阵 \(Q\)，使得 \(\tilde{H} = Q^T H Q\) 也是上Hessenberg矩阵，并且执行了位移 \(\mu\) 的QR步。根据QR算法，\(Q\) 应该是矩阵 \(H - \mu I\) 的QR分解中的 \(Q\) 矩阵。关键是，我们不需要计算整个 \(Q\)。

    • 设 \(x = (H - \mu I) e_1\)，即 \(H - \mu I\) 的第一列。我们只需要构造一个 \(n \times n\) 的正交矩阵 \(Q\)，它的第一列与 \(x\) 平行即可。这通常用一个 Householder反射器 \(P_0\) 就能做到。

  2. 相似变换：用这个 \(P_0\) 对 \(H\) 做相似变换：\(H_1 = P_0 H P_0^T\)。这会把 \(H\) 的Hessenberg形状破坏，在 (3,1), (2,1) 等位置引入“凸起”（bulge），我们称这个非零元为“bulge”。

  3. Bulge追逐（Bulge Chasing）：我们的目标是把 \(H_1\) 变回上Hessenberg形式 \(\tilde{H}\)，同时保持与 \(H\) 的正交相似。

    • 针对 \(H_1\) 中破坏形状的bulge，我们设计另一个Householder（或Givens）变换 \(P_1\)，从左边乘上去消去bulge，但从右边乘上 \(P_1^T\) 又会在稍后的位置引入新的bulge。
    • 通过一系列精心设计的正交变换 \(P_1, P_2, \dots, P_{n-2}\)，我们可以将这个bulge从矩阵的左上角“追逐”到右下角，最终消失，得到一个新的上Hessenberg矩阵 \(\tilde{H}\)。
    • 整个过程是：\(\tilde{H} = P_{n-2}^T \dots P_1^T P_0^T H P_0 P_1 \dots P_{n-2}\)。

  4. 结果：根据隐式Q定理，这样得到的 \(\tilde{H}\) 恰好就是通过显式位移QR迭代（\(H - \mu I = QR, \tilde{H} = Q^T H Q\)）得到的矩阵。我们实现了同样的效果，但避免了显式形成 \(H - \mu I\) 和进行QR分解，数值上更稳定，计算也更高效。

第5步：完整的隐式QR算法流程

  1. 输入：一个一般的 \(n \times n\) 实矩阵 \(A\)。
  2. Hessenberg化：通过Householder变换将 \(A\) 相似变换为上Hessenberg矩阵 \(H\)。令 \(H_0 = H\)。
  3. 迭代直到收敛：
    • 检查 \(H\) 的次对角线元素 \(h_{i+1, i}\)。如果某个 \(|h_{i+1, i}|\) 足够小（如小于机器精度乘以 \((|h_{ii}| + |h_{i+1, i+1}|)\)），则将其设为0。
    • 这样矩阵被解耦成更小的上Hessenberg块。对尚未收敛的最后一个 \(m \times m\) 块（通常是右下角）进行操作。
    • 位移选择：
      • 如果该块的右下角2×2矩阵有一对复共轭特征值，则采用双位移策略（Francis位移）。这等价于执行两次连续的实位移QR步，但通过一个巧妙的技巧，我们可以构造一个实的变换矩阵 \(Q\)，其第一列由该2×2矩阵的特征值决定。
      • 否则，采用单位移（如取右下角元素）。
    • 隐式QR步：对当前活跃的 \(m \times m\) 块，使用上述的“bulge chasing”过程，应用一个（单位移）或一串（双位移）正交相似变换，将其变为新的上Hessenberg矩阵，同时该块的右下角元素会更接近其特征值。
  4. 输出：最终 \(H\) 被化为实Schur形式：一个块上三角矩阵，对角块为1×1（实特征值）或2×2（共轭复特征值对）。从这些对角块中可提取出全部特征值。

第6步：总结与意义

Hessenberg矩阵的隐式位移QR迭代算法 是数值线性代数中特征值计算的基石。它的成功源于三个核心思想的结合：

  1. Hessenberg化：将问题简化到计算量更小的矩阵类。
  2. 位移策略：通过位移加速收敛，并利用双位移策略处理复特征值，同时保持实数运算。
  3. 隐式QR技术：基于隐式Q定理，通过“bulge chasing”稳定高效地实现位移QR步，避免了显式QR分解和潜在的数值问题。

这个算法被广泛应用于MATLAB的 eig 函数（对于稠密非对称矩阵）的核心步骤中，是工程和科学计算中求解特征值问题的标准工具。

Hessenberg矩阵的隐式位移QR迭代算法 我来为你讲解Hessenberg矩阵的隐式位移QR迭代算法。这个算法是计算一般矩阵所有特征值的核心方法，结合了Hessenberg化减少计算量、位移技术加速收敛，以及隐式QR技巧保证数值稳定性。 题目描述 给定一个 \( n \times n \) 的实矩阵 \( A \)，我们需要计算它的所有特征值（可能是实数或共轭复数对）。直接对原矩阵进行QR迭代计算量巨大（\(O(n^3)\) 每步）。高效的算法是： 先将 \( A \) 通过相似变换化为上Hessenberg形式 \( H \)（次对角线以下全为0），这只需 \(O(n^3)\) 一次。 然后对上Hessenberg矩阵 \( H \) 应用带位移的隐式QR迭代，逐步将其化为实Schur形式（对角块为1×1或2×2的块上三角矩阵），从而求出所有特征值。 循序渐进讲解 第1步：为什么要用Hessenberg矩阵？ 对于一个一般的稠密矩阵 \( A \)，完整的QR分解每步需要 \(O(n^3)\) 运算量。如果我们能把 \( A \) 相似变换（不改变特征值）成一个特殊形状的矩阵，就能大幅减少每步QR迭代的计算量。 上Hessenberg矩阵 \( H \) 的形式如下： \[ H = \begin{bmatrix} & * & * & \dots & * \\ & * & * & \dots & * \\ 0 & * & * & \dots & * \\ \vdots & \ddots & \ddots & \ddots & \vdots \\ 0 & \dots & 0 & * & * \end{bmatrix} \] 即所有元素 \( h_ {ij} = 0 \) 当 \( i > j+1 \)。简单说，只有主对角线、上次对角线和上三角区域有非零元。 关键好处 ：对上Hessenberg矩阵进行QR迭代，每步的计算量仅为 \(O(n^2)\)，而不是 \(O(n^3)\)。 将任意矩阵约化为上Hessenberg形式，可以使用一系列Householder变换（或Givens旋转），这是一个稳定的 \(O(n^3)\) 预处理步骤。 第2步：基本的QR迭代算法思想 从 \( A_ 0 = A \) 开始。 对 \( k = 0, 1, 2, \dots \)，重复： 计算 \( A_ k = Q_ k R_ k \) （QR分解） 令 \( A_ {k+1} = R_ k Q_ k \) 在满足一定条件下，\( A_ k \) 会收敛到一个上三角矩阵（或实Schur形式），其对角线（或2×2块的对角线）元素就是特征值。 为什么有效？ 因为 \( A_ {k+1} = Q_ k^T A_ k Q_ k \)，所以所有 \( A_ k \) 都相似，特征值不变。这个过程本质上是幂法的推广。 问题 ：基本QR迭代收敛可能很慢。 第3步：位移技术加速收敛 为了加速收敛，我们引入 位移（shift） 。其思想是：如果 \( \mu \) 是矩阵特征值的一个好的近似，那么对 \( (A - \mu I) \) 进行QR迭代，收敛会快得多。 带单位移的QR迭代步骤 （假设 \( A_ k \) 已经是上Hessenberg矩阵）： 选择一个位移量 \( \mu_ k \)（例如，取 \( A_ k \) 的右下角元素，即Rayleigh商位移）。 计算 \( A_ k - \mu_ k I = Q_ k R_ k \)。 计算 \( A_ {k+1} = R_ k Q_ k + \mu_ k I \)。 可以证明 \( A_ {k+1} = Q_ k^T A_ k Q_ k \)，仍是上Hessenberg矩阵。 通过巧妙选择位移（如双位移策略处理复共轭特征值对），可以使收敛速度达到二次甚至三次。 第4步：隐式QR技巧——算法的核心 直接按公式 \( A_ {k+1} = Q_ k^T A_ k Q_ k \) 计算有两个问题：1) 需要显式计算QR分解；2) 会引入不必要的舍入误差。 隐式QR定理（Implicit Q Theorem） 是关键。它告诉我们：如果一个正交矩阵 \( Q \) 能把一个上Hessenberg矩阵 \( H \) 变换成另一个上Hessenberg矩阵 \( \tilde{H} \)（假设两者都有正上次对角元），并且 \( Q \) 的第一列已知，那么 \( Q \) 本质上唯一确定，\( \tilde{H} \) 也唯一确定（只差对角块的排列）。 隐式QR迭代的步骤 （以单位移为例）： 计算第一列 ：我们想要一个正交矩阵 \( Q \)，使得 \( \tilde{H} = Q^T H Q \) 也是上Hessenberg矩阵，并且执行了位移 \( \mu \) 的QR步。根据QR算法，\( Q \) 应该是矩阵 \( H - \mu I \) 的QR分解中的 \( Q \) 矩阵。关键是，我们 不需要 计算整个 \( Q \)。 设 \( x = (H - \mu I) e_ 1 \)，即 \( H - \mu I \) 的第一列。我们只需要构造一个 \( n \times n \) 的正交矩阵 \( Q \)，它的第一列与 \( x \) 平行即可。这通常用一个 Householder反射器 \( P_ 0 \) 就能做到。 相似变换 ：用这个 \( P_ 0 \) 对 \( H \) 做相似变换：\( H_ 1 = P_ 0 H P_ 0^T \)。这会把 \( H \) 的Hessenberg形状 破坏 ，在 (3,1), (2,1) 等位置引入“凸起”（bulge），我们称这个非零元为“bulge”。 Bulge追逐（Bulge Chasing） ：我们的目标是把 \( H_ 1 \) 变回上Hessenberg形式 \( \tilde{H} \)，同时保持与 \( H \) 的正交相似。 针对 \( H_ 1 \) 中破坏形状的bulge，我们设计另一个Householder（或Givens）变换 \( P_ 1 \)，从左边乘上去消去bulge，但从右边乘上 \( P_ 1^T \) 又会在稍后的位置引入新的bulge。 通过一系列精心设计的正交变换 \( P_ 1, P_ 2, \dots, P_ {n-2} \)，我们可以将这个bulge从矩阵的左上角“追逐”到右下角，最终消失，得到一个新的上Hessenberg矩阵 \( \tilde{H} \)。 整个过程是：\( \tilde{H} = P_ {n-2}^T \dots P_ 1^T P_ 0^T H P_ 0 P_ 1 \dots P_ {n-2} \)。 结果 ：根据隐式Q定理，这样得到的 \( \tilde{H} \) 恰好就是 通过显式位移QR迭代（\( H - \mu I = QR, \tilde{H} = Q^T H Q \)）得到的矩阵。我们实现了同样的效果，但避免了显式形成 \( H - \mu I \) 和进行QR分解，数值上更稳定，计算也更高效。 第5步：完整的隐式QR算法流程 输入 ：一个一般的 \( n \times n \) 实矩阵 \( A \)。 Hessenberg化 ：通过Householder变换将 \( A \) 相似变换为上Hessenberg矩阵 \( H \)。令 \( H_ 0 = H \)。 迭代直到收敛 ： 检查 \( H \) 的次对角线元素 \( h_ {i+1, i} \)。如果某个 \( |h_ {i+1, i}| \) 足够小（如小于机器精度乘以 \( (|h_ {ii}| + |h_ {i+1, i+1}|) \)），则将其设为0。 这样矩阵被 解耦 成更小的上Hessenberg块。对尚未收敛的最后一个 \( m \times m \) 块（通常是右下角）进行操作。 位移选择 ： 如果该块的右下角2×2矩阵有一对复共轭特征值，则采用 双位移策略 （Francis位移）。这等价于执行两次连续的实位移QR步，但通过一个巧妙的技巧，我们可以构造一个 实 的变换矩阵 \( Q \)，其第一列由该2×2矩阵的特征值决定。 否则，采用 单位移 （如取右下角元素）。 隐式QR步 ：对当前活跃的 \( m \times m \) 块，使用上述的“bulge chasing”过程，应用一个（单位移）或一串（双位移）正交相似变换，将其变为新的上Hessenberg矩阵，同时该块的右下角元素会更接近其特征值。 输出 ：最终 \( H \) 被化为 实Schur形式 ：一个块上三角矩阵，对角块为1×1（实特征值）或2×2（共轭复特征值对）。从这些对角块中可提取出全部特征值。 第6步：总结与意义 Hessenberg矩阵的隐式位移QR迭代算法 是数值线性代数中特征值计算的基石。它的成功源于三个核心思想的结合： Hessenberg化 ：将问题简化到计算量更小的矩阵类。 位移策略 ：通过位移加速收敛，并利用双位移策略处理复特征值，同时保持实数运算。 隐式QR技术 ：基于隐式Q定理，通过“bulge chasing”稳定高效地实现位移QR步，避免了显式QR分解和潜在的数值问题。 这个算法被广泛应用于MATLAB的 eig 函数（对于稠密非对称矩阵）的核心步骤中，是工程和科学计算中求解特征值问题的标准工具。