# Generating Latin hypercube samples with numpy

I wrote some code to generate Latin hypercube samples for interpolation over high-dimensional parameter spaces. Latin hypercubes are essentially collections of points on a hypercube that are placed on a cubic/rectangular grid, which possess the property that no two points share any individual coordinate, and every row/column/higher-dimensional-axis is sampled once. e.g., for 2 dimensions and 4 total samples, this is a Latin hypercube:

| | | |x|
| |x| | |
|x| | | |
| | |x| |


One algorithm for generating this is to generate D random permutations of the integers 0 through N-1, where D is the number of dimensions and N is the desired number of samples. Concatenating these together into a DxN matrix gives a list of coordinates which will form a Latin hypercube.

I am also interested in generating a special type of Latin hypercube called a symmetric Latin hypercube. This property essentially means that the hypercube is invariant under spatial inversion. One way of expressing this condition, is that (if the samples are zero-indexed) if the hypercube has the sample (i, j, ..., k), for some integer indices i,j,k, then it also has the sample (n-1-i, n-1-j, ..., n-1-k), where n is the number of samples. For 2D and n=4, this is an example of a symmetric Latin hypercube:

| |x| | |
| | | |x|
|x| | | |
| | |x| |


The condensed version of why I want to do this is that I want to evaluate a function which is expensive to compute, and depends on many parameters. I pre-generate a list of values for this function at a set of points and use interpolation over the parameters to extend this list to any desired point within a set range. Latin hypercube samples are used for the pre-computation instead of a grid, with the idea that they will more efficiently sample the behavior of the function and result in lower interpolation errors.

I have tested the code below both using automated tests and by visual inspection of the results for the trivial 2D case. Statistical tests of higher dimensional cases have convinced me it works successfully there as well. The code can produce either random Latin hypercubes or a rather nongeneral set of symmetric ones (I am convinced my second algorithm cannot in principle generate any possible symmetric Latin hypercube, just a specific type). Also, the symmetric algorithm doesn't really work properly for odd numbers of samples, so I have outright blocked such samples from being generated.

My code can be seen along with specification and performance tests here, although not the visual tests at the moment.

I am aware that my project structure is probably overkill for the amount of code here. I am mostly looking for any insights regarding my algorithms, as while they are very fast even for large dimensions and sample sizes, I don't doubt improvements exist. However any comments are welcome.

import numpy as np

class LatinSampler:
def _check_num_is_even(self, num):
if num % 2 != 0:
raise ValueError("Number of samples must be even")

def get_lh_sample(self, param_mins, param_maxes, num_samples):
dim = param_mins.size

latin_points = np.array([np.random.permutation(num_samples) for i in range(dim)]).T

lengths = (param_maxes - param_mins)[None, :]
return lengths*(latin_points + 0.5)/num_samples + param_mins[None, :]

def get_sym_sample(self, param_mins, param_maxes, num_samples):
self._check_num_is_even(num_samples)

dim = param_mins.size

even_nums = np.arange(0, num_samples, 2)
permutations = np.array([np.random.permutation(even_nums) for i in range(dim)])
inverses = (num_samples - 1) - permutations

latin_points = np.concatenate((permutations,inverses), axis=1).T

lengths = (param_maxes - param_mins)[None, :]
return  lengths*(latin_points + 0.5)/num_samples + param_mins[None, :]


Edit: I glossed over this, but this code doesn't return Latin hypercubes per se, but rather Latin hypercubes stretched/squeezed and shifted so that the points can be over arbitrary ranges of values. However the core code will generate a Latin hypercube if the last manipulation (multiply by lengths and add param_mins) is omitted.

• I don't know if you've checked this already, but you could look at the Eight Queen Problem, which is similar to what you're trying to achieve, and is a pretty popular problem, so there are surely many solutions to the problem. Jul 5, 2019 at 15:29
• @IEatBagels: I cannot speak for the general case and I could be totally wrong, but at least for 2D this looks more like the Eight Rooks Problem since diagonally adjacent samples seem to be allowed. Jul 5, 2019 at 16:28
• @AlexV This is correct, diagonal adjacency is not disallowed. Also here the concerns are a bit different from 8 Rooks. I have no interest in counting how many LHs are possible for a given N and D; I just want to generate them. The biggest concerns are speed of this generation, memory usage, and most importantly, if the "space-filling-ness" of them can be improved. One such criterion for space-filling-ness is called orthogonal sampling and basically says the density of the samples should be optimally uniform, but this is a far harder problem than just generating a LH. Jul 5, 2019 at 17:08
• I'm not convinced get_lh_samples generates latin hypercubes. It looks more like latin hypercuboids. Are you just interested in sampling from a multivariate distribution uniformly, or do you actually want to generate latin hypercubes / symmetrical latin hypercubes? Jul 6, 2019 at 13:20
• @FirefoxMetzger Do you mean in the sense that I am warping the lengths and shifting the minimum values? If so yes, indeed it doesn't actually generate Latin hypercubes. It would generate a cube when all of the parameters span the same range. This is desireable behavior for me; I want a Latin hypercube stretched to fit over the actual parameter ranges I am interested in. I am sorry for glossing over that in the question. If you mean something else I am not sure what, could you clarify? Jul 6, 2019 at 17:27

Once I understood what you wanted to accomplish, your code became a lot more clear to me and I can provide some thoughts.

Regarding the symmetrical version, you are indeed only sampling from a subset of the available symmetrical latin hypercubes. Half of the parameter space will never be chosen as samples. Consider the 2D example with n=4: even_nums = [0, 2] and, hence, A,B,C,D are the 4 possible samples that can be generated with this method. In fact there's two possible Latin Hypercubes: [A,D] and [B,C]. A',B',C',D' denote their respective inverse positions that are then added. As you can see, half of the fields are empty and will never be chosen as a sample.

Instead of permuting indexes, I would sample pairs of symmetrical points directly, as it allows you to control for the constraint while still achieving uniform sampling:

import random
import numpy as np

num_samples = 5
dim = 2

available_indices = [set(range(num_samples)) for _ in range(dim)]
samples = []

# if num_samples is odd, we have to chose the midpoint as a sample
if num_samples % 2 != 0:
k = num_samples//2
samples.append([k] * dim)
for idx in available_indices:
idx.remove(k)

# sample symmetrical pairs
for _ in range(num_samples//2):
sample1 = list()
sample2 = list()

for idx in available_indices:
k = random.sample(idx, 1)
sample1.append(k)
sample2.append(num_samples-1-k)
idx.remove(k)
idx.remove(num_samples-1-k)

samples.append(sample1)
samples.append(sample2)

samples = np.array(samples)


You can use the LatinHypercube class which I added to SciPy in version 1.7. This can replace most of the code here.

To generate a LHS, one simply can generate a uniform sample U(0, 1) and shift it with random permutations of integers, like this:

import numpy as np

n, d = ...
rng = np.random.default_rng()
samples = rng.uniform(size=(n, d))
perms = np.tile(np.arange(1, n + 1), (d, 1))
for i in range(d):
rng.shuffle(perms[i, :])
perms = perms.T

samples = (perms - samples) / n

• Hej @tupui ... the world is a small place. This is a cool new feature and is certainly related; however, I believe the question here is/was primarily looking for symmetric LatinHypercubes. I checked the docs, but I couldn't figure out if this was an option. If there is the option to constrain samples to this subset, perhaps adding it here would address above comments? (as: the provided code can now be refactored into a one-liner) Aug 1, 2021 at 13:20
• @FirefoxMetzger indeed! There is no symmetric option yet. Reading about the problem, it's not really a symmetry but a rotation invariance. If you are interested to add it in SciPy, I would be happy to review a PR ;) TobySpeight sorry I should have put a comment I guess. I will edit to add some code at least. Aug 1, 2021 at 13:34