QR分解のためのハウスホルダー変換Matlabコード

Householder Transformation

このコードを次のコードと組み合わせると、HT最小二乗問題を実現できます。

% %===================================================== % function x=lshouse(A,b) % %---------------------------- General M-file comment % % Version : V1.0 % % Coded by : song.yz %% Date: 2012-10-29 23:30:17 *Monday* % %----------------------------------------------------- %% Purpose: QR decomposition calculation least squares problem % % Post Script : %% *. QR decomposition method adopts householder transformation % % *. % %----------------------------------------------------- % % Parameters : %% *. A----------Coefficient matrix %% *. b----------right vector %% *. x----------Solution result % %---------------------------------------------------- % [M,N]=size(A) % %Get the dimension % [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 *Monday* %----------------------------------------------------- % Purpose: Householder transformation % Reference : %----------------------------------------------------- % Parameters : % *. Q-----orthogonal matrix % *. R-----Upper triangular matrix % *. %---------------------------------------------------- [M,N]=size(A) % Get matrix dimension A1=A H1=zeros(M,M) for j=1:M H1(j,j)=1 end %k means for all columns for k=1:N %Set the initial value of the H matrix, here is set to the identity matrix H0=zeros(M,M) % Is set to the identity matrix 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 Norm of the vector calculated s=sqrt(s) u=zeros(N,1) %--------------------------------- % This paragraph is very important and is related to numerical stability % The purpose is to make the 2 norm of u as large as possible The principle of% is that if the first element is greater than zero, the first element of u is positive + positive % If the first element is less than zero, the first element of u is negative-negative 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 % Unit u length squared du=du+u(i)*u(i) end %Calculate the large H matrix 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 %Be sure to pay attention to the order of matrix multiplication % Update matrix first A2=H0*A1 A1=A2 The initial value of %H1 is the identity matrix, which will be updated gradually H1=H1*H0 %Update the transformed matrix 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 *Monday* %----------------------------------------------------- % Parameters : % *. A---coefficient matrix % *. b---right vector %---------------------------------------------------- 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