title: Why do scipy.linalg's methods run slower on 9x9 matrices?
body: <p>I'm comparing runtimes for various ways to solve linear systems, and I've found an odd pattern. The solution methods I'm testing are <code>la.solve()</code>, <code>la.inv()</code>, and <code>la.lu_factor_solve()</code>.</p>
  
  <pre><code>import scipy.linalg as la
  import numpy as np
  from time import time
  from matplotlib import pyplot as plt
  
  N = 20  # up to NxN matrix
  T = 100  # run T times
  inv_time, solve_time = [[] for _ in range(20)], [[] for _ in range(20)]
  lu_factor_solve, lu_just_solve = [[] for _ in range(20)], [[] for _ in range(20)]
  for _ in range(T):
      for n in range(1, N + 1):
          A = np.random.rand(n, n)
          b = np.random.rand(n, 1)
          np.dot(la.inv(A), b)  # the first time through is always slow,
          la.solve(A, b)        # so we run it once to get it out of the way
          start = time()
          np.dot(la.inv(A), b)
          end = time()
          inv_time[n - 1].append(end - start)
  
          start = time()
          la.solve(A, b)
          end = time()
          solve_time[n - 1].append(end - start)
  
          start = time()
          la.lu_solve(la.lu_factor(A), b)
          end = time()
          lu_factor_solve[n - 1].append(end - start)
  
          temp = la.lu_factor(A)
          start = time()
          la.lu_solve(temp, b)
          end = time()
          lu_just_solve[n - 1].append(end - start)
  inv_time = np.mean(np.array(inv_time), axis=1)
  solve_time = np.mean(np.array(solve_time), axis=1)
  lu_factor_solve = np.mean(np.array(lu_factor_solve), axis=1)
  lu_just_solve = np.mean(np.array(lu_just_solve), axis=1)
  # do some plots
  plt.plot(range(1, N + 1), inv_time, '-o', label='by inverse')
  plt.plot(range(1, N + 1), solve_time, '-o', label='by la.solve()')
  plt.plot(range(1, N + 1), lu_factor_solve, '-o', label='by lu factor solve')
  plt.yscale('log')
  plt.plot(range(1, N + 1), lu_just_solve, '-o', label='just la.lu_solve()')
  plt.legend()
  plt.show()
  </code></pre>
  
  <p>These are run on random matrices <code>A</code> and column vectors <code>b</code> produced by <code>A = np.random.rand(n, n)</code> and <code>b = np.random.rand(n, 1)</code> for values of <code>n</code> from 1 to 20. I've found that every time I run the program, it finds solutions to 9x9 matrices much slower than for matrices of other sizes, as pictured below. The red line shows the time required just to perform <code>la.lu_solve()</code>. Here is a graph of the results with <code>T = 1</code>:
  <a href="https://i.stack.imgur.com/AY67R.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/AY67R.png" alt="Graph of times with T=1"></a></p>
  
  <p>And here are the results with <code>T = 100</code>:
  <a href="https://i.stack.imgur.com/2FeRT.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/2FeRT.png" alt="Graph of times with T=100"></a></p>
  
  <p>Does this have to do with something inherent about 9x9 matrices, some optimization not available for this size of matrix, or something different?</p>