Reducing a matrix to upper Hessenberg form using Householder transformations in Matlab

Question

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

Answer 1

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