In a recent post on his personal blog, my colleague, Reid Atcheson, mentioned that the fastest numerical linear algebra routines usually exploit some kind of matrix structure. Reid went on to discuss Toeplitz matrices in particular. Reid's post got me wondering, just how much of a difference can using specialised solvers make?
In one of my talks, High Performance Computing: There's plenty of Room at the Bottom, I cite the paper A survey of the practice of computational science by Prabhu et al who say "Across disciplines, an order of magnitude performance improvement was cited as a requirement for significant changes in research quality." So, ideally anything I do to speed things up will give us at least a factor of 10 improvement.
To explore this, I made use of the NAG Library for Python within a Jupyter notebook which is what you are looking at now. The original notebook is available in the linear algebra section of NAG's Python Examples GitHub repo.
Solving linear systems of Toeplitz matrices - An example of exploiting stucture
Toeplitz matrices occur in various applications of which a couple of examples are given at http://jack.valmadre.net/notes/2015/03/28/symmetric-positive-toeplitz/
Here we solve a linear system where the coefficient matrix is a Toeplitz matrix. First, we make no use of the underlying structure and just use NAG's interface to the general purpose dgesv function
from naginterfaces.library.lapacklin import dgesv # A general solver from naginterfaces.library.linsys import real_toeplitz_solve # A toeplitz solver from pytictoc import TicToc import numpy as np import scipy.linalg timer = TicToc () # Construct a real, symmetric, positive definite toeplitz matrix matrix_size = 5000 t = np.arange(0,matrix_size); a = np.exp(-np.abs(t)/10); # The toeplitz matrix is defined by it's diagonals. # We can construct the full matrix from the diagonals using scipy A = scipy.linalg.toeplitz(a, a) # Construct and a random Right hand side np.random.seed(2) b = np.random.rand(matrix_size,1) timer.tic() [asol,ipiv,x_nag_gen] = dgesv(A, b) timer.toc()
Elapsed time is 5.394560 seconds.
We will now solve the same problem but make use of the function real_toeplitz_solve
# NAG's toeplitz solver requires that b be of dimension N instead of N,1
# So we reconstruct it with the correct dimension
np.random.seed(2)
b = np.random.rand(matrix_size)
timer.tic()
[x_nag_toe,p] = real_toeplitz_solve(matrix_size,a,b,wantp=False)
timer.toc()
Elapsed time is 0.035298 seconds.
# Let's see how many times faster it is to use the toeplitz solver 5.39/ 0.035
153.99999999999997
It is well over 100x faster to use the Toeplitz solver. We should check that we get the same answers?
Annoyingly, the solution matrices are different shapes!
#Solution from the general solver x_nag_gen.shape
(5000, 1)
#Solution from the toeplitz solver x_nag_toe.shape
(5000,)
but the values agree to high accuracy
np.max(abs(x_nag_gen - np.reshape(x_nag_toe,(matrix_size,1))))
1.1191048088221578e-13
Other linear solvers in the NAG Library for Python that make use of various matrix structures
Your application may not be able to make use of Toeplitz solvers but it may well be able to take advantage of one of the other specialised linear solvers in the NAG Library. In addition to the mathematical structure of the matrix, other specialisations include the use of packed storage and mixed precision arithmetic. If you have multiple right-hand sides to compute with the same coefficient matrix A, there are some specialised solvers for that situation too. Here are some of the available solvers with links to their Python documentation:
Real matrices
- library.linsys.real_toeplitz_solve Solution of real symmetric positive definite Toeplitz system of linear equations.
- library.lapacklin.dgbsv Computes the solution to a real banded system of linear equations.
- library.lapacklin.dsgesv Computes the solution to a real system of linear equations using mixed precision arithmetic.
- library.lapacklin.dsysv Computes the solution to a real symmetric system of linear equations.
- library.lapacklin.dspsv Computes the solution to a real symmetric system of linear equations, packed storage.
- library.lapacklin.dpbsv Computes the solution to a real symmetric positive definite banded system of linear equations.
- library.lapacklin.dposv Computes the solution to a real symmetric positive definite system of linear equations.
- library.lapacklin.dppsv Computes the solution to a real symmetric positive definite system of linear equations, packed storage.
- library.lapacklin.dsposv Computes the solution to a real symmetric positive definite system of linear equations using mixed precision arithmetic.
- library.lapacklin.dptsv Computes the solution to a real symmetric positive definite tridiagonal system of linear equations.
- library.lapacklin.dtbtrs Solution of real band triangular system of linear equations, multiple right-hand sides.
- library.lapacklin.dtrtrs Solution of real triangular system of linear equations, multiple right-hand sides.
- library.lapacklin.dtptrs Solution of real triangular system of linear equations, multiple right-hand sides, packed storage.
- library.lapacklin.dgtsv Computes the solution to a real tridiagonal system of linear equations.
Complex matrices
- library.lapacklin.zgbsv Computes the solution to a complex banded system of linear equations.
- library.lapacklin.zhesv Computes the solution to a complex Hermitian system of linear equations.
- library.lapacklin.zhpsv Computes the solution to complex Hermitian of linear equations, packed storage.
- library.lapacklin.zpbsv Computes the solution to a complex Hermitian positive definite banded system of linear equations.
- library.lapacklin.zposv Computes the solution to a complex Hermitian positive definite system of linear equations.
- library.lapacklin.zppsv Computes the solution to a complex Hermitian positive definite system of linear equations, packed storage.
- library.lapacklin.zcposv Computes the solution to a complex Hermitian positive definite system of linear equations using mixed precision arithmetic.
- library.lapacklin.zptsv Computes the solution to a complex Hermitian positive definite tridiagonal system of linear equations.
- library.lapacklin.zgesv Computes the solution to a complex system of linear equations.
- library.lapacklin.zcgesv Computes the solution to a complex system of linear equations using mixed precision arithmetic.
- library.lapacklin.zsysv Computes the solution to a complex symmetric system of linear equations.
- library.lapacklin.zspsv Computes the solution to a complex symmetric system of linear equations, packed storage.
- ibrary.lapacklin.ztbtrs Solution of complex band triangular system of linear equations, multiple right-hand sides.
- library.lapacklin.ztrtrs Solution of complex triangular system of linear equations, multiple right-hand sides.
- library.lapacklin.ztptrs Solution of complex triangular system of linear equations, multiple right-hand sides, packed storage.
- library.lapacklin.zgtsv Computes the solution to a complex tridiagonal system of linear equations.