function [R] = householder(A)
% Householder QR decomposition.
[n, m] = size(A);
for k = 1:m
    x = A(k:n, k);
    v = x;
    v(1) = x(1) + sign(x(1)) * norm(x);
    v = v / norm(v);
    A(k:n, k:m) = A(k:n, k:m) - 2 * v * (v' * A(k:n, k:m));
end
R = A;
end