11

np.shape(M) = (N, L, L) を使用して、numpy 配列 M の行列式を計算しようとしています。次のようなものです。

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

これは 3D 配列だと思うので、apply_over_axes を機能させることができませんでした。とにかく、コードのプロファイリングは、プログラムがループでほとんど時間を費やしていないことを示しています。

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 行列で同じテストを実行すると、行列が小さいにもかかわらず、合計時間は .844 秒、linalg.py:1642(det) では .821 秒になります。すると、det()小さな行列に対してオーバーヘッドが大きい関数のようです。

30x30 でそれを行うと、合計時間は 1.198 秒、時間det()は 1.172 秒になります。

70x70 の場合、合計は 3.122 で、時間det()は 3.094 で、ループ内の 1% 未満です。

いずれにしても、python ループを最適化する価値はないようです。

于 2012-11-15T12:22:56.660 に答える