function [Q,R]=householderQR(A)
% Householder method for QR factorization
% A is n x k matrix, k <= n full rank
% Q is n x n orthogonal
% R is n x k upper triangular
% A = Q * R

[n,k] = size(A); R = A; Q = eye(n);
for p = 1:k
    x = R(p:n,p);
    if x(1) ~= 0
        x(1) = x(1)+sign(x(1))*norm(x);
        x = x./norm(x);
    else
        x(1) = x(1)+norm(x);
        x = x./norm(x);
    end
    R(p:n,p:k) = R(p:n,p:k)-2*x*(x'*R(p:n,p:k));
    Q(p:n,1:n) = Q(p:n,1:n)-2*x*(x'*Q(p:n,1:n));
end
Q = Q';
