我有维度(100、100、16、16(的numpy数组'test',它为我提供了一个不同的16x16数组,用于100x100网格上的点。 我还有一些特征值和向量,其中 vals 的维度为 (100, 100, 16( 和 vecs(100, 100, 16, 16(,其中 vecs[x, y, :, i] 将是矩阵在对应于第 i 个特征值 vals[x, y, i] 的点 (x, y( 处的第 i 个特征向量。
现在我想在网格上的所有点处获取数组的第一个特征向量,使用测试矩阵做一个矩阵积,然后在网格上的所有点上用数组的所有其他特征向量做结果向量的标量乘积并将它们求和。 生成的数组应具有维度 (100, 100(。在此之后,我想取数组的第二个特征向量,矩阵将其与测试相乘,然后取结果的标量乘积与所有不是第二个的特征向量,依此类推,这样最终我就有 16 (100, 100( 或者更确切地说是一个 (100, 100, 16( 数组。到目前为止,我只成功了很多我想避免的 for 循环,但使用 tensordot 给了我错误的维度,我不明白如何为 np.dot 函数选择矢量化的轴。 我听说 einsum 可能适合这个任务,但所有不依赖 python 循环的东西对我来说都很好。
import numpy as np
from numpy import linalg as la
test = np.arange(16*16*100*100).reshape((100, 100, 16, 16))
vals, vecs = la.eig(test + 1)
np.tensordot(vecs, test, axes=[2, 3]).shape
>>> (10, 10, 16, 10, 10, 16)
编辑:好的,所以我使用np.einsum来获得正确的中间结果。
np.einsum('ijkl, ijkm -> ijlm', vecs, test)
但在下一步中,我只想对 vec 的所有其他条目进行标量积。我可以在这个 einsum 形式主义中实现一些逆克罗内克三角洲吗?还是我现在应该切换回通常的numpy?
好的,我玩了一下,用 np.einsum 找到了一种方法来做上面描述的事情。einsum的一个很好的功能是,如果你在"输出"中重复双重出现的索引(所以在"->"-事物的右边(,你可以沿着一些轴进行逐元素乘法,沿着其他一些轴收缩(这是你在手写张量代数符号中没有的(。
result = np.einsum('ijkl, ijlm -> ijkm', np.einsum('ijkl, ijkm -> ijlm', vecs, test), vecs)
这几乎可以解决问题。现在只需要删除对角线项。我们可以通过像这样减去对角线项来做到这一点:
result = result - result * np.eye(np.shape(test)[-1])[None, None, ...]