Particle method for a nonlinear multimarginal optimal transport problem

1. Problem Statement

The paper addresses a nonlinear multimarginal optimal transport (OT) problem arising in risk management. The core objective is to find a joint probability distribution (coupling) $\gamma$ of $D$ random variables (risk factors) with fixed marginal distributions $\rho_1, \dots, \rho_D$ that maximizes a specific risk measure.

• The Objective: Maximize a Spectral Risk Measure of the pushforward of the coupling by a cost function $c$.
  $$\max_{\gamma \in \Gamma(\rho_1, \dots, \rho_D)} R_\alpha(c_\# \gamma)$$
  where $R_\alpha(\mu) = \int_0^1 F_\mu^{-1}(t) \alpha(t) dt$. Here, $\alpha$ is a non-decreasing, non-negative weighting function (spectral function), and $F_\mu^{-1}$ is the quantile function.
• Context: Unlike standard OT where the objective is linear ($\int c d\gamma$), the spectral risk measure makes the objective nonlinear because it biases the average towards high values (tail risk).
• Linear Reformulation: The authors leverage a known result (Ennaji et al., [6]) showing that this nonlinear problem is equivalent to a linear multimarginal OT problem with an additional marginal $\rho_0 = \alpha_\# \mathcal{L}(0,1)$. The new problem involves finding a coupling $\eta$ between $\rho_0$ and the original marginals to maximize $\int \omega c(x) d\eta$.

2. Methodology: Lagrangian Particle Discretization

The authors propose a numerical scheme based on Lagrangian particle discretization, avoiding the inefficiency of Eulerian grid-based methods for high-dimensional OT problems.

• Discretization: The admissible couplings are approximated by uniformly weighted point clouds (Dirac masses) of $N$ points:
  $$\delta_Y = \frac{1}{N} \sum_{i=1}^N \delta_{y_i}, \quad Y = (y_1, \dots, y_N) \in \mathcal{X}^N$$
• Relaxation of Constraints: Instead of enforcing the marginal constraints $\pi_j \# \gamma = \rho_j$ strictly (which is difficult for discrete measures), the authors use quadratic Wasserstein penalization.
• The Discrete Problem ($P_N$): The problem is reformulated as a finite-dimensional minimization:
  $$\min_{Y \in \mathcal{X}^N} \left( -R_\alpha(c_\# \delta_Y) + \lambda_N \sum_{j=1}^D W_p(\rho_j, \pi_j \# \delta_Y)^p \right)$$
  • The first term is the negative spectral risk (to be minimized).
  • The second term penalizes the deviation of the projected point clouds from the target marginals $\rho_j$.
  • $\lambda_N$ is a penalty parameter that depends on $N$.
• Computation:
  • The spectral risk term is computed efficiently by sorting the values $c(y_i)$ and applying weights derived from $\alpha$.
  • The penalty terms involve solving semi-discrete optimal transport problems (or partial OT for CVaR cases), computed using Laguerre tessellations and gradient-based optimization (Limited-memory BFGS).

3. Key Contributions

1. Quantitative Convergence Analysis: The paper provides rigorous convergence rates for the discretized problem ($P_N$) to the continuous problem ($P$) as $N \to \infty$.
2. Link to Quantization Error: The convergence rate is shown to be governed by the uniform quantization error of the optimal solution.
3. Dimensional Dependence: The authors derive explicit bounds based on the box dimension (Minkowski dimension) of the support of the optimal solution.
4. Supermodular Case: For univariate marginals and supermodular cost functions, the authors prove that the optimal coupling is comonotone. This allows for sharper convergence rates dependent directly on the quantization errors of the marginals rather than the joint support.
5. Partial Transport Extension: The method is extended to the Conditional Value at Risk (CVaR) case, which reduces the problem to a multimarginal partial transport formulation (transporting only a fraction $m$ of the mass).

4. Main Results

A. General Convergence Theorem (Theorem 1 & Corollary 17)

Under assumptions of Hölder continuity for the cost function $c$ (exponent $\beta$) and integrability for the spectral function $\alpha$:

• Let $d$ be the box dimension of the support of an optimal solution.
• The error in the objective value scales as:
  $$| \min F_N - \min F | = O(u_N^\beta)$$
  where $u_N$ is the optimal $N$-point uniform quantization error.
• Specific Rates:
  • If $d > p$: Rate is $O(N^{-\beta/d})$.
  • If $d = p$: Rate is $O((\log N)^{\beta/d} N^{-\beta/d})$.
  • If $d < p$: Rate is $O(N^{-\beta/p})$.
• This implies that if the optimal solution lies on a low-dimensional manifold (common in OT), the method converges significantly faster than the curse of dimensionality would suggest.

B. Supermodular Case (Theorem 2 & Corollary 24)

If marginals are univariate and $c$ is supermodular and increasing:

• The optimal plan is the comonotone coupling (totally ordered support).
• The convergence rate depends on the quantization errors of the marginals ($\rho_j$) rather than the joint support.
• If marginals have connected support and densities bounded away from zero, the rate improves to $O(N^{-\beta})$ (effectively $d=1$), which is optimal for univariate problems.

C. Partial Transport (CVaR)

For the CVaR case (mass $m < 1$):

• The problem is formulated using the partial Wasserstein distance $W_{p, \max}$.
• Similar convergence rates are established, linking the error to the quantization of the active submeasures of the marginals.

5. Numerical Experiments & Significance

The authors validate their method across several domains:

• Risk Management (Flood Modeling): Applied to a hydrological model with 8 input variables. The method successfully identifies the "worst-case" dependence structure (coupling) that maximizes the flood overflow risk, demonstrating the ability to handle high-dimensional, non-linear risk aggregation.
• Partial Barycenters: The method is tested on partial Wasserstein barycenters. While it captures the structure, the authors note challenges with non-local minima in the partial transport setting, suggesting initialization strategies need refinement.
• Repulsive Costs: Simulations with Coulomb-like repulsive costs show the method can capture complex, low-dimensional structures (e.g., cyclical structures in density functional theory contexts).

Significance:

• Bridging Theory and Practice: The paper provides a practical algorithm for a class of problems (nonlinear multimarginal OT) that were previously difficult to solve numerically due to the lack of linear structure and high dimensionality.
• Efficiency: By utilizing particle methods and penalization, it bypasses the need for high-dimensional grids.
• Theoretical Insight: The connection between the convergence rate and the box dimension of the solution's support offers a theoretical justification for why particle methods work well even in high dimensions (since optimal OT plans often have low-dimensional support).
• Risk Management: It offers a robust tool for stress-testing financial or engineering systems by identifying the most dangerous dependence structures among risk factors, moving beyond the assumption of independence.