参考地址 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
版权声明：本文为博主原创文章，转载请附上博文链接！