Python 大数据矩阵操作答案

【问题标题】：Python Big Data Matrix manipulationPython 大数据矩阵操作
【发布时间】：2014-11-17 08:36:40
【问题描述】：

我想我有一个大数据（N = 1e6 和维度 = 3）场景。我需要在我的代码中多次执行一些矩阵操作，例如 einsum、矩阵求逆等。为了给出一个想法，我想做如下的事情。

import numpy.random as rd

ndata, kdata = 1e6, 1e5

x = rd.normal(0,1,(ndata, kdata,3,3))

y = rd.normal(0,1,(ndata, kdata,3,3))

对于小的ndata，kdata跟随是一种高效便捷的方法，

xy =  einsum('pqrs, pqsu -> pqru', x, y )

由于我有很大的 ndata 和 kdata 上述方法成为内存绑定问题，因此下一个赌注将是在 ndata 和 kdata 上嵌套 for 循环的点积，如下所示：

xyloop1 = np.empty((ndata, kdata, 3, 3))

for j in xrange(ndata):

    for k in xrange(kdata):

        xyloop1[j,k] =  np.dot(x[j,k], y[j,k] )

鉴于我所学的 for 循环在 python 中是令人讨厌的。另外，我想利用 numpy 的好处，所以认为块矩阵方法更可取，如下所示：

nstep = 200
ndiv  = ndata/nstep   

kstep = 200
kdiv  = kdata/kstep   

xyloop2 = np.empty((ndata, kdata, 3, 3))

for j in xrange(ndiv):

    ji, jf = j*nstep, (j+1)*nstep     

    for k in xrange(kdiv):

        ki, kf = k*kstep, (k+1)*kstep     

        xyloop2[ji:jf,ki:kf] =  einsum('pqrs, pqsu -> pqru', x[ji:jf,ki:kf], y[ji:jf,ki:kf] )

另外，我需要这些 xy 或 xyloop1 或 xyloop2 来进行进一步的计算。所以我必须在每次计算后编写和阅读它。考虑到系统 I/O 的带宽，您认为最好的方法是方法 3，因为与方法 2 相比，这意味着更少的 I/O 和少量的迭代？如果您有任何其他想法或需要更多信息，请告诉我。

我是堆栈的新手，所以请对我温柔:)。任何帮助将不胜感激。顺便说一句，我正在尝试解决大数据的混合建模问题。谢谢！

【问题讨论】：

这是一个很好的问题，但最终，只有你能正确回答。 :) 分析事物并查看哪个对于您的特定问题最快。两者都是很好的方法。究竟哪种方法最快将在很大程度上取决于输入数据的大小以及您正在做什么。您还可以考虑更改您的第一个版本以沿一个轴矢量化计算并在另一个轴上循环。矢量化是内存使用和速度之间的权衡。如果内存使用存在问题，通常最好只删除嵌套循环的一部分。
感谢@JoeKington 的评论。为了清楚起见，我只提到了内循环中的点积。在这种情况下，像你所说的矢量化内循环绝对是要走的路。然而，实际上我正在做几个矩阵代数，包括逆、转置、几个时间点积等。所以我猜对内循环进行矢量化或使用列表理解不是一个好主意。关于，您评论的第一部分，我将进行分析并发布。我想补充一点，使用 cvxopt LAPACK 和 BLAS 包，方法 2 等于方法 3 的速度。

标签： python numpy matrix

【解决方案1】：

虽然我同意这样的评论，即确定了解的唯一方法是自己分析事物，但有一些指导原则可以帮助您在第一次尝试时编写高效的numpy 代码。以下是针对您的问题的一些建议：

创建一个新的 numpy 数组的开销大约是加/乘成本的 1000 倍，因此方法 2 应该是低效的，因为每次调用 np.dot 都会创建一个数组，但只执行 27 次加乘。李>
如果您要在 python 中使用缓慢的 for 循环，请尽可能在最左边的轴上执行（对于 C 有序数组）。
很难高效地编写非常通用的 N 维代码，所以我的猜测是，一系列更简单的numpy 调用将比np.einsum 更有效。试试C = np.sum(A[...,:,None] * B[...,:,:], axis=-2)（虽然这很投机）。

所以我会尝试以下方法：

xyloop2 = np.empty((ndata, kdata, 3, 3))

for i in xrange(ndata):
    xyloop2[i] = np.sum(x[i,:,:,:,None] * y[i,:,None,:,:], axis=-2)

类似于方法 2，但更简单（且更高效）的 for 循环。还把矩阵乘法换成了我认为可能更快的东西。

【讨论】：

【解决方案2】：

显然，有时 einsum 是有效的。 p,q,r,s 为 100, 50, 3, 3 的示例

示例一：

%timeit tt=np.einsum('pqrs, pqsu->pqru',x,y)
100 loops, best of 3: 3.45 ms per loop

%timeit zz= np.sum(x[:,:,:,None,:]*y[:,:,:,None],axis=-2)
10000 loops, best of 3: 153 µs per loop

示例二：

%timeit zz= np.sum(x[:,:,:,None,:]*y[:,:,:,None,None],axis=-2)
1000 loops, best of 3: 274 µs per loop

%timeit tt=np.einsum('pqrs, pqs->pqr',x,y)
10000 loops, best of 3: 151 µs per loop

np.allclose(zz,tt)
True

【讨论】：