title: How to implement cholesky decomposition using numpy efficiently?
body: <p>I want to implement efficient realization of cholesky decomposition. Naive code looks like</p>
  <pre><code>import numpy as np
  def cholesky(A):
      n = A.shape[0]
      L = np.zeros_like(A)
      for i in range(n):
          for j in range(i+1):
              s = 0
              for k in range(j):
                  s += L[i][k] * L[j][k]
  
              if (i == j):
                  L[i][j] = (A[i][i] - s) ** 0.5
              else:
                  L[i][j] = (1.0 / L[j][j] * (A[i][j] - s))
      return L
  </code></pre>
  <p>I wonder if there is a way to make it more efficient. E.g vectorize it?</p>