Backward stable projection and normalization of a vector

Question

Given a machine precision unit vector $n$, and an arbitrary vector $v$, I want an unconditionally backward stable method to compute $$f(v) = \frac{v-nn'v}{\left|v-nn'v\right|}$$ In other words, project $v$ orthogonal to $n$ and then normalize. Specifically, the result should satisfy (to machine precision if possible):

1. $\left|f(v)\right| \approx 1$
2. $n \cdot f(v) \approx 0$
3. The approximation to $f(v)$ is exactly $f(\tilde{v})$ for some $\tilde{v}$ close to v except for normalization error (backward stability).

These should hold even if $n \cdot v = 0$ (in which case $f(v)$ can be any unit vector orthogonal to $n$). The trick is to find some way of stably projecting and normalizing that doesn't lose accuracy as $n \cdot v$ approaches zero.

This formula arises in the analytic computation of extrapolated normals to cylinders. For example, here is one wrong way to do it.

Answer 1

To carry out a projection operation, the best paper out there that I know of discussing the advantages and disadvantages is:

Stewart, G. W. On the Numerical Analysis of Oblique Projectors, SIMAX (2011), Vol. 32, Issue 1, p. 309-348.

For orthogonal projectors, naïvely constructing and applying the projection matrix is fine because the condition number of that operation is the norm of the projection matrix (which, for an orthogonal projector, is 1).

The most numerically stable way of constructing such a projector that I know of is to take $n$ and perform an SVD of it. The left singular vectors corresponding to zero singular values (I suppose, technically, there are no singular values corresponding to those vectors, since there will only be one singular value, and it will be nonzero...) will be the basis for the orthogonal complement of the span of $n$; let the matrix of these singular values be $U$. Then the projector you're looking for is $P = UU^{H}$, and you can calculate the vector you want one of two ways:

1. Let $w = Pv$, and then $f(v) = w / \|w\|$, unless $w = 0$, in which case $f(v) = 0$.
2. Calculate $x = U^{H}v$, then $w = Ux$, and $f(v) = w / \|w\|$, unless $w = 0$, in which case $f(v) = 0$.

Based on Stewart's comments in the paper on oblique projectors, method 1 or method 2 should give you roughly the same amount of accuracy.

You can also calculate your projector using QR instead of SVD, but SVD is more numerically stable.

EDIT: The below method also works:

3) Calculate $x = U^{H}v$, then set $w = x / \|x\|$ (unless $x = 0$, then $w = 0$). Finally, $f(v) = Uw$.

Essentially, the argument is that $\|w\|$ and $\|f(v)\|$ are the same because $U$ is a partial isometry. $f(v)$ is guaranteed to be in the span of the columns of $U$ (to within numerical error), and division by $\|x\|$ could amplify error in this subspace, as opposed to methods 1 and 2, where division by $\|w\|$ could amplify error in all of $\mathbb{R}^{n}$, including the normal vector, $n$.

Answer 2

A backward stable algorithm is not possible

The reason is discussed in Example 15.2 (Outer Product) of Trefethen and Bau, in which they state

As a rule, for problems where the dimension of the solution space $Y$ is greater than that of the problem space $X$, backward stability is rare.

Recall the definition: we say that an algorithm $\tilde f$ for computing $f : X \to Y$ is backward stable if for each $x \in X$,

$$ \tilde f(x) = f(\tilde x) \text{ for some $\tilde x$ with } \frac{\lVert \tilde x - x \rVert}{\lVert x \rvert} \in O(\epsilon) . $$

The "expand" step in your proposed algorithm is a function from $R^2$, yielding $\tilde f \in R^3$. Since in general, $\tilde f$ will not be exactly orthogonal to $n$, there will not exist a vector $\tilde v_2 \in R^2$ such that its exact "expansion" is $\tilde f$. Of course this step lacking backward stability ruins the chances for the overall algorithm.

More generally, no algorithm can for your problem can be backward stable because this problem is rank-deficient with an image that does not embed exactly in floating point. Even with an oracle that told you the exact answer, projecting it into a floating point representation would produce an approximation lying outside the image of the exact function. You can, however, have stability (which relaxes equality in the definition).

Answer 3

In 3D, at least, the most direct way is to compute a orthogonal basis $u_0, u_1$ of vectors orthogonal to $n$, project $v$ into this basis, normalize in 2D, and expand back to 3D.

As Jed notes in his answer, the property that this algorithm satisfies is not true backward stability, but it does satisfy the property that I should have asked for, namely stability.

(I didn't remember this method until after I wrote the question.)