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