See my approach here (I think that’s the version in that repo that’s most similar to your use case). Note that the likelihood could be accelerated even further using the Gaussian sufficient stats trick, but I haven’t gotten around to adding that yet.
But remember: best to get the model working correctly with smaller simulated data before delving into optimization for larger data.