# Is Latin hypercube sampling effective in multiple dimensions?

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.

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.

• Not sure what it means for $f=X_1-X_2$ to have an "almost iid random distribution, no LHS characteristics". In this case $f$ is still additive so you can expect good variance reduction by using LHS, just as with the additive function $f=100X_1+X_2$. You can verify this by simulation. – S. Catterall Reinstate Monica Aug 25 '17 at 10:24