numpy相关系数:大型阵列上的np.dot(A,AT)导致seg故障

wbg*_*wbg 5 numpy blas dot-product

注意:

Speed is not as important as getting a final result.
However, some speed up over worst case is required as well. 
Run Code Online (Sandbox Code Playgroud)

我有一个大阵列A:

A.shape=(20000,265) # or possibly larger like 50,000 x 265
Run Code Online (Sandbox Code Playgroud)

我需要计算相关系数.

np.corrcoeff # internally casts the results as doubles
Run Code Online (Sandbox Code Playgroud)

我只是借用了他们的代码并写了我自己的cov/corr而不是投入双打,因为我真的只需要32位浮点数.而且我抛弃了conj(),因为我的数据总是真实的.

cov = A.dot(A.T)/n  #where A is an array of 32 bit floats
diag = np.diag(cov)
corr = cov / np.sqrt(np.mutliply.outer(d,d))
Run Code Online (Sandbox Code Playgroud)

我仍然没有内存,我正在使用264GB的大型内存机器

我被告知,快速的C库,可能正在使用一个将点积分解成碎片的例程,并且为了优化它,元素的数量被填充到2的幂.

我真的不需要计算相关系数矩阵的对称半部分.但是,我没有看到在合理的时间内使用python循环"手动"执行此操作的方法.

有没有人知道如何让numpy获得一个像样的点产品例程,这可以平衡内存使用与速度......?

干杯

更新:

有趣的是,写这些问题往往有助于我找到更好的谷歌查询的语言.

发现这个:

http://wiki.scipy.org/PerformanceTips
Run Code Online (Sandbox Code Playgroud)

不确定我是否按照它...所以,请评论或提供有关此解决方案,您自己的想法或仅对此类问题的一般评论的答案.

TIA

编辑:我很抱歉,因为我的阵列确实比我想象的要大得多.数组大小实际上是151,000 x 265我在264 GB的机器上运行内存不足,至少230 GB可用.

令我感到惊讶的是,对于使用dgemm并且小心使用C阶阵列的那种笨拙的调用没有做下蹲.

Jai*_*ime 2

您可以尝试看看是否np.einsum比 dot 更适合您的情况:

cov = np.einsum('ij,kj->ik', A, A) / n
Run Code Online (Sandbox Code Playgroud)

的内部工作原理dot有点晦涩,因为它尝试使用 BLAS 优化例程,有时需要按 Fortran 顺序复制数组,不确定这里是否是这种情况。einsum将缓冲其输入,并在可能的情况下使用矢量化 SIMD 运算,但除此之外,它基本上将运行简单的三个嵌套循环来计算矩阵乘积。