import matplotlib
if not hasattr(matplotlib.RcParams, "_get"):
matplotlib.RcParams._get = dict.get
6.3. Solving systems of equations#
In numerical methods, problems are often cast in the form of a linear system of equations
Where \(\mathbf{K}\) is a matrix, \(\mathbf{u}\) is a vector with unknowns and \(\mathbf{f}\) is a known right hand side vector. The system of equations then needs to be solved for \(\mathbf{u}\). Mathematically, it is clear that the solution can be obtained as:
However, computing the inverse of a large matrix is a very expensive operation. Generally, the number of operations that need to be performed to compute the inverse of an \(n\times n\)-matrix is of the order \(n^3\), which implies that it blows up quickly when the matrix is large. Faster algorithms are available that compute the solution \(\mathbf{u}\) to the linear system of equations without computing the inverse of the matrix. Many different algorithms are available that do the job and which is to be preferred depends on the size and structure of the matrix.
For matrices of medium size, numpy.linalg.solve is a better choice than numpy.linalg.inv. The algorithms in numpy.linalg.solve still scale with the order \(n^3\) but they get to a result faster, typically 2-3 times faster. Besides, inv followed by matrix-vector multiplication can give rise to larger round-off errors than solve, which may be significant depending on the condition number of the matrix. Alternatively, one can use scipy instead of numpy by calling scipy.linalg.solve, which gives access to more optimized implementations. In scipy there are dedicated algorithms for when the matrix has a banded structure as one would expect from a numerical simulation on a structured grid.
For sparse matrices, better scaling with the size can be obtained with special algorithms. The function scipy.sparse.linalg.spsolve generally gives much better performance than its dense counterpart.
For very large systems, the direct solver of scipy.sparse.linalg.solve can be outperformed by iterative solvers that search for the \(\mathbf{u}\) that satisfies \(\mathbf{Ku}=\mathbf{f}\) iteratively. Examples are the Conjugate Gradient method scipy.sparse.linalg.cg for symmetric positive definite matrices and the Generalized Minimal Residual method scipy.sparse.linalg.gmres for general non-symmetric matrices. In order to get good performance out of iterative solvers, they need to be combined with a preconditioner, for instance incomplete LU (ILU) decomposition.
Example: tridiagonal matrix#
For illustration, the different methods are compared for solving a linear system of equations with a tridiagonal matrix of variable size. The tridiagonal matrix is of the following form, which is a typical matrix for a 1D elliptic equation.
The tridiagonal matrix is one for which scipy.linalg.solve has a dedicated implementation, so here it is much faster than numpy.linalg.solve. The performance of scipy.linalg.inv is compared to that of scipy.linalg.solve and scipy.sparse.linalg.spsolve, where the size of the \(n\times n\)-matrix is increased up to \(n = 10^6\), so a matrix of size 1 million by 1 million.
Fig. 6.18 Time required to solve linear system of equations with three different methods for an \(n\times n\) tri-diagonal matrix#
Fig. 6.18 shows the runtime for different sizes of the matrix. The near-linear scaling that is obtained with spsolve is only possible because of the very simple structure in the matrix. However, the observation that spsolve is many orders of magnitude faster for large matrices than dense matrix inversion is more generally true. You can estimate from the trend how much time it would take to do matrix inversion on a \(10^6\times10^6\)-matrix. In fact memory requirements would make it impossible to store such a large matrix as dense, but notice the log-scale and imagine how long it would take.
In the code block below, you can try for yourself how the different approaches scale. Here scipy.sparse.linalg.cg with preconditioning is also included. For a more thorough comparison on the timing, they would need to be run multiple times. Note that the relative performance of the different methods may differ between systems and python versions.
import numpy as np
import scipy
import time
import scipy.sparse.linalg as spla
# function to build the tridiagonal matrix of size n x n
def build_matrix(n):
row_indices = np.append(np.append(range(n), range(n-1)), range(1, n))
col_indices = np.append(np.append(range(n), range(1, n)), range(n-1))
data = np.append(np.append(2 * np.ones(n), -1 * np.ones(n-1)), -1 * np.ones(n-1))
matrix = scipy.sparse.csr_matrix((data, (row_indices, col_indices)), shape=(n, n))
return matrix
# define the size of the system of equations
n = 500
# build matrix K and right hand side vector f
Ksparse = build_matrix(n)
Kdense = Ksparse.todense()
f = np.zeros(n)
f[int(n/2)] = 1
# solve using different methods and time them
start_time = time.perf_counter()
u = scipy.linalg.inv(Kdense) @ f
t_inv = time.perf_counter() - start_time
start_time = time.perf_counter()
u = scipy.linalg.solve(Kdense, f)
t_scipy_solve = time.perf_counter() - start_time
start_time = time.perf_counter()
u = spla.spsolve(Ksparse, f)
t_spsolve = time.perf_counter() - start_time
Kcsc = build_matrix(n).tocsc()
start_time = time.perf_counter()
ilu = spla.spilu(Kcsc) # factorization
M2 = spla.LinearOperator(Kcsc.shape, matvec=ilu.solve)
u = spla.cg(Ksparse, f, M=M2)
t_cg = time.perf_counter() - start_time
print(f"Time taken for matrix inversion and multiplication: {t_inv:.4f} seconds")
print(f"Time taken for scipy.linalg.solve: {t_scipy_solve:.4f} seconds")
print(f"Time taken for scipy.sparse.linalg.spsolve: {t_spsolve:.4f} seconds")
print(f"Time taken for scipy.sparse.linalg.cg with ILU preconditioner: {t_cg:.4f} seconds")
Attribution
This chapter is written by Frans van der Meer. Find out more here.