没有临时装饰的外层产品的小数目
我希望解决的实际问题是，给定一组N个单位矢量，并且另一组M矢量针对每个单位矢量计算其与每个M矢量的点积的绝对值的平均值。 基本上，这是计算两个矩阵的外积，并求中值之间的绝对值求和和求平均值。
对于N和M不太大，这并不难，有很多方法可以继续（见下文）。 问题是，当N和M很大时，创建的临时对象是巨大的，并为所提供的方法提供了实际的限制。 这个计算可以完成而不需要创建临时对象吗？ 我所遇到的主要困难是由于绝对值的存在。 是否有“线程化”这些计算的通用技术？
作为例子，请考虑以下代码
N = 7
M = 5
# Create the unit vectors, just so we have some examples,
# this is not meant to be elegant
phi = np.random.rand(N)*2*np.pi
ctheta = np.random.rand(N)*2 - 1
stheta = np.sqrt(1-ctheta**2)
nhat = np.array([stheta*np.cos(phi), stheta*np.sin(phi), ctheta]).T
# Create the other vectors
m = np.random.rand(M,3)
# Calculate the quantity we desire, here using broadcasting.
S = np.average(np.abs(np.sum(nhat*m[:,np.newaxis,:], axis=-1)), axis=0)
这很棒，S现在是一个长度为N的数组，并包含所需的结果。 不幸的是，在这个过程中，我们创建了一些潜在巨大的数组 的结果
np.sum(nhat*m[:,np.newaxis,:], axis=-1)
是一个MXN数组。 当然，最终的结果只有N的大小。开始增加N和M的大小，我们很快就会遇到内存错误。
如上所述，如果绝对值不是必需的，那么我们可以按照以下步骤进行操作，现在使用einsum()
T = np.einsum('ik,jk,j', nhat, m, np.ones(M)) / M
即使对于相当大的N和M，这种方法也可以快速工作。 对于具体的问题，我需要包括abs()但更通用的解决方案（也许更通用的ufunc）也会引起人们的兴趣。
根据一些评论，似乎使用cython是最好的选择。 我愚蠢地从未考虑过使用cython。 事实证明，生成工作代码相对容易。
经过一番搜索，我把以下的代码放在一起。 这不是最通用的代码，可能不是编写它的最佳方式，而且可能会使效率更高。 即使如此，它只比原始问题中的einsum()代码慢25％左右，所以它并不算太坏！ 它已被写入明确使用在原始问题中创建的数组（因此是输入数组的假定模式）。
尽管有这样的警告，它确实为原始问题提供了一个合理有效的解决方案，并且可以作为类似情况的起点。
import numpy as np
cimport numpy as np
import cython
DTYPE = np.float64
ctypedef np.float64_t DTYPE_t
cdef inline double d_abs (double a) : return a if a >= 0 else -a
@cython.boundscheck(False)
@cython.wraparound(False)
def process_vectors (np.ndarray[DTYPE_t, ndim=2, mode="fortran"] nhat not None,
np.ndarray[DTYPE_t, ndim=2, mode="c"] m not None) :
if nhat.shape[1] != m.shape[1] :
raise ValueError ("Arrays must contain vectors of the same dimension")
cdef Py_ssize_t imax = nhat.shape[0]
cdef Py_ssize_t jmax = m.shape[0]
cdef Py_ssize_t kmax = nhat.shape[1] # same as m.shape[1]
cdef np.ndarray[DTYPE_t, ndim=1] S = np.zeros(imax, dtype=DTYPE)
cdef Py_ssize_t i, j, k
cdef DTYPE_t val, tmp
for i in range(imax) :
val = 0
for j in range(jmax) :
tmp = 0
for k in range(kmax) :
tmp += nhat[i,k] * m[j,k]
val += d_abs(tmp)
S[i] = val / jmax
return S
我不认为有什么简单的方法（Cython之类的）来加速你的确切操作。 但是你可能想要考虑你是否真的需要计算你正在计算的东西。 因为如果不是绝对值的平均值可以使用均方根，那么您仍然可以以某种方式对内部产品的大小进行平均，但是您可以将其作为一个单独的数字来计算：
rms = np.sqrt(np.einsum('ij,il,kj,kl,k->i', nhat, nhat, m, m, np.ones(M)/M))
这和做一样：
rms_2 = np.sqrt(np.average(np.einsum('ij,kj->ik', nhat, m)**2, axis=-1))
是的，这不完全是你要求的，但恐怕它与你用矢量化方法得到的结果差不多。 如果你决定走这条路，看看np.einsum对于大的N和M表现如何：当通过太多的参数和指数时，它倾向于np.einsum 。
这比较慢，但不会创建大的中间矩阵。
vals = np.zeros(N)
for i in xrange(N):
u = nhat[i]
for v in m:
vals[i]+=abs(np.dot(u,v))
vals[i]=vals[i]/M
编辑：移动除以for循环外的M.
编辑2：新观点，为后人保留旧观点和相关评论。
m2 = np.average(m,0)
vals = np.zeros(N)
for i in xrange(N):
u=nhat[i]
vals[i]=abs(np.dot(u,m2))
这很快，但有时会给出不同的值，我正在研究为什么，但也许它可以帮助在同一时间。
编辑3：啊，这是绝对有价值的事情。 嗯
>>> S
array([ 0.28620962, 0.65337876, 0.37470707, 0.46500913, 0.49579837,
0.29348924, 0.27444208, 0.74586928, 0.35789315, 0.3079964 ,
0.298353 , 0.42571445, 0.32535728, 0.87505053, 0.25547394,
0.23964505, 0.44773271, 0.25235646, 0.4722281 , 0.33003338])
>>> vals
array([ 0.2099343 , 0.6532155 , 0.33039334, 0.45366889, 0.48921527,
0.20467291, 0.16585856, 0.74586928, 0.31234917, 0.22198642,
0.21013519, 0.41422894, 0.26020981, 0.87505053, 0.1199069 ,
0.06542492, 0.44145805, 0.08455833, 0.46824704, 0.28483342])
time to compute S: 0.000342130661011 seconds
time to compute vals: 7.29560852051e-05 seconds
编辑4：好吧，如果你的单位向量大多数都是正值，那么它应该运行得更快，假设m中的向量总是正值，就像它们在虚拟数据中一样。
m2 = np.average(m,0)
vals = np.zeros(N)
for i in xrange(N):
u=nhat[i]
if u[0] >= 0 and u[1] >= 0 and u[2] >= 0:
vals[i] = abs(np.dot(u,m2))
else:
for j in xrange(M):
vals[i]+=abs(np.dot(u,m[j]))
vals[i]/=M
上一篇: trivial sums of outer products without temporaries in numpy