Python NUMPY HUGE 矩阵乘法

Question

我需要将两个大矩阵相乘并对它们的列进行排序。

 import numpy
 a= numpy.random.rand(1000000, 100)
 b= numpy.random.rand(300000,100)
 c= numpy.dot(b,a.T)
 sorted = [argsort(j)[:10] for j in c.T]

此过程需要大量时间和内存。有没有办法加快这个过程？如果不是，我如何计算执行此操作所需的 RAM？我目前有一个带有 4GB RAM 且没有交换的 EC2 盒子。

我想知道这个操作是否可以序列化，我不必将所有内容都存储在内存中。

Answer 1

您可以做的一件事是使用优化的 BLAS 库编译 numpy，例如ATLAS、GOTO blas 或英特尔专有的 MKL。

要计算所需的内存，您需要监控 Python 的驻留集大小（“RSS”）。以下命令在 UNIX 系统上运行（准确地说是 FreeBSD，在 64 位机器上）。

> ipython

In [1]: import numpy as np

In [2]: a = np.random.rand(1000, 1000)

In [3]: a.dtype
Out[3]: dtype('float64')

In [4]: del(a)

要获取我运行的 RSS：

ps -xao comm,rss | grep python

[编辑：请参阅ps manual page 以获得对选项的完整说明，但基本上这些ps 选项使其仅显示所有进程的命令和驻留集大小。我相信 Linux 的 ps 的等效格式是 ps -xao c,r。]

结果是；

• 启动解释器后：24880 kiB
• 导入 numpy 后：34364 kiB
• 创建后a: 42200 kiB
• 删除a后：34368 kiB

计算大小；

In [4]: (42200 - 34364) * 1024
Out[4]: 8024064

In [5]: 8024064/(1000*1000)
Out[5]: 8.024064

如您所见，计算出的大小与默认数据类型 float64 的 8 个字节非常匹配。区别在于内部开销。

MiB 中原始数组的大小约为；

In [11]: 8*1000000*100/1024**2
Out[11]: 762.939453125

In [12]: 8*300000*100/1024**2
Out[12]: 228.8818359375

这还不算太糟糕。但是，点积会太大：

In [19]: 8*1000000*300000/1024**3
Out[19]: 2235.1741790771484

这是 2235 GiB！

您可以做的是将问题拆分并执行dot操作；

• 将b 加载为ndarray
• 依次从a 加载每一行作为ndarray。
• 将行乘以b 的每一列并将结果写入文件。
• del() 该行并加载下一行。

这不会让它更快，但它会减少它使用的内存！

编辑： 在这种情况下，我建议以二进制格式编写输出文件（例如，使用 struct 或 ndarray.tofile）。这将使从文件中读取列变得更加容易，例如numpy.memmap。

Answer 2

无论如何你都有问题。正如Roland Smith 在他的回答中向您展示的那样，数据量和计算量是巨大的。你可能对线性代数不是很熟悉，所以几句解释可能有助于理解（然后希望解决）这个问题。

您的数组都是长度为 100 的向量的集合。其中一个数组有 300 000 个向量，另一个有 1 000 000 个向量。这些数组之间的点积意味着您计算每个可能的向量对的点积。有 300 000 000 000 个这样的对，因此生成的矩阵是 1.2 TB 或 2.4 TB，具体取决于您使用 32 位还是 64 位浮点数。

在我的计算机上，将 (300,100) 数组与 (100,1000) 数组相乘大约需要 1 毫秒。由此推断，您正在查看 1000 秒的计算时间（取决于内核数量）。

采用点积的好处是您可以分段进行。保持输出是另一个问题。

---

如果您在自己的计算机上运行它，可以通过以下方式计算结果矩阵：

• 在磁盘上创建一个输出数组作为np.memmap 数组
• 一次计算一行结果（如Roland Smith 所述）

这将导致使用较大 (2.4 TB) 文件的线性文件写入。

这不需要太多的代码行。但是，请确保所有内容都以合适的方式转置；转置输入数组很便宜，转置输出非常昂贵。如果您可以访问彼此靠近的元素，则访问生成的巨大数组很便宜，如果您访问彼此远离的元素，则成本很高。

必须小心地对一个巨大的内存映射数组进行排序。您应该使用对连续数据块进行操作的就地排序算法。数据存储在 4 KiB 块中（512 或 1024 个浮点数），需要读取的块越少越好。

---

现在您不是在我们自己的机器上而是在云平台上运行代码，情况发生了很大变化。通常云 SSD 存储在随机访问时非常快，但 IO 很昂贵（在金钱方面也是如此）。可能最便宜的选择是计算合适的数据块并将它们发送到 S3 存储以供进一步使用。 “合适的块”部分取决于您打算如何使用数据。如果您需要处理单个列，则一次将一列或几列发送到云对象存储。

---

但是，很大程度上取决于您的排序需求。您的代码看起来好像您最终只查看每列的几个第一项。如果是这种情况，那么您应该只计算前几个项目而不是完整的输出矩阵。这样你就可以在内存中做所有事情。

也许如果你多说一点你的排序需求，可能会有一种可行的方法来做你想做的事。

哦，一件重要的事情：你的矩阵是密集的还是稀疏的？ （稀疏意味着它们大多包含 0。）如果您希望输出矩阵大部分为零，那可能会彻底改变游戏规则。

Answer 3

DrV 和 Roland Smith 说的是很好的答案；他们应该被倾听。我的回答只不过是提供了一个使您的数据稀疏的选项，一个彻底改变游戏规则的选项。

稀疏性可以非常强大。它会将您的 O(100 * 300000 * 1000000) 操作转换为具有 k 个非零元素的 O(k) 操作（稀疏仅意味着矩阵基本上为零）。我知道 DrV 提到了稀疏性并认为它不适用，但我猜它是。

所有需要做的就是找到一个sparse representation 来计算这个变换（解释结果是另一场球赛）。简单（快速）的方法包括Fourier transform 或wavelet transform（都依赖于矩阵元素之间的相似性），但这个问题可以通过several different algorithms 推广。

拥有处理此类问题的经验，这听起来像是一个相对常见的问题，通常可以通过一些巧妙的技巧来解决。在机器学习这样的领域，这些类型的问题被归类为“简单”时，情况往往如此。