如何在 C 中高效计算矩阵的上 Hessenberg 形式？

2026-11-04数据小陈

结论

可以。Householder 矩阵是单位矩阵的秩一修正，因此无需显式构造 P。直接用 Householder 向量更新 H 的行和列，就能完成相同的相似变换。每轮计算成本可从通用矩阵乘法的 O(n^3) 降到 O(n(n-i))。整个 Hessenberg 化过程的复杂度仍为 O(n^3)，但常数更小。

用向量更新替代矩阵乘法

设当前反射作用于从下标 i 开始的子空间，Householder 向量为 v，并令 beta = 2 / (v^T v)。将 v 的前 i 个位置补零后，嵌入整个矩阵的反射为 P = I - beta vv^T。相似变换可以分两步完成：

左乘：H ← P H。对每一列计算 v^T H，再减去对应的 beta v 倍数。
右乘：H ← H P。对每一行计算 H v，再减去对应的 beta v^T 倍数。

C 中可以按下面的方式更新。这里 m = d - i，vec 是本轮的 Householder 向量。左乘和右乘都会遍历完整的行或列，因此不需要对结构零作额外假设。

// vec[0..m-1] 已构造好
 double vv = 0.0;
for (int a = 0; a < m; ++a) vv += vec[a] * vec[a];

// 零向量对应恒等变换，避免除以零
if (vv > 0.0) {
    double beta = 2.0 / vv;
    double w[d];

    // 左乘：H <- P H
    for (int j = 0; j < d; ++j) {
        double dot = 0.0;
        for (int a = 0; a < m; ++a)
            dot += vec[a] * H[i + a][j];
        w[j] = beta * dot;
    }
    for (int a = 0; a < m; ++a)
        for (int j = 0; j < d; ++j)
            H[i + a][j] -= vec[a] * w[j];

    // 右乘：H <- H P
    for (int r = 0; r < d; ++r) {
        double dot = 0.0;
        for (int a = 0; a < m; ++a)
            dot += H[r][i + a] * vec[a];
        w[r] = beta * dot;
    }
    for (int r = 0; r < d; ++r)
        for (int a = 0; a < m; ++a)
            H[r][i + a] -= w[r] * vec[a];
}

第一步先把向量点积保存到 w，再更新矩阵，避免修改后续还要读取的数据。第二步也按同样的顺序处理。先左乘、再右乘，最终结果等价于 P H P。

实现时需要检查的细节

若 vv == 0，当前向量不需要反射。跳过本轮更新可以避免 1 / vv 产生无效值。
代码中的 A 在初始化后没有参与后续计算，迭代实际使用的是 H。
选择 Householder 向量的符号时，应避免相减导致的消去误差。改写更新逻辑时，向量构造方式要与反射公式保持一致。
double A[d][d] 等变长数组属于 C 的数组语法，能否使用取决于编译器和所用的 C 标准。矩阵较大时，可以改为动态分配工作区。
对一般实矩阵，上述反射是对称的，因此左右乘使用同一个 P。处理复矩阵时，应使用相应的共轭转置公式。

处理实际规模的矩阵时，可以调用 LAPACK 的 dgehrd 完成 Hessenberg 化。它采用专门的数值线性代数实现，省去在应用代码中手写完整的反射更新流程。