title: Performance gap between np.linalg.solve and scipy.linalg.cho_solve
body: <p>I am trying to find <code>alpha</code> by solving two systems of linear equations that involve Cholesky decomposition. <code>scipy</code> has a special function to do so. There is a significant performance gap between <code>scipy</code> and <code>numpy</code>. Can I achieve performance as good as <code>scipy</code> in <code>numpy</code> with any other way? (assume that I am not allowed to use scipy).</p>
  <pre class="lang-py prettyprint-override"><code>import numpy as np
  import scipy
  
  def numpy_cho_solve(N,M):
      for seed in range(N):
          np.random.seed(seed)
          x = np.random.rand(M,1)
          y = np.random.rand(M,1)
          k = x@x.T + np.eye(M)# M*M
          L = np.linalg.cholesky(k)
          alpha = np.linalg.solve(L.T, np.linalg.solve(L, y))
          
  def scipy_cho_solve(N,M):
      for seed in range(N):
          np.random.seed(seed)
          x = np.random.rand(M,1)
          y = np.random.rand(M,1)
          k = x@x.T + np.eye(M)# M*M
          L = np.linalg.cholesky(k)
          alpha = scipy.linalg.cho_solve((L,True), y)
  
  %timeit numpy_cho_solve(100,100)
  %timeit scipy_cho_solve(100,100)
  </code></pre>
  <p>Output</p>
  <pre><code>317 ms ± 12.3 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
  76.9 ms ± 3.12 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
  </code></pre>