# What algorithm(s) do numpy and scipy use to calculate matrix inverses?

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?

• By way of comparison, there are R libraries which use the Crout method and the Doolittle method of LU decomposition. Commented Dec 10, 2023 at 22:13

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.

• 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. Commented Dec 9, 2023 at 20:14
• This is also the default algorithm choice for numpy and scipy's solve function, so that can be used instead for solving linear systems. Commented Dec 9, 2023 at 23:04

First of all, they don't take matrix inverses. They perform linear solves. Typically this is done via LU factorization.

• Sure they do compute inverses, if you call inv. Commented Dec 9, 2023 at 20:12
• Sure, but don't do that. Commented Dec 10, 2023 at 2:11
• I was sure that the documentation page for numpy.linalg.inv had a comment recommending to use numpy.linalg.solve instead whenever possible, but apparently I remember wrong.
– Stef
Commented Dec 10, 2023 at 18:27
• They use the exact same algorithm under the hood but obviously less work is done in using solve compared to inv because you just do fewer operations. Commented Dec 11, 2023 at 17:25