jancies2301 opened a new pull request, #2545:
URL: https://github.com/apache/systemds/pull/2545
### 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:
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
<img width="600" height="371" alt="Mean Log Likelihood"
src="https://github.com/user-attachments/assets/c47ee62a-045f-4763-90c4-93e918dba695";
/>
<img width="600" height="371" alt="KL-Divergence"
src="https://github.com/user-attachments/assets/14711109-4949-4fa7-9dd7-1e2e46b0b2f4";
/>
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.
--
This is an automated message from the Apache Git Service.
To respond to the message, please log on to GitHub and use the
URL above to go to the specific comment.
To unsubscribe, e-mail: [email protected]
For queries about this service, please contact Infrastructure at:
[email protected]
