How to solve linear algebra problems with scipy.linalg in Python

When faced with a system of linear equations, represented in matrix form as Ax = b, it’s tempting to reach for the mathematical definition of the solution: x = A⁻¹b. This approach suggests a two-step process: first, compute the inverse of matrix A, and second, multiply this inverse by the vector b. While this is mathematically sound, it is almost always the wrong way to solve the problem computationally. The SciPy library provides scipy.linalg.inv to compute a matrix inverse, but you should rarely use it to solve a linear system.

Consider this implementation of the “inverse-then-multiply” strategy:

import numpy as np
from scipy import linalg

# A is a 3x3 matrix, b is a vector
A = np.array([[2, 1, -1],
              [-3, -1, 2],
              [-2, 1, 2]])

b = np.array([8, -11, -3])

# Step 1: Compute the inverse of A
A_inv = linalg.inv(A)

# Step 2: Multiply A_inv by b to find x
x = A_inv @ b

# x is now approximately [2., 3., -1.]

This code works, and for a small, well-behaved matrix like this one, it yields the correct result. However, a more direct, numerically stable, and computationally efficient function exists for this exact purpose: scipy.linalg.solve. The preferred approach is to pass the original matrices A and b directly to this function.

import numpy as np
from scipy import linalg

A = np.array([[2, 1, -1],
              [-3, -1, 2],
              [-2, 1, 2]])

b = np.array([8, -11, -3])

# Solve the system Ax = b directly
x = linalg.solve(A, b)

# x is now approximately [2., 3., -1.]

The reasons for preferring solve over inv are twofold and critical. First, numerical stability. Computing a matrix inverse is a delicate operation susceptible to floating-point rounding errors. Small inaccuracies in the computed inverse can be significantly magnified when you subsequently multiply it by the vector b, leading to a less accurate final solution for x. The solve function, by contrast, typically uses an LU decomposition with partial pivoting under the hood. This algorithm is specifically designed to minimize the accumulation of rounding errors by working with the original matrix A and vector b as much as possible.

Second, computational efficiency. Inverting an n x n matrix is a more computationally expensive task than solving a single system of equations. Both operations are generally O(n³), but the constant factor for inversion is larger. Calculating the full inverse matrix provides you with a lot of information, but if your goal is simply to find the vector x, you are performing unnecessary calculations. The LU decomposition used by solve is computationally cheaper than a full inversion. Why compute all the columns of the inverse matrix when all you need is the single result of its product with b?

The performance difference becomes stark as the matrix size increases. For a moderately sized matrix, the difference is not just academic; it is tangible.

import numpy as np
from scipy import linalg
import time

# Set up a larger, well-conditioned system
np.random.seed(42)
n = 1000
A = np.random.rand(n, n) + np.eye(n) # Add identity to improve conditioning
b = np.random.rand(n)

# Time the explicit inversion method
start_time_inv = time.time()
A_inv = linalg.inv(A)
x_inv = A_inv @ b
end_time_inv = time.time()
print(f"Time using linalg.inv: {end_time_inv - start_time_inv:.4f}s")

# Time the direct solve method
start_time_solve = time.time()
x_solve = linalg.solve(A, b)
end_time_solve = time.time()
print(f"Time using linalg.solve: {end_time_solve - start_time_solve:.4f}s")

# Verify the results are equivalent
# print(f"Solutions are close: {np.allclose(x_inv, x_solve)}")

Running this code will consistently show that linalg.solve is significantly faster, often by a factor of two or more. The function linalg.inv has its uses, for instance, when the entries of the inverse matrix themselves are of theoretical interest, such as in statistical contexts for deriving variance estimates. But for the common task of solving Ax = b, it is the wrong tool. Using linalg.solve is not merely a matter of convenience; it represents a fundamentally superior computational approach. It is more robust, more accurate, and faster.

Understand how matrix decompositions power scipy.linalg

The superior performance and accuracy of functions like linalg.solve are not accidental; they are the direct result of using matrix decompositions. A matrix decomposition, or factorization, is a way of rewriting a matrix as a product of two or more other matrices that have simpler, more desirable properties. The core strategy is to transform one difficult problem-solving a general linear system-into a sequence of much easier ones. Many of the most important algorithms in numerical linear algebra are built upon this foundational concept.

For solving a general square linear system Ax = b, the most common decomposition is the LU decomposition. This factorization breaks the matrix A into the product of a lower triangular matrix L and an upper triangular matrix U, such that A = LU. In practice, to ensure numerical stability and avoid division by zero or very small numbers, the algorithm incorporates row interchanges. This is known as partial pivoting, and the resulting factorization is more accurately represented as PA = LU, where P is a permutation matrix that encodes the row swaps.

Once this decomposition is found, solving Ax = b becomes a two-step process. We substitute PA = LU into the original equation, which gives LUx = Pb. We can then define an intermediate vector y = Ux. This breaks the problem into two simpler systems:

1. Solve Ly = Pb for y.
2. Solve Ux = y for x.

The beauty of this approach is that both systems are trivial to solve. Because L is lower triangular, the first equation can be solved with a simple process called forward substitution. Similarly, because U is upper triangular, the second equation is solved using backward substitution. Both substitutions are computationally inexpensive, having O(n²) complexity, whereas the initial decomposition is the most expensive step at O(n³).

You can see this process explicitly using SciPy’s lower-level functions. The function scipy.linalg.lu performs the factorization, and scipy.linalg.solve_triangular handles the substitution steps.

import numpy as np
from scipy import linalg

# A square matrix and a vector b
A = np.array([[2., 5., 8., 7.],
              [5., 2., 2., 8.],
              [7., 5., 6., 6.],
              [5., 4., 4., 8.]])
b = np.array([1., 1., 1., 1.])

# Step 1: Perform the PA = LU factorization
P, L, U = linalg.lu(A)

# P is a permutation matrix, so we apply it to b
Pb = P @ b

# Step 2: Solve Ly = Pb for y using forward substitution
# The 'lower=True' flag tells the solver L is lower triangular
y = linalg.solve_triangular(L, Pb, lower=True)

# Step 3: Solve Ux = y for x using backward substitution
# By default, solve_triangular assumes an upper triangular matrix
x = linalg.solve_triangular(U, y)

# The result x is the solution to Ax = b
# print(x)
# Compare with the direct solver
# print(linalg.solve(A, b))

While LU decomposition is the workhorse for general square systems, it is just one of many factorizations available in scipy.linalg. The choice of decomposition depends on the properties of the matrix and the problem you are trying to solve. For a matrix that is symmetric and positive-definite (a common case in statistics, optimization, and physics), the Cholesky decomposition is a better choice. It factors A into LLᵀ, where L is a lower triangular matrix. This is roughly twice as fast as LU decomposition and is also more numerically stable, as it avoids the need for pivoting.

Another fundamental tool is the QR decomposition, which factors a matrix A into the product of an orthogonal matrix Q and an upper triangular matrix R. Orthogonal matrices are exceptionally well-behaved numerically because their inverse is simply their transpose (Q⁻¹ = Qᵀ). The QR decomposition is the standard method for solving linear least-squares problems, which arise when a system is overdetermined (more equations than unknowns).

Perhaps the most powerful and general factorization is the Singular Value Decomposition (SVD). It can be applied to any m x n matrix A, factoring it into UΣVᵀ, where U and V are orthogonal matrices and Σ is a rectangular diagonal matrix containing the singular values. SVD is a cornerstone of numerical analysis and is used for a vast range of applications, including computing the pseudoinverse for solving ill-conditioned linear systems, determining the rank of a matrix, and principal component analysis (PCA) in data science. Although it is the most computationally expensive of these common decompositions, its robustness and generality make it invaluable.

The high-level functions in scipy.linalg are designed to be smart about these choices. For instance, when you call linalg.solve, it doesn’t blindly apply one single algorithm. It first inspects the matrix A to check for special structures. You can provide hints using function arguments like sym_pos=True to assert that A is symmetric positive-definite, which directs solve to use the more efficient Cholesky-based solver. By understanding the existence and purpose of these underlying decompositions, you can better appreciate why certain functions are preferred over others and how to use the library most effectively. It’s the difference between merely using a tool and understanding how it works.

Distinguish between general and symmetric eigenvalue problems

The eigenvalue problem is another cornerstone of linear algebra, concerned with finding the special scalars (eigenvalues) and vectors (eigenvectors) for a square matrix A that satisfy the equation Av = λv. In this relationship, multiplying the eigenvector v by the matrix A has the same effect as scaling v by the scalar eigenvalue λ. These values are fundamental to understanding matrix transformations, analyzing dynamical systems, and many data analysis techniques. For a general, non-symmetric square matrix, SciPy provides the function scipy.linalg.eig to solve this problem.

This function takes a single matrix A as input and returns two arrays: one containing the eigenvalues and another whose columns are the corresponding right eigenvectors.

import numpy as np
from scipy import linalg

# A general non-symmetric matrix
A = np.array([[0, 1, -1],
              [1, 1, 0],
              [1, 0, 1]])

# Solve the general eigenvalue problem
eigenvalues, eigenvectors = linalg.eig(A)

# eigenvalues is an array of (potentially complex) scalars
# eigenvectors is a matrix where column i corresponds to eigenvalues[i]
# print("Eigenvalues:")
# print(eigenvalues)
# print("nEigenvectors:")
# print(eigenvectors)

For a non-symmetric matrix, the eigenvalues can be complex numbers, and the eigenvectors are not guaranteed to be orthogonal. The function linalg.eig is the general-purpose tool designed to handle this complexity. However, a very important special case exists when the matrix A is symmetric (or Hermitian, if it contains complex numbers, meaning A = Aᴴ, where H denotes the conjugate transpose). Symmetric matrices have two remarkable properties regarding their eigensystem: their eigenvalues are always real numbers, and their eigenvectors can be chosen to form an orthonormal basis.

These properties are not just mathematically elegant; they enable the use of far more efficient and numerically stable algorithms for finding the eigenvalues and eigenvectors. Recognizing this, SciPy provides a separate, specialized function: scipy.linalg.eigh (the ‘h’ stands for Hermitian). You should always prefer eigh over eig when you know your matrix is symmetric. Using the general-purpose eig on a symmetric matrix is not an error, but it is inefficient. It’s like using a general-purpose sorting algorithm on an already nearly-sorted list; you’re ignoring a special structure that could be exploited for a significant performance gain.

The eigh function leverages algorithms that take advantage of the matrix’s symmetry, making it substantially faster and less prone to numerical error for this class of problems. It also returns the eigenvalues as a purely real array, which is mathematically guaranteed.

import numpy as np
from scipy import linalg

# A symmetric matrix
B = np.array([[6, 3, 1, 5],
              [3, 0, 5, 1],
              [1, 5, 6, 2],
              [5, 1, 2, 4]])

# Use the specialized solver for symmetric matrices
eigenvalues_h, eigenvectors_h = linalg.eigh(B)

# eigenvalues_h is a real-valued array
# The columns of eigenvectors_h form an orthonormal set
# print("Eigenvalues (symmetric):")
# print(eigenvalues_h)
# print("nEigenvectors (symmetric):")
# print(eigenvectors_h)

The performance difference between eig and eigh is not trivial. For a symmetric matrix, eigh can be several times faster because it does less work. The general algorithm used by eig first reduces the matrix to Hessenberg form, then iteratively applies the QR algorithm. The specialized algorithm for eigh reduces the symmetric matrix to a simpler tridiagonal form, which makes the subsequent iterative steps much faster.

import numpy as np
from scipy import linalg
import time

# Create a large symmetric matrix
n = 500
np.random.seed(0)
A_rand = np.random.rand(n, n)
S = (A_rand + A_rand.T) / 2 # Ensure S is symmetric

# Time the general eigenvalue solver
start_time_eig = time.time()
vals_eig, vecs_eig = linalg.eig(S)
end_time_eig = time.time()
print(f"Time using linalg.eig on a symmetric matrix: {end_time_eig - start_time_eig:.4f}s")

# Time the specialized symmetric eigenvalue solver
start_time_eigh = time.time()
vals_eigh, vecs_eigh = linalg.eigh(S)
end_time_eigh = time.time()
print(f"Time using linalg.eigh on a symmetric matrix: {end_time_eigh - start_time_eigh:.4f}s")

Beyond speed, eigh provides a more accurate result in the context of the problem’s mathematical guarantees. If you use eig on a symmetric matrix, floating-point inaccuracies may result in eigenvalues with tiny, spurious imaginary components. The eigh function, by its design, ensures the output is purely real, aligning the computation with the underlying theory. Furthermore, if you only need the eigenvalues and not the eigenvectors, eigh offers an additional optimization. By passing the argument eigvals_only=True, you can instruct the function to compute only the eigenvalues, which is even faster.

This distinction between eig and eigh is a perfect illustration of a core principle in numerical computing: use the most specific tool available for your problem. The library’s designers provide these specialized functions because the underlying mathematical structures permit vastly superior algorithms. Ignoring them means leaving performance and accuracy on the table. Choosing correctly between the general and the specialized solver is a hallmark of an effective user of scipy.linalg.