Complexité de l'inversion de matrice dans numpy

Je résous des équations différentielles qui nécessitent d'inverser des matrices carrées denses. Cette inversion de matrice consomme la plupart de mon temps de calcul, donc je me demandais si j'utilisais l'algorithme le plus rapide disponible.

Mon choix actuel est numpy.linalg.inv . D'après mes données numériques, je vois qu'il évolue comme où n est le nombre de lignes, donc la méthode semble être l'élimination gaussienne.$O(n^3)$

Selon Wikipedia , il existe des algorithmes plus rapides disponibles. Est-ce que quelqu'un sait s'il existe une bibliothèque qui les implémente?

Je me demande, pourquoi numpy n'utilise-t-il pas ces algorithmes plus rapides?

(Cela devient trop long pour les commentaires ...)

Je suppose que vous devez réellement calculer un inverse dans votre algorithme. ¹ Premièrement, il est important de noter que ces algorithmes alternatifs ne sont pas réellement revendiqués comme étant plus rapides , juste qu'ils ont une meilleure complexité asymptotique (ce qui signifie que le nombre requis d'opérations élémentaires croît plus lentement). En fait, dans la pratique, elles sont en réalité (beaucoup) plus lentes que l'approche standard (pour donné ), pour les raisons suivantes:$n$

1. La notation cache une constante devant la puissance de , qui peut être astronomiquement grande - si grande que peut être beaucoup plus petite que pour tout qui peut être géré par n'importe quel ordinateur dans un avenir prévisible. (C'est le cas de l'algorithme Coppersmith – Winograd, par exemple.)$\mathcal{O}$$n$$C_1 n^3$$C_2 n^{2.x}$$n$

2. La complexité suppose que chaque opération (arithmétique) prend le même temps - mais cela est loin d'être vrai dans la pratique: multiplier un groupe de nombres avec le même nombre est beaucoup plus rapide que multiplier la même quantité de nombres différents . Cela est dû au fait que le goulot d'étranglement majeur de l'informatique actuelle consiste à mettre les données en cache, et non les opérations arithmétiques réelles sur ces données. Ainsi, un algorithme qui peut être réorganisé pour avoir la première situation (appelé compatible avec le cache ) sera beaucoup plus rapide que celui où cela n'est pas possible. (C'est le cas de l'algorithme Strassen, par exemple.)

De plus, la stabilité numérique est au moins aussi importante que les performances; et là encore, l'approche standard l'emporte généralement.

Pour cette raison, les bibliothèques hautes performances standard (BLAS / LAPACK, que Numpy appelle lorsque vous lui demandez de calculer une inverse) implémentent généralement cette approche. Bien sûr, il existe des implémentations Numpy, par exemple, de l'algorithme de Strassen, mais un algorithme réglé manuellement au niveau de l'assemblage battra fortement un algorithme écrit dans un langage de haut niveau pour toute taille de matrice raisonnable.$\mathcal{O}(n^3)$$\mathcal{O}(n^{2.x})$

---

¹ Mais je m'en voudrais de ne pas souligner que cela est très rarement vraiment nécessaire: chaque fois que vous avez besoin de calculer un produit , vous devriez plutôt résoudre le système linéaire (par exemple, en utilisant ) et utilisez place - c'est beaucoup plus stable, et cela peut être fait (en fonction de la structure de la matrice ) beaucoup plus rapidement. Si vous devez utiliser plusieurs fois, vous pouvez précalculer une factorisation de (qui est généralement la partie la plus coûteuse de la résolution) et la réutiliser plus tard.$A^{-1}b$$Ax=b$numpy.linalg.solve$x$$A$$A^{-1}$$A$

Vous devriez probablement noter que, enfouie profondément dans le code source numpy (voir https://github.com/numpy/numpy/blob/master/numpy/linalg/umath_linalg.c.src ), la routine inv tente d'appeler la fonction dgetrf à partir de votre package système LAPACK, qui effectue ensuite une décomposition LU de votre matrice d'origine. Ceci est moralement équivalent à l'élimination gaussienne, mais peut être réglé sur une complexité légèrement inférieure en utilisant des algorithmes de multiplication matricielle plus rapides dans un BLAS haute performance.

Si vous suivez cette route, vous devez être averti que forcer toute la chaîne de bibliothèques à utiliser la nouvelle bibliothèque plutôt que celle du système fournie avec votre distribution est assez complexe. Une alternative sur les systèmes informatiques modernes est de regarder les méthodes parallélisées en utilisant des packages comme scaLAPACK ou (dans le monde python) petsc4py. Cependant, ils sont généralement plus heureux d'être utilisés comme solveurs itératifs pour les systèmes d'algèbre linéaire que appliqués aux méthodes directes et le PETSc en particulier cible les systèmes clairsemés plus que les systèmes denses.