Scans are slow, so that’s a candidate for where to look. Is the cholesky function just doing a cholesky decomposition? If so, you should be using not pt.linalg.cholesky. That will be a lot faster because the forward computation is farmed out to LAPACK and the backward computation is done using matrix adjoints, eliminating iterative computation.