%% Transformaciones de Householder paso a paso

% Matriz sacada de la seccion 5.5 en el libro de Moler. La semana que viene
% vamos a trabajar con este ejemplo.

clear
clc

s = (0:0.2:1)'; % vector 6x1
A = [s.^2, s, ones(size(s))] % matriz 6x3: estamos usando un polinomio cuadratico como modelo

sigma = norm(A(:,1)) % A(1,1) = 0, no hay que preocuparse por el signo de sigma

u = A(:,1) + [sigma; 0; 0; 0; 0; 0]

rho = 2/(u'*u)

H1 = eye(6) - rho*u*u' % primera transformacion de Householder

A1 = H1*A % primer paso

pause

clc

A2_tilde = A1(2:end,2)

sigma = norm(A2_tilde(:,1))*sign(A2_tilde(1,1))

u = A2_tilde(:,1) + [sigma; 0; 0; 0; 0]

rho = 2/(u'*u);

H2_tilde = eye(5) - rho*u*u'

H2 = [1 0 0 0 0 0; zeros(5,1) H2_tilde] % segunda transformacion de Householder

A2 = H2*A1 % segundo paso

pause

clc

A3_tilde = A2(3:end,3)

sigma = norm(A3_tilde(:,1))*sign(A3_tilde(1,1))

u = A3_tilde(:,1) + [sigma; 0; 0; 0]

rho = 2/(u'*u);

H3_tilde = eye(4) - rho*u*u'

H3 = [1 0 0 0 0 0; 0 1 0 0 0 0; zeros(4,2) H3_tilde] % tercera transformacion de Householder

A3 = H3*A2 % tercer paso: esta ya es la R que buscamos!