LHS sampling can return perfectly correlated set of points

[dk-teknologisk-mon]

The present LHS algorithm creates the design by creating a list of N equally spaced values along each dimension, followed by a randomization step. This does not safeguard against obtaining values in a sorted order, which is not ideal. The problem is biggest for spaces with few dimensions and a low number of points. Here the risk of obtaining pearls on a string along the diagonal is greatest.

There are many ways of potentially improving the algorithm, one could be to create a relatively large number of candidate designs and pick the one that is optimized in terms of maximin, e.g. the maximal minimum pairwise distance between points in the space. An alternative could be optimizing for the lowest covariance / orthogonality of the design.

Furthermore, right now the internal LHS algorithm is hard-seeded in the construction of the Optimizer. I'm not sure this is ideal and maybe we should instead pass an RNG to the LHS call?

Here is one example of the problem:

from ProcessOptimizer import Optimizer
import matplotlib.pyplot as plt

space = [(0.0, 1.0), (0.0, 1.0)]
n_points = 5  # Number of points in the design
seed = 41

# Build optimizer
opt = Optimizer(
    dimensions=space,
    lhs=True,
    n_initial_points=n_points,
    random_state=seed,
)

x = opt.space.lhs(n=n_points, seed=seed)
# Plot the starting version
plt.figure()
for x in x:
    plt.scatter(x[0], x[1])
plt.xlim(0, 1)
plt.ylim(0, 1)
plt.title("Badly correlated initial points")

[dk-teknologisk-mon]

And after a quick round of ChatGPT, here is a maximin optimizer that works. It would need a bit of formatting to conform to the PEP style of the repo, but at least the functionality is there.

from scipy.spatial.distance import cdist

def maximin_lhs(n_points, optimizer, iterations=10000):
    """
    Generate a Latin Hypercube Sampling (LHS) design optimized using the maximin criterion.
    
    Parameters:
    - n_points: Number of sample points.
    - optimizer: The optimizer for the system we want LHS points from.
    - iterations: Number of optimization iterations.
    
    Returns:
    - Optimized LHS design as a numpy array of shape (n_points, n_dimensions).
    """
    def calculate_min_distance(points):
        """Calculate the minimum pairwise distance in the design."""
        distances = cdist(points, points)
        np.fill_diagonal(distances, np.inf)  # Ignore self-distances
        return np.min(distances)

    # Generate an initial LHS design
    best_lhs = opt.space.lhs(n=n_points, seed=42)
    best_min_dist = calculate_min_distance(best_lhs)

    # Optimize the LHS design using maximin criterion
    for _ in range(iterations):
        candidate_lhs = opt.space.lhs(n=n_points, seed=None)
        candidate_min_dist = calculate_min_distance(candidate_lhs)
        
        # Update the best design if the candidate is better
        if candidate_min_dist > best_min_dist:
            best_lhs = candidate_lhs
            best_min_dist = candidate_min_dist

    return best_lhs

Comparing before and after with 15 points, starting with the present algorithm:

Random:

Maximin optimized:

[dk-teknologisk-mon]

The function takes maybe a second to run if you increase the number of factors to e.g. 10, but not really in a problematic way since you'll likely only call this function once in an optimization.

[RuneChristensen-NN]

The GPT suggestion is quite similar to _lhsmaxmin here: https://github.com/clicumu/pyDOE2/blob/master/pyDOE2/doe_lhs.py

It seems pdist would be faster than cdist (although I have not tested it) and you would not need the np.fill_diagonal call afterwards.

Also, running 10000 iterations is perhaps quite a lot more than you should need. It is the default the Sigurd made for Scikit-optimize, but other implementations of the same algorithm has default iterations at 100.