Householder三对角化算法在求解对称矩阵特征值问题中的核心作用与步骤详解

我将为您详细讲解Householder三对角化算法——这是将实对称矩阵通过正交相似变换化为三对角矩阵的关键算法，为后续高效求解特征值奠定基础。

一、算法背景与问题描述

对于n×n实对称矩阵A，我们需要计算其特征值。直接对A应用QR迭代等算法计算量大（O(n³)每步）。但若先将A化为三对角矩阵T（仅主对角线及相邻两条对角线上有非零元素），则后续QR迭代等算法的复杂度降至O(n²)每步，且能保持对称性。

核心问题：寻找正交矩阵Q，使得QᵀAQ = T，其中T是三对角矩阵。

二、Householder反射器的数学原理

1. 基本定义：Householder反射矩阵H = I - 2uuᵀ/(uᵀu)，其中u是非零向量

   • H是正交矩阵（HᵀH = I）
   • H是对称矩阵（Hᵀ = H）
   • H是反射矩阵（H² = I）

2. 几何意义：H将任意向量x反射到关于u的正交超平面镜像

3. 关键性质：对任意向量x，可构造u使Hx = ±‖x‖₂e₁

   • 取u = x ∓ ‖x‖₂e₁
   • 这使我们能将向量后n-k个分量"清零"

三、逐步三对角化过程

第1步：初始化
设A₀ = A，维数n×n
我们将进行n-2步变换，每步消除一列和一行的非三对角元素

第2步：第k次变换（k=1到n-2）
此时矩阵A₍k-1₎已部分三对角化，其前k-1行/列已为三对角形式

1. 构造Householder向量：

   • 取当前矩阵第k列的后n-k个元素：x = a₍k+1:n, k₎
   • 计算σ = sign(x₁)‖x‖₂（符号选择为避免数值问题）
   • 计算u = x + σe₁（其中e₁是n-k维单位向量）
   • 令v = u/‖u‖₂（规范化）

2. 构造Householder矩阵：

   • Pₖ = I₍n-k₎ - 2vvᵀ

3. 嵌入到n维空间：

   • Qₖ = diag(Iₖ, Pₖ)（分块对角矩阵）

4. 应用相似变换：

   • Aₖ = QₖᵀA₍k-1₎Qₖ
   • 这会将第k列的a₍k+2:n, k₎清零，同时对称地将第k行的a₍k, k+2:n₎清零

第3步：矩阵结构变化详解

以n=5，k=1为例：
原始对称矩阵A：

[× × × × ×]
[× × × × ×]
[× × × × ×]
[× × × × ×]
[× × × × ×]

第1步后A₁：

[× β 0 0 0]
[β × × × ×]
[0 × × × ×]
[0 × × × ×]
[0 × × × ×]

其中β是产生的次对角元

四、算法伪代码

输入：n×n对称矩阵A
输出：三对角矩阵T，正交变换矩阵Q（可选）

T = A
Q = Iₙ（若需要累积变换）

for k = 1 to n-2 do
    // 提取第k列的后n-k个元素
    x = T[k+1:n, k]
    
    // 计算Householder向量
    σ = sign(x₁) * ‖x‖₂
    u = x.copy()
    u[1] = u[1] + σ
    
    // 避免除零
    if ‖u‖₂ ≠ 0 then
        u = u / ‖u‖₂
    
    // 更新T的右下子矩阵
    // 计算p = 2 * (T[k+1:n, k:n] * u) * uᵀ
    p = 2 * (T[k+1:n, k:n] * u)
    p = p - 2 * u * (uᵀ * p)  // 更稳定的计算
    
    T[k+1:n, k:n] = T[k+1:n, k:n] - u * pᵀ
    T[k:n, k+1:n] = (T[k+1:n, k:n])ᵀ  // 保持对称性
    
    // 更新Q（若需要）
    if 需要Q then
        Q[:, k+1:n] = Q[:, k+1:n] - 2 * (Q[:, k+1:n] * u) * uᵀ
    end if
    
    // 显式设置零元素（数值上实际已有）
    T[k+2:n, k] = 0
    T[k, k+2:n] = 0
end for

五、数值稳定性与计算复杂度

1. 稳定性：Householder变换是数值稳定的，不会放大误差
2. 复杂度：
   • 浮点运算：约(4/3)n³ + O(n²)
   • 比Givens旋转的(8/3)n³更高效
3. 存储：可原地存储，三对角元素存于两个一维数组

六、与特征值求解的衔接

得到三对角矩阵T后：

1. 可使用对称QR算法或分治法高效计算T的特征值
2. 若需要特征向量：
   • 累积变换：Q = Q₁Q₂...Q₍n-2₎
   • 计算T的特征向量Y
   • A的特征向量X = QY

七、关键技巧与注意事项

1. 符号选择：sign(x₁)的选择避免u接近零向量
2. 子矩阵更新：避免显式形成Pₖ，直接更新子矩阵
3. 对称性保持：利用对称性只需计算一半元素
4. 缓存优化：块化操作提高缓存利用率

八、示例演示（n=4矩阵）

设A =

[4  1 -2  2]
[1  2  0  1]
[-2 0  3 -2]
[2  1 -2  1]

步骤1（k=1）：
x = [1, -2, 2]ᵀ
‖x‖₂ = √(1+4+4) = 3
σ = 3（x₁=1>0）
u = [1+3, -2, 2]ᵀ = [4, -2, 2]ᵀ
‖u‖₂ = √(16+4+4) = √24
v = u/√24

变换后得三对角矩阵的第一行/列...

九、算法优势

1. 保持对称性：正交相似变换保持矩阵对称性
2. 向后稳定：计算的特征值对应一个微小扰动后的矩阵精确特征值
3. 为后续算法奠基：三对角形式使后续特征值算法效率提升一个数量级
4. 广泛应用：LAPACK、MATLAB等软件库的核心算法

这是对称矩阵特征值计算的标准预处理步骤，将一般稠密对称矩阵的特征值问题转化为三对角矩阵特征值问题，为高效稳定的数值计算奠定基础。