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

# 2603.18036), which introduced MST-Direct in the bivariate, unconditional, small-grid setting; this paper contributes the scalable, multivariate, and conditional extensions.
Source: [https://arxiv.org/html/2606.07578](https://arxiv.org/html/2606.07578)
## MST\-Direct at Scale: Multivariate and Conditional Geostatistical Simulation via Sinkhorn Optimal Transport††thanks:Follow\-up to the author’s prior work \(arXiv:2603\.18036\), which introduced MST\-Direct in the bivariate, unconditional, small\-grid setting; this paper contributes the scalable, multivariate, and conditional extensions\.

###### Abstract

This paper extends MST\-Direct, the Matching\-via\-Sinkhorn\-Transport approach to multivariate geostatistical simulation, from the bivariate, unconditional, small\-grid setting of the original formulation to the*multivariate*,*conditional*, and*large\-grid*regime\. We address the three open limitations identified in that work: \(i\) scalability beyond a few thousand nodes, through a sparse, candidate\-restricted Sinkhorn matcher withO​\(n​C\)O\(n\\,C\)memory; \(ii\) extension to many variables, by matching the target value tuples onto an independent FFT\-MA Gaussian backbone that carries a prescribed variogram; and \(iii\) hard\-data conditioning, by pinning the data tuples to their locations while conditioning the backbone by kriging\. Because the transport plan remains a permutation of the target tuples, the multivariate joint distribution is preserved*exactly*\. We validate the method on the same 6\-variate, heteroscedastic, strongly non\-linear reference distribution used by Direct Multivariate Simulation \(DMS\), in both unconditional \(200×200200\\times 200\) and conditional \(100×100100\\times 100, 200 hard data\) settings, and we benchmark it against the Projection Pursuit Multivariate Transform \(PPMT\)\. MST\-Direct reproduces the joint distribution with zero histogram error, honours the hard data exactly, and reproduces the prescribed spatial correlation, whereas PPMT remains an approximation\.

## IIntroduction

Geostatistical simulation generates multiple realizations of geological models that replicate the spatial features of the data and support uncertainty quantification\[[5](https://arxiv.org/html/2606.07578#bib.bib14),[2](https://arxiv.org/html/2606.07578#bib.bib13)\]\. Many geological relationships are strongly non\-linear and heteroscedastic and cannot be characterized by the linear correlation coefficient on which classical multivariate methods rely\. The Stepwise Conditional Transformation \(SCT\)\[[7](https://arxiv.org/html/2606.07578#bib.bib5)\], the Projection Pursuit Multivariate Transform \(PPMT\)\[[1](https://arxiv.org/html/2606.07578#bib.bib4)\], and Direct Multivariate Simulation \(DMS\)\[[4](https://arxiv.org/html/2606.07578#bib.bib1)\]address this by transforming the variables to a Gaussian space and back; the transform is exact only asymptotically and its error grows with the number of variables\.

In previous work we introduced*MST\-Direct*\(Matching via Sinkhorn Transport\)\[[10](https://arxiv.org/html/2606.07578#bib.bib2)\], which casts simulation as an entropy\-regularized optimal\-transport \(OT\) problem\[[12](https://arxiv.org/html/2606.07578#bib.bib23),[8](https://arxiv.org/html/2606.07578#bib.bib22),[3](https://arxiv.org/html/2606.07578#bib.bib21)\]: the target value tuples are*matched*onto a spatial template through the Sinkhorn algorithm\[[11](https://arxiv.org/html/2606.07578#bib.bib20)\]augmented with a relationalkk\-nearest\-neighbour term, and the resulting permutation reproduces the spatial structure while preserving the joint distribution exactly\. That formulation achieved100%100\\%shape preservation on five complex*bivariate*relationships, but was demonstrated only on a25×2525\\times 25grid against Gaussian\-copula and LU\-decomposition baselines, and explicitly left three problems open\[[10](https://arxiv.org/html/2606.07578#bib.bib2)\]: scalability beyond∼\\sim10 000 nodes, extension to more than two variables, and conditional simulation with hard data\.

This paper resolves those three limitations and benchmarks the method on the DMS validation problem\. Our contributions are:

1. 1\.a*scalable*MST\-Direct matcher — sparse, candidate\-restricted Sinkhorn with greedy bijection completion — that runs the200×200200\\times 200\(40 00040\\,000\-node\) and100×100100\\times 100grids in well under a minute;
2. 2\.a*multivariate*formulation that matches the target cloud onto an independent FFT\-MA Gaussian backbone\[[6](https://arxiv.org/html/2606.07578#bib.bib9)\]carrying the prescribed variogram, demonstrated on six variables;
3. 3\.*conditional*simulation that honours hard data exactly, by pinning the data tuples and conditioning the backbone by simple kriging;
4. 4\.a head\-to\-head validation against PPMT on the same 6\-variate, heteroscedastic, non\-linear reference distribution as DMS\[[4](https://arxiv.org/html/2606.07578#bib.bib1)\]\.

## IIMethod

We first summarize the MST\-Direct principle introduced in\[[10](https://arxiv.org/html/2606.07578#bib.bib2)\]; the rest of this section presents the three extensions that constitute the present contribution \(scalable matching, the multivariate backbone, and conditioning\)\. Let\{𝐳\(k\)\}k=1N\\\{\\mathbf\{z\}^\{\(k\)\}\\\}\_\{k=1\}^\{N\},𝐳\(k\)∈ℝd\\mathbf\{z\}^\{\(k\)\}\\in\\mathbb\{R\}^\{d\}, be a sample ofNNtuples from a targetdd\-variate non\-parametric distributionp​\(z1,…,zd\)p\(z\_\{1\},\\dots,z\_\{d\}\)\(the data, or a representative training set\), to be placed on a grid ofn=Nn=Nlocations so that the realization reproduces bothppand a prescribed spatial covariance\. Writing the target tuples and a spatial template\{𝐠j\}\\\{\\mathbf\{g\}\_\{j\}\\\}as discrete measures, a realization is a coupling that assigns one tuple to each location; MST\-Direct selects it by the entropy\-regularized optimal\-transport problem\[[3](https://arxiv.org/html/2606.07578#bib.bib21)\]

minM∈Π⁡⟨C,M⟩−1β​H​\(M\),\\min\_\{M\\in\\Pi\}\\;\\langle C,M\\rangle\-\\tfrac\{1\}\{\\beta\}H\(M\),\(1\)over the transport polytopeΠ=\{M≥0:M​𝟏=1n​𝟏,M⊤​𝟏=1n​𝟏\}\\Pi=\\\{M\\\!\\geq\\\!0:\\,M\\mathbf\{1\}=\\tfrac\{1\}\{n\}\\mathbf\{1\},\\,M^\{\\\!\\top\}\\mathbf\{1\}=\\tfrac\{1\}\{n\}\\mathbf\{1\}\\\}, withCi​j=∥𝐠~i−𝐳~\(j\)∥2C\_\{ij\}=\\lVert\\tilde\{\\mathbf\{g\}\}\_\{i\}\-\\tilde\{\\mathbf\{z\}\}^\{\(j\)\}\\rVert^\{2\}the squared distance in the per\-variable standardized space,HHthe entropy, andβ\\betathe regularization\. The solutionM=diag​\(𝐫\)​K​diag​\(𝐜\)M=\\mathrm\{diag\}\(\\mathbf\{r\}\)K\\,\\mathrm\{diag\}\(\\mathbf\{c\}\),K=exp⁡\(−β​C\)K=\\exp\(\-\\beta C\), is found by the Sinkhorn fixed point\[[11](https://arxiv.org/html/2606.07578#bib.bib20)\], run in the log domain\.

A relationalkk\-nearest\-neighbour reward then augments the cost so that spatially adjacent locations are matched to mutually similar tuples\[[10](https://arxiv.org/html/2606.07578#bib.bib2)\]; withAAthe row\-normalized grid adjacency,

Mi​j∝exp⁡\(β​\(−Ci​j\+λ​\[A​M​A⊤\]i​j\)\),M\_\{ij\}\\propto\\exp\\\!\\big\(\\beta\(\-C\_\{ij\}\+\\lambda\\,\[AMA^\{\\\!\\top\}\]\_\{ij\}\)\\big\),\(2\)solved by alternating the rewardA​M​A⊤AMA^\{\\\!\\top\}with Sinkhorn normalization\. The soft coupling is rounded to a permutationπ\\piby greedy assignment in order of decreasing confidence\. Sinceπ\\piis a bijection onto the target tuples, the realization𝐳i∗=𝐳\(π​\(i\)\)\\mathbf\{z\}^\{\*\}\_\{i\}=\\mathbf\{z\}^\{\(\\pi\(i\)\)\}is the*same multiset*of tuples: every marginal and every non\-linear cross\-dependence is reproduced exactly\. OT decides only*where*each tuple goes\.

The template carries the prescribed spatial structure\. As in DMS we generatedd*independent*standard\-Gaussian fields with the target \(spherical\) variogram by FFT\-MA\[[6](https://arxiv.org/html/2606.07578#bib.bib9)\],𝐠\(c\)=ℱ−1​\{S​\(𝝎\)​ℱ​\{W\(c\)\}\}\\mathbf\{g\}^\{\(c\)\}=\\mathcal\{F\}^\{\-1\}\\\{\\sqrt\{S\(\\boldsymbol\{\\omega\}\)\}\\,\\mathcal\{F\}\\\{W^\{\(c\)\}\\\}\\\}\. Neighbouring locations receive similar feature vectors and, after matching, similar tuples, transferring the variogram to every variable\. The matching itself imposes*no*isotropy or stationarity requirement: any covariance admissible by the chosen Gaussian simulator — including anisotropic models — can be used in the backbone; we adopt an isotropic spherical model only for comparability with DMS\[[4](https://arxiv.org/html/2606.07578#bib.bib1)\]\.

A dense coupling isO​\(n2\)O\(n^\{2\}\)in memory and∼O​\(n3\)\\sim O\(n^\{3\}\)per relational iteration, infeasible atn=40 000n=40\\,000\[[10](https://arxiv.org/html/2606.07578#bib.bib2)\]\. We instead restrict each location to itsCCnearest target tuples \(via akk\-d tree\), so the coupling hasO​\(n​C\)O\(nC\)nonzeros and the Sinkhorn updates become segment reductions over the candidate edge list\. Greedy rounding yields the permutation; the few locations left without a free candidate are repaired by nearest\-unused assignment, guaranteeing a complete bijection\. The relational term \([2](https://arxiv.org/html/2606.07578#S2.E2)\) is applied as a few refinement passes pulling each query towardZ¯𝒩​\(i\)\\overline\{Z\}\_\{\\mathcal\{N\}\(i\)\}, the mean tuple currently assigned to the spatialkk\-nearest neighbours𝒩​\(i\)\\mathcal\{N\}\(i\)of locationii\(Algorithm 1\)\.

Algorithm 1: Scalable MST\-Direct

1. 1\.inputtarget tuplesZ∈ℝn×dZ\\in\\mathbb\{R\}^\{n\\times d\}; coordinates; variogram; hard data \(optional\)\.
2. 2\.standardizeZZ; build backboneGG\(kriging\-conditioned if hard data\)\.
3. 3\.pin hard\-data locations to their tuples; exclude from the pool\.
4. 4\.candidates:CCnearest rows ofZZto eachGiG\_\{i\}\(kk\-d tree\)\.
5. 5\.sparse log\-domain Sinkhorn \([1](https://arxiv.org/html/2606.07578#S2.E1)\) on the candidate edges\.
6. 6\.rrrelational passes \([2](https://arxiv.org/html/2606.07578#S2.E2)\):Gi←\(1−λ\)​Gi\+λ​Z¯𝒩​\(i\)G\_\{i\}\\leftarrow\(1\-\\lambda\)G\_\{i\}\+\\lambda\\,\\overline\{Z\}\_\{\\mathcal\{N\}\(i\)\}, rematch\.
7. 7\.greedy\-round to a permutationπ\\pi; repair leftovers\.
8. 8\.return𝐳i∗=Zπ​\(i\)\\mathbf\{z\}^\{\*\}\_\{i\}=Z\_\{\\pi\(i\)\}\(a permutation ofZZ\)\.

Hard data\{𝐳ℓd\}\\\{\\mathbf\{z\}^\{\\mathrm\{d\}\}\_\{\\ell\}\\\}at locations\{ℓ\}\\\{\\ell\\\}are honoured in two steps\. First, the backbone is conditioned to the data \(in the matching space\) by simple kriging\[[6](https://arxiv.org/html/2606.07578#bib.bib9)\], so it interpolates the data and keeps the prescribed covariance elsewhere\. Second, the data locations are*pinned*: their tuples are fixed and removed from the pool, and only the remaining locations enter \([1](https://arxiv.org/html/2606.07578#S2.E1)\)\. The realization thus honours the data exactly while the kriged backbone keeps the surroundings coherent\.

## IIIExperimental Design

We reproduce the two experiments of DMS\[[4](https://arxiv.org/html/2606.07578#bib.bib1)\]using*the same*6\-variate target distribution, which is strongly non\-linear and heteroscedastic \(left panels of Figs\.[1](https://arxiv.org/html/2606.07578#S4.F1)and[3](https://arxiv.org/html/2606.07578#S4.F3)\)\. In every comparison the target \(*real*\) distribution is shown with the PPMT and the MST\-Direct realizations, in that order; diagonals are marginals and off\-diagonals bivariate histograms\. We quantify the reproduction, for both methods, by the mean squared error \(MSE\) between experimental and reference histograms, and we assess the spatial reproduction by the experimental variograms\. We benchmark against PPMT\[[1](https://arxiv.org/html/2606.07578#bib.bib4)\]rather than the Gaussian\-copula or LU\-decomposition baselines used in the bivariate study\[[10](https://arxiv.org/html/2606.07578#bib.bib2)\]: those baselines are Gaussian by construction and therefore cannot reproduce complex non\-linear joints, which would make the comparison uninformative, whereas PPMT is a non\-parametric multivariate transform designed precisely for such distributions\. PPMT is implemented with projection\-pursuit Gaussianization and uses the same FFT\-MA backbone\. The unconditional grid is200×200200\\times 200with an isotropic spherical variogram of range 40; the conditional grid is100×100100\\times 100with 200 uniformly\-located hard data and a spherical variogram fitted from the data \(range≈16\\approx 16, vs\. 18 in\[[4](https://arxiv.org/html/2606.07578#bib.bib1)\]\)\.

## IVResults

For the unconditional experiment \(200×200200\\times 200grid, spherical range 40\), Fig\.[1](https://arxiv.org/html/2606.07578#S4.F1)compares the experimental distributions\. The MST\-Direct realization \(right\) is indistinguishable from the reference \(left\) — a spatial permutation of the same tuples — whereas PPMT \(centre\) recovers the gross structure but distorts the fine non\-linear joints\. Table[I](https://arxiv.org/html/2606.07578#S4.T1)reports the histogram MSE for both methods: identically zero for MST\-Direct and of order10−610^\{\-6\}–10−510^\{\-5\}for PPMT\. Fig\.[2](https://arxiv.org/html/2606.07578#S4.F2)shows that both methods reproduce the imposed range\-40 model; MST\-Direct matches the sill exactly \(its variance equals the reference by construction\), while PPMT overshoots it by1010–20%20\\%\. The scalable bijection leaves a mild short\-lag component, reduced by the relational passes\.

![Refer to caption](https://arxiv.org/html/2606.07578v1/figures/reference_scatter.png)

![Refer to caption](https://arxiv.org/html/2606.07578v1/figures/uncond_ppmt_scatter.png)

![Refer to caption](https://arxiv.org/html/2606.07578v1/figures/uncond_mst_scatter.png)

Figure 1:Unconditional experimental distributions:*real*\(left\), PPMT \(centre\), MST\-Direct \(right\)\. MST\-Direct is identical to the reference; PPMT distorts the fine non\-linear joints\.TABLE I:Unconditional histogram MSE vs\. the reference \(×10−5\\times 10^\{\-5\}\): MST\-Direct \(left\) and PPMT \(right\)\. Diagonal: marginals; off\-diagonal: bivariate histograms\.
![Refer to caption](https://arxiv.org/html/2606.07578v1/figures/uncond_variograms.png)Figure 2:Unconditional experimental variograms of MST\-Direct \(blue\) and PPMT \(red\) against the imposed spherical model with range 40 \(dashed\)\.For the conditional experiment, a reference model is built as an unconditional realization on the100×100100\\times 100grid, and 200 of its nodes become hard data\. Both methods honour the data — MST\-Direct exactly by pinning \(maximum absolute error0over the 200 locations and all six variables\), PPMT to machine precision \(∼10−14\\sim 10^\{\-14\}\)\. Fig\.[3](https://arxiv.org/html/2606.07578#S4.F3)compares the distributions \(again MST\-Direct exact, PPMT approximate\), with the MSE in Table[II](https://arxiv.org/html/2606.07578#S4.T2)\. Fig\.[4](https://arxiv.org/html/2606.07578#S4.F4)shows the three realizations as maps \(reference / PPMT / MST\-Direct\), all honouring the data; Fig\.[5](https://arxiv.org/html/2606.07578#S4.F5)verifies the exact data reproduction \(measured vs\. simulated on the45∘45^\{\\circ\}line\)\. Fig\.[6](https://arxiv.org/html/2606.07578#S4.F6)shows that both methods approximately reproduce the reference variograms, MST\-Direct tracking the sill exactly and PPMT slightly undershooting it — consistent with\[[4](https://arxiv.org/html/2606.07578#bib.bib1)\]\.

![Refer to caption](https://arxiv.org/html/2606.07578v1/figures/reference_scatter.png)

![Refer to caption](https://arxiv.org/html/2606.07578v1/figures/cond_ppmt_scatter.png)

![Refer to caption](https://arxiv.org/html/2606.07578v1/figures/cond_mst_scatter.png)

Figure 3:Conditional experimental distributions:*real*\(left\), PPMT \(centre\), MST\-Direct \(right\)\.![Refer to caption](https://arxiv.org/html/2606.07578v1/figures/cond_maps_comparison.png)Figure 4:Conditional realizations on the100×100100\\times 100grid, rows*reference / PPMT / MST\-Direct*, columnsz1z\_\{1\}–z6z\_\{6\}\(shared colour scale per column\)\. Black dots mark the 200 hard\-data locations, honoured by all methods\.![Refer to caption](https://arxiv.org/html/2606.07578v1/figures/cond_data_honoring.png)Figure 5:Conditional MST\-Direct: measured vs\. simulated values at the 200 hard\-data locations\. All points lie on the45∘45^\{\\circ\}line \(maximum absolute error0\)\.![Refer to caption](https://arxiv.org/html/2606.07578v1/figures/cond_variograms.png)Figure 6:Conditional experimental variograms: reference model \(black\), MST\-Direct \(blue\), PPMT \(red\)\.TABLE II:Conditional histogram MSE vs\. the reference \(×10−5\\times 10^\{\-5\}\): MST\-Direct \(left\) and PPMT \(right\)\.

## VDiscussion

The experiments close the three limitations left open in\[[10](https://arxiv.org/html/2606.07578#bib.bib2)\]\.*Scalability:*the sparse candidate\-restricted Sinkhorn runs the40 00040\\,000\- and10 00010\\,000\-node grids in well under a minute, where the dense matcher would require a40 000240\\,000^\{2\}coupling\.*Many variables:*matching the target cloud onto an independent FFT\-MA backbone extends the method to six variables \(and, in principle, more\) without the factorization error that grows with dimensionality in SCT/PPMT/DMS\.*Conditioning:*pinning combined with a kriged backbone honours hard data exactly while reproducing the spatial model\. As in DMS, neither method controls the cross\-variograms directly, yet both reproduce the target statistics\. The main residual artifact is a short\-lag component in the unconditional variograms, controllable through the number of candidates and relational passes; anisotropic variograms remain future work\.

## VIConclusion

We extended MST\-Direct to multivariate, conditional, large\-grid geostatistical simulation\. The method matches the target value tuples onto an FFT\-MA Gaussian backbone through a scalable Sinkhorn optimal\-transport step and honours hard data by pinning, so the multivariate joint distribution is preserved exactly while the prescribed spatial correlation is reproduced\. On the 6\-variate, heteroscedastic, non\-linear benchmark of Direct Multivariate Simulation, in both unconditional and conditional settings, MST\-Direct reproduces the target distribution with zero histogram error and honours the conditioning data exactly, outperforming PPMT on the joint distribution while matching it on the spatial statistics\. The approach is well suited to applications with complex, strongly non\-linear multivariate dependencies\.

## Code availability

An open\-source implementation \(v2\.0\.0\), together with the scripts that reproduce the experiments in this paper, is released as themst\-directPython package\[[9](https://arxiv.org/html/2606.07578#bib.bib3)\]\.

## Acknowledgment

The author thanks PX\.Center for supporting this research\.

## References

- \[1\]R\. M\. Barnett, J\. G\. Manchuk, and C\. V\. Deutsch\(2014\)Projection pursuit multivariate transform\.Mathematical Geosciences46,pp\. 337–359\.External Links:[Document](https://dx.doi.org/10.1007/s11004-013-9497-7)Cited by:[§I](https://arxiv.org/html/2606.07578#S1.p1.1),[§III](https://arxiv.org/html/2606.07578#S3.p1.3)\.
- \[2\]J\. Chilès and P\. Delfiner\(1999\)Geostatistics: modeling spatial uncertainty\.John Wiley & Sons\.Cited by:[§I](https://arxiv.org/html/2606.07578#S1.p1.1)\.
- \[3\]M\. Cuturi\(2013\)Sinkhorn distances: lightspeed computation of optimal transport\.InAdvances in Neural Information Processing Systems,Vol\.26,pp\. 2292–2300\.Cited by:[§I](https://arxiv.org/html/2606.07578#S1.p2.4),[§II](https://arxiv.org/html/2606.07578#S2.p1.8)\.
- \[4\]L\. P\. de Figueiredo, T\. Schmitz, R\. Lunelli, M\. Roisenberg, D\. S\. de Freitas, and D\. Grana\(2021\)Direct multivariate simulation – a stepwise conditional transformation for multivariate geostatistical simulation\.Computers & Geosciences147,pp\. 104659\.External Links:[Document](https://dx.doi.org/10.1016/j.cageo.2020.104659)Cited by:[item 4](https://arxiv.org/html/2606.07578#S1.I1.i4.p1.1),[§I](https://arxiv.org/html/2606.07578#S1.p1.1),[§II](https://arxiv.org/html/2606.07578#S2.p3.2),[§III](https://arxiv.org/html/2606.07578#S3.p1.3),[§IV](https://arxiv.org/html/2606.07578#S4.p2.4)\.
- \[5\]P\. Goovaerts\(1997\)Geostatistics for natural resources evaluation\.Oxford University Press\.Cited by:[§I](https://arxiv.org/html/2606.07578#S1.p1.1)\.
- \[6\]M\. Le Ravalec, B\. Noetinger, and L\. Y\. Hu\(2000\)The FFT moving average \(FFT\-MA\) generator: an efficient numerical method for generating and conditioning gaussian simulations\.Mathematical Geology32\(6\),pp\. 701–723\.External Links:[Document](https://dx.doi.org/10.1023/A%3A1007542406333)Cited by:[item 2](https://arxiv.org/html/2606.07578#S1.I1.i2.p1.1),[§II](https://arxiv.org/html/2606.07578#S2.p3.2),[§II](https://arxiv.org/html/2606.07578#S2.p8.2)\.
- \[7\]O\. Leuangthong and C\. V\. Deutsch\(2003\)Stepwise conditional transformation for simulation of multiple variables\.Mathematical Geology35\(2\),pp\. 155–173\.External Links:[Document](https://dx.doi.org/10.1023/A%3A1023235505120)Cited by:[§I](https://arxiv.org/html/2606.07578#S1.p1.1)\.
- \[8\]G\. Peyré and M\. Cuturi\(2019\)Computational optimal transport\.Foundations and Trends in Machine Learning11\(5–6\),pp\. 355–607\.External Links:[Document](https://dx.doi.org/10.1561/2200000073)Cited by:[§I](https://arxiv.org/html/2606.07578#S1.p2.4)\.
- \[9\]T\. B\. Schmitz\(2026\)mst\-direct: matching via sinkhorn transport for multivariate geostatistical simulation \(python package\)\.Note:Version 2\.0\.0[https://pypi\.org/project/mst\-direct/](https://pypi.org/project/mst-direct/)Cited by:[Code availability](https://arxiv.org/html/2606.07578#Sx1.p1.1)\.
- \[10\]T\. B\. Schmitz\(2026\)MST\-Direct: matching via sinkhorn transport for multivariate geostatistical simulation with complex non\-linear dependencies\.arXiv preprint arXiv:2603\.18036\.External Links:[Document](https://dx.doi.org/10.48550/arXiv.2603.18036)Cited by:[§I](https://arxiv.org/html/2606.07578#S1.p2.4),[§II](https://arxiv.org/html/2606.07578#S2.p1.8),[§II](https://arxiv.org/html/2606.07578#S2.p2.2),[§II](https://arxiv.org/html/2606.07578#S2.p4.10),[§III](https://arxiv.org/html/2606.07578#S3.p1.3),[§V](https://arxiv.org/html/2606.07578#S5.p1.3)\.
- \[11\]R\. Sinkhorn\(1964\)A relationship between arbitrary positive matrices and doubly stochastic matrices\.The Annals of Mathematical Statistics35\(2\),pp\. 876–879\.External Links:[Document](https://dx.doi.org/10.1214/aoms/1177703591)Cited by:[§I](https://arxiv.org/html/2606.07578#S1.p2.4),[§II](https://arxiv.org/html/2606.07578#S2.p1.14)\.
- \[12\]C\. Villani\(2009\)Optimal transport: old and new\.Grundlehren der mathematischen Wissenschaften, Vol\.338,Springer\.External Links:[Document](https://dx.doi.org/10.1007/978-3-540-71050-9)Cited by:[§I](https://arxiv.org/html/2606.07578#S1.p2.4)\.