TopolearnLeast squares, simulations and statistics rely on solving linear systems quickly and stably.
scipy.linalg contains everything in numpy.linalg plus more, and is always built against optimised LAPACK. Prefer it when you need more than the basics.
Solving `A x = b`: use linalg.solve(A, b), never inv(A) @ b; it is faster and more accurate. Tell it about structure with assume_a="pos" (symmetric positive definite) or "sym". For banded or triangular matrices use solve_banded and solve_triangular.
Factorisations you can reuse:
lu_factor/lu_solve: factor once, solve for many right-hand sides,cho_factor/cho_solve: Cholesky for positive definite matrices such as covariances,qr,svd,schur.
Eigenproblems: linalg.eigh for symmetric/Hermitian matrices (real eigenvalues, sorted ascending, more stable) and linalg.eig for general ones.
Least squares: linalg.lstsq(A, b) for over-determined systems.
Matrix functions: expm (matrix exponential, used for linear ODE systems and Markov chains), logm, sqrtm.
Check conditioning with np.linalg.cond(A): a huge condition number means the solution is very sensitive to tiny changes in the input, and no solver can fix that.
Example
import numpy as np
from scipy import linalg
A = np.array([[4.0, 2.0], [2.0, 3.0]]) # symmetric positive definite
b = np.array([2.0, 1.0])
print(linalg.solve(A, b, assume_a="pos")) # [0.5 0. ]
c = linalg.cho_factor(A) # factor once...
for rhs in ([1.0, 0.0], [0.0, 1.0]): # ...solve many times
print(linalg.cho_solve(c, rhs).round(3))
# [ 0.375 -0.25 ]
# [-0.25 0.5 ]
vals, vecs = linalg.eigh(A)
print(vals.round(3)) # [1.438 5.562]
Q = np.array([[-0.5, 0.5], [0.2, -0.2]]) # Markov generator matrix
print(linalg.expm(Q * 10).round(3)) # transition probabilities after t=10
# [[0.286 0.714]
# [0.286 0.714]]Remember
- Use
solve, notinv(A) @ b - Factor once with
lu_factor/cho_factorto solve many right-hand sides eighfor symmetric matrices;expmfor matrix exponentials
Watch out for
- Computing an explicit inverse to solve a linear system.