How do I efficiently compute the gyration tensor in numpy?

Question

The gyration tensor of a set of N points in 3d space is defined as

assuming the condition

.

How do I compute this in numpy without using an explicit for loop? I know that I can just do something like

import numpy as np

def calculate_gyration_tensor(points):
    '''
    Calculates the gyration tensor of a set of points.
    '''
    COM = centre_of_mass(points)
    gyration_tensor = np.zeros((3, 3))
    for p in points:
        gyration_tensor += np.outer(p-COM, p-COM)
    return gyration_tensor / len(points)

but this quickly becomes inefficient for large N, because of the for loop. Is there a better way to do it?

score 1 · Answer 1 · answered Jun 02 '20 at 14:28

1

You can do with np.einsum like this:

def gyration(points):
    '''
    Calculate the gyrason tensor
    points : numpy array of shape N x 3
    '''

    center = points.mean(0)

    # normalized points
    normed_points = points - center[None,:]

    return np.einsum('im,in->mn', normed_points,normed_points)/len(points)


# test
points = np.arange(36).reshape(12,3)

gyration(points)

Output:

array([[107.25, 107.25, 107.25],
       [107.25, 107.25, 107.25],
       [107.25, 107.25, 107.25]])

answered Jun 02 '20 at 14:28

Quang Hoang

146,074
10
56
74

Thanks - can you please explain what does einsum do, especially what that string means? – Ferdinando Randisi Jun 02 '20 at 16:02
You can have a read on [this excellent answer](https://stackoverflow.com/a/33641428/4238408). Basically, your operation is just `normed_points.T @ normed_points`. – Quang Hoang Jun 02 '20 at 16:07

How do I efficiently compute the gyration tensor in numpy?

1 Answers1