我想解决的实际问题是,给定一组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的数组,包含了所需的结果。不幸的是,在这个过程中,我们创建了一些潜在的巨大数组。结果
^{pr2}$是一个MXN数组。当然,最终的结果只是大小上的N。开始增加N和M的大小,我们很快就会遇到内存错误。在
如上所述,如果不需要绝对值,那么我们可以按如下方式继续,现在使用einsum()
T = np.einsum('ik,jk,j', nhat, m, np.ones(M)) / M
即使是在相当大的N和M的情况下,这种方法也能迅速工作。对于具体的问题,我需要包括abs()
,但更通用的解决方案(也许是更通用的ufunc)也会很有趣。在
这相当慢,但不会创建大的中间矩阵。在
编辑:在for循环外移动除以M。在
第二版:新观念、旧观念及相关评论。在
^{pr2}$这很快,但有时会给出不同的值,我正在研究原因,但同时可能会有所帮助。在
编辑3:啊,这是绝对值的东西。嗯
编辑4:好吧,如果你的单位向量大部分都是正值,这应该运行得更快,假设m中的向量总是正数,就像它们在你的虚拟数据中一样。在
我不认为有任何简单的方法(除了Cython之类的)来加速你的精确操作。但是你可能需要考虑你是否真的需要计算你正在计算的东西。因为如果你可以用root mean square代替绝对值的平均值,你仍然会以某种方式平均内积的大小,但是你可以在一次射击中得到:
这与执行以下操作相同:
^{pr2}$是的,这并不完全是你所要求的,但我担心这是最接近你将得到的矢量化方法。如果您决定沿着这条路走下去,请看
np.einsum
对大的N
和M
的表现:当传递太多参数和索引时,它有陷入困境的趋势。在根据一些评论,使用cython似乎是最好的方法。我愚蠢地从来没有考虑过使用赛顿。事实证明,生成工作代码相对容易。在
经过一番搜索,我把下面的cython代码放在一起。这不是最通用的代码,可能不是最好的编写方法,而且可能会变得更高效。即使如此,它只比原始问题中的
einsum()
代码慢25%,所以它不是太糟糕!它的目的是显式地处理在原始问题中创建的数组(因此是输入数组的假定模式)。尽管有警告,但它确实为原始问题提供了一个合理有效的解决方案,并且可以作为类似情况下的起点。在
相关问题 更多 >
编程相关推荐