Householder Vector algorithm in Golub and Van Loan

Question

In the 4th edition of "Matrix Computations", Golub and Van Loan present "Algorithm 5.1.1 (Householder Vector)". The first couple of lines (translated into MATLAB-syntax) read:

m = length(x); sigma = x(2:m)'*x(2:m); v = [1; x(2:m)];
if sigma == 0 && x(1) >= 0
  beta = 0;
elseif sigma == 0 && x(1) < 0
  beta = -2;
else
  ...

The else clause handles the case where sigma is nonzero and no code after the provided snippet modifies v if sigma is zero. The matrix form of the resulting Householder transformation is $I-\beta v v^T$.

The elseif clause is a little strange. It doesn't appear in the algorithm listing in the 3rd edition, so it was added for the 4th, presumably for numerical stability. However, it seems to me to generate a $(\beta, v)$ pair that doesn't map to an orthogonal matrix. For example, if x = [-1; 0; 0] then sigma == 0 and x(1) < 0 so we get beta = -2 and v = [1; 0; 0] and a Householder transformation of [3, 0, 0; 0, 1, 0; 0, 0, 1] which is not an orthogonal matrix.

So my questions are:

• What benefit is there to handling that case separately, rather than just setting beta = 0 if sigma = 0?
• Does the resulting Householder transformation need to be applied differently?

Answer 1

In this answer I'll try to provide further context and summarize the consensus I think we achieved in the comments to the accepted answer.

The full algorithm (5.1.1) reads:

m = length(x); sigma = x(2:m)'*x(2:m); v = [1; x(2:m)];
if sigma == 0 && x(1) >= 0
  beta = 0;
elseif sigma == 0 && x(1) < 0
  beta = -2;
else
  mu = sqrt(x(1)*x(1)+sigma);
  if (x(1) <= 0)
    v(1) = x(1) - mu;
  else
    v(1) = -sigma/(x(1)+mu);
  end
  beta = 2*v(1)*v(1)/(sigma+v(1)*v(1));
  v = v/v(1);
end

and produces a v ($v$) and a beta ($\beta$) intended to be converted to a matrix via $P=I-\beta v v^T$ or applied directly as $x_r = x - \beta v (v^T x)$. That these are subtractions rather than additions is confirmed by the main else case.

The line beta = -2; is incorrect and should read beta = 2;.

Any householder algorithm will need to handle the case when $\sigma=0$ and $x_1=0$ because this implies that the norm of the supplied vector $x$ is exactly zero and so a reflection is not defined. We wish to output $\beta=0$ in this case. While not strictly a householder transformation, when converted to a matrix it simply yields the identity matrix. In this algorithm that is handled as part of the if sigma == 0 && x(1) >= 0 branch.

Note that if $\sigma = 0$ then all the elements of the supplied vector below the first are zero and no transformation need be applied to zero them out. However, there are valid Householder transformations what affect the first element, leaving the rest as zeros, that could be applied if desirable.

The main else clause is designed to avoid severe cancellation when $x$ is close to a positive multiple of $e_1$ but would have division by zero error if $\sigma=0$ and $x_1 > 0$. So the case $\sigma=0$ and $x_1 = 0$ is extended to $\sigma=0$ and $x_1 >= 0$ and $\beta = 0$ returned in that case.

In previous editions, the case $\sigma=0$ and $x_1 < 0$ was handled the same way but in the 4th edition the authors choose to handle it as a special case. In that case the matrix form of the Householder transformation is the identity matrix but with a -1 in the top left position. This leaves the zeros unchanged but flips the sign of the element on the diagonal.

This does not appear to be for reasons of numerical stability but rather to create more positive entries on the diagonal of the resulting $R$ matrix. It does not appear to create numerical instability however.

Answer 2

You are correct, this is another error in this algorithm (complete code in a github archive, are there more official online sources?). It is clearly stated in the header comments that

% v is a column m-vector with v(1) = 1 and beta is a scalar such
% that P = I - beta*v*v' is orthogonal and Px = norm(x)*e1.

and in the general case beta is computed as

    beta = 2*v(1)^2/(sigma+v(1)^2);
    v = v/v(1);

according to the documented intent. When testing this there apparently were no tests for trivial cases. On the other hand, getting sigma == 0 in the usual situation of a call from the middle of some numerical algorithm has probability zero.

---

The more serious error is in the algorithm design, in the choice of the construction of the reflection normal v as a continuous vector field on the punctured sphere, excepting the "north" pole, as done with the first special case.

The question of course remains if that makes sense at all. Comparing a floating point value with $0$ in a numerical algorithm smells of bad design. If one has to avoid a singularity in one point, then it is to be expected that results close to that singularity are ill-behaved, with a dramatic loss in accuracy or the magnification of floating-point noise to the level of the intended results.

The stated aim is to get the first component of the reflected vector to be positive. The general construction is to find one of the two bisectors of the lines along the given vectors $\vec x$, and some unit vector $\vec e$, typically $\vec e=\vec e_1=(1,0,…,0)^T$, $$ \vec v = \frac{\pm\|\vec x\|\vec e+\vec x}{\bigl\|\pm\|\vec x\|\vec e+\vec x\bigr\|} =\frac{\pm\|\vec x\|\vec e+\vec x}{\sqrt{2\|\vec x\|(\|\vec x\|\pm\vec e^T\vec x)}}. $$ The recommendation for numerical stability is to take the variant that makes the denominator largest. This switches the variants at the equator (normal plane to $\vec e$), this switch is a jump by $90°$ in the direction.

With $\vec e=\vec e_1$, the reflected vector has as first component \begin{align} x_{1,\rm reflect} &=x_1-2\frac{(\pm\|\vec x\|+x_1)\|\vec x\|(\pm x_1+\|\vec x\|)}{2\|\vec x\|(\|\vec x\|\pm x_1)} \\&=\mp \|\vec x\| \end{align}

---

So to always end up with a positive number one would have to chose the lower sign variant, with the special case where that does not work. This can even be done in a numerically stable way when $\vec x\approx \|\vec x\|\vec e_1$, as by binomial extension $$ x_1-\|\vec x\|=\frac{-(\|\vec x\|^2-x_1^2)}{\|\vec x\|+x_1} =-\frac{\|\vec x_{2:n}\|^2}{\|\vec x\|+x_1} $$

Why would that be bad? Explore the reflector construction for $\vec x=\pmatrix{1\\\bar x}$ with $\|\bar x\|\sim\mu=10^{-16}$ some quasi-random accumulation of rounding errors from previous computational steps. Then using the "bad choice" the reflector unit normal will be close to $\pmatrix{-\|\bar x\|/2\\\bar x/\|\bar x\|}$. Thus a small error in the input leads to a large difference in the result, a random point on the equator. This kind of unpredictability, "catastrophic" error magnification, is usually unwanted in stable numerical algorithms.

---

The "good choice" would always result in a reflector normal close to the north unit vector $\pmatrix{1\\0}$, this kind of dependency is what is actually desired.

Now one could ask what about points on the equator under the "good choice", then both the midpoint to the north pole and the midpoint to the south pole are equally good. While it is true that under slight perturbations one of the two is selected quasi randomly, the result still stays close to one of these two possibilities, this is still numerically stable under a slight extension of this idea.