3
$\begingroup$

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
$\endgroup$
4
$\begingroup$

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
$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$ – Ian Jun 22 '14 at 4:46
  • $\begingroup$ @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$). $\endgroup$ – Algebraic Pavel Jun 23 '14 at 17:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.