Is Latin hypercube sampling effective in multiple dimensions?

Question

I am currently using a Latin Hypercube Sampling (LHS) to generate well-spaced uniform random numbers for Monte Carlo procedures. Although the variance reduction that I obtain from LHS is excellent for 1 dimension, it does not seem to be effective in 2 or more dimensions. Seeing how LHS is a well-known variance reduction technique, I am wondering whether I may be misinterpreting the algorithm or misusing it in some way.

In particular, the LHS algorithm that I use to generate $N$ spaced out uniform random variables in $D$ dimensions is:

• For each dimension $D$, generate a set of $N$ uniformly distributed random numbers $\{u^1_D,u^2_D...u^N_D\}$ such that $u^1_D \in [0,\frac{1}{N+1}]$, $u^2_D \in [\frac{1}{N+1}, \frac{2}{N+1}]$ ... $u^N_D \in [\frac{N}{N+1}, 1]$

• For each dimension $D \geq 2$, randomly reorder the elements from each set. The first $U(0,1)^D$ produced by LHS is the a $D$ dimensional vector containing the first element from each reordered set, the second $U(0,1)^D$ produced by LHS is the a $D$ dimensional vector containing the second element from each reordered set, and so on...

I have included some plots below to illustrate the variance reduction I get in $D = 1 $ and $D = 2$ for a Monte Carlo procedure. In this case, the problem involves estimating the expected value of a cost function $E[c(x)]$ where $c(x) = \phi(x)$, and $x$ is a $D$-dimensional random variable distributed between $[-5,5]$. In particular, the plots show the mean and the standard deviation of 100 sample mean estimates of $E[c(x)]$ for sample sizes of 1000 to 10000.

I get the same type of variance reduction results regardless of whether I use my own implementation or the lhsdesign function in MATLAB. Also, the variance reduction does not not change if I permute all sets of random numbers instead of just the ones corresponding to $D \geq 2$.

The results make sense since stratified sampling in $D = 2$ means that we should sample from $N^2$ squares instead of $N$ squares that are guaranteed to be well spread.

Answer 1

I have split the issues described in your post into three questions below. A good reference for results on Latin Hypercube Sampling and other variance reduction techniques is this book chapter. Also, this book chapter provides information on some of the 'basics' of variance reduction.

Q0. What is variance reduction? Before going into the details, it's helpful to recall what 'variance reduction' actually means. As explained in the 'basics' book chapter, the error variance associated with a Monte Carlo procedure is typically of the form $\sigma^2/n$ under IID sampling. To reduce the error variance, we can either increase the sample size $n$ or find a way to reduce $\sigma$. Variance reduction is concerned with ways of reducing $\sigma$, so such methods may not have any effect on the way in which the error variance changes as $n$ varies.

Q1. Has Latin Hypercube Sampling been correctly implemented? Your written description seems correct to me and is consistent with the description in the book chapter. My only comment is that the ranges of the $u^i_D$ variables don't seem to fill the whole unit interval; it seems that you actually require $u^i_D \in [\frac{i-1}{N}, \frac{i}{N}]$, but hopefully this error did not creep into your implementation. Anyway, the fact that both of the implementations gave similar results would suggest that your implementation is likely to be correct.

Q2. Are your results consistent with what you might expect from LHS? Proposition 10.4 in the book chapter states that the LHS variance can never be (much) worse than the variance obtained from IID sampling. Often, the LHS variance is much less than the IID variance. More precisely, Proposition 10.1 states that, for the LHS estimate $\hat{\mu}_{LHS}=\frac{1}{n} \sum_{i=1}^n f(X_i)$, we have $$\mathrm{Var}(\hat{\mu}_{LHS})=n^{-1}\int e(x)^2dx+o(n^{-1})$$ where $e(x)$ is the 'residual from additivity' of the function $f$ i.e. $f$ minus its best additive approximation (see p.10 of book chapter for details, $f$ is additive if we can write $f(x)=\mu+\sum_{j=1}^D f_j (x_j)$).

For $D=1$, every function is additive so $e=0$ and $\mathrm{Var}(\hat{\mu}_{LHS})=o(n^{-1})$ from Proposition 10.1. In fact, for $D=1$ LHS is equivalent to grid based stratification (Section 10.1 in book chapter) so the variance is actually $O(n^{-3})$ (equation 10.2 in book chapter; assumes $f$ is continuously differentiable). This seems not inconsistent with your first graph. The main point is that $D=1$ is a very special case!

For $D=2$, it is likely the case that $e\neq 0$ so you might expect a variance of order $O(n^{-1})$. Again, this is not inconsistent with your second graph. The actual variance reduction achieved (in comparison to IID sampling) will depend on how close your chosen function is to being additive.

In summary, LHS can be effective in low to moderate dimensions and especially for functions well approximated by additive functions.

Answer 2

http://statweb.stanford.edu/~owen/mc/Ch-var-adv.pdf

This paper discusses the variance reduction of Latin Hypercube Sampling in multiple dimensions. LHS doesn't enforce uniformity when sampling in multiple dimensions because it simply samples in each dimension independently and then combines the dimensions randomly. Stratified sampling of N² bins as you mention is also referred to as Orthogonal Sampling as discussed on the Wikipedia page: https://en.wikipedia.org/wiki/Latin_hypercube_sampling and more enforces multi-dimensional uniformity by sampling from the bins of the all dimensions combined instead.

With a few tweaks to this style of sampling the error variance can be shown to be O(N^-1-2/d) (in ref above). Although this provides large gains for small dimensions, in larger dimensions it begins to degrade back to the performance of ordinary Monte Carlo.

Answer 3

I want to comment on "additivity". LHS makes e.g. sure that X1 and X2 are distributed well (usually in (0,1)), so if a design depends only on one variable you will get a "perfect" histogram and strong variance reduction. For the integration of f=100*X1+X2 you will get good results too, but not for X1-X2! This difference has an almost i.i.d. random distribution, no LHS characteristics. In electronics, designs often exploit that 2 parameters influences will mostly cancel eachother (differential pair, current mirror, replica circuits, etc.), but the effect of mismatch X1-X2 is still present and often dominant. Thus LHS MC analysis behave not better then rnd MC in many electrical designs.