Định thức của mảng đa chiều

Question

Tôi đang cố gắng để tính toán các yếu tố quyết định của một mảng NumPy M, với np.shape (M) = (N, L, L) là một cái gì đó như:

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[_]) 

Có cách mà không lặp? "N" là một số đơn đặt hàng có cường độ lớn hơn "L". Tôi nghĩ về một cái gì đó như:

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

Có cách nào nhanh hơn để làm những gì tôi muốn không? Tôi đoán chi phí vòng lặp là một nút cổ chai hiệu suất bởi vì yếu tố quyết định của ma trận nhỏ là một hoạt động tương đối rẻ (hoặc tôi sai?).

Answer 1

Tôi không thể áp dụng phương thức apply_over_axes vì ​​đây là mảng 3D mà tôi nghĩ. Dù sao, profiling mã cho thấy rằng chương trình dành ít thời gian trong vòng lặp,

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) 

Kết quả là "0,955 giây" của thời gian chạy, với 0.930s thời gian tích lũy chi tiêu trong linalg.py:1642(det).

Nếu tôi thực hiện thử nghiệm tương tự với ma trận 2x2, tôi nhận được .844 tổng thời gian và .821s trong linalg.py:1642 (det), mặc dù ma trận nhỏ. Sau đó, có vẻ như chức năng det() có chi phí lớn cho các ma trận nhỏ.

Làm điều đó với 30x30, tôi nhận được tổng thời gian 1.198 và thời gian ở det() trong 1.172.

Với 70x70, tổng là 3.122 và thời gian ở det() là 3.094, với ít hơn 1% trong vòng lặp.

Dường như trong mọi trường hợp, nó không đáng để tối ưu hóa vòng lặp python.

Answer 2

Bạn cần sửa đổi np.linalg.det để tăng tốc độ. Ý tưởng là det() là một hàm Python, nó thực hiện rất nhiều kiểm tra đầu tiên, và gọi thường trình fortran, và có một số mảng tính toán để có được kết quả.

Đây là mã từ 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) 

Để tăng tốc chức năng này, bạn có thể bỏ qua việc kiểm tra (nó trở thành trách nhiệm của bạn để giữ cho các đầu vào bên phải), và thu thập kết quả fortran trong một mảng, và thực hiện các phép tính cuối cùng cho tất cả các mảng nhỏ mà không có vòng lặp.

Dưới đây là kết quả của tôi:

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())) 

và tốc độ là:

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 

Vì vậy, bằng cách làm này, bạn có thể tăng tốc bằng 15 lần. Đó là một kết quả tốt mà không có bất kỳ mã được biên dịch nào.

lưu ý: Tôi bỏ qua kiểm tra lỗi cho thường trình fortran.

Answer 3

Tôi vừa hoàn thành mã mã VBA để tính toán Xác định ma trận bằng phương pháp Pivot.

Hy vọng điều này hữu ích

---

Option Explicit 


    Public Function Deterpivot(Matrix() As Double) As Double 

Dim j As Integer 
Dim i As Integer 
Dim k As Integer 
Dim n As Integer 
Dim x As Integer 
Dim First As Integer 
Dim Last As Integer 
Dim lastm As Integer 
Dim Row() As Double 
Dim Dvalue As Double 
Dim Divisor As Double 
Dim Multiplier As Double 
Dim Pmatrix() As Double 
Dim Nmatrix() As Double 
'This function assumes square matrix 
'with lower bound same for each dimension 
Last = UBound(Matrix, 1) 
lastm = UBound(Matrix, 1) 

'create new matrix that store value of matrix of variable, since matrix variable is of fixed type not dynamic type 
ReDim Nmatrix(1 To Last, 1 To Last) 
    For j = 1 To Last 
    For k = 1 To Last 
     Nmatrix(j, k) = Matrix(j, k) 
    Next k 
    Next j 

    If Last = 2 Then 

    Deterpivot = Det2m(Matrix()) 

    Else 
    Dvalue = 1 'initialize the DET value of matrix 
    ReDim Row(0 To lastm - 2) 
    Row(0) = 1 

     For i = 1 To lastm - 2 

      ReDim Pmatrix(1 To Last - 1, 1 To Last - 1) 

      For j = 1 To Last - 1 
       For k = 1 To Last - 1 

        Row(i) = Nmatrix(1, 1) 
        Pmatrix(j, k) = (Nmatrix(1, 1) * Nmatrix(j + 1, k + 1) - Nmatrix(j + 1, 1) * Nmatrix(1, k + 1))/Row(i - 1) 


       Next k 
      Next j 


      ReDim Nmatrix(1 To Last - 1, 1 To Last - 1) 
      For j = 1 To Last - 1 
       For k = 1 To Last - 1 
       ' if pivot point =0 then exchange row number, DET will multiply with -1 

        If (Nmatrix(1, 1) = 0 And Nmatrix(j, 1) <> 0) Then 
        MatrixRowExchange Nmatrix(), 1, j 
        Dvalue = -Dvalue 
        End If 

        Nmatrix(j, k) = Pmatrix(j, k) 
       Next k 
      Next j 

      Last = UBound(Nmatrix, 1) 


     Next i 

     Dvalue = Dvalue * Det2m(Nmatrix()) 
     Deterpivot = Dvalue/Row(lastm - 2) 

End If 


End Function 


Public Function Det2m(Matrix() As Double) As Double 

    Dim j As Integer 
    Dim i As Integer 
    Dim k As Integer 
    Dim n As Integer 

Det2m = Matrix(1, 1) * Matrix(2, 2) - Matrix(1, 2) * Matrix(2, 1) 

End Function 

Public Sub MatrixRowExchange(Matrix() As Double, rowb As Integer, rowa As Integer) 'sub to exchange row of matrix 
Dim i As Integer 
Dim j As Integer 
Dim k As Integer 
Dim evalue As Double 
    For i = LBound(Matrix, 1) To UBound(Matrix, 2) 
     For j = LBound(Matrix, 1) To UBound(Matrix, 2) 
      If i = rowb Then 

       evalue = Matrix(rowb, j) 
       Matrix(rowb, j) = Matrix(rowa, j) 
       Matrix(rowa, j) = evalue 

      End If 

     Next j 
    Next i 

End Sub 

---