%% EFFICIENT HOUSEHOLDER

%% Setup

A = [4 -5 ; 0 6 ; -4 2 ; 2 5];
b = [1 2 3 9]';

A = [4 -10 ; 0 6 ; -4 2 ; 2 3]
b = [-6 6 -2 5]';

[m, n] = size(A);

I = eye(m);
A_orig = A;
b_orig = b;

%% Step 1 (1st column)
a = A(:,1);
v = a + sign(a(1))*norm(a,2)*I(:,1);
A = A - 2*v*(v'*A)/(v'*v)
b = b - 2*v*(v'*b)/(v'*v);

%% Step 2 (2nd column)
a = [0 ; A(2:end,2)];
v = a + sign(a(2))*norm(a,2)*I(:,2);
A(:,2:end) = A(:,2:end) - 2*v*(v'*A(:,2:end))/(v'*v)
b = b - 2*v*(v'*b)/(v'*v);

%% LS
x = A(1:n,1:n)\b(1:n)
test = A_orig\b_orig


