Algorithm to get an arbitrary perpendicular vector

Question

Is there a more efficient/direct way to find an arbitrary perpendicular vector of another vector?

def perpendicular_vector(v):
 r""" Finds an arbitrary perpendicular vector to *v*. """
 a, b = random.random(), random.random()
 if not iszero(v.z):
 x, y, z = v.x, v.y, v.z
 elif not iszero(v.y):
 x, y, z = v.x, v.z, v.y
 elif not iszero(v.x):
 x, y, z = v.y, v.z, v.x
 else:
 raise ValueError('zero-vector')
 c = (- x * a - y * b) / z
 if not iszero(v.z):
 return Vector(a, b, c)
 elif not iszero(v.y):
 return Vector(a, c, b)
 elif not iszero(v.x):
 return Vector(b, c, a)

• slepic: You better not use random. Arbitrary does not mean differrent on each call. Better use size that comes out of your algorithm. Either it be same size as input vector or 1 could be a good choice or whatever comes out of your equations. Using random to explicitly resize the output vector Is unnecesary performance drop.

Answer 1

The cross product of two vectors is perpendicular to both vectors, unless both vectors are parallel. So you could simply take the cross product of your first vector with (1, 0, 0), unless it is parallel to (1, 0, 0), in which case you could use (0, 1, 0). If you use lists rather than dedicated classes with attributes and are willing to use numpy, this gets ridiculously short:

import numpy as np
def perpendicular_vector(v):
 if v[1] == 0 and v[2] == 0:
 if v[0] == 0:
 raise ValueError('zero vector')
 else:
 return np.cross(v, [0, 1, 0]]
 return np.cross(v, [1, 0, 0])

Answer 2

Hmm. The lack of comments make it slightly non-obvious what is happening. The inverted conditions not iszero(...) don't make it any easier to understand. The same set of these conditions occurs two times, which is a bit confusing. I also don't know where the iszero(...) function is from, and will replace it by ... == 0 in the following.

Now a bit of math first. Two vectors are perpendicular if their scalar product is zero:

1st vector (x, y, z)
2nd vector (a, b, c)
0 = a·x + b·y + c·z

You are correct in avoiding the trivial solution x = y = z = 0. Note that you do not avoid the solution a = b = c = 0, because random.random() can return zero.

If one of x, y, z is zero, then the above equation can be solved by choosing a non-zero value for the corresponding direction, and setting the other variables to zero. Example: given (1, 0, 3) as a first vector, the equation can be solved by a second vector (0, 1, 0).

The above rule can also be applied in reverse if the first vector has two zero fields, but it's only required to change one of the zero fields to a non-zero value.

We can encode these simple cases in pseudocode as

match vec1 with
| (0, 0, 0) -> ValueError('zero-vector')
| (0, _, _) -> Vector(1, 0, 0)
| (_, 0, _) -> Vector(0, 1, 0)
| (_, _, 0) -> Vector(0, 0, 1)
| (x, y, z) -> a more complex case which we'll handle in a moment

If all parts of the input vector are non-zero, the calculation is a bit more complex because we have three variables a, b, c to determine but only one equation – this means we can choose two values arbitrarily. Choosing a random value might make sense when hardening your application, but it's difficult to test, and we can make the code simpler by choosing a = b = 1 (we cannot choose a = b = 0). The equation can now be used to calculate c:

c = -(x + y)/z

which you essentially have used as well.

So to wrap it all up, I would write this code:

def perpendicular_vector(v):
 r""" Finds an arbitrary perpendicular vector to *v*."""
 # for two vectors (x, y, z) and (a, b, c) to be perpendicular,
 # the following equation has to be fulfilled
 # 0 = ax + by + cz
 # x = y = z = 0 is not an acceptable solution
 if v.x == v.y == v.z == 0:
 raise ValueError('zero-vector')
 # If one dimension is zero, this can be solved by setting that to
 # non-zero and the others to zero. Example: (4, 2, 0) lies in the
 # x-y-Plane, so (0, 0, 1) is orthogonal to the plane.
 if v.x == 0:
 return Vector(1, 0, 0)
 if v.y == 0:
 return Vector(0, 1, 0)
 if v.z == 0:
 return Vector(0, 0, 1)
 # arbitrarily set a = b = 1
 # then the equation simplifies to
 # c = -(x + y)/z
 return Vector(1, 1, -1.0 * (v.x + v.y) / v.z)

Answer 3

Jaime beat me to it!

Anyway, here's an optimised version that uses the same idea and doesn't call numpy.cross:

def perpendicular_vector(v):
 if iszero(v.x) and iszero(v.y):
 if iszero(v.z):
 # v is Vector(0, 0, 0)
 raise ValueError('zero vector')
 # v is Vector(0, 0, v.z)
 return Vector(0, 1, 0)
 return Vector(-v.y, v.x, 0)

• Toby Speight: You have presented an alternative solution, but haven't reviewed the code. Please edit to show what aspects of the question code prompted you to write this version, and in what ways it's an improvement over the original. It may be worth (re-)reading How to Answer.

Answer 4

Maybe I'm too late, but if you are using numpy arrays, here is a version that works with vectors of arbitrary length and that doesn't have ugly if statements to test which coordinate is zero:

import numpy as np
def perpendicular_vector(v0):
 idx_max = np.argmax(np.abs(v0))
 v1 = np.zeros(v0.shape)
 v1[idx_max] = -v0[(idx_max+1) % len(v0)]/v0[idx_max]
 v1[(idx_max+1) % len(v0)] = 1
 
 return v1

• pacmaninbw: This answer contains a meaningful observation about the original code, but FYI alternate solutions without meaningful observations are bad answers on Code Review. It would be a good idea to read How do I write a good answer.