Multidimensional array determinant

Question

I am trying to calculate the numpy M array identifier, and np.shape (M) = (N, L, L) is something like:

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

Is there a way without a loop? "N" is several orders of magnitude greater than "L". I thought of something like:

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

Is there a faster way to do what I want? I believe that loop overhead is a performance bottleneck since the small matrix determinant is a relatively cheap operation (or am I wrong?).

Answer 1

You need to modify np.linalg.det to get speed. The idea is that det() is a Python function, it first does a lot of checking and calls the fortran procedure and computes some array to get the result.

Here is the code from 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) 

To speed up this function, you can omit the check (it becomes your responsibility to keep the input to the right) and collect the fortran results in an array and perform the final calculations for all small arrays without a loop.

Here is my result:

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

and speed:

 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 

So, having done this, you can speed it up 15 times. This is a good result without compiled code.

Note: I skipped error checking for the fortran procedure.

Answer 2

I could not get apply_over_axes to work because it is a 3D array, I think. In any case, profiling the code shows that the program spends little time in the loop,

 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) 

The result is a “0.955 second” run time, from 0.930 with the total time spent in linalg.py:1642(det).

If I perform the same test with 2x2 matrices, I get .844s of the total time and .821s in linalg.py:1642(det), even though the matrices are small. Then, it seems that the function is det() , which has large overhead for small matrices.

By doing this with 30x30, I get the total time of 1.198s and the time in det() from 1.172.

With 70x70, the total is 3.122, and the time in det() is 3.094, with less than 1% in the loop.

It seems like you shouldn't optimize the python loop anyway.

Answer 3

I just finished VBA code to calculate matrix determinant using the Pivot method.

Hope this is helpful

---

  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 

---