Computing the angle between two vectors (vectorized) for small angles and with few copies

Question

I am implementing a function that computes the angle between two vectors when given two n-dimensional arrays and an axis along which to operate. I want to do this with as few copies as possible, and as numerically stable as possible. (The latter gets problematic for small angles (see here)).

What I would like reviewed is primarily the speed of the implementation. There is a sign error hidden somewhere, so if you find it I'd love to know about it, too.

Here is a (primitive) test case for the function I implemented.

vec_a = (1, 0, 0)
vec_b = (0, 1, 0)
result = angle_between(vec_a, vec_b)
assert np.allclose(result, np.pi/2)
vec_a = (1, 0)
vec_b = (0, 1)
result = angle_between(vec_a, vec_b)
assert np.allclose(result, np.pi/2)
result = angle_between(vec_b, vec_a)
assert np.allclose(result, -np.pi/2)
vec_a = (1, 0)
vec_b = (0, 2)
result = angle_between(vec_a, vec_b)
assert np.allclose(result, -np.pi/2)
vec_a = ((1, 0), (2, 0))
vec_b = ((0, 2), (0, 1))
result = angle_between(vec_a, vec_b)
assert np.allclose(result, (-np.pi/2, np.pi/2))

Note that the test currently fails, because of the sign error, but you can remove the offending assert.

Here is the function

def angle_between(vec_a: ArrayLike, vec_b: ArrayLike, *, axis:int=-1) -> np.ndarray:
 """Computes the signed angle from vec_a to vec_b (in a right-handed frame)
 Notes
 -----
 Implementation is based on this StackOverflow post:
 https://scicomp.stackexchange.com/a/27694
 """
 vec_a = np.asarray(vec_a)[None, :]
 vec_b = np.asarray(vec_b)[None, :]
 if axis >= 0:
 axis += 1
 len_c = np.linalg.norm(vec_a - vec_b, axis=axis)
 len_a = np.linalg.norm(vec_a, axis=axis)
 len_b = np.linalg.norm(vec_b, axis=axis)
 flipped = np.where(len_a >= len_b, 1, -1)
 mask = len_a >= len_b
 tmp = np.where(mask, len_a, len_b)
 np.putmask(len_b, ~mask, len_a)
 len_a = tmp
 mask = len_c > len_b
 mu = np.where(mask, len_b - (len_a - len_c), len_c - (len_a - len_b))
 numerator = ((len_a - len_b) + len_c) * mu
 denominator = (len_a + (len_b + len_c)) * ((len_a - len_c) + len_b)
 mask = denominator != 0
 angle = np.divide(numerator, denominator, where=mask)
 np.sqrt(angle, out=angle)
 np.arctan(angle, out=angle)
 angle *= 2
 angle[~mask] = np.pi
 angle *= flipped
 return angle[0]

Answer 1

I want to do this with as few copies as possible

This should not be a concern while your function is incorrect (your "sign error"). But also: should this -

vec_a = (1, 0)
vec_b = (0, 1)

really have the opposite angle sign from this?

vec_a = (1, 0)
vec_b = (0, 2)

A vector's argument is invariant under uniform scaling, so I think the answer must be "no", regardless of whether you're running a generic cosine similarity or some exotic low-angle algorithm.

Other bits -

Don't use two calls to asarray. If you're making a library-like function that does array coercion, then you're better off issuing a single call to broadcast_vectors.

I really don't understand why your code is prepending a new axis. It shouldn't need to do that.

Your code at least superficially resembles your linked Heron's method, but if it doesn't work, then no amount of optimisation makes sense. I suggest performing a comparative test between it and a better-known method such as the cosine similarity demonstration below.

import numpy as np
def angle_between(vec_a: np.typing.ArrayLike, vec_b: np.typing.ArrayLike, *, axis: int = -1) -> np.ndarray:
 """
 Computes the angle in signed radians from vec_a to vec_b in a right-handed frame.
 In all cases the magnitude of the angle follows simple cosine similarity.
 In the two-dimensional case, the angle sign is based on a simple vector-argument comparison.
 In the three-dimensional case, consider the plane shared by vec_a and vec_b. The sign of
 the plane's normal is chosen so that it minimises the angular deflection from an
 all-ones vector to disambiguate chirality. The angle sign is then the same as the 2D case
 for the two vectors projected to the plane oriented by that normal.
 Signed rotation demands a well-defined plane of rotation to disambiguate chirality.
 Since two vectors cannot fix a hyperplane in four-dimensional or higher space, those
 higher cases are all returned with a positive sign.
 """
 vec_a, vec_b = np.broadcast_arrays(vec_a, vec_b)
 # Simple cosine similarity
 abs_angle = np.acos(
 np.vecdot(vec_a, vec_b, axis=axis)
 /np.linalg.norm(vec_a, axis=axis)
 /np.linalg.norm(vec_b, axis=axis)
 )
 dim = vec_a.shape[axis]
 if 2 <= dim <= 3:
 if dim == 3:
 # normal vector orienting the plane intersecting a and b
 normal = np.linalg.cross(vec_a, vec_b, axis=axis)
 # coerce the normal to conventional orientation
 ones_basis = np.full(shape=normal.shape, fill_value=np.sqrt(1 / dim)) # use the normal closest to this
 normal_to_ones = np.acos(
 np.vecdot(normal, ones_basis, axis=axis)
 / np.linalg.norm(normal, axis=axis)
 # the ones basis is already a unit vector, so no second norm in denominator
 )
 normal *= 1 - 2*(normal_to_ones > 0.5*np.pi) # -1 if > 180 # now the normal has conventional orientation
 # project to the plane
 u_basis = vec_a
 v_basis = np.linalg.cross(normal, vec_a, axis=axis)
 ax = np.vecdot(vec_a, u_basis, axis=axis) # u-dimension in first vector
 bx = np.vecdot(vec_b, u_basis, axis=axis) # u-dimension in second vector
 ay = np.vecdot(vec_a, v_basis, axis=axis) # v-dimension in first vector
 by = np.vecdot(vec_b, v_basis, axis=axis) # v-dimension in second vector
 else:
 indexer = [slice(None)] * len(vec_a.shape) # to start, select all axes
 xidx = list(indexer) # copy
 xidx[axis] = 0 # select first (x) position in axis of interest
 xidx = tuple(xidx) # multi-axis indexer must be tuple
 ax = vec_a[xidx] # x-dimension in first vector
 bx = vec_b[xidx] # x-dimension in second vector
 yidx = list(indexer) # copy
 yidx[axis] = 1 # select second (y) position in axis of interest
 yidx = tuple(yidx) # multi-axis indexer must be tuple
 ay = vec_a[yidx] # y-dimension in first vector
 by = vec_b[yidx] # y-dimension in second vector
 ta = np.atan2(ay, ax) # argument, first vector
 tb = np.atan2(by, bx) # argument, second vector
 sign = np.ones_like(ta) # sign vector, preserving shape and dtype
 sign[tb < ta] = -1
 else:
 return abs_angle # implied sign of 1
 return sign*abs_angle
def unit_test() -> None:
 vec_a = (1, 0, 0)
 vec_b = (0, 1, 0)
 result = angle_between(vec_a, vec_b)
 assert np.isclose(result, 0.5*np.pi, atol=0, rtol=1e-14)
 result = angle_between(vec_b, vec_a)
 assert np.isclose(result, -0.5*np.pi, atol=0, rtol=1e-14)
 vec_a = (-0.5, 0, 0)
 vec_b = ( 0, 3, 0)
 result = angle_between(vec_a, vec_b)
 assert np.isclose(result, -0.5*np.pi, atol=0, rtol=1e-14)
 result = angle_between(vec_b, vec_a)
 assert np.isclose(result, 0.5*np.pi, atol=0, rtol=1e-14)
 vec_a = (1, 0)
 vec_b = (0, 1)
 result = angle_between(vec_a, vec_b)
 assert np.isclose(result, 0.5*np.pi, atol=0, rtol=1e-14)
 result = angle_between(vec_b, vec_a)
 assert np.isclose(result, -0.5*np.pi)
 vec_a = (1, 0)
 vec_b = (0, 2)
 result = angle_between(vec_a, vec_b)
 assert np.isclose(result, 0.5*np.pi, atol=0, rtol=1e-14)
 result = angle_between(vec_b, vec_a)
 assert np.isclose(result, -0.5*np.pi, atol=0, rtol=1e-14)
 vec_a = (
 (1, 0),
 (2, 0),
 (3, 0),
 )
 vec_b = (
 (0, 2),
 (0, 1),
 (0, 0.5),
 )
 result = angle_between(vec_a, vec_b)
 assert np.allclose(result, (0.5*np.pi, 0.5*np.pi, 0.5*np.pi), atol=0, rtol=1e-14)
 result = angle_between(vec_b, vec_a)
 assert np.allclose(result, (-0.5*np.pi, -0.5*np.pi, -0.5*np.pi), atol=0, rtol=1e-14)
if __name__ == '__main__':
 unit_test()