17

What is the simplest and most efficient ways in numpy to generate two orthonormal vectors a and b such that the cross product of the two vectors equals another unit vector k, which is already known?

I know there are infinitely many such pairs, and it doesn't matter to me which pairs I get as long as the conditions axb=k and a.b=0 are satisfied.

Physicist
  • 2,848
  • 8
  • 33
  • 62
  • `a` and `b` are shape (3,) (1d with 3 elements)? How would you do this 'by hand'? – hpaulj Nov 11 '15 at 19:54
  • yes. a, b and k are all of shape(3,). I know how to do it by hand. I have 5 equations in 6 unknowns, does it mean there are no simple ways (with just a few lines) and I have to use numerical methods in scipy in solve it? – Physicist Nov 11 '15 at 19:58
  • I was thinking of algebraic versions of the cross and dot products in the accepted answer. For small arrays those would be just as fast. – hpaulj Nov 11 '15 at 20:39

3 Answers3

21

The Gram-Schmidt procedure will do exactly this. For example:

>>> k  # an arbitrary unit vector k is not array. k is must be numpy class. np.array
np.array([ 0.59500984,  0.09655469, -0.79789754])

To obtain the 1st one:

>>> x = np.random.randn(3)  # take a random vector
>>> x -= x.dot(k) * k       # make it orthogonal to k
>>> x /= np.linalg.norm(x)  # normalize it

To obtain the 2nd one:

>>> y = np.cross(k, x)      # cross product with k

and to verify:

>>> np.linalg.norm(x), np.linalg.norm(y)
(1.0, 1.0)
>>> np.cross(x, y)          # same as k
array([ 0.59500984,  0.09655469, -0.79789754])
>>> np.dot(x, y)            # and they are orthogonal
-1.3877787807814457e-17
>>> np.dot(x, k)
-1.1102230246251565e-16
>>> np.dot(y, k)
0.0
Martijn Courteaux
  • 67,591
  • 47
  • 198
  • 287
behzad.nouri
  • 74,723
  • 18
  • 126
  • 124
  • 3
    This doesn't work where k is different vector ``` >>> k = np.array([ 0.0, 0.0215, -1.334]) >>> x = np.random.randn(3) >>> x -= x.dot(k) * k >>> x /= np.linalg.norm(x) >>> np.dot(x, k) 0.38532052441276377 ``` – Miae Kim Jul 23 '19 at 03:55
  • Excellent thanks, together with René Wirnata answer below this works even if the first input vector is NOT a unit vector. – Charly Empereur-mot Aug 04 '19 at 08:38
13

Sorry, I can't put it as a comment because of a lack of reputation.

Regarding @behzad.nouri's answer, note that if k is not a unit vector the code will not give an orthogonal vector anymore!

The correct and general way to do so is to subtract the longitudinal part of the random vector. The general formula for this is here

So you simply have to replace this in the original code:

>>> x -= x.dot(k) * k / np.linalg.norm(k)**2
Anshul Rai
  • 772
  • 7
  • 21
René Wirnata
  • 131
  • 1
  • 3
1

Assume the vector that supports the orthogonal basis is u.

b1 = np.cross(u, [1, 0, 0])   # [1, 0, 0] can be replaced by other vectors, just get a vector orthogonal to u
b2 = np.cross(u, b1)
b1, b2 = b1 / np.linalg.norm(b1), b2 / np.linalg.norm(b2)

A shorter answer if you like.

Get a transformation matrix

B = np.array([b1, b2])
TransB = np.dot(B.T, B)
u2b = TransB.dot(u) # should be like [0, 0, 0]
william liang
  • 11
  • 3