Python NUMPY巨型矩阵乘法

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]

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

我想知道是否可以将此操作串行化，而不必将所有内容存储在内存中。

- pg2455

---

所以你实际上只对最小的10个标量积感兴趣？那么你可以分块迭代，并在过程中丢弃很多。 - eickenberg

那似乎是个好主意.. 我该怎么做呢？我的意思是，是否有 Pythonic 的方式，还是我必须编写自己的算法？ - pg2455

如果你只需要最小的10个，你应该看一下np.argpartition，我认为它会节省你一些时间。关于其他建议，请参考：https://dev59.com/h2w05IYBdhLWcg3w_Gqw - user2379410

看一下 Dask https://www.dask.org/ 它可以自动批处理类似于 NumPy 的操作，使它们适合内存。 - C. Yduqoli

Answer 1

你可以做的一件事是使用优化的BLAS库（如ATLAS、GOTO Blas或Intel专有的MKL）编译numpy，以加快速度。

要计算所需的内存，需要监视Python的Resident Set Size（"RSS"）。以下命令在UNIX系统上运行（准确地说，是在64位机器上的FreeBSD系统上运行）。

> 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手册页，以获取选项的完整说明，但基本上这些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这样的工具更容易从文件中读取列。

- Roland Smith

---

嘿，感谢你的帮助。你能解释一下你的bash脚本中的xao comm,rss是什么吗？我已经完全理解了内存问题。非常感谢。 - pg2455

在Python中，也可以通过调用resource.getrusage(resource.RUSAGE_SELF).ru_maxrss直接获取驻留集大小。详见https://docs.python.org/3/library/resource.html。 - vaer-k

Answer 2

What DrV和Roland Smith所说的都是好答案，应该听取。我的回答只是提供了一种使您的数据变得稀疏的选项，这是一个完全改变游戏规则的选择。

稀疏性可以非常强大。它将把您的O(100 * 300000 * 1000000)操作转换为一个O(k)操作，其中k是非零元素（稀疏性仅意味着矩阵大部分为零）。我知道DrV已经提到过稀疏性并被忽略为不适用，但我认为它可能是适用的。

需要做的就是找到计算此转换的稀疏表示（解释结果是另一回事）。简单（且快速）的方法包括傅里叶变换或小波变换（两者都依赖于矩阵元素之间的相似性），但是该问题可通过几种不同的算法进行推广。

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

- Scott

---

那真是一次大开眼界的经历。非常感谢你的分享。 - pg2455

Answer 3

您在任何情况下都会遇到问题。正如 Roland Smith 在他的答案中所示，数据量和计算次数都非常巨大。您可能对线性代数不是很熟悉，因此需要一些解释来帮助您理解（然后希望解决）这个问题。

您的数组都是长度为100的向量集合。其中一个数组有300,000个向量，另一个数组有1,000,000个向量。这两个数组之间的点积意味着您要计算每对向量的点积。这样的对数有300,000,000,000个，因此得到的矩阵无论使用32位还是64位浮点数，都是1.2 TB或2.4 TB。

在我的电脑上，将(300,100)的数组与(100,1000)的数组进行点乘需要大约1毫秒。推算一下，您需要1000秒的计算时间（取决于核心数量）。

点积的好处是可以逐步完成。然后保留输出是另一个问题。

---

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

• 将输出数组创建为np.memmap数组写入磁盘
• 按行逐个计算结果（由Roland Smith解释）

这将导致一个大文件的线性写入（2.4 TB）。

这不需要太多代码。但是，请确保以适当的方式转置一切；输入数组的转置很便宜，输出的转置非常昂贵。如果您可以访问彼此靠近的元素，则访问结果巨大的数组很容易，否则非常昂贵。

对于巨大的memmapped数组进行排序必须小心处理。您应该使用在连续数据块上操作的原地排序算法。数据存储在4 KiB块中（512或1024个浮点数），读取的块越少，效果越好。

---

现在，您不是在自己的机器上运行代码，而是在云平台上运行，情况发生了很大变化。通常，云SSD存储对随机访问非常快，但是IO代价高（也是金钱成本高）。可能最便宜的选项是计算适当大小的数据块并将其发送到S3存储以供进一步使用。 "适当大小" 的部分取决于您如何使用数据。如果您需要处理单个列，则每次将一个或几个列发送到云对象存储。

---

然而，很多事情取决于您的排序需求。您的代码看起来好像最终只查看每列的前几个项目。如果是这种情况，则应仅计算前几个项目而不是完整的输出矩阵。这样，您可以在内存中完成所有操作。

也许，如果您能更详细地说明您的排序需求，就可以找到可行的方法来实现您想要的结果。

哦，还有一件重要的事情：您的矩阵是密集的还是稀疏的？（稀疏意味着它们主要包含0。）如果您预计输出矩阵大部分都是零，那可能会完全改变游戏规则。

- DrV

---

嗨，感谢您的帮助。不幸的是，我的矩阵并不是稀疏的。我刚刚了解到np.memmap。它使用Python的内存映射对象，不允许文件>2GB。正如您所说，对于数据需要约2TB，有没有办法增加memmap的大小？ - pg2455

你能解释一下如何对数据块进行排序吗？在Python中有什么特定的库可以实现吗？据我所知，我的问题可以通过将一个向量乘以矩阵并对其进行排序，然后将其存储在文件中来解决。 - pg2455

现代版本的memmap在64位Python中没有任何实际限制，2 GB的限制是旧的，并且可能取决于操作系统。此外，如果您只对一行内的项目进行排序（如果需要，请转置矩阵），那么简单的sort可以非常高效地对该行进行排序，这样您的生活会变得更加轻松。 （我理解您需要移动矩阵中的列或行，这将更加困难。） - DrV