# 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


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

• You can make this a little bit more succinct, and also closer to the version that works in $\mathbb{C}^n$, by replacing the two lines that determine $\alpha$ by $\alpha = \text{sign}(v(1))*\text{norm}(v)$. This does require you to pick a nonzero value for $\text{sign}(0)$, however. – Ian Jun 22 '14 at 4:46
• @Ian Yes, that's exactly what the line there does (implements the "sign" function by taking care of a possibly zero entry $v_1$. It is a bit more concise than a bunch of "if's". In the complex case, $\alpha$ is usually given by $\alpha=-\phi(v_1)\|v\|$, where $\phi(v_1)$ denotes the phase of $v_1$ so that in forming $v_1\leftarrow v_1-\alpha$, we add the two complex numbers with the same phase and avoid cancellation in computing both real and imaginary parts (one can take $\phi(v_1)=1$ if $v_1=0$). – Algebraic Pavel Jun 23 '14 at 17:15