我需要执行一种特殊类型的张量收缩。我想要这样的东西：
A_{bg} = Sum_{a，a'，a''} ( B_{a} C_{a'b} D_{a''g} )
其中所有索引的值都可以为 0,1，并且对于 a+a'+a'' = 1 或 a+a'' = 2 的所有情况，将携带 a、a' 和 a'' 的总和。所以这就像爱因斯坦求和惯例的反面：我只想在三个指数中的一个与其他指数不同时才求和。
此外，我希望对未求和的索引数量有一定的灵活性：在示例中，生成的张量有 2 个索引，总和是 3 个张量的元素的乘积，一个有一个索引，另外两个有两个索引。这些索引的数量会有所不同，所以总的来说，我希望能够写这样的东西：
A_{...} = Sum_{a，a'，a''} ( B_{a...}C_{a...}D_{a''...} )
我想指出的是，索引的数量不是固定的，但它是受控的：我可以知道并指定每个张量在每个步骤中有多少个索引。
我尝试了np.einsum()，但显然我被迫对标准爱因斯坦约定中的重复索引求和，我不知道如何实现我在这里暴露的条件。
而且我不能用各种for编写所有内容，因为正如我所说，所涉及的张量的索引数量不是固定的。
有人有想法吗？
来自评论：
我会用这样的编程语言写我放在这里的内容：
tensa = np.zeros((2,2))
for be in range(2):
for ga in range(2):
for al in range(2):
for alp in range(2):
for alpp in range(res(al,alp),prod(al,alp)):
tensa[be,ga] += tensb[al] * tensc[alp,be] * tensd[alpp,ga]
其中res和prod是确保al+alp+alpp = 1 or 2的两个函数。这样做的问题是我需要指定所有涉及的索引，而我不能在所有晶格的一般计算中这样做。
首先，让我们在 Python 循环中写出您的示例，以便有一个比较的基线。如果我理解正确，这就是您想要做的：
b, g = 4, 5
B = np.random.rand(2)
C = np.random.rand(2, b)
D = np.random.rand(2, g)
out = np.zeros((b, g))
for j in (0, 1):
for k in (0, 1):
for l in (0, 1):
if j + k + l in (1, 2):
out += B[j] * C[k, :, None] * D[l, None, :]
当我运行这个时，我得到这个输出：
>>> out
array([[ 1.27679643, 2.26125361, 1.32775173, 1.5517918 , 0.47083151],
[ 0.84302586, 1.57516142, 1.1335904 , 1.14702252, 0.34226837],
[ 0.70592576, 1.34187278, 1.02080112, 0.99458563, 0.29535054],
[ 1.66907981, 3.07143067, 2.09677013, 2.20062463, 0.65961165]])
你不能直接用np.einsum来解决这个问题，但你可以运行它两次，得到你的结果作为这两者的区别：
>>> np.einsum('i,jk,lm->km', B, C, D) - np.einsum('i,ik,im->km', B, C, D)
array([[ 1.27679643, 2.26125361, 1.32775173, 1.5517918 , 0.47083151],
[ 0.84302586, 1.57516142, 1.1335904 , 1.14702252, 0.34226837],
[ 0.70592576, 1.34187278, 1.02080112, 0.99458563, 0.29535054],
[ 1.66907981, 3.07143067, 2.09677013, 2.20062463, 0.65961165]])
对np.einsum的第一个调用是将所有内容相加，而不管索引加起来是多少。第二个只将所有三个指数都相同的指数相加。所以很明显你的结果是两者的区别。
理想情况下，您现在可以继续编写以下内容：
>>>(np.einsum('i...,j...,k...->...', B, C, D) -
... np.einsum('i...,i...,i...->...', B, C, D))
无论 C 和 D 数组的尺寸如何，都能获得结果。如果您尝试第一个，您将收到以下错误消息：
ValueError: operands could not be broadcast together with remapped shapes
[original->remapped]: (2)->(2,newaxis,newaxis) (2,4)->(4,newaxis,2,newaxis)
(2,5)->(5,newaxis,newaxis,2)
问题是，由于您没有指定要对张量的b和g维执行的操作，因此它会尝试将它们一起广播，并且由于它们是不同的，因此失败了。您可以通过添加尺寸为 1 的额外尺寸来使其工作：
>>> (np.einsum('i...,j...,k...->...', B, C, D[:, None]) -
... np.einsum('i...,i...,i...->...', B, C, D[:, None]))
array([[ 1.27679643, 2.26125361, 1.32775173, 1.5517918 , 0.47083151],
[ 0.84302586, 1.57516142, 1.1335904 , 1.14702252, 0.34226837],
[ 0.70592576, 1.34187278, 1.02080112, 0.99458563, 0.29535054],
[ 1.66907981, 3.07143067, 2.09677013, 2.20062463, 0.65961165]])
如果您希望将 B 的所有轴放在 C 的所有轴之前，并且这些轴放置在 D 的所有轴之前，则以下内容似乎有效，至少就创建正确形状的输出而言，尽管您可能需要仔细检查结果是否真的是您想要的：
>>> B = np.random.rand(2, 3)
>>> C = np.random.rand(2, 4, 5)
>>> D = np.random.rand(2, 6)
>>> C_idx = (slice(None),) + (None,) * (B.ndim - 1)
>>> D_idx = C_idx + (None,) * (C.ndim - 1)
>>> (np.einsum('i...,j...,k...->...', B, C[C_idx], D[D_idx]) -
... np.einsum('i...,i...,i...->...', B, C[C_idx], D[D_idx])).shape
(3L, 4L, 5L, 6L)
编辑从注释中，如果不是每个张量的第一个轴必须缩小，而是前两个轴，那么上面的可以写成：
>>> B = np.random.rand(2, 2, 3)
>>> C = np.random.rand(2, 2, 4, 5)
>>> D = np.random.rand(2, 2, 6)
>>> C_idx = (slice(None),) * 2 + (None,) * (B.ndim - 2)
>>> D_idx = C_idx + (None,) * (C.ndim - 2)
>>> (np.einsum('ij...,kl...,mn...->...', B, C[C_idx], D[D_idx]) -
... np.einsum('ij...,ij...,ij...->...', B, C[C_idx], D[D_idx])).shape
(3L, 4L, 5L, 6L)
更一般地说，如果减少超过d指数，C_idx和D_idx将如下所示：
>>> C_idx = (slice(None),) * d + (None,) * (B.ndim - d)
>>> D_idx = C_idx + (None,) * (C.ndim - d)
对np.einsum的调用需要在索引中具有d字母，在第一次调用中是唯一的，在第二次调用中重复。
编辑 2那么C_idx和D_idx到底是怎么回事？举最后一个例子，用B、C和D形状(2, 2, 3)、(2, 2, 4, 5)和(2, 2, 6)。C_idx由两个空切片组成，加上与B减 2 的维度数一样多的Nones，因此当我们取C[C_idx]时，结果具有形状(2, 2, 1, 4, 5).类似地，D_idxC_idx加Nones 与C减去 2 的维数一样多，因此D[D_idx]的结果具有形状(2, 2, 1, 1, 1, 6)。这三个数组不会一起进行建模，但np.einsum增加了大小为 1 的额外维度，即上面错误的"重新映射"形状，因此生成的数组结果是有额外的尾随数组，并且形状如下：
(2, 2, 3, 1, 1, 1)
(2, 2, 1, 4, 5, 1)
(2, 2, 1, 1, 1, 6)
前两个轴被缩小，因此从输出中消失，在其他情况下，广播适用，其中复制大小为 1 的维度以匹配较大的维度，因此输出按我们想要(3, 4, 5, 6)。
@hpaulj提出了一种使用"Levi-Civita like"张量的方法，理论上应该更快，请参阅我对原始问题的评论。下面是一些用于比较的代码：
b, g = 5000, 2000
B = np.random.rand(2)
C = np.random.rand(2, b)
D = np.random.rand(2, g)
def calc1(b, c, d):
return (np.einsum('i,jm,kn->mn', b, c, d) -
np.einsum('i,im,in->mn', b, c, d))
def calc2(b, c, d):
return np.einsum('ijk,i,jm,kn->mn', calc2.e, b, c, d)
calc2.e = np.ones((2,2,2))
calc2.e[0, 0, 0] = 0
calc2.e[1, 1, 1] = 0
但是在运行它时：
%timeit calc1(B, C, D)
1 loops, best of 3: 361 ms per loop
%timeit calc2(B, C, D)
1 loops, best of 3: 643 ms per loop
np.allclose(calc1(B, C, D), calc2(B, C, D))
Out[48]: True
一个令人惊讶的结果，我无法解释...