Placeholder

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

In [1]:
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

In [2]:
# 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.
In [3]:
# Let's see how many times faster it is to use the toeplitz solver
5.39/ 0.035
Out[3]:
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!

In [5]:
#Solution from the general solver
x_nag_gen.shape
Out[5]:
(5000, 1)
In [6]:
#Solution from the toeplitz solver
x_nag_toe.shape
Out[6]:
(5000,)

but the values agree to high accuracy

In [8]:
np.max(abs(x_nag_gen - np.reshape(x_nag_toe,(matrix_size,1))))
Out[8]:
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

Complex matrices

NAG Library for Python more information
Author
Mike Croucher
Leave a Comment