Most efficient way to obtain samples from high-dimensional multivariate distributions?

Question

Is MultinormalDistribution[] efficient and easy to use for high dimensions?

I have a variable $n$ representing the dimension of a Monte Carlo integration I do on a multivariate Gaussian copula, where typical values of $n$ are near 100. I am using a simple correlation matrix made from ConstantArray[ρ, {n,m}] (except on the diagonal).

For now I simply wrote a function that calls RandomVariate[NormalDistribution[0,1], {n,m}], generates the correlation matrix, calculates its Cholesky decomposition, and then multiplies to obtain the correlated multivariate normal samples. It works fine and was easy to code and understand.

However, generating correlated multivariate $t$ samples is more involved, so it would be convenient to just drop in a call to MultivariateTDistribution[]. That leaves me with two issues:

1. Setting up MultinormalDistribution[] with an arbitrary variable count is hard because it seems to want variable names, which I am having trouble generating programmatically.
2. I am not sure that the internals of RandomVariate on MultinormalDistribution[] and MultivariateTDistribution[] are set up to efficiently obtain high-dimensional sample sets.

If the efficiency is not expected to be that great I will stick with my current approach. Otherwise I would appreciate advice on using MultinormalDistribution in this high-dimensional context, since it would be worth investing the time in this more elegant and Mathematica-like approach.

Answer 1

Multivariate distributions do not require any named variables. My hunch is that your confusion in this regard is due to the excessive use of {{1, ρ}, {ρ, 1}} in the documentation for MultivariateTDistribution. You might've missed the fact that in most cases, the value for ρ is being substituted via a Table or a ParallelTable.

You can directly input your covariance matrix (and any additional parameters depending on the distribution) and call RandomVariate to draw from that. For your example:

Σ = With[{ρ = 0.2}, SparseArray[{i_, i_} -> 1, {100, 100}, ρ]];
x = RandomVariate[MultivariateTDistribution[Σ, 10], {100}];

This will generate 100 samples of multivariate T distributed random vectors, each of length 100.

Answer 2

Yes.

In principle, one call to RandomVariate will incur some overhead to diagonalize the matrix and then cost a fixed amount of time per output vector. Let's see by first creating some data to describe such a distribution:

μ = RandomReal[{0, 1}, 100];
Σ = DiagonalMatrix[Exp[RandomReal[{0, 1}, 100]]];

Now request a bunch of random draws (each one of them a 100-vector):

Timing[RandomVariate[MultinormalDistribution[μ, Σ], 100000];]

The response is 0.687 seconds: that is, we can get about 150,000 100-vectors per second. Not shabby and unlikely to be bettered by other approaches. (MultinormalDistribution appears to be parallelized automatically, so the timing is the sum of (a) pre-processing overhead, including RAM allocation for the output, all performed by one core, and (b) generation of random values, which appears to be performed in parallel by all the available (licensed) cores. My timings reflect a four-core license.)

Incidentally, because there must be some matrix multiplication going on for each output vector, we ought to expect slightly sublinear scaling with dimension. E.g., generating 10,000 1000-vectors produces the same quantity of output numbers (i.e., $10^7$ of them in both cases) but requires 1.30 seconds--about twice as long.