if you want to save the inverse mass matrix from a run of the sampler given that each chain converges to a different inverse mass matrix, which mass matrix should be saved for reuse?

You want to save the inverse mass matrix from each chain. Then if you want to restart four chains, you’ll have the data to do that with.

The basic API is only one chain at a time, though. Maybe we should think about generalizing that? Especially if we could multi-thread everything for users.

Really? My intuition would opt for an appropriate average of the different inverse mass matrices. That thing is the global variance-covariance matrix, so I don’t quite see why one would want to not average it which should get the noise down in its estimate.

Again, this is my intuition which can well be wrong. Saving each matrix can also make sense, so maybe one should have the option.

It depends whether you want to be able to stop and start and get the same answers as if you’d just kept going. If you’re willing to change everything and restart with a new seed, sure.

I also figured at the interface level, the easiest thing to do would be to make them all available. At the very least, we need them for diagnostics.

How do you average covariance matrices to get the same result as if you’d thrown all the data together and estimated a covariance (assuming equal length chains)? I’m pretty sure you’d need the means (and counts?) too. We could do it at the level of sufficient statistics like the Welford accumulators.