23

我正在寻找最节省内存的方法来计算复杂的 numpy ndarray 的绝对平方值

arr = np.empty((250000, 150), dtype='complex128')  # common size

我还没有找到一个完全可以做到的 ufunc np.abs()**2。

由于这种大小和类型的数组占用大约半 GB，我正在寻找一种主要节省内存的方法。

我也希望它是可移植的，所以理想情况下是一些 ufunc 的组合。

到目前为止，我的理解是这应该是最好的

result = np.abs(arr)
result **= 2

它将不必要地计算(**0.5)**2，但应该**2就地计算。总共峰值内存需求只是原始数组大小 + 结果数组大小，应该是 1.5 * 原始数组大小，因为结果是真实的。

如果我想摆脱无用的**2电话，我必须做这样的事情

result = arr.real**2
result += arr.imag**2

但如果我没记错的话，这意味着我必须为实部和虚部计算分配内存，因此峰值内存使用量将是 2.0 * 原始数组大小。这些arr.real属性还返回一个不连续的数组（但这不太重要）。

有什么我想念的吗？有没有更好的方法来做到这一点？

编辑1：很抱歉没有说清楚，我不想覆盖arr，所以我不能用它。

4

5 回答 5

10

感谢numba.vectorize最近版本的 numba，为该任务创建一个 numpy 通用函数非常容易：

@numba.vectorize([numba.float64(numba.complex128),numba.float32(numba.complex64)])
def abs2(x):
    return x.real**2 + x.imag**2

在我的机器上，与创建中间数组的纯 numpy 版本相比，我发现速度提高了三倍：

>>> x = np.random.randn(10000).view('c16')
>>> y = abs2(x)
>>> np.all(y == x.real**2 + x.imag**2)   # exactly equal, being the same operation
True
>>> %timeit np.abs(x)**2
10000 loops, best of 3: 81.4 µs per loop
>>> %timeit x.real**2 + x.imag**2
100000 loops, best of 3: 12.7 µs per loop
>>> %timeit abs2(x)
100000 loops, best of 3: 4.6 µs per loop
于 2016-06-15T22:01:59.590 回答
4

编辑：此解决方案的最低内存要求是最低内存要求的两倍，而且速度稍快。不过，评论中的讨论很有参考价值。

这是一个更快的解决方案，结果存储在res：

import numpy as np
res = arr.conjugate()
np.multiply(arr,res,out=res)

我们利用复数的 abs 的性质，即abs(z) = sqrt(z*z.conjugate)，使得abs(z)**2 = z*z.conjugate

于 2015-05-25T13:32:36.147 回答
1

如果您的主要目标是节省内存，NumPy 的 ufunc 采用一个可选out参数，可让您将输出定向到您选择的数组。当您想就地执行操作时，它会很有用。

如果您对第一个方法进行了这个小修改，那么您可以arr完全就地执行操作：

np.abs(arr, out=arr)
arr **= 2

一种只使用一点额外内存的复杂方法可能是arr就地修改，计算新的实数值数组，然后恢复arr。

这意味着存储有关符号的信息（除非您知道您的复数都有正实部和虚部）。每个实数或虚数的符号只需要一个位，因此这使用了（除了您创建的新浮点数组之外）1/16 + 1/16 == 1/8的内存。arr

>>> signs_real = np.signbit(arr.real) # store information about the signs
>>> signs_imag = np.signbit(arr.imag)
>>> arr.real **= 2 # square the real and imaginary values
>>> arr.imag **= 2
>>> result = arr.real + arr.imag
>>> arr.real **= 0.5 # positive square roots of real and imaginary values
>>> arr.imag **= 0.5
>>> arr.real[signs_real] *= -1 # restore the signs of the real and imagary values
>>> arr.imag[signs_imag] *= -1

以存储符号位为代价，arr保持不变并result保持我们想要的值。

于 2015-05-25T12:36:11.267 回答
0

arr.real并且arr.imag只是复杂数组的视图。所以没有分配额外的内存。

于 2015-05-25T12:15:57.883 回答
0

如果你不想要sqrt（应该比乘法重得多），那么不要abs。

如果你不想要双内存，那么没有real**2 + imag**2

然后你可以试试这个（使用索引技巧）

N0 = 23
np0 = (np.random.randn(N0) + 1j*np.random.randn(N0)).astype(np.complex128)
ret_ = np.abs(np0)**2
tmp0 = np0.view(np.float64)
ret0 = np.matmul(tmp0.reshape(N0,1,2), tmp0.reshape(N0,2,1)).reshape(N0)
assert np.abs(ret_-ret0).max()<1e-7

无论如何，我更喜欢numba解决方案

于 2021-12-08T14:20:29.813 回答