Using scipy for Cholesky decomposition
From what I've seen the scipy cholesky routine is twice as fast when inputs are in np.float32 (which can be the case if !43 (closed) is finalized and merged), whereas the numpy version has the same computational time for float32 or float64. I also have the feeling the scipy cholesky can reduced the memory impact by using the same work array (sigma).