%
% this matlab program tests stability of two versions of the 
% householder transformations.
%
% written by Ming Gu for Math 221, Fall 2007
%
% for random vectors, they work equally well.
%
n = 10;
x = randn(n,1);
u1= house(x);
u2= house2(x);
h1= x - 2*u1*(u1'*x);
h2= x - 2*u2*(u2'*x);
ff = ['Random vector of length ', num2str(n)];
disp(ff);
disp([norm(h1(2:end)), norm(h2(2:end))]);
%
% In pathological cases, things can go wrong badly. 
%
x = [rand; randn(n-1,1)*sqrt(eps/2)];
u1= house(x);
u2= house2(x);
h1= x - 2*u1*(u1'*x);
h2= x - 2*u2*(u2'*x);
ff = ['Pathological vector of length ', num2str(n)];
disp(ff);
disp([norm(h1(2:end)), norm(h2(2:end))]);

function u = house(x)
%
% compute householder transformation
% this is a stable version.
% 
% written by Ming Gu for Math 221, Fall 2007
%
x  = x(:);
n  = length(x);
if (norm(x) == 0)
    u = ones(n,1)/sqrt(n);
    return;
end
x(1) = x(1) +sign(x(1))*norm(x);
u    = x/norm(x);

function u = house2(x)
%
% compute householder transformation
% this is an unstable version.
%
% written by Ming Gu for Math 221, Fall 2007
%
x  = x(:);
n  = length(x);
if (norm(x) == 0)
    u = ones(n,1)/sqrt(n);
    return;
end
x(1) = x(1) -norm(x);
u    = x/norm(x);