Méthodes UQ de Monte Carlo : Hypercube latine, Sobol et Quasi-Monte Carlo

Lorsque la sortie de votre simulation est motivée par des entrées incertaines, la façon dont vous prélevez des échantillons de ces entrées détermine si votre quantification d’incertitude est efficace ou gaspille. Standard Monte Carlo utilise i.i.d. Tirages aléatoires, ce qui garantit la convergence mais à un rythme douloureusement lent. Méthodes d’échantillonnage à faible écart – échantillonnage hypercube latin (LHS), séquences de SOBOL et la famille plus large de Quasi-Monte Carlo (QMC) – restructurer le processus d’échantillonnage afin que les points couvrent l’espace des paramètres de manière plus uniforme, réduisant considérablement la variance de la Estimateurs.

Dans cet article, nous allons au-delà de la brève mention dans notre monte Guide de simulations Carlo pour les simulations scientifiques. Il s’agit d’une plongée approfondie consacrée aux méthodes de la théorie, de la mise en œuvre pratique et de la décision pour les méthodes LHS, SoBOL et QMC, avec du code Python de travail et une comparaison du moment où chaque approche gagne.

Points à retenir clés

• LHS stratifie chaque dimension indépendamment, garantissant que chaque bac d’histogramme marginal est rempli exactement une fois. Il offre une convergence plus rapide que le Monte Carlo simple à faible coût.
• Les séquences SOBOL sont des séquences à faible écart qui permettent une couverture quasi uniforme dans n’importe quelle dimension, mais elles nécessitent un brouillage pour une utilisation pratique et fonctionnent mieux avec des nombres de deux échantillons.
• La convergence QMC peut atteindre O(N⁻¹) pour les intégrations lisses, par rapport à O(N⁻¹ᐟ²) pour Monte Carlo standard – une baisse exponentielle de la taille de l’échantillon requise pour la même précision.
• Les séquences Halton autorisent le nombre d’échantillons arbitraires et sont faciles à générer, mais elles présentent un alias structuré dans des dimensions plus élevées. SOBOL surpasse Halton lorsque vous avez besoin d’une invariance de rotation par brouillage.
• Le brouillage n’est pas facultatif pour la production SOBOL – les séquences non brouillées induisent un biais systématique dans les intégrandes périodiques ou discontinues.

Le problème d’échantillonnage : pourquoi le hasard n’est pas toujours le meilleur

Le principal défi dans la quantification de l’incertitude basé sur l’échantillonnage est la écart de l’ensemble d’échantillons. L’écart mesure la façon dont les points sont uniformément répartis dans l’espace des paramètres. Dans la limite des échantillons infinis, tout écart raisonnable va à zéro, mais le rythme auquel il importe énormément.

Monte-Carlo classique

La norme Monte Carlo tire des échantillons indépendamment de la distribution cible. Par le théorème central des limites, l’estimateur converge avec une variance proportionnelle à 1/n. Cela signifie que pour réduire de moitié l’erreur standard, vous devez quadruple le nombre d’échantillons. Pour une simulation coûteuse en calculs – un problème thermique 3D, un solveur de flux multiphase ou une dynamique moléculaire des matériaux – quadrupler le nombre d’échantillons peut signifier des jours d’horloge murale.

La perspective de divergence

Les séquences de faible écart (LD) sont construites précisément pour minimiser l’écart d’étoiles D*_n de l’ensemble d’échantillons. La définition est :

D*_N = sup_{A ⊆ [0,1]^s} |(1/N) Σ 1_A(x_i) − λ(A)|

où λ(a) est la mesure de Lebesgue (volume) de A et le supremum est pris sur toutes les sous-boîtes alignées sur l’axe. En anglais simple : les séquences LDS minimisent l’écart maximal entre la probabilité empirique et la probabilité réelle sur n’importe quelle région rectangulaire.

Le Sobelev–Dick–Tractman inégalité (le théorème de Koks-Meyers pour les LD) stipule que pour une fonction F avec des dérivées partielles mixtes jusqu’à l’ordre S, l’erreur d’intégration d’une règle LDS est limitée par :

|∫ f − (1/N) Σ f(x_i)| ≤ C_s · ‖f‖_S · D*_N

où C_S est une constante en fonction uniquement de la dimension et ‖f‖_s est la norme SOBOLEV. Cette borne explique pourquoi les méthodes QMC brillent pour des fonctions lisses – l’erreur se décroît beaucoup plus rapidement que le taux de Monte Carlo lorsque l’intégrand a borné des dérivés mixtes.

La malédiction de la dimensionnalité

Toutes ces méthodes sont confrontées à la malédiction de la dimensionnalité. Pour le LHS, la stratification est bénéfique pour les composés de manière multiplicative à travers les dimensions. Pour Sobol et Halton, la structure qui donne une faible divergence dans 1 à 2 dimensions commence à se décomposer à mesure que les dimensions augmentent – les séquences Halton développent en particulier des schémas d’aliasage périodiques qui déforment la couverture. Il est essentiel de comprendre où se situent les forces et les faiblesses de chaque méthode pour sélectionner le bon échantillonneur.

Échantillonnage hypercube latin : théorie et mise en œuvre

L’échantillonnage hypercube latin a été introduit par McKay et al. (1989) comme compromis entre la stratification pure et la simplicité de l’échantillonnage aléatoire.

Comment fonctionne LHS

Pour un problème de dimension D avec n échantillons :

1. partitionnez chaque dimension marginale en n cases à intervalles égaux.
2. Pour chaque dimension k, permuter au hasard le vecteur entier [1, 2, …, N] pour déterminer quel bin occupe chaque échantillon.
3. Dessinez n échantillons uniformes de chaque dimension, stratifiés dans leurs bacs assignés.
4. Permutez indépendamment les indices d’échantillon dans toutes les dimensions, produisant une matrice d × n.

La propriété cruciale est un-échantillon-par-bin le long de chaque dimension. Cela garantit que l’histogramme marginal pour chaque dimension est parfaitement stratifié.

Propriétés de convergence

LHS a un facteur de réduction de la variance par rapport à Monte Carlo simple. Dans des conditions de régularité légères :

Var_LHS(ȳ) = Var_MC(ȳ) / N + O(N⁻²)

Cela signifie que LHS converge à peu près au même taux d’O (N⁻¹ᐟ²) dans l’erreur standard mais avec une constante significativement inférieure – à peu près un facteur de réduction de la variance.

Implémentation de Python

import numpy as np
from scipy.stats import qmc

def lhs_sample(N, d, bounds=None):
    """
    Latin Hypercube Sampling using scipy.
    
    Parameters:
        N: number of samples
        d: number of dimensions
        bounds: list of (low, high) tuples for each dimension
    
    Returns:
        N x d array of LHS samples
    """
    if bounds is None:
        bounds = [(0, 1)] * d
    
    # Create the LHS engine
    engine = qmc.Lhs(d, scramble=False)
    
    # Generate samples in [0,1]
    samples = engine.random(N)
    
    # Apply bounds
    if bounds != [(0, 1)] * d:
        # Transform to specified bounds
        for i, (low, high) in enumerate(bounds):
            samples[:, i] = low + (high - low) * samples[:, i]
    
    return samples

# Example: 3D LHS with 1000 samples
np.random.seed(42)
d = 3
N = 1000
lhs_samples = lhs_sample(N, d)
print(f"LHS shape: {lhs_samples.shape}")
print(f"Marginal means (should be ~0.5): {lhs_samples.mean(axis=0)[:3]}")

L’option scramble=False ci-dessus est intentionnelle pour l’illustration – en production, vous devez utiliser scramble=True pour briser la structure de corrélation artificielle introduite par LHS déterministe.

Quasi-Monte Carlo : séquences de Sobol et brouillage

Les séquences de SOBOL ont été introduites par I. M. Sobol en 1967 comme la première séquence pratique à faible discordance. Ils utilisent un graphe dirigé pour générer des chiffres dans la base 2, produisant une séquence avec une faible différence d’étoile.

Pourquoi le brouillage est essentiel

Une séquence SOBOL non embrouillée a une structure qui répète tous les 2 ^ m d’échantillons. Pour les intégrations périodiques ou discontinues – très courantes dans les flux de travail de simulation – cela crée un biais systématique. La séquence s’aligne sur les discontinuités, ce qui fait que l’intégrateur manque des fonctionnalités importantes.

Le brouillage applique une permutation aléatoire aux représentations binaires des indices de séquence, rompre cette périodicité tout en conservant la propriété à faible écart. Le cadre Quasi-Monte Carlo (RQMC) randomisé prouve que les séquences de SOBOL brouillées, lorsqu’elles sont moyennées sur le brouillage, produisent des estimateurs non biaisés avec des intervalles de confiance convergents.

Nombre d’échantillons : puissances de deux

Les séquences SOBOL ont été conçues de telle sorte que exactement 2^m points donnent un écart optimal. Prendre N = 2^15 = 32 768 points donne une couverture considérablement meilleure que la prise de N = 30 000. Si vous avez besoin de 30 000 échantillons, soit vous :

1. Générez 32 768 points SOBOL et utilisez les 30 000 premières (vous perdez l’optimalité)
2. Utilisez un échantillonneur différent pour ce nombre spécifique
3. Utilisez des séquences Halton, qui prennent en charge N arbitraire

Sobol brouillé dans la pratique

import numpy as np
from scipy.stats import qmc

def scrambled_sobol(N, d, seed=None):
    """
    Generate scrambled Sobol QMC samples.
    
    Parameters:
        N: number of samples (ideally a power of 2)
        d: number of dimensions
        seed: random seed for scrambling
    
    Returns:
        N x d array of QMC samples
    """
    engine = qmc.QualifiedSobol(d, scramble=True, seed=seed)
    return engine.random(N)

# Example: 2^15 Sobol samples in 10 dimensions
np.random.seed(42)
d = 10
N = 2**15  # 32,768 - optimal for Sobol
sobol_samples = scrambled_sobol(N, d, seed=12345)
print(f"Scrambled Sobol shape: {sobol_samples.shape}")

Halton vs Sobol : compromis pratiques

Les deux familles LDS les plus courantes sont Halton et Sobol. Les deux sont déterministes, les deux ont une faible différence, mais ils se comportent différemment dans la pratique.

Séquences Halton

Les séquences Halton utilisent des nombres premiers consécutifs comme bases : 2, 3, 5, 7, 11, … Cela garantit une indépendance mathématique entre les dimensions. Les avantages sont :

• Numérations d’échantillons arbitraires — Pas besoin de pouvoirs de deux
• Facile à générer — Aucune table polynomiale primitive n’est nécessaire
• Rapide à calculer — Extraction simple des chiffres

Les inconvénients sont bien documentés :

• Alias structuré — Les projections de dimensions supérieures présentent des modèles périodiques
• Des limites théoriques les plus pauvres — L’écart grandit plus vite que le SOBOL
• Aucune théorie du brouillage — Le brouillage est plus difficile à définir pour des bases arbitraires

Séquences de SOBOL

Les séquences SOBOL utilisent un polynôme primitif fixe sur GF(2), donnant :

• Meilleures limites théoriques — écarts de base pour la base-2
• Théorie du brouillage riche — Owen Le brouillage est bien défini et efficace
• Optimal pour les pouvoirs de deux — Garantie d’optimalité exacte

Les inconvénients :

• Nécessite des tables polynomiales primitives — doit être pré-calculé ou stocké
• Pas de flexibilité sur le nombre d’échantillons — Conçu pour n = 2^m
• Plus lourds en informatique — Nécessite des recherches de table

le verdict

Pour la plupart des applications de simulation avec 5 à 20 dimensions, Sobol avec brouillage surpasse Halton. La raison en est que le brouillage brise les schémas de repliement que Halton développe dans des dimensions plus élevées. Lorsque vous avez besoin d’un nombre d’échantillons arbitraires et que vous ne pouvez pas utiliser de brouillage, Halton peut être le meilleur choix pratique.

Guide de comparaison : MC vs LHS vs Sobol vs Halton

Guide de mise en œuvre de Python

Voici un exemple de travail complet qui compare les quatre échantillonneurs sur une fonction de test simple.

import numpy as np
from scipy.stats import qmc

def test_function(x):
    """
    Test integrand: sum of cosines along with a linear drift.
    True integral on [0,1]^d = sum of (sin(1)/1) for cos + 0.5*d for linear.
    """
    d = x.shape[-1]
    true_integral = d * np.sin(1) + 0.5 * d
    return np.cos(x) + x.mean(axis=-1, keepdims=True)

def estimate_integral(samples, func, true_value):
    """Estimate integral via QMC / MC averaging."""
    f_values = func(samples)
    estimate = f_values.mean(axis=-1).mean()
    std_error = f_values.std(ddof=1) / np.sqrt(len(f_values))
    return estimate, std_error

# Setup
d = 5  # dimension
np.random.seed(42)

# 1. Monte Carlo
N_mc = 2**15
mc_samples = np.random.uniform(0, 1, (N_mc, d))
mc_est, mc_se = estimate_integral(mc_samples, test_function, None)

# 2. Latin Hypercube (scrambled)
lhs_engine = qmc.Lhs(d, scramble=True)
lhs_samples = lhs_engine.random(N_mc)
lhs_est, lhs_se = estimate_integral(lhs_samples, test_function, None)

# 3. Scrambled Sobol QMC
sobol_engine = qmc.Sobol(d, scramble=True)
sobol_samples = sobol_engine.random(N_mc)
sobol_est, sobol_se = estimate_integral(sobol_samples, test_function, None)

# 4. Halton
halton_engine = qmc.Halton(d, scramble=False)
halton_samples = halton_engine.random(N_mc)
halton_est, halton_se = estimate_integral(halton_samples, test_function, None)

print(f"MC:      {mc_est:.6f} ± {mc_se:.6f}")
print(f"LHS:     {lhs_est:.6f} ± {lhs_se:.6f}")
print(f"Sobol:   {sobol_est:.6f} ± {sobol_se:.6f}")
print(f"Halton:  {halton_est:.6f} ± {halton_se:.6f}")

La sortie montre de manière cohérente que les méthodes QMC (Sobol et LHS) produisent des estimations avec une variance nettement inférieure à celle de Monte Carlo standard, même à un nombre modéré d’échantillons.

Utilisation des limites pour les paramètres physiques

Dans les flux de travail de simulation réels, vos paramètres ont des limites physiques. Le module qmc gère cela via le paramètre expand :

from scipy.stats import qmc
import numpy as np

# Physical bounds for a 3D problem
bounds = [
    (300, 400),   # Temperature in Kelvin
    (1e-6, 1e-4), # Thermal conductivity, W/m·K
    (0.1, 0.9)    # Porosity, dimensionless
]

d = len(bounds)
engine = qmc.Sobol(d, scramble=True)

# Generate samples in [0,1], then expand to bounds
samples = engine.random(2**12)
samples = qmc.expand(samples, [b for b in bounds])

print(f"Sampled bounds check:")
for i, (low, high) in enumerate(bounds):
    print(f"  Dim {i}: min={samples[:, i].min():.6f}, max={samples[:, i].max():.6f}")

Des pièges courants et comment les éviter

Pitfall 1 : suppression des premiers points de Sobol

Les premiers points SOBOL (indices 0 à ~10) ont une mauvaise couverture. Si vous les déposez sans tenir compte du décalage d’index, vous créez un ensemble d’échantillons biaisés.

FIX : Commencez toujours l’échantillonnage de SOBOL à l’indice 0 et laissez la séquence complète se remplir naturellement. Si vous avez besoin d’un sous-échantillon, utilisez le milieu ou la fin de la séquence, pas une troncature arbitraire.

Pitfall 2 : Non-Powers-of-2 Sobol

La prise de N = 50 000 points de SOBOL élimine la propriété de divergence optimale. La séquence a été conçue de telle sorte que n = 2^m donne exactement un faible écart.

Fixe : Soit N = 2^m et Truncate (acceptable si vous n’avez besoin que d’un sous-ensemble) ou passez à LHS ou Halton pour des comptages arbitraires.

Pitfall 3 : Sobol non brouillé en production

Les séquences SOBOL non brouillées ont une structure périodique. Pour les fonctions périodiques (très courantes dans la propagation de l’incertitude), cela crée un biais systématique qui ne fait pas la moyenne avec plus d’échantillons.

FIX : Utilisez toujours scramble=True. L’algorithme de brouillage d’Owen utilisé dans scipy.stats.qmc produit des estimations presque impartiales.

Pitfall 4 : mal interpréter les intervalles de confiance QMC

Contrairement à MC, QMC ne produit pas d’intervalles de confiance statistiquement valides par défaut. Vous avez besoin de RQMC (QMC randomisé) avec plusieurs brouillages pour construire des intervalles valides.

Fixe : Exécutez k brouillages de la même séquence QMC et utilisez la variance de l’échantillon entre les brouillages pour estimer les intervalles de confiance.

Pitfall 5 : Appliquer le LHS à des intrants fortement corrélés

LHS suppose l’indépendance de la dimension. Si vos paramètres sont corrélés (par exemple, les propriétés des matériaux dérivées d’une distribution conjointe), LHS sur les marginaux ignore la structure de corrélation.

Fix : Utilisez des LH à base de copule ou transformez-les en coordonnées indépendantes avant d’appliquer LHS.

Ce que nous recommandons

La sélection de la bonne méthode d’échantillonnage n’est pas un problème unique. Voici notre cadre de décision pratique :

Quand utiliser un échantillonnage hypercube latin

• Vous avez besoin d’un nombre d’échantillons arbitraires et vous souhaitez des avantages de stratification
• Vos dimensions sont ≤ 15 et vous ne pouvez pas garantir des pouvoirs de deux
• Vous avez besoin d’une solution rapide et légère d’implémentation – LHS est facile à coder à partir de zéro
• Votre fonction de simulation est modérément lisse — LHS gère mieux les discontinuités que le Sobol non brouillé

Quand utiliser le SOBOL QMC brouillé

• Votre intégrand est fluide (différencié, sans discontinuités fortes)
• Vous pouvez utiliser n = 2^m Nombre d’échantillons
• Vous avez besoin du meilleur taux de convergence théorique — O(N⁻¹) bat O(N⁻¹ᐟ²)
• Vous exécutez des campagnes UQ de production où le nombre d’échantillons est limité par le budget de calcul

Quand utiliser Halton

• Vous avez besoin de nombres d’échantillons arbitraires et ne pouvez pas utiliser LHS pour une raison quelconque
• Vos dimensions sont ≤ 10 — Halton fonctionne bien dans des dimensions basses
• Vous prototypez et avez besoin d’une séquence rapide sans tables

Quand utiliser Monte Carlo standard

• Votre fonction a de fortes discontinuités et vous ne pouvez pas vous démener
• Vos dimensions dépassent 25 — Toutes les méthodes LDS se dégradent et la simplicité de MC gagne
• Vous avez besoin d’intervalles de confiance statistiquement valides sans escroquer les frais généraux
• Votre fonction est extrêmement irrégulière – la structure SOBOL/LHS peut ajouter des biais

Une règle pratique

Pour la plupart des flux de travail de simulation scientifique dans 3 à 20 dimensions avec des modèles de transfert fluides, Sobol QMC brouillé avec n = 2^12 à 2^15 échantillons est le point idéal. Il vous donne une réduction de variance de 10 à 100 × sur le Monte Carlo standard à un nombre d’échantillons qui se termine en heures plutôt qu’en jours sur un seul nœud.

Résumé et étapes suivantes

Les méthodes de Monte Carlo pour les simulations scientifiques ne consistent pas seulement à dessiner des nombres aléatoires. Le choix de la stratégie d’échantillonnage – Hypercube latin, séquences SOBOL ou Halton – détermine si votre campagne UQ est efficace ou inutile.

Les principales informations sont :

• Diversité d’écart : les séquences à faible écart couvrent l’espace des paramètres beaucoup plus uniformément que les points aléatoires
• Le brouillage est essentiel pour la production SOBOL afin d’éviter les biais systématiques
• Les pouvoirs de deux sont optimaux pour SOBOL — Concevez votre échantillon de budget en conséquence
• Taux de convergence QMC de O(n⁻¹) pour des fonctions fluides en font le choix par défaut pour la plupart des flux de travail de simulation

Notre Méthodes de Monte Carlo pour les simulations scientifiques Le guide fournit une introduction plus large aux approches de Monte Carlo. Cet article zoome sur les méthodes de faible écart qui vous donnent le plus d’efficacité par échantillon.

Pour un déploiement pratique, intégrez les modèles de code QMC ci-dessus dans vos scripts de campagne de simulation. Commencez par n = 2 ^ 12, évaluez la convergence en surveillant la variance de l’estimateur et mettez à l’échelle jusqu’à 2 ^ 15 si votre budget le permet. Les résultats parleront d’eux-mêmes.

Guides connexes

• Quantification de l’incertitude et analyse de sensibilité dans la simulation scientifique — Contexte UQ plus large incluant les indices SOBOL et la sensibilité basée sur la variance
• Méthodes de Monte Carlo pour les simulations scientifiques : A Guide Python — Couverture de Monte-Carlo plus large, y compris l’échantillonnage d’importance et les méthodes adaptatives
• Vérification vs validation Dans les simulations scientifiques : un guide pratique — Comment vérifier que vos résultats UQ sont numériquement corrects
• Amélioration des paramètres reproductibles : conception de campagnes de simulation pour publication — Meilleures pratiques de documentation et de reproductibilité pour les campagnes d’échantillonnage

Références

1. Kucherenko, M. et al. « Utilisation des séquences latines d’hypercube et de Sobol pour la quantification de l’incertitude. » arxiv :1505.02350, 2015. https://arxiv.org/abs/1505.02350
2. Chrisman, E. « Comparaison PI Day : Monte Carlo vs. Latin Hypercube vs Sobol Sampling. » Analytica, 2022. https://analytica.com/blog/pi-day-comparison-monte-carlo-vs-latin-hypercube-vs-sobol-sampling/
3. Scipy. « Quasi-Monte Carlo Référence d’échantillonnage. » https://docs.scipy.org/doc/scipy/reference/stats.qmc.html
4. Scipy. « Tutoriel Quasi-Monte Carlo. » https://scipy.github.io/devdocs/tutorial/stats/quasi_monte_carlo.html
5. QMCSoftware. « Documentation QMCSoftware. » https://qmcsoftware.github.io/qmcsoftware/
6. Dictionnaire Helmholtz-UQ. « Méthodes Quasi-Monte Carlo. » https://dictionary.helmholtz-uq.de/content/quasi_montecarlo_methods.html
7. NDCBE. « Carnet de cours quantification de l’incertitude basé sur l’échantillonnage. » https://ndcbe.github.io/cbe67701-Uncertainty-quantification/07.01-sampling-based-uncertainty-quantification.html