5

笔记：

Speed is not as important as getting a final result.
However, some speed up over worst case is required as well.

我有一个大数组A：

A.shape=(20000,265) # or possibly larger like 50,000 x 265

我需要计算相关系数。

np.corrcoeff # internally casts the results as doubles

我只是借用了他们的代码并编写了我自己的 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))

我仍然内存不足，我正在使用一台大内存机器，264GB

有人告诉我，快速的 C 库可能正在使用将点积分解为多个部分的例程，为了优化这一点，元素的数量被填充为 2 的幂。

我真的不需要计算相关系数矩阵的对称一半。但是，我没有看到一种方法可以在合理的时间内通过 python 循环“手动”执行此操作。

有谁知道向 numpy 询问一个体面的点积例程的方法，它可以平衡内存使用与速度......？

干杯

更新：

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

发现这个：

http://wiki.scipy.org/PerformanceTips

不确定我是否遵循它....所以，请评论或提供有关此解决方案的答案、您自己的想法，或者只是对此类问题的一般评论。

TIA

编辑：我很抱歉，因为我的数组比我想象的要大得多。数组大小实际上是 151,000 x 265 我在一台 264 GB 且至少有 230 GB 可用空间的机器上内存不足。

我很惊讶对 blas dgemm 的 numpy 调用并小心使用 C 顺序数组并没有蹲下。

4

3 回答 3

4

使用 intel 的 mkl 编译的 Python 将在大约 30 秒内使用 12GB 内存运行它：

>>> A = np.random.rand(50000,265).astype(np.float32)
>>> A.dot(A.T)
array([[ 86.54410553,  64.25226593,  67.24698639, ...,  68.5118103 ,
         64.57299805,  66.69223785],
         ...,
       [ 66.69223785,  62.01016235,  67.35866547, ...,  66.66306305,
         65.75863647,  86.3017807 ]], dtype=float32)

如果您无法访问 intel 的 MKL，请下载 python anaconda并安装具有 30 天试用版的加速包，或者对于包含 mkl 编译的学者免费。其他各种 C++ BLAS 库也应该可以工作——即使它将数组从 C 复制到 F，它也不应该占用超过 30GB 的内存。

我能想到的唯一一件事就是您的安装尝试将整个 50,000 x 50,000 x 265 阵列保存在内存中，坦率地说，这非常糟糕。作为参考，float32 50,000 x 50,000 数组只有 10GB，而上述数组为 2.6TB...

如果它是一个 gemm 问题，您可以尝试一个块 gemm 公式：

def chunk_gemm(A, B, csize):
    out = np.empty((A.shape[0],B.shape[1]), dtype=A.dtype)
    for i in xrange(0, A.shape[0], csize):
        iend = i+csize
        for j in xrange(0, B.shape[1], csize):
            jend = j+csize
            out[i:iend, j:jend] = np.dot(A[i:iend], B[:,j:jend])
    return out

这会更慢，但希望能克服你的记忆问题。

于 2014-03-04T22:20:37.763 回答
2

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

cov = np.einsum('ij,kj->ik', A, A) / n

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

于 2014-03-04T19:42:18.057 回答
0

更新：结果点积完成且没有错误，但经过仔细检查，输出数组由 151,000 列中的 95,000 到末尾的零组成。也就是说，所有行的 out[:,94999] = 非零但 out[:,95000] = 0 ...

这个超级烦。。。

另一个 Blas 描述

交流中，提到了我也想过的事情……既然 blas 是 fortran，那么输入的顺序不应该是 F 顺序吗……？正如下面的 scipy 文档页面所示，C 命令。

尝试 F 顺序导致分段错误。所以我回到第一方。

原始帖子我终于找到了我的问题，就像往常一样在细节上。

我正在使用存储为 F 顺序的 np.float32 数组。据我所知，我无法控制 F 顺序，因为数据是使用成像库从图像中加载的。

import scipy
roi = np.ascontiguousarray( roi )# see roi.flags below
out = scipy.linalg.blas.sgemm(alpha=1.0, a=roi, b=roi, trans_b=True)

这个 3 级 blas 例程可以解决问题。我的问题有两个：

roi.flags
C_CONTIGUOUS : False
F_CONTIGUOUS : True
OWNDATA : True
WRITEABLE : True
ALIGNED : True
UPDATEIFCOPY : False

而且...我使用的是 blas dgemm 而不是 sgemm。“d”代表“双”，“s”代表“单”。请参阅此 pdf：BLAS 摘要 pdf

我看了一遍，不知所措……我回去阅读了关于 blas 例程的维基百科文章，以了解 3 级与其他级别：wikipedia article on blas

现在它适用于 A = 150,000 x 265，执行：

A \dot A.T

谢谢大家的想法......知道它可以完成是最重要的。

于 2014-03-11T21:31:25.000 回答