# Householder Vector algorithm in Golub and Van Loan

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?

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.

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.

• It was observed that in constructions such as the pseudo-inverse, the error in pinv(pinv(A))-A increased by a factor of $4$ (or more?). Extreme cases like x = [1e+10, 1e-10] gave grave errors in the first Lapack versions using this method (v 3.2, 2009), but that was corrected with the binomial trick some time later (and the universal use of LARFP was rolled back, making ...P versions of the routines using it). Commented Dec 7, 2022 at 11:33

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.

• In Golub and Van Loan they define their Householder transformations as $P = I-\beta v v^T$ and $\beta = 2/(v^T v)$ with no requirement that $v$ has unit norm. In your notation, that would be $\vec x_{\rm reflect}=\vec x-\beta\vec v(\vec v^T\vec x)$. For convenience of storage, they scale $v$ such that its first element is 1, so the norm of $v$ won't be 1 or $\sqrt 2$ but something that need to be tracked separately in $\beta$. Commented Dec 1, 2022 at 19:22
• That is just a question of where you store the normalization factor, on the theoretical side it makes no difference for the argument. It allows the trick to set $β=0$ for the singular case, and that is not a good idea. Commented Dec 1, 2022 at 22:56
• It is valid in exact arithmetic, my complaints are about the floating point case. What happens if `sigma´ is very small but not zero. I'm suspicious of the idea that made that trick for the special case necessary. Commented Dec 2, 2022 at 0:45
• Do we agree or disagree on the question of whether the provided code produces a valid Householder transformation? It seems to me that the output of the function (the pair $(v, \beta)$) is valid exactly when $\beta = 2/(v^T v)$ or if $\beta=0$. I'm claiming that, in the noted case, the sign of $\beta$ is wrong. In that case $v=e_1$ and so $v^T v = 1$, hence $\beta=2$. But the code has $\beta=-2$ in that case. Commented Dec 2, 2022 at 0:52
• I have found a source for the complete code, your question is about an actual error (in a zero-probability event). Adapted the answer to this new (for me) information. I had previously only seen the Lapack Fortran implementation of this variant. Commented Dec 4, 2022 at 18:15