Determinante della matrice multidimensionale

Question

sto cercando di calcolare il determinante di una matrice NumPy M, con np.shape (M) = (N, L, L) è qualcosa di simile:

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

Esiste un modo senza loop? "N" è un ordine di grandezza maggiore di "L". Ho pensato a qualcosa tipo:

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

C'è un modo più veloce di fare quello che voglio? Immagino che il sovraccarico del loop sia un collo di bottiglia delle prestazioni perché il determinante delle piccole matrici è un'operazione relativamente economica (o sbaglio?).

Answer 1

Non ho potuto applicare apply_over_axes perché questa è una matrice 3D, penso. In ogni caso, il profiling del codice mostra che il programma spende poco tempo nel ciclo,

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) 

Il risultato è "0,955 secondi" di runtime, con 0.930s di tempo cumulativo trascorso in linalg.py:1642(det).

Se eseguo lo stesso test con le matrici 2x2, ottengo 0,844 di tempo totale e .821 in linalg.py:1642(det), nonostante le matrici siano piccole. Quindi, sembra che la funzione det() abbia un overhead elevato per le piccole matrici.

Facendolo con 30x30, ottengo un tempo totale di 1.198s e un tempo in det() di 1.172.

Con 70x70, il totale è 3.122 e il tempo in det() è 3.094, con meno dell'1% nel ciclo.

Sembra che in ogni caso, non vale la pena ottimizzare il loop python.

Answer 2

È necessario modificare np.linalg.det per ottenere la velocità. L'idea è che det() è una funzione Python, fa un sacco di controlli prima e chiama la routine fortran e calcola alcuni array per ottenere il risultato.

Ecco il codice da 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) 

per accelerare questa funzione, è possibile omettere il controllo (diventa vostra responsabilità di mantenere l'ingresso a destra), e raccogliere i risultati FORTRAN in una matrice, e eseguire i calcoli finali per tutti i piccoli array senza ciclo.

Ecco il mio risultato:

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

e la velocità è:

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 

Quindi, in questo modo, è possibile velocizzare di 15 volte. Questo è un buon risultato senza alcun codice compilato.

note: ometto il controllo degli errori per la routine fortran.

Answer 3

Ho appena terminato di codificare il codice VBA per calcolare il Determinante della matrice utilizzando il metodo Pivot.

Spero che questo utile

---

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 

---