标签: numpy-einsum

我很难理解究竟是如何einsum运作的.我看过文档和一些例子,但它似乎并不坚持.

这是我们在课堂上看到的一个例子:

C = np.einsum("ij,jk->ki", A, B)

对于两个数组A和B

我想这会采取A^T * B,但我不确定(它正在将其中一个的转置正确吗？).任何人都可以告诉我这里发生了什么(通常在使用时einsum)？

在tensorflow中,函数tf.einsum,tf.matmul和tf.tensordot都可以用于相同的任务.(我意识到tf.einsum并且tf.tensordot有更多的一般定义;我也意识到它tf.matmul具有批处理功能.)在可以使用三者中的任何一种的情况下,一个函数是否最快？还有其他推荐规则吗？

例如,假设这A是一个rank-2张量,并且b是rank-1张量,并且您想要计算产品c_j = A_ij b_j.在三个选项中:

c = tf.einsum('ij,j->i', A, b)

c = tf.matmul(A, tf.expand_dims(b,1))

c = tf.tensordot(A, b, 1)

是否通常比其他人更好？

想象一下，我有整数n,q和具有以下维度的向量/数组：

import numpy as np
n = 100
q = 102
A = np.random.normal(size=(n,n))
B = np.random.normal(size=(q, ))

C = np.einsum("i, jk -> ijk", B, A)
D = np.einsum('ijk, ikj -> k', C, C)

如果所有中间数组都适合内存，则工作正常。

现在假设我可以存储 size 的内存数组(n,n)，(q,n)但不能存储任何三维数组，例如 shape (n,n,q)。我无法存储在C上面的内存数组中。相反，要计算D，

D1 = np.einsum('i, jk, i, kj -> k', B, A, B, A, optimize='optimal')

工作正常，np.einsum通常足够聪明，可以找到一个einsum_path3d 数组。伟大的！

现在让我们稍微复杂一点：

C = np.einsum("i, jk -> ijk", B, A)    # …

我有一些嵌套循环（总共三个），我试图使用 numpy.einsum 来加速计算，但我正在努力使符号正确。我设法摆脱了一个循环，但我无法弄清楚另外两个。这是我到目前为止所得到的：

import numpy as np
import time

def myfunc(r, q, f):
    nr = r.shape[0]
    nq = q.shape[0]
    y = np.zeros(nq)
    for ri in range(nr):
        for qi in range(nq):
            y[qi] += np.einsum('i,i',f[ri,qi]*f[:,qi],np.sinc(q[qi]*r[ri,:]/np.pi))
    return y

r = np.random.random(size=(1000,1000))
q = np.linspace(0,1,1001)
f = np.random.random(size=(r.shape[0],q.shape[0]))

start = time.time()
y = myfunc(r, q, f)
end = time.time()

print(end-start)

虽然这比原来快得多，但这仍然太慢，大约需要 30 秒。请注意，没有 einsum 调用的原始内容如下（看起来需要大约 2.5 小时，迫不及待地想确定）：

def myfunc(r, q, f):
    nr = r.shape[0]
    nq = q.shape[0]
    y = np.zeros(nq)
    for ri in range(nr): …

我np.einsum用来乘以概率表,如:

np.einsum('ijk,jklm->ijklm', A, B)

问题是我总共处理超过26个随机变量(轴),所以如果我为每个随机变量分配一个字母,我就会用完字母.有没有其他方法可以指定上述操作来避免这个问题,而不是诉诸混乱np.sum和np.dot操作？

我在'np.einsum'的文档中理解的是,置换字符串将给出向量中轴的置换.以下实验证实了这一点:

>>> M = np.arange(24).reshape(2,3,4)
>>> M.shape
(2, 3, 4)
>>> np.einsum('ijk', M).shape
(2, 3, 4)
>>> np.einsum('ikj', M).shape
(2, 4, 3)
>>> np.einsum('jik', M).shape
(3, 2, 4)

但这我无法理解:

>>> np.einsum('kij', M).shape
(3, 4, 2)

我希望(4,2,3)相反......我的理解有什么问题？

我有一组使用 meshgrid() 生成的 3d 坐标。我希望能够围绕 3 个轴旋转这些。

我尝试解开网格并在每个点上旋转，但网​​格很大并且内存不足。

这个问题使用 einsum() 在 2d 中解决了这个问题，但是在将其扩展到 3d 时我无法弄清楚字符串格式。

我已经阅读了有关 einsum() 及其格式字符串的其他几页，但一直无法弄清楚。

编辑：

我称我的网格轴为 X、Y 和 Z，每个轴的形状为 (213, 48, 37)。此外，当我尝试将结果放回网格时，实际出现了内存错误。

当我试图“解开”它以逐点旋转时，我使用了以下函数：

def mg2coords(X, Y, Z):
    return np.vstack([X.ravel(), Y.ravel(), Z.ravel()]).T

我用以下内容循环了结果：

def rotz(angle, point):
    rad = np.radians(angle)
    sin = np.sin(rad)
    cos = np.cos(rad)
    rot = [[cos, -sin, 0],
           [sin,  cos, 0],
           [0, 0, 1]]

    return np.dot(rot, point)

旋转后，我将使用这些点进行插值。

我正在使用numpy einsum对一些3维和4维张量进行相当复杂的操作.

我的实际代码是

np.einsum('oij,imj,mjkn,lnk,plk->op',phi,B,Suu,B,phi)

这就是我想要的.

使用einsum_path,结果是:

>>> path = np.einsum_path('oij,imj,mjkn,lnk,plk->op',phi,B,Suu,B,phi)

>>> print(path[0])
['einsum_path', (0, 1), (0, 3), (0, 1), (0, 1)]

>>> print(path[1])
  Complete contraction:  oij,imj,mjkn,lnk,plk->op
         Naive scaling:  8
     Optimized scaling:  5
      Naive FLOP count:  2.668e+07
  Optimized FLOP count:  1.340e+05
   Theoretical speedup:  199.136
  Largest intermediate:  7.700e+02 elements
--------------------------------------------------------------------------
scaling                  current                                remaining
--------------------------------------------------------------------------
   4                imj,oij->moj                     mjkn,lnk,plk,moj->op
   5               moj,mjkn->nok                          lnk,plk,nok->op
   4                plk,lnk->npk                              nok,npk->op
   4                 npk,nok->op                                   op->op

这表明理论加速约为200倍.

如何使用此结果来加速我的代码？我如何"实现"einsum_path告诉我的内容？

我有一段代码类型：

nnt = np.real(np.einsum('xa,xb,yc,yd,abcde->exy',evec,evec,evec,evec,quartic))

其中evec是（比如说）一个 L x Lnp.float32阵列，并且quartic是一个 L x L x L x L x Tnp.complex64阵列。

我发现这个例程相当慢。

我认为由于所有evec的都是相同的，可能有更快的方法吗？

提前致谢。

我想将以下 numpy einsum 编写为特征张量操作

import numpy as np

L = np.random.rand(2, 2, 136)
U = np.random.rand(2, 2, 136)

result = np.einsum('ijl,jkl->ikl', U, L)

我可以像 C++ 中那样用 for 循环编写它

  for (int i = 0; i < 2; i++) {
    for (int j = 0; j < 2; j++) {
      for (int k = 0; k < 2; k++) {
        for (int l = 0; l < 136; l++) {
          result(i, k, l) += U(i, j, l) * L(j, k, l);
        } …