
% Demonstration of QR factorization with Householder Transforms
% H. Motteler, Feb 1999
%
% These examples are intended for transparency, rather than
% efficiency.  The semicolon at the end of a statement can
% be deleted to watch the value change with each iteration.


% The array A and vector b from Example 3.5 of the text
%
A0 = [1   -1    1;
     1   -0.5  0.25;
     1    0    0;
     1    0.5  0.25;
     1    1    1];

b0 = [1 0.5, 0, .5, 2]';


% a larger array A and vector b, for a more serious test
%
% A0 = rand(20, 18);
% b0 = rand(20, 1);

[m, n] = size(A0);
I = eye(m);


% this loop does a Householder reduction of A with explicit
% multiplication by matrices H(v), reducing A in place to R
%
A = A0;
b = b0;
Qt = eye(m);
for i = 1:n

  % calculate the zapping vector for column i
  e = I(:,i);
  a = A(:,i);
  a(1:i-1) = zeros(i-1, 1);
  v = a + sign(a(i)) * norm(a) * e 

  % explicit calculation and application of H
  H = I - 2 * (v*v') / (v'*v) 
  A = H * A 
  b = H * b ;

  % explicit calculation of Qt
  Qt = H * Qt;

end

% check residuals
%
R = A;
Q = Qt';
norm(A0 - Q*R)

% check against the Matlab QR factorization routine
%
% [Q0, R0] = qr(A0);
% norm(R0 - R)
% norm(Q0 - Q)
% norm(Q0'*b0 - b)


% This loop does a Householder reduction of A with implicit
% multiplication by matrices H(v), reducing A in place to R.
%
% Though faster in general, the inner loop on columns used to
% perform the implicit H(v)*A multiplication makes this version
% slower in Matlab.
%
A = A0;
b = b0;
Qt = eye(m);
for i = 1:n

  % calculate the zapping vector for column i
  e = I(:,i);
  a = A(:,i);
  a(1:i-1) = zeros(i-1, 1);
  v = a + sign(a(i)) * norm(a) * e ;

  % implicit calculation of H*A 
  for j = 1:n
    A(:,j) = A(:,j) - (2*v'*A(:,j) / (v'*v)) * v;
  end

  % implicit calculation of H*b
  b = b - (2*v'*b / (v'*v)) * v;

  % implicit calculation of Qt = Hm* ... *H2*H1
  for j = 1:m
    Qt(:,j) = Qt(:,j) - (2*v'*Qt(:,j) / (v'*v)) * v;
  end

end

% check residuals
%
R = A;
Q = Qt';
norm(A0 - Q*R)

% check against the Matlab QR factorization routine
%
% [Q0, R0] = qr(A0);
% norm(R0 - R)
% norm(Q0 - Q)
% norm(Q0'*b0 - b)

