基于Householder反射的约化一般矩阵为Hessenberg形式的算法详解

我将为您讲解一个在线性代数数值计算中非常重要的预处理算法：将一般矩阵（特别是稠密矩阵）通过相似变换化为上Hessenberg形式，这是QR算法等特征值计算方法的必要前置步骤。

题目描述

给定一个n×n的实矩阵A，我们希望找到一个正交相似变换QᵀAQ = H，其中H是一个上Hessenberg矩阵（即当i > j+1时，h_ij = 0）。这能极大简化后续的特征值计算（如QR算法）。

什么是上Hessenberg矩阵？
形式上，矩阵H是上Hessenberg矩阵，如果：

• 对所有i > j+1，h_ij = 0
• 即主对角线以下的第一条次对角线（subdiagonal）可以有非零元，但更下面的元素都为零

目标： 构造一个正交矩阵Q（实际上是若干个Householder反射器的乘积），使得QᵀAQ是上Hessenberg形式。

循序渐进讲解

第一步：理解为什么需要化为Hessenberg形式

1. 计算效率：QR算法迭代计算特征值时，对完全稠密矩阵的每次QR分解是O(n³)操作。而对Hessenberg矩阵的QR分解只需O(n²)运算。
2. 保持结构：QR迭代在Hessenberg矩阵上运算时，能保持Hessenberg结构不变，计算量大大减少。
3. 数值稳定性：使用正交变换（Householder反射）而不是高斯消元，能保证良好的数值稳定性。

第二步：Householder反射器基础回顾

一个Householder反射器形式为：
P = I - 2vvᵀ/(vᵀv)，其中v是一个非零向量。

性质：

• P是正交且对称的：Pᵀ = P, PᵀP = I
• 对任意向量x，选择v = x ± ‖x‖₂e₁，则Px = ∓‖x‖₂e₁（将x"反射"到第一个坐标轴上）

第三步：算法的核心思想（分步约化）

我们从左到右（列）进行约化。对于n×n矩阵A，需要n-2步。

第1步：处理第1列

• 目标：将第1列中第3行到第n行的元素化为0
• 只看A的第一列：x = [a₂₁, a₃₁, ..., aₙ₁]ᵀ
• 构造一个(n-1)维的Householder反射器P₁'，使得P₁'x = [*, 0, ..., 0]ᵀ
• 构造n×n的正交矩阵P₁ = diag(1, P₁')
• 计算A₁ = P₁AP₁ᵀ = P₁AP₁（因为P₁对称正交）

关键洞察：从右边也乘以P₁ᵀ是为了保持相似变换（A₁ = P₁AP₁ᵀ），但这样会影响第一行。因为P₁只改变了第2到n行和列，所以第一行也被改变了，但这是可以接受的。

第四步：第k步的通用模式

假设已完成k-1步，矩阵已经是部分Hessenberg形式：

• 前k-1列已满足Hessenberg结构
• 第k列中，行k+2到n还有非零元素需要清零

第k步：

1. 考虑子矩阵：只看A的(k+1):n, k列（即第k列的下部）
2. 令x = [a_{k+1,k}, a_{k+2,k}, ..., a_{n,k}]ᵀ
3. 构造一个(n-k)维的Householder反射器P_k'，使得P_k'x = [*, 0, ..., 0]ᵀ
4. 构造n×n的正交矩阵P_k = diag(I_k, P_k')，其中I_k是k×k单位阵
5. 计算A ← P_kAP_kᵀ

注意：从左边乘以P_k（A ← P_kA）会清零第k列下方的元素；从右边乘以P_kᵀ（A ← AP_kᵀ）会保持相似性，但会填充相应的行。

第五步：具体示例（4×4矩阵）

设A是4×4矩阵：

A = [a11 a12 a13 a14]
    [a21 a22 a23 a24]
    [a31 a32 a33 a34]
    [a41 a42 a43 a44]

第1步（k=1）：

• 看第1列下部：x = [a21, a31, a41]ᵀ
• 构造3×3的Householder反射器P₁'，使得P₁'x = [β₁, 0, 0]ᵀ
• P₁ = diag(1, P₁')
• 计算A₁ = P₁AP₁ᵀ

结果形式：

A₁ = [h11 h12 h13 h14]
     [h21 h22 h23 h24]
     [0   h32 h33 h34]
     [0   h42 h43 h44]

现在第1列的下部已清零（除了h21）

第2步（k=2）：

• 看第2列下部：x = [h32, h42]ᵀ
• 构造2×2的Householder反射器P₂'，使得P₂'x = [β₂, 0]ᵀ
• P₂ = diag(I_2, P₂')，其中I_2是2×2单位阵
• 计算H = P₂A₁P₂ᵀ

最终Hessenberg形式：

H = [h11 h12 h13 h14]
    [h21 h22 h23 h24]
    [0   h32 h33 h34]
    [0   0   h43 h44]

第六步：算法伪代码

输入：n×n实矩阵A
输出：上Hessenberg矩阵H，正交矩阵Q使得QᵀAQ = H

1. 初始化Q = I_n (n×n单位矩阵)
2. 对于k = 1 到 n-2:
   a. 令x = A[k+1:n, k]  (第k列的下部)
   b. 计算Householder向量v_k，使得P_k'x = β_k * e_1
      (其中P_k' = I - 2v_kv_kᵀ/(v_kᵀv_k))
   c. 从左边作用：A[k+1:n, k:n] = (I - 2v_kv_kᵀ/(v_kᵀv_k)) * A[k+1:n, k:n]
   d. 从右边作用：A[1:n, k+1:n] = A[1:n, k+1:n] * (I - 2v_kv_kᵀ/(v_kᵀv_k))
   e. 累积Q：Q[1:n, k+1:n] = Q[1:n, k+1:n] * (I - 2v_kv_kᵀ/(v_kᵀv_k))
      (或者保存v_k，后续需要时再生成Q)
3. H = A

计算量：约(10/3)n³次浮点运算，是标准QR分解的一半。

第七步：关键细节与技巧

1. 存储效率：

   • 不需要显式存储P_k，只需存储Householder向量v_k
   • v_k可以存储在A的下三角部分（被清零的位置）
   • 最终的H存储在A的上三角和第一条次对角线上

2. 数值稳定性：

   • 在构造Householder向量时，选择符号以避免减法抵消：
     v = x + sign(x₁)‖x‖₂e₁
   • 或者更稳定的：v = x + sign(x₁)‖x‖₂e₁，其中sign(0)=1

3. 处理复数矩阵：

   • 对复矩阵，使用Hermitian反射器：P = I - 2vvᴴ/(vᴴv)
   • 算法结构完全相同

第八步：完整示例（数值计算）

设：

A = [1  2  3  4]
    [5  6  7  8]
    [9  10 11 12]
    [13 14 15 16]

第一步：
x = [5, 9, 13]ᵀ, ‖x‖₂ ≈ 16.8819
v = x + sign(5)×‖x‖₂×[1,0,0]ᵀ = [5+16.8819, 9, 13]ᵀ = [21.8819, 9, 13]ᵀ
归一化：v = v/‖v‖₂
构造P₁'，计算A₁ = P₁AP₁ᵀ
得到：

A₁ ≈ [1    -2.182  -2.773  -3.364]
     [-16.88 20.91  23.73   26.55]
     [0     -0.465  -0.530  -0.595]
     [0     -0.728  -0.614  -0.500]

第二步：
x = [-0.465, -0.728]ᵀ
类似构造P₂'，计算H = P₂A₁P₂ᵀ
最终得到Hessenberg形式。

第九步：应用与意义

1. QR算法预处理：这是隐式QR算法高效计算特征值的关键第一步
2. Krylov子空间方法：Arnoldi过程本质上是计算Hessenberg形式的另一种方式
3. 特征值问题简化：从O(n³)降到O(n²)的迭代计算

总结：基于Householder反射的Hessenberg化算法通过n-2步正交相似变换，将一般稠密矩阵化为上Hessenberg形式，为后续特征值计算提供了高效的预处理步骤，同时保持了良好的数值稳定性。