While you are on this, could you please implement the two different cases if N>D and D>N.
This covariance matrix is D^2 and this is impossible to store for very large D and small N.
The other way stores a matrix that is N^2.

Or do you want to do this later? I would do it all on the fly as you are touching the code/integration tests anyways...

See my issue description at #1876