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