[SYSTEMDS-3925] Add lhs builtin function

[jancies2301]

Overview

This PR introduces a new builtin function lhs() to Apache SystemDS that implements Latin Hypercube Sampling (LHS), a statistical sampling method for generating well-distributed samples from a multidimensional uniform distribution.

Implementation

New builtin function: lhs() - performs Latin hypercube sampling with multiple algorithm choices.

Parameters

• N: Number of samples
• d: Number of dimensions
• method: Sampling method (random, build, cp_sweep, genetic)
• reps: Repetitions for optimization
• goal: Optimization criterion (maximin, avg_dist, opt, sum_inv)
• eps: Convergence threshold for cp_sweep
• genetic_generations, genetic_population_size, genetic_mutation_rate: Genetic algorithm parameters
• return_type: Output format (buckets or numerical)
• seed: which seed to use (default is -1)

Returns:

• M: the sampled latin hypercube as a Nxd matrix
• G: final value of the optimized goal (negative for maximin and avg_dist goals)

Following sampling methods are supported:

• Random - the hypercube is generated at random. In case reps>1, multiple hypercubes are generated and the best one is returned
• Build - the hypercube is built by iteratively generating a number of candidate points and picking the best one
• CP Sweep - Columnwise Pairwise Swapping optimizes a random hypercube by swapping two elements from each column in an optimal way
• Genetic - Genetic algorithm-based optimization where a population of hypercubes evolves over multiple generations, where the top half of hypercubes survive each generation and produce offspring via column exchanges with occasional random mutations (swaps within columns).

Each of the methods optimizes a given objective function, to make the hypercube as space-filling as possible. Following goals are supported:

• Maximin - Maximize the minimum distance between two points
• Opt - Minimize difference between two closest points to the theoretical optimal distance N/(N^(1/d))
• Sum Inv - Minimize sum of inverse squared Euclidean distances between points
• Avg Dist - Maximize average Euclidean distance between points

When using the Maximin and Avg Dist goals, the negative value of the goal is minimized, so that all goals are optimized in the same direction.

The return_type parameter controls how the matrix should be returned:

• buckets: the matrix contains integers in range (0,N]
• numerical: the matrix contains values from uniform distribution in range [0,1]

Method Comparison

The following table contains the results of running the following snippet:

N = 10
d = 5
reps = 10
seed = 1410
methods = list("random", "cp_sweep", "build", "genetic")
goals = list("opt","maximin","avg_dist","sum_inv")
for(j in 1:length(goals)){
    goal = as.scalar(goals[j])
    for (i in 1:length(methods)) {
        method_name = as.scalar(methods[i])
        [M,G] = lhs(N = N, d = d, reps = reps, method = method_name, goal = goal, seed = seed)
        print("Method: " + method_name + ", Goal: " + goal + ": G = " + G)
    }
}

Method	Mean Distance	Min Distance	Divergence from Opt	Sum of inverse squared distances
build	9.19523	5.19615	0.22681	0.68072
cp_sweep	9.48126	7.21110	0.01498	0.52165
genetic	9.43672	6.32456	0.01498	0.56041
random	9.32597	3.31662	2.99295	0.67605

Note: the table contains positive values for Min Distance and sum of inverse squared distance for readability, in practice the function return negative values.

In all categories the CP-Sweep method provides the best results, it is however the most computationally expensive method, since all possible swaps per column have to be compared.

Benchmarking against Random Sampling

To evaluate the effectiveness of LHS the following benchmarking experiment has been conducted:

systemds /systemds/src/test/scripts/functions/builtin/lhs_gmm_benchmark.dml --nvargs method=cp_sweep reps=10 goal=sum_inv seed=42

1. Generate a true 5-dimensional Multivariate Normal distribution with a randomized mean vector $\mu_{\text{true}}$ and a symmetric positive semi-definite covariance matrix $\Sigma_{\text{true}}$.
2. Draw training samples of sizes $N \in [10,13,15,20,25,30,40,50,75,100]$ using Simple Random Sampling from standard normal distribution and Latin Hypercube Sampling (LHS) mapped to normal space qnorm() function. Both samples are then mapped to the true distribution
3. Fit a single-component GMM (using model = "VVV" to estimate full covariance) on both training datasets. With one component, this reduces to standard full-covariance maximum-likelihood estimation.
4. Evaluation metrics:

• Mean Log-Likelihood (LL): Tracks both training LL and out-of-sample LL evaluated on an independent test set of $10^6$ samples (likelihood of the data originating from the estimated distribution).
• KL Divergence ($D_{KL}$): Computes the information loss between the true distribution ($P$) and the estimated GMM ($Q$) distribution

Results

The results show that the KL-Divergence converges to nearly zero already with N < 20, while random sampling requires a larger number of samples to converge. Similarly in case of the mean LL, for small N there is a wide gap between train and test mean LL for random sampling, while LHS converges immediately, showcasing the ability of LHS to draw samples representative of the underlying distribution, even with very small N.

References:
[1] Stocki, Rafal. "A method to improve design reliability using optimal Latin hypercube sampling." Computer Assisted Mechanics and Engineering Sciences 12.4 (2005): 393.
[2] Beachkofski, Brian, and Ramana Grandhi. "Improved distributed hypercube sampling." 43rd AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference. 2002.