Calculate the angle between the rows of two matrices in numpy

Question

I have two matrices consisting of 3d vectors (numpy 1D arrays) and I need to calculate the angle between the vectors, row-wise, and return the results in a 1d array. I know how to calculate the angle between two 1d vectors. What is the the proper way to do this?

*** The resulting angles are in degrees not radians.

By now I have this:

import numpy as np

A = np.array([[1,0,0],
              [0,1,0],
              [0,0,1]])

B = np.array([[1,0,1],
              [1,1,0],
              [0,1,0]])

def angle(V1,V2):
    """
    angle between vectors V1 and V2 in degrees using
    angle = arccos ( V1 dot V2 / norm(V1) * norm(V2) ) *180/np.pi
    """

    cos_of_angle = V1.dot(V2) / (np.linalg.norm(V1) * np.linalg.norm(V2)) 
    return np.arccos(np.clip(cos_of_angle,-1,1))  * 180/np.pi

Note the scaling term 180/np.pi for the conversion from rad to deg.

I would like to have an array:

C = [ angle(A[0],B[0]) , angle(A[1],B[1])...... and so on]

Really appreciated if someone could help.

Eller

The indentation and parenthesis looked wrong. Take a look and edit further if needed. – Mikol

Corrected. Thanks @Mikol – Eller

Answer 1

We could use einsum to replace the dot-product computations and axis param for the norm ones to have a vectorized solution, like so -

def angle_rowwise(A, B):
    p1 = np.einsum('ij,ij->i',A,B)
    p2 = np.linalg.norm(A,axis=1)
    p3 = np.linalg.norm(B,axis=1)
    p4 = p1 / (p2*p3)
    return np.arccos(np.clip(p4,-1.0,1.0))

We could optimize further and bring in more of einsum, specifically to compute norms with it. Hence, we could use it like so -

def angle_rowwise_v2(A, B):
    p1 = np.einsum('ij,ij->i',A,B)
    p2 = np.einsum('ij,ij->i',A,A)
    p3 = np.einsum('ij,ij->i',B,B)
    p4 = p1 / np.sqrt(p2*p3)
    return np.arccos(np.clip(p4,-1.0,1.0))

---

Hence, to solve our case to get output in degrees -

out = angle_rowwise(A, B)*180/np.pi

Mikol

are these angles relative to the origin (0, 0, 0)? I am trying to translate the idea to using 3D points. – Achaemenid

@Achaemenid I am not really sure about the theory there. I was just focusing on translating it to a vectorized code. – Mikol

it must be something different since you get different results than the implementation here #2827893 a = np.array([[1, 0, 0], [1, 0, 0], [1, 0, 0]]) and b = np.array([[1, 0, 0], [0, 1, 0], [0, -1, 0]]) to represent identical, same and opposite directional vectors pairing – Achaemenid

@Achaemenid That implementation looks different than the one posted by OP here. I am following this one here. – Mikol

Dear friends, thank you for the kind help.I edited the original post, as I have noticed I forgot to include the np.arccos in the return statement of the angle function. And as pointed by @Nan, I have included the correction suggested at (#2827893) where a Nan would be returned in case of parallel vectors. – Eller

@RaphaelSantos So, are you sure about the * 180/np.pi that you are proposing in this question? I don't see that in the linked question. – Mikol

@Mikol the original post returned radians, the 180/np.pi is to convert to degrees... preference for explicit reading is simply np.degrees which does the same – Achaemenid

@RaphaelSantos Edited answer. – Mikol

@Achaemenid and Divakar I can include np.degrees in the original post if you think it makes it easier for others to understand, otherwise I'll call it a wrap for now. – Eller

@RaphaelSantos Thin, you can just add in the question that the final output required is in degrees and hence that scaling with 180/np.pi. Looks good otherwise. – Mikol

To all.. deja vu codereview.stackexchange.com/questions/54347/… – Achaemenid

@Divakar, still finding my way around here. Thank you for pointing it out. – Eller

Answer 2

If you're working with 3D vectors, you can do this concisely using the toolbelt vg. It's a light layer on top of numpy and it works equally well with individual vectors and stacks of vectors.

import numpy as np
import vg

A = np.array([[1,0,0],
              [0,1,0],
              [0,0,1]])

B = np.array([[1,0,1],
              [1,1,0],
              [0,1,0]])

vg.angle(A, B)

I created the library at my last startup, where it was motivated by uses like this: simple ideas which are verbose or opaque in NumPy.

Elaterin