Comment trouver le déterminant de la matrice de grande taille

Bien, pas beaucoup d'entre nous qui travaillent dans le domaine considèrent 18x18 comme une matrice de grande taille et presque toute la technique que vous choisissez, il doit être assez rapide sur n'importe quel ordinateur moderne. Ni beaucoup d'entre nous aborder les questions de matrice avec des algorithmes récursifs, beaucoup plus enclins à utiliser itératif, mais qui pourrait être un reflet du fait que beaucoup de personnes qui travaillent sur la matrice des problèmes sont des scientifiques et des ingénieurs informaticiens.

Je vous suggère de regarder les Recettes ou Numérique en C++. Pas forcément le meilleur code que vous trouverez, mais c'est un texte pour l'étude et l'apprentissage. Pour mieux les codes, BOOST a une bonne réputation et il y a toujours BLAS et des choses comme l'Intel Math Kernel Library ou AMD de Base en Mathématiques de la Bibliothèque. Je pense que toutes ces implémentations de déterminant-trouver des routines qui va s'attaquer à un 18x18 matrice très rapidement.

OriginalL'auteur High Performance Mark

La fonction det_recursive travaille pour une matrice carrée de taille quelconque. Toutefois, depuis qu'elle est récursive à l'aide de méthode naïve de l'expansion de la formule de Laplace, il est très lent pour des matrices de grande taille.

Une autre technique consiste à transformer la matrice dans un haut de forme triangulaire à l'aide de l'élimination de gauss technique. Alors le déterminant de la matrice est simplement les produits des éléments de la diagonale de la forme triangulaire de transformer la forme de la matrice d'origine.

Fondamentalement numpy est le plus rapide mais en interne, il utilise une sorte de linéaire de la matrice de transformation méthode similaire à ce que l'élimination de gauss. Cependant, je ne suis pas sûr de ce qu'il est!

In[1]
import numpy as np
In[2]
mat = np.random.rand(9,9)
print("numpy asnwer = ", np.linalg.det(mat))
Out[2] 
numpy asnwer =  0.016770106020608373
In[3]
def det_recursive(A):
if A.shape[0] != A.shape[1]:
raise ValueError('matrix {} is not Square'.format(A))
sol = 0
if A.shape != (1,1):
for i in range(A.shape[0]):
sol = sol +  (-1)**i * A[i, 0] * det_recursive(np.delete(np.delete(A, 0, axis= 1), i, axis= 0))
return sol
else:
return A[0,0]
​
In[4]
print("recursive asnwer = ", det_recursive(mat))
Out[4]
recursive asnwer =  0.016770106020608397
In[5]
def det_gauss_elimination(a,tol=1.0e-9):
"""
calculate determinant using gauss-elimination method
"""
a = np.copy(a)
assert(a.shape[0] == a.shape[1])
n = a.shape[0]
# Set up scale factors
s = np.zeros(n)
mult = 0
for i in range(n):
s[i] = max(np.abs(a[i,:])) # find the max of each row
for k in range(0, n-1): #pivot row
# Row interchange, if needed
p = np.argmax(np.abs(a[k:n,k])/s[k:n]) + k 
if abs(a[p,k]) < tol: 
print("Matrix is singular")
return 0
if p != k: 
a[[k,p],:] = a[[p, k],:] 
s[k],s[p] = s[p],s[k]
mult = mult + 1
​
# convert a to upper triangular matrix
for i in range(k+1,n):
if a[i,k] != 0.0: # skip if a(i,k) is already zero
lam = a [i,k]/a[k,k] 
a[i,k:n] = a[i,k:n] - lam*a[k,k:n]
​
deter = np.prod(np.diag(a))* (-1)**mult   
return deter
In[6] 
print("gauss elimination asnwer = ", det_gauss_elimination(mat))
Out[6] 
gauss elimination asnwer =  0.016770106020608383
In[7] 
print("numpy time")
%timeit -n3 -r3 np.linalg.det(mat)
print("\nrecursion time")
%timeit -n3 -r3 det_recursive(mat)
print("\ngauss_elimination time")
%timeit -n3 -r3 det_gauss_elimination(mat)
Out[7]
numpy time
40.8 µs ± 17.2 µs per loop (mean ± std. dev. of 3 runs, 3 loops each)
recursion time
10.1 s ± 128 ms per loop (mean ± std. dev. of 3 runs, 3 loops each)
gauss_elimination time
472 µs ± 106 µs per loop (mean ± std. dev. of 3 runs, 3 loops each)

OriginalL'auteur Khalil Al Hooti