% The following codes demonstrate how to form Householder transformations
% and how to use them to obtain a QR decomposition of a matrix. The codes
% are written in sections, each of which can be run independently. To run
% the codes, start MATLAB in the directory where you have saved this file. 
% Then type "edit householder_examples.m" to bring the file up in the 
% edit window. To run a section, first click somewhere within the section
% to highlight it; then click on "Run Section" in the menu bar above. 

%% Example 1: Householder transformation #1
% This demonstrates constructing a Householder transformation that takes 
% a given vector x to a given vector y. Note that x and y must have the 
% same (Euclidean) norm for this to be possible. 

% Create vectors x and y having the same norm. 
n = 4;
x = (1:n)'
y = (n:-1:1)'
% If x and y don't already have the same norm, uncomment the next line 
% to achieve this. 
% y = (norm(y)\norm(x))*y;

% Form the unit vector u for the Householder transformation H - 2uu^T. 
u = norm(y-x)\(y-x)

% The Housholder transformation is H = I - 2(uu^T). However, there is 
% rarely a need to actually form it. Instead, just store u and, to 
% apply H to a vector v, form 2(u'v) and then Hv = v - [2(u'v)]u. 

% Follow this procedure to check whether this u does the job. 

fprintf('\n Checking.\n');

Hx = x - (2*(u'*x))*u

fprintf('\n norm(y-Hx) = %g\n', norm(y-Hx));

%% Example 2: Householder transformation #2
% This demonstrates constructing a Housholder transformation that takes 
% a given vector x to +/- ||x|| e_1 using the efficient, numerically 
% sound procedure given in class. This uses the formula H = I - c*w*w', 
% where w = x + sign(x_1)||x|| need not have norm 1 and 
% c = (||w||^2)\2) = (||x||(||x||+|x_1|))\1. 

% Create x. 
n = 4;
x = (1:n)'

% Create the constant and (non-normalized) w for the Householder 
% transformation. 
xnorm = norm(x)
c = (xnorm*(xnorm+abs(x(1))))\1;
w = x; w(1) = w(1)+sign(x(1))*xnorm;

% Check to make sure this does the job. 
Hx = x - c*(w'*x)*w


%% Example 3: Householder QR
% This demonstrates performing QR decomposition on an mxn matrix A  
% using Householder transformations. 
% Notes: 
% * To keep things simple, this code does NOT do column pivoting. 
% * This code is meant to illustrate the process; it is not meant 
%   to be efficient. 
% * There is a "pause" statement in the loop that is currently set for 
%   1 second. You can be change to different pause lengths or remove it
%   altogether, as you wish. 

% Choose m and n with m >= n. 
m = 6; n = 3;

% Generate a random mxn matrix (with randn, each entry is normally
% distributed with mean zero and variance one). 
A = randn(m,n)

% Compute the Householder vectors w_1, ..., w_(n) and constants 
% c_1,...,c_(n), storing the vectors as columns of a matrix W and 
% the constants as components of a vector c. Concurrently compute the 
% R factor. 

% Initialize R = A; preallocate storage for c and W (always a good idea). 

R = A; c = zeros(n,1); W = zeros(m,n); 

fprintf('\n Reduce A to a triangular matrix R. Initally, set R = A.\n');
%fprintf('\n Initally, set R = A.\n');

% Generate the Householder vectors and reduce R to right-triangular form. 
for k = 1:n
    % Form the kth Householder vector. 
    % Note that the first (k-1) components are zero.
    rnorm = norm(R(k:m,k));
    c(k) = (rnorm*(rnorm + abs(R(k,k))))\1;
    W(k:m,k) = R(k:m,k);
    W(k,k) = W(k,k) + sign(R(k,k))*rnorm;
    % Update the kth column of R. 
    R(k:m,k) = [-sign(R(k,k))*rnorm; zeros(m-k,1)];
    % Update columns (k+1) to m of R. 
    for j = k+1:n
        R(k:m,j) = R(k:m,j) - (c(k)*(W(k:m,k)'*R(k:m,j)))*W(k:m,k);
    end
    fprintf('\n After step %d:\n', k);
    R 
    %fprintf('Strike any key to continue.\n');
    pause(1)
end

% Form Q for illustration. There is no need to do this in practice. 
Q = eye(m,m);
for k = n:-1:1
    Q = Q - W(:,k)*(c(k)*(W(:,k)'*Q));
end

fprintf('\n Final results:\n');
Q, R

QR = Q*R;

fprintf('\n norm(A-QR) = %g\n',norm(A-QR));





