numpy：计算大矩阵的x.T*x

Question

在 numpy 中，计算 x.T * x 的最有效方法是什么？这里的 x 是一个很大的矩阵，大小是 200,000 行和 1000 列，数据类型是 float32，而 .T 是转置操作符。

为了避免误解，结果的大小是 1000 行和 1000 列。

补充说明：在我最初的问题中，我提到 np.dot(x.T, x) 需要几个小时来计算。后来发现我的矩阵中有一些 NaNs（表示无效或缺失的值），这让 np.dot 的性能大幅下降（有人知道为什么吗？）现在这个问题已经解决，但我最初的问题依然存在。

Answer 1

嗯，x 大约是 800 兆字节，如果结果也需要差不多的内存，你确定你的电脑有足够的物理内存吗？是不是在用虚拟内存了？

除此之外，numpy 应该会使用一个叫 BLAS 的函数。虽然 numpy 默认使用的库可能比较慢，但对于这个大小的数据来说，应该还是能正常工作的。

编辑

import numpy as npy
import time

def mm_timing():
  print "   n   Gflops/s"
  print "==============="
  m = 1000
  n = 200000
  a = npy.random.rand(n, m)
  flops = (2 * float(n) - 1) * float(m)**2
  t1 = time.time()
  c = npy.dot(a.T, a)
  t2 = time.time()
  perf = flops / (t2 - t1) / 1.e9
  print "%4i" % n + "     " + "%6.3f" % perf

mm_timing()

Answer 2

首先，确保你使用的是优化过的blas/lapack，这样可以带来很大的性能提升（可能快十倍）。比如，如果你使用的是多线程的ATLAS，它能比较高效地利用你电脑的所有核心（不过你需要使用一个比较新的ATLAS版本，而且编译ATLAS的过程会比较麻烦）。

至于为什么NaN会让一切变慢：这几乎是不可避免的，处理NaN的速度比“正常”的浮点数要慢很多，尤其是在CPU层面上。具体的速度影响还和你的CPU型号、使用的指令集，以及你所用的算法和实现方式有关。你可以查看这个链接了解更多信息：http://www.cygnus-software.com/papers/x86andinfinity.html。

Answer 3

这可能不是你想要的答案，但有一种方法可以大大加快速度，那就是用显卡（gpu）代替中央处理器（cpu）。如果你有一块性能不错的显卡，它的速度会比你的cpu快很多，即使你的系统调试得很好。

如果你想和numpy很好地结合使用，可以试试theano（前提是你的显卡是nvidia的）。下面这段代码的计算在我这里只用了几秒钟（不过我有一块非常强大的显卡）：

$ THEANO_FLAGS=device=gpu0 python
Python 2.6.5 (r265:79063, Apr 16 2010, 13:57:41) 
[GCC 4.4.3] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> import theano
Using gpu device 0: GeForce GTX 480
>>> from theano import tensor as T
>>> import numpy
>>> x = numpy.ones((200000, 1000), dtype=numpy.float32)
>>> m = T.matrix() 
>>> mTm = T.dot(m.T, m)
>>> f = theano.function([m], mTm)
>>> f(x)
array([[ 200000.,  200000.,  200000., ...,  200000.,  200000.,  200000.],
       [ 200000.,  200000.,  200000., ...,  200000.,  200000.,  200000.],
       [ 200000.,  200000.,  200000., ...,  200000.,  200000.,  200000.],
       ..., 
       [ 200000.,  200000.,  200000., ...,  200000.,  200000.,  200000.],
       [ 200000.,  200000.,  200000., ...,  200000.,  200000.,  200000.],
       [ 200000.,  200000.,  200000., ...,  200000.,  200000.,  200000.]], dtype=float32)
>>> r = f(x)
>>> r.shape
(1000, 1000)

我本来想等着看看 >>> numpy.dot(x.T, x) 需要多长时间来做比较，但我等得有点无聊了……

如果你没有nvidia显卡，也可以尝试PyCuda或PyOpenCL，虽然我不太确定它们对numpy的支持是否那么简单。