Householder变换

整理自:《数值线性代数（徐树方）》

Householder变换是一种能将n维向量x变换到任一n维向量y的正交变换，由于从几何上看Householder变换通过x和y之间的垂直平分面将x“反射”到y，因此Householder变换又叫镜面变换；

Householder的主要应用在于它能够将x变换成任意一个等长的若干个分量为0的向量（这种向量具有某些良好的性质，尤其是在最小二乘法的正交化解法的应用），只需要对变换后的向量再进行一次Householder变换，就能变回x；

本篇先介绍Householder变换的定义及其性质，再推导一种用于求Householder变换的数值化方法

一、Householder变换及其性质

    定义：

    Householder变换：设ω∈Rⁿ, ||ω||₂=1,定义：

                   H=I-2ωω^T(H∈R^n×n)                                                                          公式1

    称H为Householder变换(矩阵)

    性质：

    1.对称性：H^T=H

    2.正交性：H^TH=I

    3.对合性：HH=I

    4.反射性：对任意x∈Rⁿ，Hx是x关于ω的垂直超平面(即span{ω^⊥})的镜面反射。

                                                            

    性质1,2,3不难证，这里仅证性质4：

    设x∈Rⁿ，则可以将x表示为x=u+αω，其中u∈span{ω^⊥}（即ω的正交补空间），α∈Rⁿ，即有：Hx=H(u+αω)=(I-2ωω^T)u+(I-2ωω^T)αω=u-αω，得证。

    从以上证明过程可以看出，H将x沿ω的分量映射到超平面的反方向，而没有改变垂直ω（即沿超平面方向）的分量方向，因此导致x经过H变换以后变为了关于ω的垂直超平面的镜面反射，实际上，以上证明的本质可以概括为H的以下两个性质，即：Hu=u,Hω=-ω。

    （由于Householder变换的反射性，Householder变换又被称为初等反射矩阵或镜像变换）

    定理1：

            给定任何两个向量x和y(x，y∈Rⁿ且||x||₂=||y||₂)，都可以找到一个Householder变换H，使得y=Hx。

    采用构造性的方法证明：令ω=(x-y)/||x-y||₂，H=I-2ωω^T，即有y=Hx，得证。

    由定理1自然得到定理2：

    定理2：

           设0≠x∈Rⁿ， 则可构造单位向量ω∈Rⁿ，使得由公式1定义的Householder变换H满足Hx=αe₁，其中α=±||x||₂。

二、Householder算法

    正如定理2显示的那样，Householder变换的主要用途在于，它能和Guass变换一样，通过选取指定的单位向量，把一个给定向量的若干个分量置为0，Householder算法就是用来寻找满足定理2的H，即对任意一个0≠x∈Rⁿ，找到一个H满足Hx=αe₁。虽然ω和H的求法定理1已揭示出，但是Householder算法从数值计算的角度，考虑到计算误差和时间、空间复杂度的问题，对求解过程做了一定的修改，使得求解算法更加高效准确。

    接下来我们推导求解将x变为||x||₂e₁的Householder变换的算法：

    首先，从定理1和定理2可推出，ω和H的基本构造方法:

    (1) 计算v=x±||x||₂e₁     

    (2)计算ω=v/||v||₂

    (3)计算H=I-2ωω^T

    以上方法从数学的角度非常美观，但是从计算机的角度，存在着计算误差以及时间、空间复杂度的问题，下面就对其缺点以及解决方法作一定的分析，在最后贴出解决了这些问题的最终版算法。

    在步骤(1)中，通常选取v=x-||x||₂e₁，但这样在计算时可能会遇到一个问题：如果x₁为正 向量且和||x||₂大小上比较接近，计算x₁-||x||₂时，会严重地损失有效数字，甚至造成下溢。解决地方法就是对该式做一定的等价变形，即：

                   v₁=x₁-||x||₂=(x₁²-||x||₂²)/(x₁+||x||₂)=-(x₂²+x₃²+...+x_n²)/(x₁+||x||₂)       公式2

    （只有在x₁为正时才需要做这种变换，当x₁为负时x₁-||x||₂不存在精度损失的问题）

     在步骤(2)中，需要计算||v||₂，其中包含开方运算，开方运算的效率较低，要尽量避免，将(2)式直接代入(3)式，恰好可以直接避免开方运算：

                  H=I-2ωω^T=I-2vv^T/v^Tv                                                                       公式3

     整理公式3，令β=2/v^Tv，即：

                  H=I-βvv^T                                                                                           公式4

    此外，为了避免x²过大造成的上溢出问题，我们在步骤(1)之前令x=x/||x||_∞，利用规格化后的x来求β和v，这样相当于在原来的v之前乘了1/||x||_∞，注意，这样做对最终的H没有影响，因为1/||x||_∞v与v的单位化向量相同，即vv^T/v^Tv 不变。

    此外，我们可以在步骤(3)之后，令v=v/v₁，这样做可以使v₁=1，v的后n-1个分量正好存在x的后n-1个为0的向量之中（v₁=1不需保存），但是注意要将β做相应调整。

    综合以上改进，最终的算法为：

                                            

    最后补充两点：

    利用Householder变换在一个向量中引入零元素，并不局限于Hx=αe₁的形式，例如，我们需要将x的第k+1至j个元素置为0，那么可以构造y=(x₁,  x₂, ...x_k-1, α, 0, 0...0 ,x_j+1,...x_n)，（a=Σ^j_i=kx_i²）注意到||x||₂=||y||₂，

再构造v=x-y即v=(0,...0,x_k-a,x_k+1,...x_j,0,...0)即可。

    注意到算法1中，并没有直接求出H，而是给出了β和v，这样做的原因使求vv^T的成本太大，实际上，不需要求出H，而利用H=I-βvv^T来进行变换更有效率，例如对矩阵A做Householder变换，可以通过以下公式进行：

                  HA=(I-βvv^T )A=A-βv(A^Tv)²=A-vw                                                      公式5

    其中w=βA^Tv，总的运算量为4mn