Householder QR Decomposition

qr.m

A = [
  -1 -1 1;
  1 3 3;
  -1 -1 5
];

[QQ,RR] = qr(A);

[m, n] = size(A);

for k=1:n
  y = A(1+k-1:end,k);

  e1 = zeros(m-(k-1), 1);
  e1(1) = 1;

  w = y + sign(y(1)) * norm(y) * e1;

  v = w / norm(w);

  if k == 1
    H = (eye(size(v,1)) - 2*(v*v'));
    A = H*A;
  else
    H = H * [
        eye(k-1) zeros(k-1,m-(k-1));
        zeros(m-(k-1),k-1) eye(size(v,1)) - 2*(v*v')
        ];
    A = [
      eye(k-1) zeros(k-1,m-(k-1));
      zeros(m-(k-1),k-1) eye(size(v,1)) - 2*(v*v')
    ]*A;
  end
end

R = A;
Q = H;