linear_algebra/real_toeplitz_solve.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exploiting Matrix structure in the solution of Linear Systems"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In a [recent post on his personal blog](http://www.reidatcheson.com/fun/linear%20algebra/2019/02/27/linear-algebra-secrets.html), NAG employee Reid Atcheson mentioned that the fastest numerical linear algebra routines usually exploit some kind of matrix structure. Reid went on to discuss [Toeplitz matrices](http://en.wikipedia.org/wiki/Toeplitz_matrix) in particular. Reid's post got me wondering, just how much of a difference can using specialised solvers make?\n",
    "\n",
    "In one of my talks, [High Performance Computing: There's plenty of Room at the Bottom](http://mikecroucher.github.io/RedOak/), I cite the paper [A survey of the practice of computational science](http://files.software-carpentry.org/training-course/2012/08/prabhu-survey-2011.pdf) 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.\n",
    "\n",
    "To explore this, I made use of the [NAG library for Python](http://www.nag.co.uk/nag-library-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](http://github.com/numericalalgorithmsgroup/NAGPythonExamples).\n",
    "\n",
    "\n",
    "## Solving linear systems of Toeplitz matrices - An example of exploiting stucture\n",
    "\n",
    "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/\n",
    "\n",
    "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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Elapsed time is 5.394560 seconds.\n"
     ]
    }
   ],
   "source": [
    "from naginterfaces.library.lapacklin import dgesv   # A general solver\n",
    "from naginterfaces.library.linsys import real_toeplitz_solve # A toeplitz solver\n",
    "from pytictoc import TicToc\n",
    "import numpy as np\n",
    "import scipy.linalg\n",
    "\n",
    "timer = TicToc()\n",
    "\n",
    "# Construct a real, symmetric, positive definite toeplitz matrix \n",
    "matrix_size = 5000\n",
    "t = np.arange(0,matrix_size);\n",
    "a = np.exp(-np.abs(t)/10);\n",
    "# The toeplitz matrix is defined by it's diagonals.  We can construct the full matrix from the diagonals using scipy\n",
    "A = scipy.linalg.toeplitz(a, a)\n",
    "# Construct and a random Righ hand side\n",
    "np.random.seed(2)\n",
    "b = np.random.rand(matrix_size,1)\n",
    "\n",
    "timer.tic()\n",
    "[asol,ipiv,x_nag_gen] = dgesv(A, b)\n",
    "timer.toc()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now solve the same problem but make use of the function **real_toeplitz_solve**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Elapsed time is 0.035298 seconds.\n"
     ]
    }
   ],
   "source": [
    "# NAG's toeplitz solver requires that b be of dimension N instead of N,1\n",
    "# So we reconstruct it with the correct dimension\n",
    "np.random.seed(2)\n",
    "b = np.random.rand(matrix_size)\n",
    "\n",
    "timer.tic()\n",
    "[x_nag_toe,p] = real_toeplitz_solve(matrix_size,a,b,wantp=False)\n",
    "timer.toc()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "153.99999999999997"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Let's see how many times faster it is to use the toeplitz solver\n",
    "5.39/0.035"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It is **well over 100x faster** to use the Toeplitz solver.  We should check that we get the same answers?  \n",
    "Annoyingly, the solution matrices are are different shapes!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(5000, 1)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Solution from the general solver\n",
    "x_nag_gen.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(5000,)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Solution from the toeplitz solver\n",
    "x_nag_toe.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "but the values agree to high accuracy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.1191048088221578e-13"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.max(abs(x_nag_gen - np.reshape(x_nag_toe,(matrix_size,1))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Other linear solvers in the NAG library for Python that make use of various matrix structures \n",
    "\n",
    "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:\n",
    "\n",
    "\n",
    "### Real matrices\n",
    "\n",
    "* [library.linsys.real_toeplitz_solve](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.linsys.html#naginterfaces.library.linsys.real_toeplitz_solve) Solution of **real symmetric positive definite Toeplitz** system of linear equations.\n",
    "* [library.lapacklin.dgbsv](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.dgbsv) Computes the solution to a **real banded** system of linear equations.\n",
    "* [library.lapacklin.dsgesv](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.dsgesv) Computes the solution to a **real** system of linear equations using **mixed precision arithmetic**.\n",
    "* [library.lapacklin.dsysv](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.dsysv) Computes the solution to a **real symmetric** system of linear equations.\n",
    "* [library.lapacklin.dspsv](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.dspsv) Computes the solution to a **real symmetric** system of linear equations, **packed storage**.\n",
    "* [library.lapacklin.dpbsv](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.dpbsv) Computes the solution to a **real symmetric positive definite banded** system of linear equations.\n",
    "* [library.lapacklin.dposv](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.dposv) Computes the solution to a **real symmetric positive definite** system of linear equations.\n",
    "* [library.lapacklin.dppsv](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.dppsv) Computes the solution to a **real symmetric** positive definite system of linear equations, **packed storage**.\n",
    "* [library.lapacklin.dsposv](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.dsposv) Computes the solution to a **real symmetric positive definite** system of linear equations using **mixed precision arithmetic**.\n",
    "* [library.lapacklin.dptsv](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.dptsv) Computes the solution to a **real symmetric positive definite tridiagonal** system of linear equations.\n",
    "* [library.lapacklin.dtbtrs](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.dtbtrs) Solution of **real band triangular** system of linear equations, **multiple right-hand sides**.\n",
    "* [library.lapacklin.dtrtrs](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.dtrtrs) Solution of **real triangular** system of **linear equations, multiple right-hand sides**.\n",
    "* [library.lapacklin.dtptrs](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.dtptrs) Solution of **real triangular** system of linear equations, **multiple right-hand sides, packed storage**.\n",
    "* [library.lapacklin.dgtsv](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.dgtsv) Computes the solution to a **real tridiagonal** system of linear equations.\n",
    "\n",
    "### Complex matrices\n",
    "\n",
    "* [library.lapacklin.zgbsv](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.zgbsv) Computes the solution to a **complex banded** system of linear equations.\n",
    "* [library.lapacklin.zhesv](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.zhesv) Computes the solution to a **complex Hermitian** system of linear equations.\n",
    "* [library.lapacklin.zhpsv](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.zhpsv) Computes the solution to  **complex Hermitian** of linear equations, **packed storage.**\n",
    "* [library.lapacklin.zpbsv](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.zpbsv) Computes the solution to a **complex Hermitian positive definite banded** system of linear equations.\n",
    "* [library.lapacklin.zposv](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.zposv) Computes the solution to a **complex Hermitian positive definite** system of linear equations.\n",
    "* [library.lapacklin.zppsv](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.zppsv) Computes the solution to a **complex Hermitian positive definite** system of linear equations, **packed storage**.\n",
    "* [library.lapacklin.zcposv](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.zcposv) Computes the solution to a **complex Hermitian positive definite** system of linear equations using **mixed precision** arithmetic.\n",
    "* [library.lapacklin.zptsv](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.zptsv) Computes the solution to a **complex Hermitian positive definite tridiagonal** system of linear equations.\n",
    "* [library.lapacklin.zgesv](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.zgesv) Computes the solution to a **complex** system of linear equations.\n",
    "* [library.lapacklin.zcgesv](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.zcgesv) Computes the solution to a **complex** system of linear equations using **mixed precision** arithmetic.\n",
    "* [library.lapacklin.zsysv](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.zsysv) Computes the solution to a **complex symmetric** system of linear equations. \n",
    "* [library.lapacklin.zspsv](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.zspsv) Computes the solution to a **complex symmetric** system of linear equations, **packed storage**.\n",
    "* [ibrary.lapacklin.ztbtrs](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.ztbtrs) Solution of **complex band triangular** system of linear equations, **multiple right-hand sides**.\n",
    "* [library.lapacklin.ztrtrs](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.ztrtrs) Solution of **complex triangular** system of linear equations, **multiple right-hand sides**.\n",
    "* [library.lapacklin.ztptrs](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.ztptrs) Solution of **complex triangular** system of linear equations, **multiple right-hand sides**, **packed storage**.\n",
    "* [library.lapacklin.zgtsv](https://www.nag.co.uk/numeric/py/nagdoc_latest/naginterfaces.library.lapacklin.html#naginterfaces.library.lapacklin.zgtsv) Computes the solution to a **complex tridiagonal** system of linear equations.\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}