MST-Direct at Scale: Multivariate and Conditional Geostatistical Simulation via Sinkhorn Optimal Transport

Technische Samenvatting: MST-Direct op Schaal

Probleemstelling
Geostatistische simulatie beoogt meerdere realisaties van geologische modellen te genereren die ruimtelijke kenmerken reproduceren en onzekerheidskwantificering ondersteunen. Een aanzienlijke uitdaging ontstaat wanneer geologische relaties sterk niet-lineair en heteroscedastisch zijn, waardoor klassieke multivariate methoden die vertrouwen op lineaire correlatiecoëfficiënten ontoereikend zijn. Bestaande benaderingen, zoals Stepwise Conditional Transformation (SCT), Projection Pursuit Multivariate Transform (PPMT) en Direct Multivariate Simulation (DMS), proberen dit op te lossen door variabelen te transformeren naar een Gaussische ruimte en weer terug. Deze transformaties zijn echter alleen asymptotisch exact, en hun fouten nemen toe naarmate het aantal variabelen stijgt.

Vorig werk introduceerde MST-Direct (Matching via Sinkhorn Transport), wat simulatie frameert als een entropie-geregulariseerd optimaal transportprobleem (OT). Hoewel deze benadering 100% vormbehoud bereikte voor bivariate relaties, was de oorspronkelijke formulering beperkt tot kleine roosters (bijv. 25 × 25), onvoorwaardelijke instellingen en bivariate data. Het liet expliciet drie onopgeloste beperkingen open: schaalbaarheid voorbij ~10.000 knopen, uitbreiding naar meer dan twee variabelen, en de mogelijkheid om conditionele simulatie met harde data uit te voeren.

Methodologie
Dit artikel breidt MST-Direct uit om deze drie geïdentificeerde beperkingen aan te pakken via een schaalbaar, multivariat framework. Het kernprincipe blijft het framen van simulatie als een optimaal transportprobleem waarbij doelwaarde-tuples worden gematcht op een ruimtelijke template.

1. Schaalbare Matching: Om grote roosters (bijv. 40.000 knopen) aan te kunnen, vervangen de auteurs de dichte $O(n^2)$ koppeling door een ijle, kandidaat-beperkte Sinkhorn-matcher. Elke roosterlocatie wordt beperkt tot zijn $C$ dichtstbijzijnde doel-tuples (geïdentificeerd via een k-d tree), waardoor de geheugencomplexiteit wordt gereduceerd tot $O(nC)$. De oplossing wordt gevonden met behulp van het Sinkhorn-algoritme in het log-domein. Een greedily assignment rondt de zachte koppeling af tot een permutatie, met een herstelstap voor locaties die zonder vrije kandidaat zijn achtergebleven.
2. Multivariaat Ruggegraat: Om de methode uit te breiden naar $d$ variabelen, wordt de doelwolk gematcht op een onafhankelijke FFT-MA (Fast Fourier Transform Moving Average) Gaussische ruggegraat. Deze ruggegraat draagt de voorgeschreven variogram (bijv. sferisch). Het matchingproces draagt de ruimtelijke structuur van de ruggegraat over naar elke variabele zonder eisen te stellen aan isotropie of stationariteit, hoewel een isotroop sferisch model hier wordt gebruikt voor vergelijkbaarheid.
3. Conditionele Simulatie: Harde data worden gehonoreerd via een tweestaps-proces. Eerst wordt de Gaussische ruggegraat geconditioneerd aan de locaties van de harde data met behulp van eenvoudige kriging, waarbij wordt gewaarborgd dat de voorgeschreven covariantie elders behouden blijft. Ten tweede worden de locaties van de harde data "gepin": hun specifieke tuples worden vastgezet en verwijderd uit de matching-pool, terwijl de resterende locaties worden gesimuleerd. Dit zorgt voor een exacte reproductie van de waarden van de harde data.

De relationele term, die de kostenfunctie aanvult om ervoor te zorgen dat ruimtelijk aangrenzende locaties wederzijds vergelijkbare tuples ontvangen, wordt toegepast via verfijningsrondes die query-punten naar het gemiddelde van hun ruimtelijke buren trekken.

Belangrijkste Bijdragen
Het artikel presenteert vier primaire bijdragen:

• Een schaalbare MST-Direct matcher die gebruikmaakt van ijle, kandidaat-beperkte Sinkhorn met greedy bijwerking, in staat om 200 × 200 (40.000-knopen) en 100 × 100 roosters te verwerken in minder dan een minuut.
• Een multivariaat formulering die doelwolken matcht op een onafhankelijke FFT-MA ruggegraat, succesvol gedemonstreerd op zes variabelen.
• Een conditioneel simulatiemechanisme dat harde data exact respecteert door data-tuples te pinnen en de ruggegraat te conditioneren via kriging.
• Een head-to-head validatie tegen de Projection Pursuit Multivariate Transform (PPMT) op een 6-variabele, heteroscedastische, niet-lineaire referentieverdeling, die zowel onvoorwaardelijke als conditionele instellingen dekt.

Resultaten
De methode werd gevalideerd op dezelfde 6-variabele, heteroscedastische, sterk niet-lineaire referentieverdeling als gebruikt door DMS.

• Onvoorwaardelijke Simulatie (200 × 200 rooster): MST-Direct produceerde realisaties die ononderscheidbaar waren van de referentieverdeling, waarbij een nul-histogramfout (Mean Squared Error, MSE) werd bereikt. In contrast hiermee herstelde PPMT de grove structuur, maar vervormde fijne niet-lineaire verbindingen, resulterend in MSE's in de orde van grootte $10^{-6}$ tot $10^{-5}$. Wat betreft de ruimtelijke statistiek mat MST-Direct de variogram-sill exact, terwijl PPMT deze met 10–20% overschoot.
• Conditionele Simulatie (100 × 100 rooster, 200 harde data): MST-Direct honoreerde de harde data exact (maximale absolute fout van 0), terwijl PPMT deze honoreerde tot machineprecisie ($\sim 10^{-14}$). MST-Direct reproduceerde opnieuw de doel-jointverdeling met een nul-histogramfout, terwijl PPMT een benadering bleef. Beide methoden reproduceerden de ruimtelijke variogrammen, waarbij MST-Direct de sill exact volgde en PPMT deze licht onderschatte.
• Efficiëntie: De ijle matcher handteerde het 40.000-knopen rooster succesvol, een schaal waarbij de dichte matcher onhaalbaar zou zijn.

Betekenis en Claims
Het artikel claimt de drie openstaande beperkingen van de oorspronkelijke MST-Direct formulering op te lossen. De auteurs stellen dat omdat het transportplan een permutatie van de doel-tuples blijft, de multivariate joint verdeling exact wordt behouden, waardoor de asymptotische fouten en dimensionaliteitsafhankelijke factorisatiefouten geassocieerd met SCT, PPMT en DMS worden vermeden.

De methode wordt gepresenteerd als zeer geschikt voor toepassingen waarbij complexe, sterk niet-lineaire multivariate afhankelijkheden aanwezig zijn. Hoewel de aanpak erin slaagt de doelstatistieken en ruimtelijke modellen te reproduceren, merken de auteurs een residu-artefact op: een milde short-lag component in onvoorwaardelijke variogrammen, die controleerbaar is via het aantal kandidaten en relationele passes. Het paper concludeert dat anisotrope variogrammen toekomstig werk blijven. Een open-source implementatie (v2.0.0) is beschikbaar als het mst-direct Python-pakket.