Any suggestions on speeding up the model?

Haha. Yes, I wrote the Cholesky decomposition by hand but finally did not use it. At beginning, I was think the pm.MvNormal on a 62-dimention may cost a lot of time and want to simplify it by cholesky decomposition. Then I found cholesky decomposition did not improved a lot (May be I am wrong, I will try the built in decomposition according to you suggestion). I keep that function just because I wrote it by myself and reluctant to delete, like a watermark for myself.