ENH: spatial: Matrix multiplication of Rotation

[moradpur]

Is your feature request related to a problem? Please describe.

In scipy.spatial.transform.Rotation.mul , If p and q are two rotations, and both p and q contain N rotations, each rotation p[i] is composed with the corresponding rotation q[i], and the output contains N rotations which is not suitable for matrix multiplication and machine learning purposes.

Describe the solution you'd like.

For machine learning purposes, the matrix multiplication of multiple rotations should be possible. I suggest that it gives N*N rotations instead of N Rotation.

Describe alternatives you've considered.

Instead of changing __mul__ method you can add matmul method which is available in numpy too. However, using similarity of rotations as (rotation_1*rotation_2.inv()).magnitude() and giving a covariance N*N matrix will be a blockbuster feature.

Additional context (e.g. screenshots, GIFs)

No response

[scottshambaugh]

I don't think we should make the change to __mul__, since it would break standard numpy broadcasting. And putting this behavior in __matmul__ doesn't pattern match with how matmul works in numpy, e.g. a @ b == np.dot(a,b)

However, you can get this to work by storing the rotations in numpy arrays and using its broadcasting rules. Here's an example of how to get what you're looking for:

import numpy as np
from scipy.spatial.transform import Rotation as R

rotations_array_4 = np.array(R.random(4))[:, np.newaxis]  # Shape will be (4, 1)
rotations_array_3 = np.array(R.random(3))  # Shape will be (3,)
resultant_rotations43 = rotations_array_4 * rotations_array_3 

print(resultant_rotations43.shape)
# prints (4, 3), all Rotation objects

You can also keep this going to higher dimensions:

resultant_rotations433 = resultant_rotations43[:, :, np.newaxis] * rotations_array_3
print(resultant_rotations433.shape)
# prints (4, 3, 3), all Rotation objects

[scottshambaugh]

Or, to do your blockbuster feature 🙂, you can do the following:

R3 = R.random(3, random_state=0)
rotations_array_3 = np.array(R3)[:, np.newaxis]
rotations_array_3_inv = np.array(R3.inv())
vectorized_magnitude = np.vectorize(lambda x: x.magnitude())
angles = vectorized_magnitude(rotations_array_3 * rotations_array_3_inv)

# In one line:
# angles = np.vectorize(lambda x: x.magnitude())(np.array(R3)[:, np.newaxis] * np.array(R3.inv()))

print(angles)
''' prints out:
[[1.14439170e-16, 2.10407917e+00, 1.50255252e+00],
 [2.10407917e+00, 0.00000000e+00, 2.75493397e+00],
 [1.50255252e+00, 2.75493397e+00, 0.00000000e+00]]
'''

[scottshambaugh]

@nmayorov thoughts?

[nmayorov]

In scipy.spatial.transform.Rotation.mul , If p and q are two rotations, and both p and q contain N rotations, each rotation p[i] is composed with the corresponding rotation q[i], and the output contains N rotations which is not suitable for matrix multiplication and machine learning purposes.

I think the current rules are most suitable for targeted engineering tasks (describe attitude of rotating frames, project or rotate vectors between the frames).

For machine learning purposes, the matrix multiplication of multiple rotations should be possible. I suggest that it gives N*N rotations instead of N Rotation.

Just changing the __mul__ rule is clearly not possible. As @scottshambaugh said the required broadcasting can be achieved by working with matrix (or possibly quaterion) representation of rotations.

Instead of changing mul method you can add matmul method which is available in numpy too. However, using similarity of rotations as (rotation_1rotation_2.inv()).magnitude() and giving a covariance NN matrix will be a blockbuster feature.

I believe this is the case where we can confidently recommend to implement the required logic on top what is already provided. I don't see how we can logically fit "outer product" broadcasting rules into Rotation.

[moradpur]

Or, to do your blockbuster feature 🙂, you can do the following:

R3 = R.random(3, random_state=0)
rotations_array_3 = np.array(R3)[:, np.newaxis]
rotations_array_3_inv = np.array(R3.inv())
vectorized_magnitude = np.vectorize(lambda x: x.magnitude())
angles = vectorized_magnitude(rotations_array_3 * rotations_array_3_inv)

# In one line:
# angles = np.vectorize(lambda x: x.magnitude())(np.array(R3)[:, np.newaxis] * np.array(R3.inv()))

print(angles)
''' prints out:
[[1.14439170e-16, 2.10407917e+00, 1.50255252e+00],
 [2.10407917e+00, 0.00000000e+00, 2.75493397e+00],
 [1.50255252e+00, 2.75493397e+00, 0.00000000e+00]]
'''

[moradpur]

@scottshambaugh
Actually, I have tried both of these solutions before. The problem is that they are much slower than just doing a for loop over the regular Rotation(single)*Rotation(Multiple).