矩阵海森伯格形式转换：实数矩阵正常，复数矩阵失效求助

我用Python实现了一个生成方阵海森伯格形式的函数（仅临时使用linalg.norm，后续会替换），处理实数矩阵完全正常，但复数矩阵时失效。严格遵循算法实现，按说不该出问题。

诡异的是，转换后矩阵的特征值和预期值非常接近。矩阵Q看起来是酉矩阵，但用Q和原始矩阵相乘后得不到正确的海森伯格形式，怀疑问题出在u向量的计算环节，尝试过换各种范数，结果反而更糟。

当前函数代码：

import numpy as np

def Hessenberg(A):
    n = np.shape(A)[0]
    H = np.copy(A)
    Q = np.eye(n,dtype=complex)
    for k in range(n - 2):
        u = np.copy(H[k+1:, k])
        sign_x = np.array((-np.sign(u[0].real), -np.sign(u[0].imag)))
        sign_x /= np.linalg.norm(sign_x)
        u[0] = u[0] - (sign_x[0]+ 1j*sign_x[1]) * np.linalg.norm(u)
        u /= np.linalg.norm(u)
        
        H[k+1:, k:] -= 2 * np.outer(u, np.dot(u.conj().T, H[k+1:, k:]))
        H[:, k+1:] -= 2 * np.outer(np.dot(H[:, k+1:], u), u.conj().T)
        Q[k+1:, :] -= 2 * np.outer(u, np.dot(u.conj().T, Q[k+1:, :]))
    return H, Q.T 

问题根源分析 #

问题出在复数Householder变换的符号选择上。当前对复数的符号处理是分别取实部和虚部的符号再归一化，不符合复数Householder变换的标准定义。

标准复数Householder变换中，向量u的构造逻辑是：

• 计算目标向量的范数norm_x = ||x||
• 构造符号因子sign = x[0] / |x[0]|（若x[0]为0则取1）
• 构造向量x = x + sign * norm_x * e1（注意是加法，符号方向需对应）
• 最后归一化x得到u

代码里的符号计算错误，导致Householder反射矩阵不正确，进而无法正确消去H矩阵次对角线下方的元素。

修复后的代码 #

import numpy as np

def Hessenberg(A):
    n = np.shape(A)[0]
    H = np.copy(A)
    Q = np.eye(n, dtype=complex)
    for k in range(n - 2):
        x = H[k+1:, k].copy()
        norm_x = np.linalg.norm(x)
        if norm_x == 0:
            continue  # 向量已为0，无需变换
        # 计算复数符号因子：x[0]/|x[0]|
        sign = 1.0 + 0.0j if x[0] == 0 else x[0] / np.abs(x[0])
        # 构造Householder向量
        x[0] += sign * norm_x
        u = x / np.linalg.norm(x)
        
        # 更新H矩阵右下块
        H[k+1:, k:] -= 2 * np.outer(u, u.conj().T @ H[k+1:, k:])
        # 更新H矩阵右上部
        H[:, k+1:] -= 2 * np.outer(H[:, k+1:] @ u, u.conj().T)
        # 更新Q矩阵
        Q[k+1:, :] -= 2 * np.outer(u, u.conj().T @ Q[k+1:, :])
    return H, Q.T

验证说明 #

• 修复后的代码正确实现了复数Householder变换的符号选择，确保次对角线下方元素被消去
• 验证时可检查H矩阵次对角线下方元素是否接近0，以及Q.conj().T @ A @ Q是否与H一致
• 特征值接近是因为Householder变换属于相似变换，特征值保持不变，但之前的错误变换导致H矩阵结构不符合海森伯格形式

内容的提问来源于stack exchange，提问作者Gauwain Savary