11

我正在尝试计算 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)

有没有更快的方法来做我想做的事？我猜循环开销是一个性能瓶颈，因为小矩阵的行列式是一个相对便宜的操作（或者我错了吗？）。

4

2 回答 2

8

您需要修改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 例程的错误检查。

于 2012-11-15T12:54:15.013 回答
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，时间为 3.094，det()循环中的时间不到 1%。

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

于 2012-11-15T12:22:56.660 回答