---
title: Optimal Transport Global Sensitivity Analysis
url: https://www.emergentmind.com/topics/optimal-transport-based-global-sensitivity-analysis
type: topic
---

# Optimal Transport Global Sensitivity Analysis

Searching arXiv for the specified papers and closely related work on optimal transport-based sensitivity analysis.
Optimal Transport-based Global Sensitivity Analysis denotes a family of sensitivity-analysis frameworks in which discrepancies between outputs generated under different inputs, parameter values, or intervention rules are quantified by optimal-transport functionals, most commonly Wasserstein distances. Across the literature, the approach appears in several distinct but related settings: sensitivity of PDE solutions such as the LWR traffic-flow model on networks [1608.00126]; global sensitivity analysis for computer codes whose outputs are probability measures or cumulative distribution functions, and for stochastic computer codes [2007.12378]; sensitivity analysis for Fokker–Planck equations through quantitative Wasserstein bounds [2602.03174]; sensitivity analysis for stochastic interventions in causal inference via Monge–Kantorovich formulations [2411.14285]; and model-agnostic post-processing tools for OT-based sensitivity indices implemented in software [2507.18588]. In all of these formulations, the central idea is to replace scalar dispersion summaries by a transport-based comparison of full output laws, thereby measuring “mass shift” rather than only pointwise or moment-wise error [1608.00126].

## 1. Conceptual scope and problem classes

Optimal Transport-based Global Sensitivity Analysis is not a single index but a collection of methodologies adapted to the geometry of distribution-valued, stochastic, or law-dependent outputs. In the formulation for black-box codes, one considers a code
$$
Z = f(X_1,\ldots,X_p)
$$
with independent inputs $X_1,\ldots,X_p$ and output $Z$ taking values in a metric space $(\mathcal X,d)$; a principal special case is $\mathcal X=\mathcal W_q(\mathbb R)$, the space of real probability measures with finite $q$th moment, endowed with the $q$-Wasserstein distance [2007.12378]. In this setting, sensitivity indices are defined directly on the space of output laws rather than on scalar summaries.

A second problem class concerns stochastic computer codes. There the code has the form
$$
f_s: E\times D\to\mathbb R,\quad (x,d)\mapsto y,
$$
where $d$ is an unobserved “random seed.” For each fixed $x$, the law of $f_s(x,\cdot)$ is a probability measure $\mu_x$, which induces the “ideal deterministic” code
$$
f: E\to\mathcal W_q(\mathbb R),\qquad x\mapsto \mu_x
$$
[2007.12378]. This construction embeds stochastic simulation outputs into Wasserstein space and permits a standard global-sensitivity analysis at the level of output distributions.

A third class is PDE- and transport-model sensitivity. For traffic networks, the LWR density on each directed edge satisfies
$$
\partial_t\rho_e(x,t)+\partial_x f(\rho_e(x,t))=0,\qquad \rho_e(x,0)=\rho_e^0(x),
$$
with junction dynamics encoded through distribution matrices $A^v=(\alpha^v_{r,r'})$ and sensitivity measured by Wasserstein distance between density fields generated by different inputs [1608.00126]. For diffusion-type models, the Fokker–Planck equation
$$
\partial_t \rho + \nabla\cdot(b(x,\theta)\rho)=\nabla\cdot(A(x,\theta)\nabla\rho)
$$
is analyzed with respect to parameter perturbations using explicit upper bounds on $W_p(\rho(t;\theta),\rho(t;\theta'))$ [2602.03174].

A fourth class arises in causal sensitivity analysis under stochastic or generalized interventions. With observed data $O=(X,A,Y)$, generalized policies induce couplings between observational and target treatment distributions, and worst-case sensitivity bounds can be written as Monge–Kantorovich problems over couplings $J\in\mathcal J(\Pi,Q)$ [2411.14285]. This suggests that OT-based sensitivity analysis is equally a framework for uncertainty propagation and for partial-identification geometry.

## 2. Mathematical formulations in Wasserstein space

The unifying object is an optimal-transport cost between probability measures. In the Kantorovich formulation, for measures $\mathbb P,\mathbb Q$ on $\mathcal Y$ and nonnegative cost $c:\mathcal Y\times\mathcal Y\to[0,\infty)$,
$$
K(\mathbb P,\mathbb Q)=\inf_{\pi\in\Pi(\mathbb P,\mathbb Q)} \int_{\mathcal Y\times\mathcal Y} c(y,y')\,d\pi(y,y').
$$
When $c(y,y')=\|y-y'\|^p$, $K$ is the $p$-Wasserstein cost [2507.18588]. In one dimension,
$$
W_q^q(F,G)=\int_0^1 |F^{-1}(u)-G^{-1}(u)|^q\,du
$$
for c.d.f.s $F,G$ [2007.12378].

For OT-based sensitivity indices in model-agnostic GSA, if $\mathbf Y$ is the model output and $X_i$ an input, a “separation” functional $\zeta$ satisfies $\zeta(\mathbb P,\mathbb P)=0$, and the local separation is
$$
\zeta(x_i)=\zeta(\mathbb P_{\mathbf Y},\mathbb P_{\mathbf Y\mid X_i=x_i}).
$$
The associated global index is
$$
\xi^\zeta(\mathbf Y,X_i)=\mathbb E_{X_i}[\zeta(X_i)].
$$
Using OT as the separation, one sets
$$
\zeta^K(x_i)=K(\mathbb P_{\mathbf Y},\mathbb P_{\mathbf Y\mid X_i=x_i}),
$$
with numerator
$$
\xi^K(\mathbf Y,X_i)=\mathbb E\!\left[K(\mathbb P_{\mathbf Y},\mathbb P_{\mathbf Y\mid X_i})\right].
$$
Normalization uses
$$
\mathbb M^K[\mathbf Y]=\mathbb E[c(\mathbf Y,\mathbf Y')],
$$
where $\mathbf Y'$ is an independent copy, and yields the first-order OT sensitivity index
$$
\iota^K(\mathbf Y,X_i)=\frac{\xi^K(\mathbf Y,X_i)}{\mathbb M^K[\mathbf Y]}\in[0,1]
$$
[2507.18588].

In the Wasserstein-space Sobol-style construction, if $\mathcal Y$ is a random element in $\mathcal W_p$ with Fréchet mean
$$
m=\arg\min_{\eta\in\mathcal W_p}\mathbb E[W_p^p(\mathcal Y,\eta)],
$$
then the Fréchet variance is
$$
V=\mathbb E[W_p^p(\mathcal Y,m)].
$$
For a subset $u\subseteq\{1,\ldots,p\}$, the conditional Fréchet mean $m_u(X_u)$ defines the sensitivity index
$$
S_u=1-\frac{\mathbb E[W_p^p(\mathcal Y,m_u)]}{V}
=\frac{\mathbb E[W_p^p(m_u,m)]}{V}.
$$
In one dimension for $p=2$,
$$
S_u=
\frac{\int_0^1 \mathrm{Var}[\mathbb E[\mathcal Y^{-1}(v)\mid X_u]]\,dv}
{\int_0^1 \mathrm{Var}[\mathcal Y^{-1}(v)]\,dv}
$$
[2007.12378].

For networked conservation laws, if $\rho^s,\rho^d$ are nonnegative densities on a network $\mathcal N$ with equal total mass $M$, then
$$
W_p(\rho^s,\rho^d)^p
=\inf_{\gamma\in\Gamma(\rho^s,\rho^d)}\int_{\mathcal N\times\mathcal N} D(x,y)^p\,d\gamma(x,y),
$$
where $D(x,y)$ is the network geodesic distance [1608.00126]. For treatment-policy sensitivity, the worst-case deviation from the identified component $t_Q(X)$ is
$$
\sup_{J\in\mathcal J_Q}\int \Gamma(a,a',X)\,dJ(a,a'|X),
$$
and minimizing worst-case width over $J$ yields the Monge–Kantorovich problem
$$
W_\Gamma(X)=\min_{J\in\mathcal J_Q(X)}\int c(a,a',X)\,dJ(a,a')
$$
with dual
$$
W_\Gamma(X)=\max_{f,g}\int f(a)\,d\Pi(a)+\int g(a')\,dQ(a')
\quad\text{subject to } f(a)+g(a')\le c(a,a',X)
$$
[2411.14285].

## 3. Principal index constructions and structural properties

Several OT-based indices have been proposed, reflecting different output objects and inferential goals. For distribution-valued outputs, the Wasserstein–Fréchet approach uses the Fréchet mean and Fréchet variance in $\mathcal W_p$ and yields Sobol-style indices $S_u$ with $0\le S_u\le 1$, $\sum_u S_u=1$, and invariance under translation, isometries or nondegenerate scaling of the output [2007.12378]. A second construction in the same work is based on Hoeffding decomposition of indicator functions of Wasserstein balls. For fixed $r\ge 0$ and target distributions $F,G$,
$$
I_r^{F,G}(Y)=\mathbf 1\{W_p(Y,F)\le r\}.
$$
Integrating the associated partial variances over reference pairs leads to
$$
S_{2,W_p}^u=
\frac{\int \mathrm{Var}\!\left[\mathbb E[T_{(F_1,F_2)}(Y)\mid X_u]\right]\,d\mathbb P(F_1)\,d\mathbb P(F_2)}
{\int \mathrm{Var}[T_{(F_1,F_2)}(Y)]\,d\mathbb P(F_1)\,d\mathbb P(F_2)},
$$
where
$$
T_{(F_1,F_2)}(Y)=\mathbf 1\{W_p(Y,F_1)\le W_p(F_1,F_2)\}
$$
[2007.12378].

The model-agnostic OT index implemented in software is based on conditional-vs-marginal distributional separation. Its normalization is designed so that $\iota^K=0\iff X_i\perp\!\!\!\perp \mathbf Y$ and $\iota^K=1\iff \mathbf Y=f(X_i)$ almost surely, with monotonicity and normalization properties [2507.18588]. This framework is especially adapted to multivariate outputs $\mathbf Y\in\mathbb R^k$ and does not require independence of the inputs $\{X_i\}$ [2507.18588].

When the cost is squared Euclidean, the index admits a decomposition using the Gelbrich formula:
$$
K(\mathbb P,\mathbb Q)=\|\mu-\mu'\|^2
+\mathrm{Tr}\Bigl(\Sigma+\Sigma'-2(\Sigma^{1/2}\Sigma'\Sigma^{1/2})^{1/2}\Bigr)
+\Gamma(\mathbb P,\mathbb Q),
$$
leading to
$$
\iota^K=\iota^V+\iota^\Sigma+\iota^\Gamma.
$$
The Wasserstein-Bures semi-metric index is $\iota^{WB}=\iota^V+\iota^\Sigma$ [2507.18588]. This decomposition connects OT-based sensitivity to mean effects, covariance effects, and higher-order moment effects.

In PDE settings, the role of the sensitivity index is often played by a normalized transport discrepancy rather than a variance ratio. For the LWR model, given a baseline input vector $\theta^s$ and perturbed input $\theta^d$, one computes
$$
\varepsilon_i(t)=\frac{1}{M}H(\rho^s(\cdot,t),\rho^d(\cdot,t)),
$$
where $H$ is the discrete transport optimum, and records either $S_i=\varepsilon_i(T)$ or the time-average
$$
\langle \varepsilon_i\rangle=\frac{1}{T}\int_0^T \varepsilon_i(t)\,dt
$$
[1608.00126]. For Fokker–Planck equations, Morange’s framework defines a transport-based sensitivity index
$$
S_p(t):=\frac{1}{\Delta_\Theta}\,\esssup_{\theta\in\Theta}W_p(\rho(t;\theta),\rho(t;\bar\theta))
$$
relative to a nominal $\bar\theta$, and proves $S_p(t)\le C_p(t)$ [2602.03174].

## 4. Computational procedures and algorithms

A core practical feature of OT-based sensitivity analysis is the reduction of sensitivity computation to optimal-transport solvers. For the LWR network model, each edge is discretized into uniform cells $C_{e,j}$ of length $\Delta x$, producing an undirected graph $H=(V_H,E_H)$ whose nodes correspond to cell centers. Masses are accumulated as
$$
s_i=\rho^s(x_i)\Delta x,\qquad d_i=\rho^d(x_i)\Delta x,
$$
with equal total mass $M$, and the cost matrix is built from shortest-path distances $c_{ij}=D_H(i,j)$, for example by Dijkstra. The resulting Hitchcock linear program is
$$
\min H=\sum_{i=1}^J\sum_{j=1}^J c_{ij}x_{ij}
$$
subject to
$$
\sum_{j=1}^J x_{ij}=s_i,\qquad
\sum_{i=1}^J x_{ij}=d_j,\qquad
x_{ij}\ge 0.
$$
Then $W_1(\rho^s,\rho^d)\approx H$, with error bound $|W_1-H|\le M\Delta x$ [1608.00126].

For generic input-output samples, the implemented OT-GSA workflow partitions the domain of an input $X_i$ into $H$ bins and, for each bin, computes an empirical OT cost between the marginal output sample and the conditional output sample restricted to that bin. The numerator is averaged over bins, and the result is normalized by an empirical upper bound
$$
\hat{\mathbb M}^K=\frac{1}{N(N-1)}\sum_{n\ne m} c(y^{(n)},y^{(m)})
$$
to obtain $\hat\iota^K$ [2507.18588]. The special cases of one-dimensional outputs and Gaussian/Bures geometry admit closed forms handled by `ot_indices_1d` and `ot_indices_wb` [2507.18588].

Two solvers are emphasized in software. Network Simplex solves the classical OT linear program with worst-case complexity roughly $O((N+N_h)^3)$, while Sinkhorn–Knopp solves the entropic-regularized problem
$$
K_\varepsilon(\mathbb P,\mathbb Q)
=\min_{\pi\in\Pi(\mathbb P,\mathbb Q)} \int c\,d\pi
+\varepsilon\,\mathrm{KL}(\pi\|\mathbb P\otimes\mathbb Q),
$$
with per-iteration cost $O(N^2)$ and linear convergence in the sense that $\|\pi^{(t)}-\pi^*\|\to 0$ exponentially fast; the “sinkhorn-stable” variant avoids underflow/overflow for small $\varepsilon$, and $K_\varepsilon\to K$ as $\varepsilon\to 0$ [2507.18588].

For stochastic computer codes in Wasserstein spaces, practical estimation proceeds by generating $N$ outer samples of the inputs and, for each outer sample, $n$ inner runs corresponding to different random seeds to form empirical measures
$$
\mu_{x,n}=\frac1n\sum_{k=1}^n \delta_{f_s(x,D_k)}.
$$
These empirical measures are then plugged into pick-freeze estimators, U-statistic estimators, or rank-based formulas for $S_u$ or $S_{2,W_p}^u$ [2007.12378].

For Fokker–Planck sensitivity, implementation is analytic rather than solver-centered. One first estimates or bounds $L_b$, $L_A$, and the uniform ellipticity $m$, then assembles
$$
C_{1,d,p}=(L_b+L_A)(2p-1)+\frac{L_A(p-1)^2}{2m},
\qquad
C_{2,d,p}=L_b+L_A+\frac{(p-1)L_A}{2m},
$$
and uses the resulting bound to control output-law perturbations over a parameter-uncertainty radius [2602.03174].

## 5. Quantitative bounds, analytical results, and canonical couplings

The Fokker–Planck setting provides an explicit global-in-time Lipschitz-in-parameter transport estimate. Under the uniform-ellipticity and global-Lipschitz assumptions, for all $t\ge 0$,
$$
W_p(\rho(t;\theta),\rho(t;\theta'))
\le
W_p(\rho_0(\cdot;\theta),\rho_0(\cdot;\theta'))e^{C_{1,d,p}t}
+
C_{1,d,p}C_{2,d,p}(e^{C_{1,d,p}t}-1)\|\theta-\theta'\|^p.
$$
If the initial law is unchanged, this becomes
$$
W_p(\rho(t;\theta),\rho(t;\theta'))
\le
C_{2,d,p}^{1/p}(e^{C_{1,d,p}t}-1)^{1/p}\|\theta-\theta'\|
$$
[2602.03174]. Two proof strategies yield the same constants: synchronous coupling of the underlying SDEs and differentiation of the Kantorovich dual [2602.03174].

In the overdamped Langevin case,
$$
\partial_t\rho=\beta^{-1}\Delta\rho+\nabla\cdot(\rho\nabla V(x,\theta)),
$$
with $\nabla^2V(x,\theta)\ge k\,\mathrm{Id}$ and $|\nabla V(x,\theta)-\nabla V(x,\theta')|\le L_3|\theta-\theta'|$, one has the sharper estimate
$$
W_p(\rho(t;\theta,\beta),\rho(t;\theta',\beta'))
\le
W_p(\rho_0(\theta,\beta),\rho_0(\theta',\beta'))e^{-\lambda t}
+K_{1,d,p}|\theta-\theta'|^p
+K_{2,d,p}|\beta^{-1}-\beta'^{-1}|^p(1-e^{-\lambda t}),
$$
where $\lambda=pk$ [2602.03174]. The separation of contraction and parameter-shift effects is a distinctive feature of this result.

In causal sensitivity analysis, the sharp worst-case bounds depend on the cost structure. Under Model 2, where
$$
c(a,a',X)=\Gamma\mathbf 1\{a\ne a'\},
$$
the optimal cost is
$$
W_\Gamma(X)=\Gamma\cdot \mathrm{TV}(\Pi,Q),
$$
so that
$$
E[Y(d)\mid X]\in
[t_Q(X)-\Gamma\,\mathrm{TV}(\Pi,Q),\,
 t_Q(X)+\Gamma\,\mathrm{TV}(\Pi,Q)].
$$
As $Q\to \Pi$, the bounds collapse to $t_Q(X)$ [2411.14285]. Under Model 3, where $c(a,a',X)=\Gamma|a-a'|^p$, the minimal cost is the $p$-Wasserstein distance,
$$
W_p^p(\Pi,Q)=\min_{J\in\mathcal J_Q}\int |a-a'|^p\,dJ(a,a'),
$$
and
$$
E[Y(d)\mid X]\in
[t_Q(X)-\Gamma W_p^p(\Pi,Q),\,
 t_Q(X)+\Gamma W_p^p(\Pi,Q)]
$$
[2411.14285].

The optimal couplings in these causal models are analytically characterized. For Model 2, the optimal coupling is the maximal coupling $J^*$ with
$$
J^*[a=a']=1-\mathrm{TV}(\Pi,Q),
$$
implemented by a two-stage policy that keeps $\bar A=A$ with probability $p(a|X)=\min\{1,q(a|X)/\pi(a|X)\}$ and otherwise redraws $\bar A\sim Q(\cdot|X)$ [2411.14285]. For Model 3 with continuous $\Pi$ and $Q$, the optimal transport for convex $h(|a-a'|)$ is the increasing rearrangement
$$
\bar a=d_Q^{mon}(X,A)=Q^{-1}(\Pi(A|X)|X)
$$
[2411.14285]. These results show that OT-based sensitivity analysis can identify not only bound widths but also extremal generalized policies.

## 6. Applications, numerical findings, and limitations

The LWR traffic-flow study is a canonical application of OT-based sensitivity analysis to large networks. The baseline input vector $\theta^s$ collects the initial condition $\rho^0$, fundamental-diagram parameters $(\sigma,f_{\max})$, the distribution matrices $\{A^v\}$, and the network topology $\mathcal N$; perturbations are introduced factor by factor, and the transport discrepancy $\varepsilon_i(t)$ is evaluated over time [1608.00126]. The reported numerical findings are specific: initial-condition errors produced by small-scale shifts yield a $W_1$ that decays rapidly in time, so $S_{IC}$ is small for large $T$; sensitivity to $\sigma$ and $f_{\max}$ grows roughly linearly with $|\Delta\sigma|$ and $|\Delta f_{\max}|$, with a proportionality constant increasing with network size; a small bias in a junction distribution coefficient causes significant $W_1$ growth in time and saturation; perturbing all junctions amplifies sensitivity roughly proportionally to network diameter; and topology changes such as closing one central road produce sensitivity of order one, with long-time $W_1$ scaling linearly with network size [1608.00126]. The paper further states that the LP-based computation is robust but scales as $O(J^3)$ in the worst case, so moderate grid resolution and network size are recommended in practice [1608.00126].

For distribution-valued and stochastic computer codes, numerical studies illustrate both first-level and second-level GSA. The toy cdf model permits closed-form expressions for $S_u$ and $S_{2,W_2}^u$, while empirical comparisons show that rank-based estimators often have substantially smaller MSE than pick-freeze estimators at equal code-call budgets [2007.12378]. In the stochastic-code version of the same toy model, empirical output measures are built from i.i.d. draws and the rank-based method again outperforms pick-freeze for first-order indices [2007.12378]. Second-level GSA is formulated by treating the distributions of the inputs themselves as uncertain and computing OT-based indices with respect to input-law parameters [2007.12378].

The software package `gsaot` operationalizes the model-agnostic OT-GSA framework as a post-processing step requiring only an input matrix `x` and an output matrix `y`; it supports `ot_indices`, `ot_indices_wb`, `ot_indices_1d`, `ot_indices_smap`, `irrelevance_threshold`, and associated plotting methods [2507.18588]. Practical examples include a linear Gaussian test, a spruce budworm ODE system, and a climate module with a custom $L^3$ cost matrix [2507.18588]. The package emphasizes multivariate outputs, correlated inputs, and visualization of local separation curves [2507.18588].

Several limitations recur across the literature. In the Fokker–Planck theory, the current results cover $p\ge 2$; the $W_1$ case might be handled via reflection coupling, while the $W_\infty$ case remains open [2602.03174]. The drift and diffusion must be globally Lipschitz, and $A(x,\theta)$ uniformly elliptic; singular or degenerate diffusions require new ideas [2602.03174]. In generic OT solvers, worst-case complexity is cubic, entropic regularization introduces bias for finite $\varepsilon$, and binning smooths local structure [2507.18588]. Rank-based methods for Wasserstein-space GSA are not yet available for higher-order indices or vector-valued inputs, and U-statistic kernels can become high-dimensional [2007.12378].

A common misconception is that transport-based methods merely repackage variance-based sensitivity. The cited works indicate otherwise. OT-based indices are formulated on full conditional and marginal laws, can be applied to multivariate outputs and correlated inputs, and in the causal setting yield sharp nonparametric bound widths that vanish as the target policy approaches the observational regime [2411.14285; 2507.18588]. A plausible implication is that optimal transport serves not only as a discrepancy metric but as a structural language for global sensitivity whenever the output object is itself distributional, stochastic, or intervention-dependent.

Source: https://www.emergentmind.com/topics/optimal-transport-based-global-sensitivity-analysis