clear
format compact

%%
% We will use Householder reflections to produce a QR factorization of the
% matrix
A = magic(7);
A = A(:,1:3);
[m, n] = size(A);
Aorig = A;
Q = eye(m);

%%
% Our first step is to introduce zeros below the diagonal in column 1.
% Define the vector 
z = A(:,1);

%%
% Applying the Householder definitions gives us
e1 = eye(m, 1);
v = norm(z)*e1 - z;
P = eye(m) - 2*(v*v')/(v'*v);   % reflector

%%
% Now we let 
A = P*A;
Q = P*Q;

% We are set to put zeros into column 2. We must not use row 1 in any way,
% lest it destroy the zeros we just introduced. To put it another way, we
% can repeat the process we just did on the smaller submatrix
%
j = 2;   % 'loop' counter
z = A(j:m,j);
e1 = eye(m-j+1,1);
v = norm(z)*e1 - z;
P = eye(m-j+1) - 2*(v*v')/(v'*v);   % reflector

%%
% We now apply the reflector to the submatrix.
A(j:m,j:n) = P*A(j:m,j:n);
Q(j:m,:) = P*Q(j:m,:);
A;

%%
% We need one more iterations of this process.
j = 3;  % 'loop' counter
z = A(j:m,j);
e1 = eye(m-j+1,1);
v = norm(z)*e1 - z;
P = eye(m-j+1) - 2*(v*v')/(v'*v);   % reflector
A(j:m,j:n) = P*A(j:m,j:n);
Q(j:m,:) = P*Q(j:m,:);

%%
% We have now reduced the original A to an upper triangular matrix
% using four orthogonal Householder reflections:
A
R = triu(A)   % enforce exact triangularity
Q = Q'
norm(Aorig - Q*R)