三阶矩计算 - numpy答案

【问题标题】：Third order moment calculation - numpy三阶矩计算 - numpy
【发布时间】：2016-07-04 10:10:56
【问题描述】：

在 python 中，我有一个包含 N 行（示例数）和 n 列（特征数）的数组 X。

如果要计算二阶矩矩阵C

C[i,j] = E(x_i x_j)

那么我有两种可能：

首先，做循环：
对于范围内的 i (N)：对于范围（n）中的j：对于范围（n）中的k： C[j,k] = C[j,k] + X[i,j]*X[i,k]/N
第二，更简单，使用numpy乘积矩阵：
导入 numpy np C = np.transpose(X).dot(X)/N
实践中的第二个版本非常快。

如果现在我要计算三阶矩矩阵T

T[i,j,k] = E(x_i x_j x_k)

那么循环替代很简单：

对于范围内的 i (N)：对于范围（n）中的j：对于范围（n）中的k：对于范围内的 m (n)： T[j,k,m] = T[j,k,m] + X[i,j]*X[i,k]*X[i,m]/N

有没有一种使用 numpy 库来计算最后一个张量的快速方法，比如二阶矩？

【问题讨论】：

标签： python algorithm numpy matrix

【解决方案1】：

您可以使用NumPy's einsum notation 来解决您的二阶和三阶案例。

二阶：

np.einsum('ij,ik->jk',X,X)/N

三阶：

np.einsum('ij,ik,il->jkl',X,X,X)/N

可以看出，将其扩展到更高阶的情况会更容易/直观。

【讨论】：

将np.einsum('ij,ik->jk',X,X)/N 应用到我的数组（图像）时，我得到一个ValueError: operand has more dimensions than subscripts given in einstein sum。
@calocedrus 您的图像很可能是 3D 数组。所以，这行不通。请发布一个新问题。
确实，感谢您指出这一点！在这种情况下，正确的语法似乎应该是 np.einsum('ij...,ik...->jk...',X,X)/N 并且等效于三阶矩。

【解决方案2】：

我知道它在速度方面并不完美，但为什么不使用np.power(x, 3).sum() / N？它比点积慢，但比循环快。

In [1]: import numpy as np
In [2]: x = np.random.rand(10000)
In [3]: x.dot(x.T)
Out[3]: 3373.6189738897856
In [4]: %timeit(x.dot(x.T))
The slowest run took 48.74 times longer than the fastest. This could   mean that an intermediate result is being cached.
100000 loops, best of 3: 4.39 µs per loop
In [5]: %timeit(np.power(x, 2).sum())
The slowest run took 4.14 times longer than the fastest. This could mean that an intermediate result is being cached.
10000 loops, best of 3: 140 µs per loop
In [6]: np.power(x, 2).sum()
Out[6]: 3373.6189738897865

顺便说一句，这就是我计算时刻的方式......

【讨论】：