Linear algebra with scipy.linalg

Solve systems with the right factorisation, decompose a matrix on purpose, and know when scipy.linalg beats numpy.linalg.

Solving without inverting

import numpy as np
from scipy import linalg

A = np.array([[4.0, 1.0], [1.0, 3.0]])
b = np.array([1.0, 2.0])

x = linalg.solve(A, b)                    # one square system
X = linalg.solve(A, np.eye(2))            # many right-hand sides at once

L = np.tril(A)
y = linalg.solve_triangular(L, b, lower=True)

M = np.random.default_rng(0).normal(size=(100, 3))
res = linalg.lstsq(M, np.random.default_rng(1).normal(size=100))
print(res[0].shape, res[1].shape)          # coefficients, residual info

⚠️

Never compute inv(A) @ b. It costs about three times as much, loses precision, and hides a singular or ill-conditioned matrix that solve would have flagged with a warning or a LinAlgError.

Factorisations

Factorisation	Use when	Call
LU with pivoting	General square system, reused right-hand sides	`lu_factor` + `lu_solve`
Cholesky	Symmetric positive definite (covariances, kernels)	`cho_factor` + `cho_solve`
QR	Least squares, orthogonalisation, rank detection	`qr` with `pivoting=True`
SVD	Rank, conditioning, dimensionality reduction	`svd`
Eigen	Symmetric or general eigenproblem	`eigh` or `eig`

lu, piv = linalg.lu_factor(A)
x1 = linalg.lu_solve((lu, piv), b)        # cheap for each new b

c, low = linalg.cho_factor(A)             # A must be SPD
x2 = linalg.cho_solve((c, low), b)

U, s, Vt = linalg.svd(M, full_matrices=False)
print("condition number:", s[0] / s[-1])

w, V = linalg.eigh(A)                     # symmetric: real, sorted eigenvalues
print(w)

For a symmetric matrix always reach for eigh: it is faster, returns real sorted eigenvalues, and gives orthonormal eigenvectors. Using general eig there is a common mistake that produces tiny imaginary parts you then have to strip by hand.

scipy.linalg vs numpy.linalg

scipy.linalg assumes a full LAPACK/BLAS link and adds the missing pieces: Cholesky, LU reuse, solve_triangular, banded and Toeplitz solvers.
numpy.linalg is part of the base install and is fine for one-off solve, svd or eig calls.
SciPy adds overwrite_a=True and check_finite=False, which avoid a copy and a NaN scan in hot loops. Only use them when you control the array.
lstsq in SciPy can run several LAPACK drivers (gelsd, gelss, gelsy); the default is accurate, and gelsy is faster when the matrix is full rank.

FAQ

Why is my solution full of NaN?

The matrix is singular or nearly so. Check np.linalg.cond(A); a condition number near 1/eps means the result carries no reliable digits. Regularise, reformulate, or solve in a least-squares sense with lstsq.

Can I solve many systems with the same matrix quickly?

Factor once and reuse: lu_factor followed by lu_solve per right-hand side is far cheaper than a fresh solve each time. If the matrices are unrelated, a plain loop over solve is fine.

SciPy arrays and NumPy interop Sparse matrices in depth

Last refreshed 2026-09-18.

Linear algebra with scipy.linalg

Solving without inverting

Factorisations

scipy.linalg vs numpy.linalg

FAQ

Related