Any suggestions on continous error "ValueError: shape mismatch: value array of shape (62,) could not be broadcast to indexing result of shape (0,)"

I have seen in some instances that supplying chol instead of cov speeds up mvnormal sampling. And that seems to be the recommended usage:

https://www.pymc.io/projects/docs/en/stable/api/distributions/generated/pymc.MvNormal.html

Also have you tried NUTs instead of Metropolis? Finally you can also try numpyro: