# Backward stable projection and normalization of a vector

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.

-
By your definition, $f(v)$ is nor a continuous function if you consider vectors $v_j\rightarrow v^\ast$ so that $v_j \cdot n \neq 0$ but $v^\ast\cdot n=0$. It's not quite clear to me how you intend to get a stable method out of this. –  Wolfgang Bangerth Jan 16 at 14:57
Lack of continuity makes forward stability impossible, but backward stability is fine. –  Geoffrey Irving Jan 16 at 17:46

### 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).

-
True. To be specific, what I want is the weak combination of forward and backward stability: find $y$ s.t. $y \approx f(\tilde{x})$ and $\tilde{x} \approx x$. Calling this just "stability" seems too ambiguous; is there a dedicated name for it? –  Geoffrey Irving Jan 17 at 4:30
That is exactly the definition given in Trefethen and Bau (Equations 14.3 and 14.4). –  Jed Brown Jan 17 at 4:41
Via wikipedia, "mixed stability" also works. –  Geoffrey Irving Jan 17 at 4:48

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$.

-
The projection involved is orthogonal, so I don't think the oblique projector paper is necessary. Also, both (1) and (2) suffer catastrophic precision loss as $Pv \to 0$, so they don't solve the problem. As I mention in my answer, (2) works if you normalize $x$ instead of $w$. –  Geoffrey Irving Jan 17 at 22:54
$Pv \rightarrow 0$ is equivalent to saying as $v \rightarrow n$ (modulo normalization). The nullspace of $P$ is going to be the orthogonal complement of the range of $U$, so in your scenario, as $Pv \rightarrow 0$, $U^{H}v \rightarrow 0$ because both columns of $U$ are orthogonal to $n$. You'll suffer from the same catastrophic precision loss. –  Geoff Oxberry Jan 17 at 23:03
If you normalize $x$ instead of $w$, the result will be orthogonal to $n$ to machine precision, which is one of the requirements. If $w$ only has a few bits of precision, it and its normalization will only be orthogonal to $n$ to a few bits, which is insufficient. For example, if I'm projecting a point onto the surface of an infinite cylinder, I want to guarantee that the result is on the cylinder regardless of singularities. –  Geoffrey Irving Jan 18 at 0:58
I'm under the impression that the loss of precision was due to dividing by what could potentially be a small number, which is why I'm not following your argument. Both $\|x\|$ and $\|w\|$ will become small as $v \rightarrow n$. Not only that, but $\|x\| = \|w\|$, because you can make a unitary matrix $R$ from $U$ such that $Rw = (x_{1}, x_{2}, \ldots, x_{n-1}, 0)$, and unitary matrices are norm-preserving. What am I missing here? –  Geoff Oxberry Jan 18 at 1:24
@GeoffreyIrving Are you overlooking the (excellent) stability of floating point multiplication? We have that $U (\alpha x) \approx \alpha U x$ for any scalar $\alpha$ (assuming no denormals). It should not matter whether you normalize $x$ or $w$ because the direction is effectively the same. –  Jed Brown Jan 18 at 3:44
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.