3

我有一个基于矩阵的问题，我认为可以使用 numpy（也许是 einsum？）在一行代码中解决（计算成本低），但无法解决。

我想知道是否有人可以提出任何建议？问题如下：

A = 4d array of shape: (5, 5, 5, 5)
M = 2d array of shape: (100, 5)

Target result, R = 2d array of shape (100, 5) where...
R[z, p] = sum for p having values of [i, j, k, l] over: A[i,j,k,l] * M[z,i] * M[z, j] * sgn[p]

where sgn[p] = -ve if p = i or j; +ve if p = k or l

举个例子：

# Enter two values in A array
A = np.zeros([5, 5, 5, 5])
A[2, 1, 3, 0] = 1
A[2, 2, 4, 1] = 1

M = np.zeros([100, 5])
M[None, :] = np.arange(5)  # fill with dummy data

R = SOLUTION_FUNC(A, M), computed as:

# For R[:, 0] there's only one term to add
# since 0 only appears once as an index for non-zero values in the A array.
# As 0 appears in second half of the index set, we keep the value positive
R[z, 0] = A[2, 1, 3, 0] * M[z, 2] * M[z, 1]

# For R[:, 1], we add both terms in A since 1 appears as an index in both non-zero values
# As 1 is in first half of the list the first time, we make this value negative
R[z, 1] = -1 * A[2, 1, 3, 0] * M[z, 2] * M[z, 1] + A[2, 2, 4, 1] * M[z, 2] * M[z,2]

# For R[z:, 2], there are three terms to add since 2 appears as a 
# non-zero index three times in the A matrix.
# The 2 always appears in first half of list, so we make all terms negative
R[z, 2] = -1* A[2, 1, 3, 0] * M[z, 2] * M[z, 1] + -1*A[2, 2, 4, 1] * M[z, 2] * M[z,2] + -1 * A[2, 2, 4, 1] * M[z, 2] * M[z,2]

# etc....
R[z, 3] = A[2, 1, 3, 0] * M[z, 2] * M[z, 1]
R[z, 4] = A[2, 2, 4, 1] * M[z, 2] * M[z, 2]

如果有帮助，A 相对稀疏（可能包含 < 20 个非零项 - 但这些可能位于任何随机位置）。

编辑：添加循环

我可以使用嵌套循环来实现这一点，但这似乎是一个糟糕的解决方案：

length = 5
A = np.zeros([length, length, length, length])
A[2, 1, 3, 0] = 1
A[2, 2, 4, 1] = 2

M = np.zeros([100, length])
M[None, :] = np.arange(length) + 1 # fill with dummy data

result = np.zeros_like(M, dtype=float)

# Many loop solution
non_zero_indexes = np.array(np.nonzero(A)).T

for z in range(len(M)):
    for non_zero_index_set in non_zero_indexes:
        for i in range(length):
            # The set indices are inspected in two halves: if i in first half, return negative value
            if i in non_zero_index_set[:2]:
                # If a number appears twice in the list, we need to add twice, so include a factor
                factor = -1 * list(non_zero_index_set).count(i)  # negative factor
                result[z, i] += factor * A[non_zero_index_set[0], non_zero_index_set[1], non_zero_index_set[2], non_zero_index_set[3]]   \
                                 * M[z, non_zero_index_set[0]] * M[z, non_zero_index_set[1]]
            elif i in non_zero_index_set[-2:]:
                factor = +1 * list(non_zero_index_set).count(i)  # positive factor
                result[z, i] += factor * A[non_zero_index_set[0], non_zero_index_set[1], non_zero_index_set[2], non_zero_index_set[3]]   \
                                 * M[z, non_zero_index_set[0]] * M[z, non_zero_index_set[1]]

我已经尝试了许多排列einsum，我认为这可能能够解决这个问题。例如形式np.einsum('ijkl,mi->mi')，但不知何故包括所有 ijkl 的总和，如上所述。

希望这是有道理的，并且任何指针（或使用 einsum / 其他数组处理路线的进一步建议）都会非常高兴。

4

0 回答 0