Reducing a matrix to upper Hessenberg form using Householder transformations in Matlab I have the below Matlab code based on what my professor gave me in class.  The last line of this code is giving me an incompatible dimensions error.  When I checked the size, I got that the matrix $Q2D$ is $(n-1) \times n$, but that $vv^T$ is $n \times n$.  What should I be doing to modify the last line so the algorithm works?  It should be equivalent to calculating $Q_nQ_{n-1}...Q_1AQ^T_1Q^T_2...Q_T^N$.  Thank you!
for k=1:n
    x = Q2D(k:n,k);
    e = zeros(n-k+1,1);     % Initialize e as the standard basis vector e_1
    e(1) = 0;

   if sign(x(1)) == 0       % Calculate v
     v = norm(x)*e + x;
   else
     v = sign(x(1))*norm(x)*e + x;
   end

   v = v/norm(v);          % Orthonormalize v

   Q2D(k:n,k:n) = Q2D(k:n,k:n) - 2*v*(v.'*Q2D(k:n,k:n));

   % size(v*v.')
   % size(Q2D(1:n,k+1:n))

   Q2D(1:n,k+1:n) = Q2D(1:n,k+1:n)-2*(Q2D(1:n,k+1:n)*v)*v.';  % This line is where 
                                                              % I'm having issues
end

 A: Step by step:
1) Transforming a matrix to the upper Hessenberg form means we want to introduce some zeros in the columns $1,\ldots,n-2$. So why a loop over $1,\ldots,n$? Replace
for k = 1 : n

by
for k = 1 : n - 2

2) Transforming a matrix to the upper Hessenberg form also means that we want to zero out components $k+2,\ldots, n$ in the given column $k$ (and not $k+1,\ldots,n$). Also, there is no need to introduce two vectors. Hence replace
x = Q2D(k:n, k);

by
v = Q2D(k+1:n, k);

3) I do not understand the meaning of
e = zeros(n-k+1,1);
e(1) = 0;

It results in the zero vector.
4) There is no need to use any vector $e$ as it contains only one nonzero component and adding/subtracting it to/from a vector does mostly nothing except changing one single component. Replace
if sign(x(1)) == 0
    v = norm(x)*e + x;
else
    v = sign(x(1))*norm(x)*e + x;
end

by
alpha = -norm(v);
if (v(1) < 0) alpha = -alpha; end
v(1) = v(1) - alpha;

Note that $\alpha$ is the number we should obtain in the entry $(k+1,k)$.
5) Normalisation is OK. Leave 
v = v / norm(v);

untouched :-).
6) Replace
Q2D(k:n,k:n) = Q2D(k:n,k:n) - 2 * v * (v.' * Q2D(k:n,k:n));

by
Q2D(k+1:n,k:n) = Q2D(k+1:n,k:n) - 2 * v * (v.' * Q2D(k+1:n,k:n));

According to point (2), we do not triangularise the matrix!
Also, you can use the fact that actually you know what exactly should be in the column $k$ and use this instead:
Q2D(k+1:n,k+1:n) = Q2D(k+1:n,k+1:n) - 2 * v * (v.' * Q2D(k+1:n,k+1:n));
Q2D(k+1,k) = alpha;
Q2D(k+2:n,k) = 0;

This in particular avoids creating tiny nonzeros (due to roundoff) in entries which should be exactly zero.
7) The line
Q2D(1:n,k+1:n) = Q2D(1:n,k+1:n)-2 * (Q2D(1:n,k+1:n) * v) * v.';  % This line is where 
                                                                 % I'm having issues

is actually correct!
So this is what remains at the end:
for k = 1 : n - 2
    v = Q2D(k+1:n,k);
    alpha = -norm(v);
    if (v(1) < 0) alpha = -alpha; end
    v(1) = v(1) - alpha; 
    v = v / norm(v);
    Q2D(k+1:n,k+1:n) = Q2D(k+1:n,k+1:n) - 2 * v * (v.' * Q2D(k+1:n,k+1:n));
    Q2D(k+1,k) = alpha;
    Q2D(k+2:n,k) = 0;
    Q2D(1:n,k+1:n) = Q2D(1:n,k+1:n) - 2 * (Q2D(1:n,k+1:n) * v) * v.';
end

