function [Q,R] = slow_householder_qr(A)
% Compute the QR factorization of real A using Householder reflectors
% ** implementation designed for transparency, not efficiency **
 
  [m,n] = size(A);
  Q = eye(m);
  for k=1:min(m-1,n)
     ak  = A(k:end,k);                                  % vector to be zeroed out
     vk  = ak + sign(ak(1))*norm(ak)*[1;zeros(m-k,1)];  % construct v_k that defines the reflector
     Hk = eye(m-k+1) - 2*vk*vk'/(vk'*vk);               % construct reflector
     Qk = [eye(k-1) zeros(k-1,m-k+1); zeros(m-k+1,k-1) Hk];
     A = Qk*A;                                          % update A
     Q = Q*Qk;                                          % accumulate Q
  end
  R = A;