计算以下表达式的最佳（最快）方法是什么：
\ sum_ {i \ in I} \ alpha_i \ sum_ {j \ in J} \ beta_j M [:, i，j]对于给定的numpy数组：
I，J包含索引；
形状（1，| I |）和形状（1，| J |）的系数的alpha和beta；
ndim的M = 3。
在更一般的情况下，我需要针对具有相同M数组的许多alpha，beta，I，J输入进行计算。因此，可以认为alphas nad Is具有形状（N，4），betas和Js具有形状（N，3），我需要为range（N）中的每个n计算此表达式。
先感谢您。
根据一些意见，为了使问题更清楚并添加一些上下文，以下是采用实际大小的天真的方法：
M 有形状 (500, 200000, 20)
I 有形状 (10^6, 4)
J 有形状 (10^6, 3)
alpha 有形状 (10^6, 4)
beta 有形状 (10^6, 3)
N = 10**6
M_new = np.zeros(M.shape[0], N)
for n in range(N):
for i in range(4):
for j in range(3):
M_new[:, n] += alpha[n, i] * beta[n, j] * M[:, I[n, i], J[n, j]]
因此，问题在于如何尽快计算M_new。
解决方案
到目前为止，最快的解决方案是@jdehesa使用Numba提出的解决方案。
@ Han-KwangNienhuys用其他方法介绍了速度比较。
编辑：
我首先误解了这个问题。我在下面留下了原始答案，但它没有做什么问题。
冒着听起来明显的危险，您可以随时求助于Numba：
import numpy as np
import numba as nb
# Original loop implementation
def comb_loop(m, ii, jj, alpha, beta):
n = ii.shape[0]
m_new = np.zeros((m.shape[0], n))
for col in range(n):
for i in range(4):
for j in range(3):
m_new[:, col] += alpha[col, i] * beta[col, j] * m[:, ii[col, i], jj[col, j]]
return m_new
# Numba implementation
@nb.njit(parallel=True)
def comb_nb(m, ii, jj, alpha, beta):
n = ii.shape[0]
m_new = np.empty((m.shape[0], n), m.dtype)
for col in nb.prange(n):
for row in range(m.shape[0]):
val = 0
for i in range(4):
for j in range(3):
val += alpha[col, i] * beta[col, j] * m[row, ii[col, i], jj[col, j]]
m_new[row, col] = val
return m_new
# Test
np.random.seed(0)
N = 1_000 # Reduced for testing
m = np.random.rand(500, 200_000, 20)
ii = np.random.randint(m.shape[1], size=(N, 4))
jj = np.random.randint(m.shape[2], size=(N, 3))
alpha = np.random.rand(N, 4)
beta = np.random.rand(N, 3)
# Check results match
print(np.allclose(comb_loop(m, ii, jj, alpha, beta), comb_nb(m, ii, jj, alpha, beta)))
# True
# Timings
%timeit comb_loop(m, ii, jj, alpha, beta)
# 181 ms ± 1.77 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
%timeit comb_nb(m, ii, jj, alpha, beta)
# 31.1 ms ± 2.03 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
原始错误答案
您可以使用np.einsum：
import numpy as np
def comb(alpha, beta, m):
return np.einsum('i,j,nij->n', alpha, beta, m)
# Test
np.random.seed(0)
alpha = np.random.rand(10)
beta = np.random.rand(20)
m = np.random.rand(30, 10, 20)
result = comb(alpha, beta, m)
print(result.shape)
# (30,)
本文收集自互联网，转载请注明来源。
如有侵权，请联系 [email protected] 删除。
我来说两句