Test case for slowdown bug in SciPy 0.13 expm.

test-mat.npy

time_expm.ipynb

{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This problem manifests itself in SciPy 0.13 and later."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import scipy\n",
      "print scipy.__version__"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "0.13.2\n"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "from scipy.linalg import expm, expm2, expm3"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "For this case, we chose a matrix typical to our problem, a complex and approximately skew-Hermetian matrix of shape $(36, 36)$."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "G = np.load('test-mat.npy')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print G.shape"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "(36, 36)\n"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print np.linalg.norm(G + G.conj().T) / np.linalg.norm(G)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "0.00588108750055\n"
       ]
      }
     ],
     "prompt_number": 5
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "In SciPy versions less than 0.13, ``expm(G)`` is fastest, as we expect, but in later versions, ``expm(G)`` is significantly slower."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit expm(G)\n",
      "%timeit expm2(G)\n",
      "%timeit expm3(G)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100 loops, best of 3: 8.12 ms per loop\n",
        "1000 loops, best of 3: 1.4 ms per loop"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "1000 loops, best of 3: 1.73 ms per loop"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Profiling, about one third of the cumulative time is taken in ``onenormest`` (on slower computers, this is more severe still):"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%prun -r\n",
      "for idx in xrange(1000):\n",
      "    expm(G)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 7,
       "text": [
        "<pstats.Stats instance at 0x10c7a0b48>"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "stats = _"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import sys\n",
      "stats.stream = sys.stdout"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "stats.print_stats(10)\n",
      "stats.print_callees('expm')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "         2416054 function calls in 9.087 seconds\n",
        "\n",
        "   Ordered by: internal time\n",
        "   List reduced from 87 to 10 due to restriction <10>\n",
        "\n",
        "   ncalls  tottime  percall  cumtime  percall filename:lineno(function)\n",
        "   363000    2.278    0.000    2.278    0.000 {method 'reduce' of 'numpy.ufunc' objects}\n",
        "   131000    1.242    0.000    1.242    0.000 {method 'dot' of 'numpy.ndarray' objects}\n",
        "   324000    1.107    0.000    3.674    0.000 _onenormest.py:161(norm_1d_inf)\n",
        "     5000    0.925    0.000    7.568    0.002 _onenormest.py:272(_onenormest_core)\n",
        "   334000    0.392    0.000    2.754    0.000 fromnumeric.py:2048(amax)\n",
        "   337000    0.336    0.000    2.390    0.000 _methods.py:15(_amax)\n",
        "    70228    0.216    0.000    0.415    0.000 _onenormest.py:121(vectors_are_parallel)\n",
        "    70228    0.199    0.000    0.199    0.000 {numpy.core._dotblas.dot}\n",
        "     1000    0.172    0.000    0.297    0.000 matfuncs.py:581(_pade13)\n",
        "    28228    0.148    0.000    0.607    0.000 _onenormest.py:138(column_needs_resampling)\n",
        "\n",
        "\n",
        "   Ordered by: internal time\n",
        "   List reduced from 87 to 2 due to restriction <'expm'>\n",
        "\n",
        "Function               called...\n",
        "                           ncalls  tottime  cumtime\n",
        "matfuncs.py:267(expm)  ->    2000    0.005    0.115  matfuncs.py:64(_exact_1_norm)\n",
        "                             1000    0.003    0.013  matfuncs.py:73(_ident_like)\n",
        "                             1000    0.005    0.088  matfuncs.py:93(_is_upper_triangular)\n",
        "                             3000    0.012    4.773  matfuncs.py:183(_onenormest_matrix_power)\n",
        "                             1000    0.004    1.569  matfuncs.py:226(_onenormest_product)\n",
        "                             1000    0.020    0.229  matfuncs.py:352(_solve_P_Q)\n",
        "                             1000    0.032    1.457  matfuncs.py:488(_ell)\n",
        "                             1000    0.172    0.297  matfuncs.py:581(_pade13)\n",
        "                             5000    0.003    0.003  {max}\n",
        "                            11000    0.435    0.435  {method 'dot' of 'numpy.ndarray' objects}\n",
        "                             1000    0.000    0.000  {min}\n",
        "                             1000    0.002    0.002  {range}\n",
        "matfuncs.py:81(expm)   ->    1000    0.097    9.077  matfuncs.py:267(expm)\n",
        "\n",
        "\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 11,
       "text": [
        "<pstats.Stats instance at 0x10c7a0b48>"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Finally, we check that the accuracy meets with our expectations."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "A = np.empty((3,3))\n",
      "for i1, f1 in enumerate([expm, expm2, expm3]):\n",
      "    for i2, f2 in enumerate([expm, expm2, expm3]):\n",
      "        A[i1,i2] = np.linalg.norm(f1(G)-f2(G))\n",
      "print A"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[  0.00000000e+00   1.12766022e-12   1.84546262e+38]\n",
        " [  1.12766022e-12   0.00000000e+00   1.84546262e+38]\n",
        " [  1.84546262e+38   1.84546262e+38   0.00000000e+00]]\n"
       ]
      }
     ],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}