gpt4 book ai didi

python - numpy 中没有临时项的外积的非平凡总和

转载 作者:太空狗 更新时间:2023-10-29 21:53:14 26 4
gpt4 key购买 nike

我想解决的实际问题是,给定一组 N 单位向量和另一组 M 向量,计算每个单位向量的平均值它与每个 M 向量的点积的绝对值。本质上,这是计算两个矩阵的外积,然后求和并取中间的绝对值取平均值。

对于 NM 不是太大,这并不难,有很多方法可以继续(见下文)。问题是当 NM 很大时,创建的临时文件很大,并为所提供的方法提供了实际限制。可以在不创建临时对象的情况下完成此计算吗?我遇到的主要困难是由于绝对值的存在。是否存在用于“线程化”此类计算的通用技术?

例如考虑以下代码

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)

是一个 M X N 数组。当然,最终结果只有 N 大小。开始增加 NM 的大小,我们很快就会遇到内存错误。

如上所述,如果不需要绝对值,那么我们可以按如下方式进行,现在使用 einsum()

T = np.einsum('ik,jk,j', nhat, m, np.ones(M)) / M

即使对于非常大的 NM ,这也能快速工作。对于我需要包含 abs() 的特定问题,但更通用的解决方案(可能是更通用的 ufunc)也很有趣。

最佳答案

根据一些评论,似乎使用 cython 是最好的方法。我愚蠢地从未研究过使用 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

关于python - numpy 中没有临时项的外积的非平凡总和,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/17624166/

26 4 0
Copyright 2021 - 2024 cfsdn All Rights Reserved 蜀ICP备2022000587号
广告合作:1813099741@qq.com 6ren.com