如何有效地计算 numpy 中的回转张量？

Question

3d空间中一组N个点的gyration tensor定义为

假设条件

.

如何在不使用显式 for 循环的情况下在 numpy 中计算它？我知道我可以做类似的事情

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)

但这对于大 N 来说很快就会变得低效，因为 for 循环。有没有更好的方法？

Answer 1

您可以像这样使用np.einsum：

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)    

输出：

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