豪斯霍尔德变换

豪斯霍尔德变换（Householder transformation）是关于过原点超平面的正交反射。对非零实向量 $v$，矩阵为 $$H=I-2\frac{v{v}^{𝖳}}{v^{𝖳}v}.$$ 它反转法向量 $v$ 的分量，保留与 $v$ 垂直的分量，常用于QR分解中一次消去一整段列元素。

从分量分解证明反射

任意向量写为 $x=x_{\perp }+\alpha v$，其中 $v^{𝖳}{x}_{\perp }=0$。代入公式得到 $Hx=x_{\perp }-\alpha v$。因此 $H^{2}x=x$ 对所有 $x$ 成立，且 $H^{𝖳}=H$，从而 $H^{𝖳}H=I$。长度保持、两次回到原位、反射面内不动，均由同一分解得到。法向一维上的特征值是负一，其正交补上的特征值为一，所以行列式为负一。

把3,4变成一个非零分量

设 $x=(3,4{)}^{𝖳}$，希望映到与 $e_1$ 同向或反向、长度仍为五的向量。选 $y=(-5,0{)}^{𝖳}$，令 $v=x-y=(8,4{)}^{𝖳}$。有 $v^{𝖳}v=80$，因此 $$H=I-\frac{1}{40}\left(\begin{array}{cc}64& 32\\ 32& 16\end{array}\right)=\left(\begin{array}{cc}-3/5& -4/5\\ -4/5& 3/5\end{array}\right).$$ 相乘得 $Hx=(-5,0{)}^{𝖳}$。两点的中点 $(-1,2)$ 满足 $2x_1+x_2=0$，正落在反射面上。

反射法线方向是2,1，与输入和输出的差8,4平行。

一般地，若 $\Vert x\Vert =\Vert y\Vert$ 且 $x\ne y$，取 $v=x-y$。由 $$v^{𝖳}x=\Vert x{\Vert }^2-y^{𝖳}x=\frac{1}{2}{v}^{𝖳}v$$ 可知 $Hx=x-v=y$。这证明构造，不是仅在主例上验证。

稳定符号与不显式形成H

为消去向量 $x\in {ℝ}^n$ 的后 $n-1$ 个分量，通常选 $\alpha =-\mathrm{sign}(x_1)\Vert x\Vert$，零处可取符号一，再令 $v=x-\alpha e_1$。相反符号避免当 $x_1\approx \Vert x\Vert$ 时发生两个接近数相减的消减误差。

若归一化成 $u=v/\Vert v\Vert$，应用只需 $Hz=z-2u(u^{𝖳}z)$。这需要一次点积与一次向量更新，时间和额外存储均为 $O(n)$，不必创建 $n\times n$ 矩阵。

以下函数输入非空有限实数序列，不改写输入；返回单位反射法向 $u$ 和首分量目标 $\alpha$。全零向量返回零法向，约定它代表恒等变换。先除以最大绝对分量，再计算归一化方向，既避免 $x_1-\alpha$ 溢出，也避免用已经舍入的次正规范数作分母。目标 $\alpha$ 单独按原尺度计算；范数超出浮点数可表示范围时明确抛出 OverflowError。

Python 3

import math

def householder_vector(values):
    x = [float(a) for a in values]
    if not x or not all(math.isfinite(a) for a in x):
        raise ValueError("require a nonempty finite real vector")
    scale = max(abs(a) for a in x)
    if scale == 0:
        return [0.0] * len(x), 0.0
    scaled = [a / scale for a in x]
    scaled_norm = math.hypot(*scaled)
    norm = math.hypot(*x)
    if not math.isfinite(norm):
        raise OverflowError("vector norm is not representable")
    sign = math.copysign(1.0, x[0])
    w = [a / scaled_norm for a in scaled]
    w[0] += sign
    length = math.hypot(*w)
    return [a / length for a in w], -sign * norm

def apply_reflector(u, z):
    if len(u) != len(z):
        raise ValueError("dimension mismatch")
    u, z = list(map(float, u)), list(map(float, z))
    if not z or not all(math.isfinite(a) for a in u + z):
        raise ValueError("require nonempty finite vectors")
    if not any(u):
        return z.copy()  # zero normal represents the identity
    scale = max(abs(b) for b in z)
    if scale == 0:
        return z.copy()
    scaled = [b / scale for b in z]
    dot = math.fsum(a * b for a, b in zip(u, scaled))
    result = [scale * (b - 2 * a * dot) for a, b in zip(u, scaled)]
    if not all(math.isfinite(b) for b in result):
        raise OverflowError("reflected component is not representable")
    return result

u, alpha = householder_vector([3, 4])
print(alpha)                       # -5.0
print(apply_reflector(u, [3, 4]))   # approximately [-5.0, 0.0]

apply_reflector 的 $u$ 应来自前一函数或已知单位向量，不能传任意未归一化法向。应用时也先缩放待反射向量，避免 $2u(u^{𝖳}z)$ 的中间乘积溢出；结果分量不能表示时抛出 OverflowError。浮点数下末尾零通常只是接近零；在次正规范围，输出还受浮点刻度限制。例如两个分量都是 math.ulp(0.0) 时，法向仍应为约 $(0.92388,0.38268)$，不能先用已舍入的原尺度范数反除来求方向。

变式输入 $(-2,0)$，构造把它映到 $(2,0)$；输入零向量则保持零。反射的行列式负一与吉文斯旋转的行列式一不同，二者都保长度，但几何动作并不相同。

参考文献

Sheldon Axler，Linear Algebra Done Right，第四版，§7D（QR分解）；Higham的Householder作者说明；本文例子与图解独立构造。