2012-11-23 24 views
3

我有一个m * m * n的numpy数组(称之为A),我想在这个数组中找到每个子矩阵A[:,:,n]的特征值。我可以在linalg.eig()的环境中做到这一点,但是确实应该有一种方法来实现这一点。就像ufunc,但它可以处理子向量而不是单个元素。这可能吗?vectorising linalg.eig()in numpy

+2

在一个侧面说明,如果你碰巧处理对称矩阵,您可以用'np.linalg.eigh',这是显著快(再次,对于对称输入的具体但常见的情况 - 例如协方差矩阵)。当然,这并不能解决你的矢量化问题(@jorgeca是正确的,没有办法直接对它进行矢量化,如果'm'非常小,'n'非常大,它只会显着加快。 –

+1

假设'm'很小,'n'不是,假设你正在多核机器上运行它,一种可能的加速方法是使用'multiprocessing'模块:只为每个进程分配几个矩阵。 –

回答

1

特征值和特征向量的计算不能被矢量化,因为通常没有办法共享不同矩阵的工作。 np.linalg.eig(对于实际输入)只是dgeev的包装,它根据文档只接受每个调用的单个矩阵,并且计算相当昂贵,所以对于不小的矩阵,python循环的开销可以忽略不计。

虽然,如果你为许多非常小的矩阵做这个,它可能会变得太慢。有several questions与此相关,解决方案通常最终是一个编译扩展。正如一位神秘的物理学家在评论中所说的那样,以ufuncs相同的方式处理子向量和子矩阵的想法一般是有用的。这些被称为泛化ufuncs,已经在numpy的development version。我发现它更快的约8倍的形状(1000, 3, 3)的矩阵:

In [2]: np.__version__ 
Out[2]: '1.8.0.dev-dcf7cac' 

In [3]: A = np.random.rand(1000, 3, 3) 

In [4]: timeit np.linalg.eig(A) 
P100 loops, best of 3: 9.65 ms per loop 

In [5]: timeit [np.linalg.eig(Ai) for Ai in A] 
10 loops, best of 3: 74.6 ms per loop 

In [6]: a1 = np.linalg.eig(A) 

In [7]: a2 = [np.linalg.eig(Ai) for Ai in A] 

In [8]: all(np.allclose(a1[i][j], a2[j][i]) for j in xrange(1000) for i in xrange(2)) 
Out[8]: True 
+0

好吧,碰巧m 非常小。 (只有4个)否则我不会问这个问题。相信我哈哈。我正在考虑更直接地使用特征值解算器的想法。我有一些代码工作,我想我可以像我希望的那样快速地工作。但是,似乎ufuncs根本就存在,并且它们确实存在,即使是以同样的方式处理子向量和子矩阵的想法在一般情况下也是有用的。 – enigmaticPhysicist

+0

@enigmaticPhysicist你想要的行为(通用ufuncs处理子数组)已经在numpy的[开发版本](https://github.com/numpy/numpy/pull/3220)中,如果你愿意建立和使用它。 – jorgeca

+0

嘿,这很酷,豪尔赫!如果我有时间,我一定要看看。 – enigmaticPhysicist