Householder变换计算最小二乘MATLAB源代码

建不了的模。

H变换计算最小二乘是最小二乘算法中最重要的方法之一，由于论坛里不好贴公式。这里只能给出源代码。
关于其算法部分，可以参考《C#科学计算讲义》，里面有算法详细的步骤。

代码如下：
%=====================================================
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]=house(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
复制代码