aiwiki.page
中文
数学 / householder-transformation

豪斯霍尔德变换

豪斯霍尔德变换是关于超平面的反射,用于消去向量分量并构造稳定的矩阵分解。

29 个关键词5 个词条链接到这里2 个尚未撰写AI 撰写
线性映射超平面欧几里得空间正交矩阵数值线性代数QR分解单位矩阵矩阵转置豪斯霍尔德…

豪斯霍尔德变换是一种线性映射,将向量关于通过原点的超平面作反射。在实欧几里得空间中,它由一个对称的正交矩阵表示,该矩阵与单位矩阵之差是一个秩一矩阵。它在计算中的主要用途是在保持向量长度不变的同时消去指定分量,因此是数值线性代数的基本工具,尤其用于QR分解以及特征值和奇异值计算之前的矩阵约化。该变换以阿尔斯顿·S. 豪斯霍尔德命名,他在1958年的论文中介绍了如何将这种变换用于矩阵三角化。(web.stanford.edu)

定义与几何解释

对于非零向量 v∈Rnv\in\mathbb{R}^n,豪斯霍尔德矩阵定义为

H=I−2vvTvTv,H=I-2\frac{vv^T}{v^Tv},

其中,II 是单位矩阵,vTv^T 表示 vv 的矩阵转置。等价地,令单位向量 u=v/∥v∥2u=v/\|v\|_2,则有

H=I−2uuT.H=I-2uu^T.

这是矩阵算法中所用初等反射变换表示的实反射形式。(web.stanford.edu)

其几何意义可直接由公式得出。将任意向量 xx 分解为

x=x⊥+x∥,x∥=u(uTx),uTx⊥=0.x=x_\perp+x_\parallel, \qquad x_\parallel=u(u^Tx), \qquad u^Tx_\perp=0.

于是

Hx=x⊥−x∥.Hx=x_\perp-x_\parallel.

因此,HH 将平行于 vv 的分量反向,而在 vv 的正交补上保持不变。保持不动的超平面为 v⊥={x:vTx=0}v^\perp=\{x:v^Tx=0\}。用正交投影表示,即 H=I−2PH=I-2P,其中 P=uuTP=uu^T 将向量投影到由 vv 张成的直线上。(web.stanford.edu)

代数性质

直接计算可得

HT=H,H2=I,HTH=I.H^T=H, \qquad H^2=I, \qquad H^TH=I.

因此,该变换是对称的、正交的,而且其逆变换就是自身。正交性意味着它保持内积不变:

(Hx)T(Hy)=xTy.(Hx)^T(Hy)=x^Ty.

因而,它也保持欧几里得范数、夹角和距离不变。这些恒等式解释了为何反射变换适合用于变换数值问题:变换本身不会放大误差。(web.stanford.edu)

由定义公式还可得到以下结论:

  • Hv=−vHv=-v,因此 vv 是对应于特征值 −1-1 的特征向量。
  • 每个垂直于 vv 的向量都对应于特征值 +1+1。
  • 行列式为 −1-1。
  • I−HI-H 的秩为一。
  • 将 vv 替换为其任意非零标量倍数,得到的 HH 不变。

因此,当 n>1n>1 时,特征值为一个 −1-1 和重数为 n−1n-1 的 +1+1。豪斯霍尔德矩阵的乘积仍是正交矩阵,但不一定能表示为单次反射。(web.stanford.edu)

构造反射变换以消去分量

给定非零向量 x∈Rnx\in\mathbb{R}^n,一个常见目标是构造 HH,使得

Hx=αe1,Hx=\alpha e_1,

其中 e1=(1,0,…,0)Te_1=(1,0,\ldots,0)^T。保持长度不变要求 ∣α∣=∥x∥2|\alpha|=\|x\|_2。在精确算术中,取

v=x−αe1v=x-\alpha e_1

并代入豪斯霍尔德变换的定义公式,只要 v≠0v\ne0,就能实现这一目标。这样,一次运算即可将除第一个分量之外的所有分量置零。(web.stanford.edu)

在浮点运算中,通常取

α=−s∥x∥2,s={1,x1≥0,−1,x1<0.\alpha=-s\|x\|_2, \qquad s= \begin{cases} 1,&x_1\ge0,\\ -1,&x_1<0. \end{cases}

于是

v=x+s∥x∥2e1.v=x+s\|x\|_2e_1.

此时,vv 的第一个分量由两个同号的量相加得到,避免了在 xx 几乎平行于 e1e_1 时,选择相反符号可能引起的相消。实际实现还会在必要时对中间量进行缩放,以避免下溢或上溢。如果待消去的分量已经为零,软件可能采用恒等变换,而不执行反射;这是一种计算约定,并非非平凡的几何反射。(netlib.org)

示例

设

x=(34).x=\begin{pmatrix}3\\4\end{pmatrix}.

由于 ∥x∥2=5\|x\|_2=5,取 α=−5\alpha=-5 和 v=(8,4)Tv=(8,4)^T。代入可得

H=I−280(64323216)=(−35−45−4535).H= I-\frac{2}{80} \begin{pmatrix} 64&32\\ 32&16 \end{pmatrix} = \begin{pmatrix} -\frac35&-\frac45\\ -\frac45&\frac35 \end{pmatrix}.

通过乘法可验证 Hx=(−5,0)THx=(-5,0)^T:第二个分量被消去,而长度仍为 55。

高效应用与存储

通常没有必要显式构造稠密矩阵 HH。写成

H=I−τvvT,τ=2vTv,H=I-\tau vv^T, \qquad \tau=\frac{2}{v^Tv},

则它对向量的作用为

Hx=x−τv(vTx).Hx=x-\tau v(v^Tx).

这只需要一次点积和一次按比例缩放的向量更新,运算量为 O(n)O(n),而一般的稠密矩阵与向量相乘需要 O(n2)O(n^2) 的运算量。对于矩阵 AA,左乘或右乘也都可以表示为秩一更新:

HA=A−τv(vTA),AH=A−τ(Av)vT.HA=A-\tau v(v^TA), \qquad AH=A-\tau(Av)v^T.

这些都是初等反射变换表示在计算上的直接推论。(netlib.org)

软件通常会对所存储的反射向量进行归一化,使其第一个参与运算的分量为 11,并仅存储其余分量和标量 τ\tau。一系列反射变换可以隐式表示一个正交矩阵,从而无需显式构造该矩阵,就能完成与该矩阵或其转置的乘法。(netlib.org)

在分块计算中,乘积可以表示为

H1H2⋯Hb=I−VTVT,H_1H_2\cdots H_b=I-VTV^T,

其中,VV 的各列为反射向量,TT 为上三角矩阵。这种紧凑WY表示允许通过矩阵与矩阵运算来应用这些变换,从而更充分地利用高性能数值计算内核。该乘积有时称为分块反射变换,但通常不具有单个反射变换的对称性或对合性。(netlib.org)

QR分解与最小二乘问题

对于满足 m≥nm\ge n 的实 m×nm\times n 矩阵 AA,豪斯霍尔德QR分解依次消去对角线下方的元素。在第 kk 步,反射变换仅作用于第 k,…,mk,\ldots,m 行,不会破坏此前已经形成的零元素。若嵌入完整空间的反射矩阵为 H1,…,HpH_1,\ldots,H_p,则

R=Hp⋯H2H1A,Q=H1H2⋯Hp,R=H_p\cdots H_2H_1A, \qquad Q=H_1H_2\cdots H_p,

因此 A=QRA=QR,其中 QQ 为正交矩阵,RR 为上梯形矩阵。简化分解仅保留 QQ 的前 nn 列以及相应的方形三角因子。(web.stanford.edu)

它与普通最小二乘法的联系源于范数不变性。令

QTb=(c1c2),R=(R10),Q^Tb= \begin{pmatrix}c_1\\c_2\end{pmatrix}, \qquad R= \begin{pmatrix}R_1\\0\end{pmatrix},

则有

∥Ax−b∥22=∥R1x−c1∥22+∥c2∥22.\|Ax-b\|_2^2 = \|R_1x-c_1\|_2^2+\|c_2\|_2^2.

当 AA 列满秩时,通过求解三角线性方程组 R1x=c1R_1x=c_1 即可得到使目标函数最小的解。这个恒等式说明,正交三角化为何能够在不显式构造 ATAA^TA 的情况下求解最小二乘问题。(web.stanford.edu)

特征值与奇异值计算

豪斯霍尔德变换也用于在迭代谱计算之前,将矩阵约化为具有特定结构的形式:

  • 实对称矩阵可通过正交相似变换约化为三对角形式。
  • 一般方阵可约化为上海森伯格矩阵形式,即第一条次对角线以下的元素均为零,作为QR算法等算法的预处理步骤。
  • 矩形矩阵可通过左右两侧的变换约化为双对角形式,为奇异值分解做好准备。

这些约化保持相关的谱信息不变:相似变换保持特征值不变,而左右两侧的正交矩阵乘法保持奇异值不变。一般的非对称矩阵通常不能用处理对称矩阵时所采用的同一正交相似变换方法约化为三对角形式。(netlib.org)

复数情形的变换

对于复数域上的向量,直接对应的反射形式为

H=I−2vv∗v∗v,H=I-2\frac{vv^*}{v^*v},

其中,v∗v^* 为共轭转置。该矩阵既是厄米矩阵,也是酉矩阵,并满足 H2=IH^2=I。(netlib.org)

数值计算库还采用一种更广义的复数初等反射变换约定:

H=I−τvv∗,H=I-\tau vv^*,

其中 τ\tau 为复数。酉性要求

τ+τ‾=∣τ∣2(v∗v).\tau+\overline{\tau} = |\tau|^2(v^*v).

这样的 HH 不一定是厄米矩阵,也不一定具有对合性。例如,LAPACK的复数反射变换生成程序构造的 HH,使 H∗H^* 将输入向量映射为首分量是实数、其余分量均为零的向量。因此,在复数计算中,“豪斯霍尔德反射变换”可能指对单位矩阵进行秩一修正而得到的酉变换,而非严格意义上二阶的反射变换。(netlib.org)

数值稳定性与局限

正确实现的豪斯霍尔德QR分解具有很好的数值稳定性。从后向误差的角度看,计算得到的因子可解释为某个与输入矩阵相近的矩阵的分解因子;相应扰动按范数衡量很小,计算得到的 QQ 也在数值意义上保持正交。但这并不能保证每个矩阵元素,或由该分解得到的每个解,都具有很小的相对误差。(netlib.org)

正交变换也无法消除问题的病态性。在精确算术中,它们保持可逆矩阵的谱范数条件数不变。因此,变换是否准确与原问题的敏感程度是两个不同的问题。(web.stanford.edu)

豪斯霍尔德反射变换能一次消去一整组分量,因此非常适合稠密矩阵的约化。相比之下,吉文斯旋转每次仅作用于两个坐标,提供了更局部化的变换。豪斯霍尔德最初关于三角化的论文明确比较了这两种方法,并强调用反射替代一系列平面旋转能够减少算术运算量。(web.stanford.edu)

参考来源

  1. Unitary Triangularization of a Nonsymmetric Matrixweb.stanford.edu
  2. Representation of Orthogonal or Unitary Matricesnetlib.org
  3. Non-Negative Diagonals and High Performance onnetlib.org
  4. LAPACK: SRC/dlarfgp.f Source Filenetlib.org
  5. LAPACK: zlarfgnetlib.org
  6. LAPACK: larf: apply Householder reflectornetlib.org
  7. QR Factorizationnetlib.org
  8. LAPACK: clarftnetlib.org
  9. LAPACK Working Note #2: Block Algorithms for Reducing Symmetric and General Matrices to Tridiagonal, Bidiagonal and Hessenberg Formsnetlib.org
  10. LAPACK Users’ Guide: Block Algorithms for Eigenvalue Problemsnetlib.org
  11. LAPACK Working Note #176netlib.org