双步位移QR算法在拟上三角矩阵特征值计算中的应用

算法描述

双步位移QR算法是计算实矩阵全部特征值的一种高效数值方法。当我们将一个一般实矩阵通过正交相似变换（例如，先通过Householder约化）转化为上Hessenberg矩阵后，标准的QR算法虽然有效，但在处理实矩阵时，若矩阵有复特征值，则单步位移QR算法必须进入复数域计算，这会显著增加计算量并降低数值稳定性。双步位移QR算法的核心思想是：通过一次使用两次实数位移（通常选取矩阵末尾2x2子矩阵的两个特征值），使得算法在实数域内完成，避免复数运算，从而高效且稳定地计算实Hessenberg矩阵（即拟上三角矩阵，指除了次对角线下方的元素可能非零外，其余结构同于上三角矩阵）的全部特征值。

具体问题：给定一个n×n的实上Hessenberg矩阵H，设计并应用双步位移QR算法，计算H的所有特征值。最终目标是通过一系列正交相似变换，将H化为实的拟上三角矩阵（实Schur型），其对角线上的1x1或2x2块给出了矩阵的全部特征值。

解题过程

第一步：理解基础与动机

1. 特征值问题：对于一个方阵A，求标量λ和非零向量x，使得Ax = λx。λ是特征值。
2. QR算法基本框架：给定矩阵A₀ = A，对于k = 0, 1, 2, ... 执行：
   • QR分解：Aₘ = QₘRₘ，其中Qₘ是正交矩阵，Rₘ是上三角矩阵。
   • 矩阵重构：Aₘ₊₁ = RₘQₘ。
   • 在满足一定条件下，Aₘ会收敛到一个块上三角矩阵（实Schur型），其特征值可由对角线块（1x1或2x2）的特征值给出。
3. 位移加速：为了加速收敛，引入位移σ：Aₘ - σI = QₘRₘ，然后Aₘ₊₁ = RₘQₘ + σI。位移σ通常选择Aₘ的右下角元素或其2x2子矩阵的特征值。
4. 复特征值的困境：如果矩阵有复共轭特征值对，在实数域中，最优的单位移σ也应是复数。这迫使QR分解和随后的乘法在复数域中进行，计算代价高且可能引入复数舍入误差。

第二步：双步位移策略

1. 位移对的选择：设H是当前实上Hessenberg矩阵。我们观察其右下角的2x2子矩阵：
   
   计算这个2x2矩阵的特征值σ₁和σ₂。它们是方程λ² - (hₙ₋₁,ₙ₋₁ + hₙ,ₙ)λ + (hₙ₋₁,ₙ₋₁hₙ,ₙ - hₙ,ₙ₋₁hₙ₋₁,ₙ) = 0的根。这两个特征值要么都是实数，要么是一对共轭复数。
2. 双步思想：即使σ₁和σ₂是复数，它们的和与积都是实数（因为它们是共轭复根二次方程的系数）。双步位移QR算法连续执行两步带有复数位移的QR迭代，但通过巧妙的安排，最终等价于一步完全在实数域内执行的运算。
   • 等效操作：首先执行 (H - σ₁I)(H - σ₂I) = H² - sH + tI，其中 s = σ₁ + σ₂, t = σ₁σ₂。由于σ₁, σ₂是共轭复数或都是实数，s和t是实数。因此，矩阵M = H² - sH + tI是一个实矩阵。
   • 核心等价：如果我们对M进行QR分解：M = QR（其中Q是实的正交矩阵），然后用Q对H进行相似变换：H' = QᵀHQ，那么H'等同于连续执行两次单步位移QR迭代（位移分别为σ₁和σ₂）后的矩阵，但整个过程完全避免了复数运算。

第三步：隐式双步QR迭代的详细步骤

为了避免显式计算M = H² - sH + tI（计算量大且可能破坏稀疏性），双步位移QR算法采用“隐式Q定理”来实现上述变换。

步骤1：计算第一列向量

1. 计算M的第一列向量v₁。由于H是上Hessenberg矩阵，M = H² - sH + tI的第一列仅有前三个元素可能非零：
   • m₁₁ = h₁₁² + h₁₂h₂₁ - s h₁₁ + t
   • m₂₁ = h₂₁(h₁₁ + h₂₂ - s)
   • m₃₁ = h₂₁h₃₂
     其他元素为0。

步骤2：构造初等反射变换
2. 构造一个3x3的Householder反射矩阵P₀，使得P₀ * v₁（v₁是M的第一列向量的前三个元素构成的向量）与第一个坐标轴（e₁ = [1, 0, 0]ᵀ）对齐，即P₀v₁ = γ e₁。
* 这个P₀是一个实的正交矩阵。
* 我们将P₀嵌入到一个n×n的正交矩阵中，形式为：

步骤3：执行相似变换
3. 将P应用于H：H₁ = Pᵀ H P。这个操作只影响H的前三行和前两列。由于P₀只作用在前三行前三列，H₁仍然是上Hessenberg矩阵，但可能在前三条次对角线（即(i, i-2)位置，其中i=3,4,...）引入非零元素，破坏了严格的Hessenberg形。我们称此时的H₁为“上Hessenberg形加上一个小凸起(bulge)”。

步骤4：消去凸起（Bulge Chasing）
4. 接下来的目标是通过一系列额外的正交相似变换（Householder反射或Givens旋转），将H₁恢复为标准的上Hessenberg形H'，而不改变由双步位移QR迭代所决定的特征值序列。这个过程称为“凸起追逐(Bulge Chasing)”。
* 从第二行开始，使用一个Householder反射（或Givens旋转）作用于第2到第4行，消去新引入的、在(4,2)位置的凸起元素。这个变换会影响到第2到第4行和第1到第3列。
* 然后向右下方移动一个位置，对第3到第5行执行类似操作，消去下一个凸起，依此类推。
* 这个追逐过程一直持续到矩阵的右下角。每一次小的正交变换都保持了矩阵与初始矩阵的正交相似性，并逐步将凸起向下向右“推”，直到它被推出矩阵的边界，从而矩阵恢复为上Hessenberg形。
* 隐式Q定理 保证了，如果我们从同一个初始Hessenberg矩阵H出发，并通过正交相似变换得到了两个具有相同第一列的正交矩阵Q和Q̃，且它们都将H化为上Hessenberg形，那么在相差一个对角块（±1）的意义下，Q和Q̃是相同的。在我们的算法中，初始的P₀恰好保证了变换后矩阵与显式计算双步QR迭代有相同的第一列，因此整个追逐过程最终产生的矩阵H'就是我们要找的双步位移QR迭代结果。

步骤5：判断收敛与收缩（Deflation）
5. 经过一次双步位移QR迭代后，检查新的上Hessenberg矩阵H'。
* 如果某个次对角元素h_{i+1,i}的绝对值小于一个预设的容差（例如，ε * (|h_{ii}| + |h_{i+1,i+1}|)），我们可以将其视为零。
* 这允许我们将原特征值问题“收缩”为两个更小的子问题：处理左上角的i×i块和右下角的(n-i)×(n-i)块。它们可以独立求解。
* 收缩极大地提高了计算效率，尤其是当特征值已被分离出来时。

第四步：整体算法流程

1. 输入：一个n×n的实矩阵A。
2. 初始化：通过Householder反射将A相似变换为实上Hessenberg矩阵H。（这是QR算法的标准预处理步骤）
3. 迭代循环：
   • while (未处理的矩阵块尺寸 > 2)
     • 检查当前活动子矩阵（从第l行到第u行）的次对角线元素。若某个h_{i+1,i}可视为零，则设置其为0，并将问题分解为两个更小的子问题。
     • 若当前活动子矩阵尺寸m = u-l+1 <= 2，则直接计算其特征值（对于1x1块就是元素本身，对于2x2块解二次方程），并收缩（l, u更新）。
     • 否则，对当前活动子矩阵（它是上Hessenberg形）执行一次隐式双步位移QR迭代（即上述第三步）。
   • end while
4. 输出：最终得到的拟上三角矩阵（实Schur型）T。T对角线上的1x1块给出实特征值，2x2块给出复共轭特征值对。对应的特征向量可以通过逆迭代等方法从T中进一步求出。

总结

双步位移QR算法通过巧妙利用实数位移对，使得即使处理具有复特征值的实矩阵，也能完全在实数域内进行计算。隐式实现的凸起追逐过程避免了显式计算H²，保持了数值稳定性和计算效率。通过收缩策略，算法可以聚焦于尚未收敛的子矩阵，显著减少计算量。它是现代特征值求解软件包（如LAPACK中的xHSEQR例程）的核心算法，用于可靠且高效地计算实矩阵的全部特征值。