Householder Transformation 进行QR分解Matlab代码

Householder Transformation 进行QR分解Matlab代码

该代码结合下面的代码可实现HT求最小二乘法问题

% %=====================================================
% function  x=lshouse(A,b)
% %---------------------------- General  M-file  comment
% %  Version   :  V1.0    
% %  Coded by  :  song.yz
% %  Date      :  2012-10-29 23:30:17          *星期一*
% %-----------------------------------------------------
% %  Purpose     :  QR分解计算最小二乘问题
% %  Post Script :
% %       *.    QR分解方法采用householder变换
% %       *.
% %-----------------------------------------------------
% %  Parameters :
% %       *.  A----------系数矩阵
% %       *.  b----------右向量
% %       *.  x----------解算结果
% %----------------------------------------------------
%      [M,N]=size(A);
%          %获得维数     
%          [Q,R]=house(A);         
%          b=Q'*b;
%          
%          R1=R(1:N,1:N);
%          
%          b1=b(1:N);
%          
%          x=uptri(R1,b1);
%          
% end

%=====================================================
function  [Q,R]=houseQR(A)
%---------------------------- General  M-file  comment
%  Version   :  V1.0    
%  Coded by  :  song.yz
%  Date      :  2012-10-29 21:14:07          *星期一*
%-----------------------------------------------------
%  Purpose     :  Householder变换
%  Reference   :     
%-----------------------------------------------------
%  Parameters :
%       *.  Q-----正交矩阵
%       *.  R-----上三角矩阵
%       *.
%----------------------------------------------------

[M,N]=size(A);
%获得矩阵维数
A1=A;
H1=zeros(M,M);
for j=1:M
 H1(j,j)=1;
end 
%k表示对所有的列
for k=1:N
 %设置H矩阵初值，这里设置为单位矩阵
  H0=zeros(M,M);
  %设置为单位矩阵
  for i=1:M
  H0(i,i)=1;
  end 
  s=0;
  for i=k:M
   s=s+A1(i,k)*A1(i,k);
  end 

%算的向量的2范数  
  s=sqrt(s);
  u=zeros(N,1);
%---------------------------------
% 这段甚为重要，关系到数值稳定性问题
% 目的是使得u的2范数尽可能大
% 原则是，如果首元素大于零，则u的第一个元素是正+正
%         如果首元素小于零，则u的第一个元素是负-负
  if (A1(k,k)>=0) 
  u(k)=A1(k,k)+s;
  else 
  u(k)=A1(k,k)-s;
  end 
%-------------------------------  
  for i=k+1:M
  u(i)=A1(i,k);
  end   
  du=0;
  for i=k:M
  %求的单位u 长度平方
  du=du+u(i)*u(i);
  end 

  %计算得到大的H矩阵
  for i=k:M
    for j=k:M
        H0(i,j)=-2*u(i)*u(j)/du;            
        if i==j 
        H0(i,j)=1+H0(i,j);
        end

    end 
  end 
%千万要注意矩阵相乘的次序
%先更新矩阵
A2=H0*A1;
A1=A2;
%H1初值为单位矩阵，后逐步更新
H1=H1*H0;
%更新变换后的矩阵
end 

Q=H1;

R=A1;
end

%=====================================================
function  x=uptri(A,b)
%---------------------------- General  M-file  comment
%  Version   :  V1.0
%  Coded by  :  song.yz
%  Date      :  2012-10-29 23:36:29          *星期一*
%-----------------------------------------------------
%  Parameters :
%       *.  A---系数矩阵
%       *.  b---右向量
%----------------------------------------------------
  N=length(b);
  x(N)=b(N)/A(N,N);

  for i=N-1:-1:1
    x(i)=b(i);
        for j=i+1:N
        x(i)=x(i)-A(i,j)*x(j);
        end
        x(i)=x(i)/A(i,i);
  end
  x=x';
end