在不到一秒内计算巨大矩阵的范数：NUMPY，PYTHON

Question

我有一个很大的矩阵，大小是（1,000,000，100），也就是有一百万行和一百列。我想计算每一行的“范数”，也就是每一行的长度。目前我使用的方法是：

import numpy
a =  numpy.random.rand((1000000,100))
b = numpy.linalg.norm(a, axis =1) # this takes 2.78 seconds.

有没有什么办法可以让计算这些范数的时间更短呢？

Answer 1

简单来说，普通的规范化是一个受内存限制的问题，用numpy处理起来效率不高。要提高效率，使用numpy的办法就是进行一些数据分块，避免多次读取同样的数据：

import numpy as np
a =  np.random.rand(1000000,100)
print np.linalg.norm(a, axis =1).shape

def g(d, out=None):
    bs = 2000
    if out is None:
        r = np.empty(d.shape[0])
    else:
        r = out
    for i in range(0, d.shape[0], bs):
        u = min(i + bs, d.shape[0])
        r[i:u] = np.linalg.norm(d[i:u], axis=1)
    return r


print  (g(a) == numpy.linalg.norm(a, axis =1)).all()
print "blocked"
%timeit -n 10 g(a)
print "normal"
%timeit -n 10 numpy.linalg.norm(a, axis =1) 

在我的机器上（DDR2内存）使用numpy 1.9，这样做能带来一些小的提升：

blocked
10 loops, best of 3: 294 ms per loop
normal
10 loops, best of 3: 561 ms per loop

理论上，增加一层使用线程的分块应该能进一步提高性能，但在我的机器上（glibc 2.19），这会导致过多的页面错误，因为glibc在多线程模式下频繁调整内存，这样反而没有帮助。

strace -f  -e madvise ipython test_threaded.ipy 2>&1  | grep MADV_DONTNEED -c
14228

Answer 2

使用更小的 dtype 类型;

In [8]: a = numpy.random.rand(1000,1000)

In [9]: %timeit numpy.linalg.norm(a, axis=1)
100 loops, best of 3: 10.4 ms per loop

In [11]: b = numpy.asarray(a, numpy.float32)

In [12]: %timeit numpy.linalg.norm(b, axis=1)
100 loops, best of 3: 3.5 ms per loop

根据你使用的numpy版本，使用一个优化过的或并行化的LAPACK库可能也会有所帮助。

在一台运行FreeBSD 10 amd64的Intel Core2 Quad (Q9300)上使用numpy，我得到了：

In [14]: a =  numpy.random.rand(1000000,100)

In [15]: %timeit -n 10 numpy.linalg.norm(a, axis =1)
10 loops, best of 3: 1.08 s per loop

还有一种加速的方法，前提是你有一台多核的机器;

• 把数组写入一个文件
• 写一个函数，从这个文件中读取某一行，计算所有数字的平方和，并返回结果。
• 使用 multiprocessing.Pool.map() 在多个进程中运行这个函数。
• 当 map() 返回时，把结果列表中的数字加起来，然后对这个和取平方根。

加速效果大约等于可用核心的数量。