多维数组的决定因素

Question

我试图计算numpy数组M的行列式,np.shape(M)=(N,L,L)是这样的:

import numpy as np

M = np.random.rand(1000*10*10).reshape(1000, 10, 10)
dm = np.zeros(1000)
for _ in xrange(len(dm)):
    dm[_] = np.linalg.det(M[_])

有没有循环的方法？"N"比"L"大一些数量级.我想到了类似的东西:

np.apply_over_axes(np.linalg.det(M), axis=0) 

做我想要的更快的方式吗？我想循环开销是一个性能瓶颈,因为小矩阵的行列式是一个相对便宜的操作(或者我错了？).

Answer 1

您需要修改np.linalg.det以获得速度.这个想法是det()一个Python函数,它先做了很多检查,然后调用fortran例程,然后做一些数组计算得到结果.

这是来自numpy的代码:

def slogdet(a):
    a = asarray(a)
    _assertRank2(a)
    _assertSquareness(a)
    t, result_t = _commonType(a)
    a = _fastCopyAndTranspose(t, a)
    a = _to_native_byte_order(a)
    n = a.shape[0]
    if isComplexType(t):
        lapack_routine = lapack_lite.zgetrf
    else:
        lapack_routine = lapack_lite.dgetrf
    pivots = zeros((n,), fortran_int)
    results = lapack_routine(n, n, a, n, pivots, 0)
    info = results['info']
    if (info < 0):
        raise TypeError, "Illegal input to Fortran routine"
    elif (info > 0):
        return (t(0.0), _realType(t)(-Inf))
    sign = 1. - 2. * (add.reduce(pivots != arange(1, n + 1)) % 2)
    d = diagonal(a)
    absd = absolute(d)
    sign *= multiply.reduce(d / absd)
    log(absd, absd)
    logdet = add.reduce(absd, axis=-1)
    return sign, logdet

def det(a):
    sign, logdet = slogdet(a)
    return sign * exp(logdet)

要加速此功能,您可以省略检查(保持输入正确的责任),并在数组中收集fortran结果,并在没有for循环的情况下对所有小数组进行最终计算.

这是我的结果:

import numpy as np
from numpy.core import intc
from numpy.linalg import lapack_lite

N = 1000
M = np.random.rand(N*10*10).reshape(N, 10, 10)

def dets(a):
    length = a.shape[0]
    dm = np.zeros(length)
    for i in xrange(length):
        dm[i] = np.linalg.det(M[i])
    return dm

def dets_fast(a):
    m = a.shape[0]
    n = a.shape[1]
    lapack_routine = lapack_lite.dgetrf
    pivots = np.zeros((m, n), intc)
    flags = np.arange(1, n + 1).reshape(1, -1)
    for i in xrange(m):
        tmp = a[i]
        lapack_routine(n, n, tmp, n, pivots[i], 0)
    sign = 1. - 2. * (np.add.reduce(pivots != flags, axis=1) % 2)
    idx = np.arange(n)
    d = a[:, idx, idx]
    absd = np.absolute(d)
    sign *= np.multiply.reduce(d / absd, axis=1)
    np.log(absd, absd)
    logdet = np.add.reduce(absd, axis=-1)
    return sign * np.exp(logdet)

print np.allclose(dets(M), dets_fast(M.copy()))

而且速度是:

timeit dets(M)
10 loops, best of 3: 159 ms per loop

timeit dets_fast(M)
100 loops, best of 3: 10.7 ms per loop

因此,通过这样做,您可以加速15倍.没有任何编译代码,这是一个很好的结果.

注意:我省略了fortran例程的错误检查.

Answer 2

我无法让 apply_over_axes 工作，因为我认为这是一个 3D 数组。不管怎样，分析代码表明程序在循环中花费的时间很少，

import cProfile
import pstats
N=10000
M = np.random.rand(N*10*10).reshape(N, 10, 10)
def f(M):
    dm = np.zeros(N)
    for _ in xrange(len(dm)):
        dm[_] = np.linalg.det(M[_])
    return dm
cProfile.run('f(M)','foo')
p = pstats.Stats('foo')
res = p.sort_stats('cumulative').print_stats(10)

结果是运行时间为“0.955 秒”，linalg.py:1642(det) 中花费了 0.930 秒的累积时间。

如果我使用 2x2 矩阵执行相同的测试，我会得到 0.844 秒的总时间，以及 linalg.py:1642(det) 中的 0.821 秒，尽管矩阵很小。然后，似乎这个det()函数对于小矩阵来说有很大的开销。

使用 30x30 进行此操作，我得到的总时间为 1.198 秒，时间为det()1.172。

对于 70x70，总计为 3.122，时间det()为 3.094，循环中的时间少于 1%。

看来无论如何，都不值得优化python循环。