CODEWITHSUNDEEP
pythonpython-3.xnumpyprecisiontrigonometry
Bryson of Heraclea
Bryson of Heraclea

Reputation: 63

What is a robust method for calculating the cosine of the angle of two vectors in python?

I want to calculate the cosine of the angle of two 2D-vectors that lie on the same plane in Python. I use the classic formula cosθ = v1∙v2/(|v1||v2|). However, depending on how I implement it, the rounding errors can give values greater than 1 or less than -1 for parallel vectors. Three examples of the implementation follow:

import numpy as np

def cos2Dvec_1(v1,v2):
    v1norm2 = v1[0]**2+v1[1]**2
    v2norm2 = v2[0]**2+v2[1]**2
    v1v2 = v1[0]*v2[0]+v1[1]*v2[1]
    if v1v2 > 0:
        return np.sqrt(v1v2**2/(v1norm2*v2norm2))
    else:
        return -np.sqrt(v1v2**2/(v1norm2*v2norm2))

def cos2Dvec_2(v1,v2):
    r1 = np.linalg.norm(v1) 
    r2 = np.linalg.norm(v2)
    return np.dot(v1,v2)/(r1*r2)

def cos2Dvec_3(v1,v2):
    v1 = v1/np.linalg.norm(v1)
    v2 = v2/np.linalg.norm(v2)
    return np.dot(v1,v2)

v2 = np.array([2.0, 3.0])
v1 = -0.01 *v2
print(cos2Dvec_1(v2,v1), cos2Dvec_2(v1,v2), cos2Dvec_3(v1,v2))

When I execute this code on my machine (Python 3.7.6, Numpy 1.18.1), I get:

-1.0 -1.0000000000000002 -1.0

The difference between the 2nd and 3rd version of the function is especially weird to me...

If I use the second function to calculate the angle via the inverse cosine, I will get a nan, which I want to avoid (and I don't want to check every time, if the result lies in [-1,1]). I was wondering though whether versions 1 and 3 (which give the correct result) are numerically stable in all cases.

Upvotes: 2

Views: 288

Answers (2)

Bryson of Heraclea
Bryson of Heraclea

Reputation: 63

After timing numpy.clip(x,-1,1), I saw that a much faster alternative is just plain old y = max(-1, min(x,1)). Timing results for the following code snippet:

t1 = time.time()
vals = [-1.001, -0.99, 0.99, 1.001]
for ii in range(1):
    for x in vals:
        y = np.clip(x, -1, 1)
#       y = max(-1, min(x,1))
        print(y)
t2 = time.time()
print(t2-t1)

Numpy clip was almost 23 times slower on my laptop.

Upvotes: 0

Pac0
Pac0

Reputation: 23174

You could use the clip function to ensure values will be coerced (or "clamped" or "clipped") between -1 and 1.

return np.clip(YourOriginalResult, -1, 1)

Rounding errors are inevitable when manipulating floats. Unfortunately, I don't have enough expertise on your calculations to recommend any of your methods.

I'd suggest you to to check which ones have correct performance, after adding the clip part to mitigate those "out of range" issues.

Upvotes: 1

Related Questions