%% EFFICIENT HOUSEHOLDER

% In loop

%% Setup

A = [4 -5 ; 0 6 ; -4 2 ; 2 5];
b = [1 2 3 9]';
[m, n] = size(A);

%% Solve
for k = 1:n
    a = [zeros(k-1,1) ; A(k:end,k)];
    v = a + sign(a(k))*norm(a,2)*I(:,k);
    A(:,k:end) = A(:,k:end) - 2*v*(v'*A(:,k:end))/(v'*v);
    b = b - 2*v*(v'*b)/(v'*v);
end
x = A(1:n,1:n)\b(1:n)
