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);

%%
% 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

%%
% By design we can use the reflector to get the zero structure we seek:
P*z

%%
% Now we let 
A = P*A

%%
% 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
A(2:m,2:n)

%%
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);
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);

%%
% We have now reduced the original A to an upper triangular matrix
% using three orthogonal Householder reflections:
A