参考地址 https://www.cnblogs.com/reasno/p/9606224.html
Householder变换是一种能将n维向量x变换到任一n维向量y的正交变换,由于从几何上看Householder变换通过x和y之间的垂直平分面将x“反射”到y,因此Householder变换又叫镜面变换;
Householder的主要应用在于它能够将x变换成任意一个等长的若干个分量为0的向量(这种向量具有某些良好的性质,尤其是在最小二乘法的正交化解法的应用),只需要对变换后的向量再进行一次Householder变换,就能变回x;
本篇先介绍Householder变换的定义及其性质,再推导一种用于求Householder变换的数值化方法
一、Householder变换及其性质
定义:
Householder变换:设ω∈Rn, ||ω||2=1,定义:
H=I-2ωωT(H∈Rn×n) 公式1
称H为Householder变换(矩阵)
性质:
1.对称性:HT=H
2.正交性:HTH=I
3.对合性:HH=I
4.反射性:对任意x∈Rn,Hx是x关于ω的垂直超平面(即span{ω⊥})的镜面反射。
![]() | |
---|---|
性质1,2,3不难证,这里仅证性质4:
设x∈Rn,则可以将x表示为x=u+αω,其中u∈span{ω⊥}(即ω的正交补空间),α∈Rn,即有: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∈Rn且||x||2=||y||2),都可以找到一个Householder变换H,使得y=Hx。
采用构造性的方法证明:令ω=(x-y)/||x-y||2,H=I-2ωωT,即有y=Hx,得证。
由定理1自然得到定理2:
定理2:
设0≠x∈Rn, 则可构造单位向量ω∈Rn,使得由公式1定义的Householder变换H满足Hx=αe1,其中α=±||x||2。
二、Householder算法
正如定理2显示的那样,Householder变换的主要用途在于,它能和Guass变换一样,通过选取指定的单位向量,把一个给定向量的若干个分量置为0,Householder算法就是用来寻找满足定理2的H,即对任意一个0≠x∈Rn,找到一个H满足Hx=αe1。虽然ω和H的求法定理1已揭示出,但是Householder算法从数值计算的角度,考虑到计算误差和时间、空间复杂度的问题,对求解过程做了一定的修改,使得求解算法更加高效准确。
接下来我们推导求解将x变为||x||2e1的Householder变换的算法:
首先,从定理1和定理2可推出,ω和H的基本构造方法:
(1) 计算v=x±||x||2e1
(2)计算ω=v/||v||2
(3)计算H=I-2ωωT
以上方法从数学的角度非常美观,但是从计算机的角度,存在着计算误差以及时间、空间复杂度的问题,下面就对其缺点以及解决方法作一定的分析,在最后贴出解决了这些问题的最终版算法。
在步骤(1)中,通常选取v=x-||x||2e1,但这样在计算时可能会遇到一个问题:如果x1为正相反且和||x||2大小上比较接近,计算x1-||x||2时,会严重地损失有效数字,甚至造成下溢。解决地方法就是对改式做一定的等价变形,即:
v1=x1-||x||2=(x12-||x||22)/(x1+||x||2)=-(x22+x32+…+xn2)/(x1+||x||2) 公式2
(只有在x1为正时才需要做这种变换,当x1为负时x1-||x||2不存在精度损失的问题)
在步骤(2)中,需要计算||v||2,其中包含开方运算,开方运算的效率较低,要尽量避免,将(2)式直接代入(3)式,恰好可以直接避免开方运算:
H=I-2ωωT=I-2vvT/vTv 公式3
整理公式3,令β=2/vTv,即:
H=I-βvvT 公式4
此外,为了避免x2过大造成的上溢出问题,我们在步骤(1)之前令x=x/||x||∞,利用规格化后的x来求β和v,这样相当于在原来的v之前乘了1/||x||∞,注意,这样做对最终的H没有影响,因为1/||x||∞v与v的单位化向量相同,即vvT/vTv 不变。
此外,我们可以在步骤(3)之后,令v=v/v1,这样做可以使v1=1,v的后n-1个分量正好存在x的后n-1个为0的向量之中(v1=1不需保存),但是注意要将β做相应调整。
综合以上改进,最终的算法为:
最后补充两点:
利用Householder变换在一个向量中引入零元素,并不局限于Hx=αe1的形式,例如,我们需要将x的第k+1至j个元素置为0,那么可以构造y=(x1, x2, …xk-1, α, 0, 0…0 ,xj+1,…xn),(a=Σji=kxi2)注意到||x||2=||y||2,
再构造v=x-y即v=(0,…0,xk-a,xk+1,…xj,0,…0)即可。
注意到算法1中,并没有直接求出H,而是给出了β和v,这样做的原因使求vvT的成本太大,实际上,不需要求出H,而利用H=I-βvvT来进行变换更有效率,例如对矩阵A做Householder变换,可以通过以下公式进行:
HA=(I-βvvT )A=A-βv(ATv)2=A-vw 公式5
其中w=βATv,总的运算量为4mn
使用matlab,基于householder变化写了QR的实现过程
1、Householder变化算法
function [ H, v, beta ] = householder( x )
% x : inout param. x is a vector which size is n*1
% v and beta : is param which construct H matrix
% H is hoseholder Matrix. H = I - beta*v*v'
%derive parameter v and beta
v = zeros(size(x));
beta = zeros(size(x));
%Houserholder Algorithm begin
x_len = length(x);
x_max = max(abs(x));
x = x./x_max;
zgama = x(2:end)'*x(2:end);
v(1) = 1;
v(2:end) = x(2:end);
if zgama == 0
beta = 0;
else
alpha = sqrt( x(1)^2 + zgama);
if x(1) <=0
v(1) = x(1) - alpha;
else
v(1) = -zgama./( x(1) + alpha);
end
beta = 2*v(1)^2./( zgama + v(1)^2 );
v = v./v(1);
end
%beta = 2./(v'*v);
H = eye(x_len,x_len) - beta*v*v';
end
QR分解
A = rand(300,20);
tic
A_hang = size(A,1);
A_lie = size(A,2);
H = cell( A_lie, 1 );
Hs = eye(A_hang, A_hang);
R2 = A;
Q2 = eye(A_hang,A_hang);
for i = 1:A_lie
[Hi, vi, betai] = householder(R2(i:end,i));
%H{i} = blkdiag(eye(i-1), Hi);
R2(i:end,i:end) = Hi*R2(i:end,i:end);
Q2 = Q2*blkdiag(eye(i-1), Hi);
end
Q2(find(abs(Q2) < 1e-10)) = 0;
R2(find(abs(R2) < 1e-10)) = 0;
toc
[Q1, R1] = qr(A); %matlab内部算法(非常非常的快)
Q1*R1 - A
Q2*R2 - A
erroQ = Q1 + Q2
erroR = R1 + R2
原文:https://blog.csdn.net/xiaoxiao133/article/details/80066972
版权声明:本文为博主原创文章,转载请附上博文链接!