% Householder.m
% orthogonalization by means of reflections. The problem is to find the
% proper subspace in which to mirror.
clear

N=5;
A=randn(N,3);

% the first subspace is a line

v1=A(:,1);
h1=zeros(N,1);
h1(1)=sign(v1(1));	% The desired mirror image is norm(v1)*h1
s1=v1/norm(v1)+h1;	% The angle of this vector is midway between v1 and desired image
s1=s1/norm(s1);
S1=2*s1*s1'-eye(N);	% S=2P-I, P is s1*s1' as s1 is normalized so that s1'*s1=1
u1=S1*v1;		% u1 is the desired vector with zeros in all positions but the first
A1=S1*A;

% To produce zeros in a second vector without
% destroying what is already done, we need to mirror
% in a plane that contains a unity vector along the first coordinate

v2=A1(:,2);
h2=zeros(N,1);
h2(2)=sign(v2(2));		% Desired image if v2(1) were zero
temp=v2; temp(1)=0;		% To produce "v2(1)=0"
s1=temp/norm(temp)+h2; s1=s1/norm(s1);
s2=zeros(N,1); s2(1)=1;
S2=2*[s1 s2]*[s1 s2]'-eye(N);	% As s1 and s2 are orthonormal by construction
u2=S2*v2;
A2=S2*A1;

% In the next step, we must mirror in a 3D space containing the first two coordinates.

v3=A2(:,3);
h3=zeros(N,1);
h3(3)=sign(v3(3));
temp=v3; temp(1)=0; temp(2)=0;
s1=temp/norm(temp)+h3; s1=s1/norm(s1);
s2=zeros(N,1); s2(1)=1;
s3=zeros(N,1); s3(2)=1;
S3=2*[s1 s2 s3]*[s1 s2 s3]'-eye(N);
u3=S3*v3;
A3=S3*A2;

% By now, the principle is clear.

U=S3*S2*S1;	% U is an orthogonal matrix by construction, A3=U*A
		% It follows that A=U'*A3 is a qr-factorization of A
norm(A-U'*A3,'fro'), A3