When writing eigen solvers, it’s useful to be able to generate an orthonormal basis containing a particular vector. Another way of saying this is: given some fixed vector, find other vectors that are perpendicular to it and each other. I am under the impression from some websites that this is also useful in computer graphics (but I know nothing about that - let me know if you do!).
I’ve seen people ask how to do this on the internet (for example: here and here) and in general the recommended process is:
- If only one perpendicular vector is required, use a heuristic to come up with something perpendicular (see previous links).
- If two perpendicular vectors are required, create one perpendicular vector using the above heuristic then use the cross-product to come up with the final vector.
- For more than two perpendicular vectors, generate random vectors and use Gram Schmidt to (hopefully) orthogonalise them.
- For 3D graphics, see the “Update” heading below…
The first two usually require branching in the heuristic. The second may not work (if you’re really really unlucky) and requires multiple normalisation operations. Follows is an approach which I haven’t seen documented elsewhere which requires one square root for an arbitrary dimension vector set.
What we want is some technique which creates an orthogonal matrix from a vector. This is exactly what the Householder transformation matrix does:
Where
Explicitly, this looks like:
If
Here is a very inefficient C implementation (the square root is not actually necessary - this can be seen in the fact that two
/* The first n values in m is a normalised vector which is desired to be
* part of the basis set. Following this vector, n-1 new vectors of length n
* will be added which are orthogonal to it. */
void house_ortho(double *m, unsigned n)
{
const double x = sqrt(1.0 - m[0]);
/* As m[0] approaches 1, x approaches 0. This detects that and will end
* up forcing the result to be the identity matrix. */
const double s = (x < 1e-20) ? 0.0 : (-1.0 / x);
unsigned i, j;
/* We first convert the m vector into the v vector. Note that we have
* incorporated the multiply by 2 present in the Householder operation
* into v by multiplying every element by sqrt(2). */
for (m[0] = x, i = 1; i < n; i++)
[i] *= s;
m/* Create all of the orthogonal vectors using the definition of the
* Householder transformation matrix: I - 2 v v*. */
for (i = 1; i < n; i++) {
for (j = 0; j < n; j++) {
[i*n+j] = ((i == j) ? 1.0 : 0.0) - m[i] * m[j];
m}
}
/* Convert the v vector back into m. Probably could have written this
* code in a way where this didn't need to happen. */
for (m[0] = 1.0 - x * x, j = 1; j < n; j++)
[j] = -m[j] * x;
m}
The termination condition for the for loop (over i) in the middle of the code sample could be modified to generate a smaller number of vectors.
This process could be modified to work with a complex valued
Update
I posted a link to this blog on Twitter and reached out to see if anyone else had used this method. Warren Moore helpfully responded and pointed me at this paper “Building an Orthonormal Basis, Revisited” (Pixar) which describes numerical and implementation improvements to an algorithm proposed in this paper “Building an Orthonormal Basis from a 3D Unit Vector Without Normalization” (Frisvad, J. R.).
The algorithm which I’ve described in terms of finding the vector
Note: we don’t actually need
The Pixar paper addresses the singularity which occurs as
In this case, the vector
Which moves the singularity to occur as
Efficient implementation of 3-dimension basis finder
As shown in the Pixar paper, there is no need to test the sign of
Modifying some signs in the explicit matrices given previously, we can re-write the matrices as:
Or more simply:
Where:
And here’s the C implementation:
/* The first three values in v are a three point unit vector. The
* function produces two more three point vectors of the same coordinate
* space which are orthogonal. */
void ortho3(float *v)
{
float q1 = v[0];
float q2 = v[1];
float q3 = v[2];
float sign = copysignf(1.0f, q1);
float scale = 1.0f / (q1 + sign);
float q2scale = q2 * scale;
float q3scale = q3 * scale;
float q3q2scale = q2scale * q3;
[3] = q2;
v[4] = q2 * q2scale - sign;
v[5] = q3q2scale;
v[6] = q3;
v[7] = q3q2scale;
v[8] = q3 * q3scale - sign;
v}