I am solving differential equations that require inverting dense square matrices, and I wanted to know what algorithm(s) do numpy and scipy use to calculate matrix inverses?
$\begingroup$
$\endgroup$
1
2 Answers
$\begingroup$
$\endgroup$
2
Documentation to numpy.linalg.inv and scipy.linalg.inv does not mention the algorithm used.
Judging from the source, numpy.linalg.inv calls a C++ wrapper to LAPACK gesv subroutine, which uses LU decomposition (more specifically, LU with partial pivoting, $A = P L U$), and scipy.linalg.inv uses slightly more low-level LAPACK computational routines, performing the same decomposition.
answered Dec 9, 2023 at 20:02
-
3$\begingroup$ Note that there are various different ways to obtain an inverse from the PLU decomposition of a matrix, so even the step following PLU is not completely trivial. Higham's "Accuracy and stability of numerical algorithms" has an entire chapter on it. $\endgroup$ Commented Dec 9, 2023 at 20:14
-
1$\begingroup$ This is also the default algorithm choice for numpy and scipy's
solvefunction, so that can be used instead for solving linear systems. $\endgroup$ Commented Dec 9, 2023 at 23:04
$\begingroup$
$\endgroup$
4
First of all, they don't take matrix inverses. They perform linear solves. Typically this is done via LU factorization.
-
5$\begingroup$ Sure they do compute inverses, if you call
inv. $\endgroup$ Commented Dec 9, 2023 at 20:12 -
1
-
$\begingroup$ I was sure that the documentation page for
numpy.linalg.invhad a comment recommending to usenumpy.linalg.solveinstead whenever possible, but apparently I remember wrong. $\endgroup$– StefCommented Dec 10, 2023 at 18:27 -
$\begingroup$ They use the exact same algorithm under the hood but obviously less work is done in using
solvecompared toinvbecause you just do fewer operations. $\endgroup$ Commented Dec 11, 2023 at 17:25
Rlibraries which use the Crout method and the Doolittle method of LU decomposition. $\endgroup$