%% Householder-based QR computation
% Martin Stoll, 2018
% generate test matrix
m = 100;
n = 99;
A = randn(m,n);
R = A;
%% loop for QR, not efficient
Q = speye(m); % This matrix collects our orthogonal transformations
for i=1:n
    Pk=speye(m);   % The current orthogonal matrix
    v = R(i:m,i); v(1) = v(1)-norm(R(i:m,i));    % The vector for the Householde reflection
    Pk(i:m,i:m) = speye(m-i+1)-2*(v*v')/(v'*v); % Embed the matrix into the larger matrix
    R = Pk*R; % Apply the transformation
    Q = Pk*Q; % Collect all transformations
end

Q = Q'; % transpose the matrix  so everything is pretty

%% sanity check
norm(A-Q*R)