Measure-Preserving Martingale Sinkhorn Algorithm

• The MPMS algorithm is an iterative entropic method that generalizes the classical Sinkhorn approach by incorporating martingale (conditional mean preservation) constraints.
• It employs iterative Bregman projections with dual variable updates to enforce row-mass, column-mass, and martingale fits, ensuring exponential convergence under mild conditions.
• MPMS is applied in stochastic volatility model calibration and quantitative finance, delivering robust, arbitrage-free solutions with high precision.

The Measure-Preserving Martingale Sinkhorn (MPMS) algorithm is an entropic iterative method for solving martingale optimal transport (MOT) and related constrained transport problems. MPMS generalizes the classical Sinkhorn algorithm—widely used for entropic optimal transport—by incorporating the martingale constraint, which enforces conditional mean preservation at every source location. The algorithm thus yields plans that simultaneously match prescribed marginals and preserve the martingale property, with guarantees of exponential convergence under broad conditions. MPMS has become a central computational tool in several domains, including stochastic volatility model calibration and robust finance.

1. The Entropic Martingale Optimal Transport Problem

The primary setting is the entropic martingale optimal transport (EMOT) problem. Given probability measures $\mu,\nu \in \mathcal P(\mathbb{R})$ with finite first moment and cost function $c(x,y,z)$, the objective is to find a martingale coupling $\pi$ on $\mathbb{R}^3$ (with prescribed marginals and mean-preserving constraint) that minimizes the relative entropy to a reference law $Q$ with density

$$\frac{dQ}{d(\mu \otimes \nu \otimes \rho)} = \exp(-c(x,y,z)).$$

The admissible set is

$$\mathcal{M}(\mu,\nu) = \left\{ \pi \in \mathcal{P}(\mathbb{R}^3) : \pi_x = \mu,\, \pi_y = \nu,\, \mathbb{E}^{\pi}[Y-X \mid X] = 0 \;\text{a.s.} \right\}.$$

The EMOT problem is

$$\inf_{\pi \in \mathcal{M}(\mu,\nu)} H(\pi \mid Q) = \inf_{\pi \in \mathcal{M}(\mu,\nu)} \int \log \left( \frac{d\pi}{dQ} \right) \, d\pi,$$

equivalently,

$$\inf_{\pi \in \mathcal{M}(\mu,\nu)} H(\pi\mid \mu\otimes\nu\otimes\rho) + \int c\,d\pi,$$

where $H(\pi\mid \cdot)$ denotes relative entropy.

Entropic regularization via $c(x,y,z)$0 induces strict convexity and smoothness, rendering the problem well-posed and numerically tractable. In the zero-temperature limit, the formulation recovers the classical (non-entropic) MOT problem (Chen et al., 2024).

2. Duality and Variational Structure

The dual formulation introduces potentials $c(x,y,z)$1, $c(x,y,z)$2, and $c(x,y,z)$3, yielding the Lagrangian

$c(x,y,z)$4

The dual problem is

$c(x,y,z)$5

with

$c(x,y,z)$6

where

$c(x,y,z)$7

Under compactness and Lipschitz continuity of the cost, there is strong duality and no primal-dual gap (Chen et al., 2024, Eckstein et al., 27 Feb 2025).

3. Algorithmic Structure: Iterative Bregman Projections

MPMS alternates Bregman (KL-divergence) projections onto affine constraint sets, exactly generalizing the Sinkhorn iterative proportional fitting procedure (IPFP). At each iteration, updates enforce the three constraint families:

1. Row-mass fit: Enforce $c(x,y,z)$8 (source marginal).
2. Martingale fit: Enforce $c(x,y,z)$9 (conditional mean preservation).
3. Column-mass fit: Enforce $\pi$0 (target marginal).

The steps are:

• Update exponentials or dual variables (row/col scaling factors and martingale multipliers).
• For the martingale projection, solve small-dimensional nonlinear equations per row (typically by Newton or bisection).
• Repeat until convergence in all constraints.

In the most general setting (with cost $\pi$1), one has (Langbridge et al., 25 Jun 2026):

$\pi$2

where $\pi$3 and $\pi$4 are scaling vectors, and $\pi$5 handles the martingale constraint.

The discrete MPMS pseudocode aligns with the above structure (see (Chen et al., 2024, Langbridge et al., 25 Jun 2026, Tang et al., 25 Aug 2025)):

Step	Constraint	Update Structure
Row-mass	$\pi$6	$\pi$7
Martingale	$\pi$8	Solve $\pi$9 or $\mathbb{R}^3$0 nonlinearly
Column-mass	$\mathbb{R}^3$1	$\mathbb{R}^3$2

Each projection is exact in the KL geometry (Langbridge et al., 25 Jun 2026).

4. Convergence and Theoretical Guarantees

The central convergence result asserts exponential (linear) convergence in dual value and L² error under mild and checkable conditions:

• Compact support of $\mathbb{R}^3$3
• Bounded, Lipschitz cost
• Existence of an interior admissible coupling

Quantitatively,

$\mathbb{R}^3$4

for some $\mathbb{R}^3$5.

In the abstract setting, the contraction rate is governed by the sum of linear subspaces (for marginal and martingale constraints) and their respective Friedrichs angles. Specifically, the explicit MPMS contraction rate is (Eckstein et al., 27 Feb 2025):

$\mathbb{R}^3$6

where $\mathbb{R}^3$7 is the cosine of the Friedrichs angle between marginal and martingale constraint spaces, and $\mathbb{R}^3$8, $\mathbb{R}^3$9 are strong convexity and Lipschitz constants of the dual functional. Convergence remains exponential except in the degenerate limit where the angle vanishes.

For practical finite-dimensional problems, 200–500 cycles suffice for high-precision results ($Q$0 accuracy) (Langbridge et al., 25 Jun 2026).

5. Variant Algorithms and Newton Acceleration

Several variations have been developed for computational acceleration and broader structural constraints:

• In (Tang et al., 25 Aug 2025), MPMS is paired with sparse-Newton acceleration: after warm-starting with standard Sinkhorn cycles, a second-order local step exploits approximate sparsity in the Hessian for rapid convergence. The Hessian is thresholded to maintain computational feasibility.
• MPMS generalizes to supermartingale constraints (by imposing row-wise inequalities).
• For one-step martingale Benamou–Brenier problems (Bass martingale computation), a continuous variant of MPMS solves backward and forward heat equations with alternating pushforwards to match marginals, remaining well-posed in higher dimensions via $Q$1-dimensional parabolic equations (Joseph et al., 2023).

Key elements for stable implementation include log-domain storage of scaling vectors, root-finding for martingale multipliers, and (optionally) Anderson or Aitken acceleration (Langbridge et al., 25 Jun 2026).

6. Applications: Stochastic Volatility Model Calibration and Quantitative Finance

MPMS provides the computational backbone for several robust calibration tasks:

• In multi-period stochastic volatility calibration, known marginals at option maturities and a prior model (e.g., Heston) yield a chain of EMOT subproblems. MPMS ensures preservation of calibration and arbitrage-free, no-loss-of-martingale structure at each link in the chain (Chen et al., 2024).
• One-step (mBB) and multi-step discrete martingale problems for statistical or financial applications are now tractable at scale.
• Numerical experiments in (Joseph et al., 2023) confirm high-precision convergence across synthetic and real market data, with typical mean-squared CDF error orders $Q$2–$Q$3 in tens of cycles.

7. Geometric and Analytical Insights

The convergence rate and stability of MPMS are governed by the geometry of the affine constraint subspaces. The key determinant is the Friedrichs angle between the aggregated marginal constraints and the martingale constraint space. The closer (less orthogonal) these spaces, the slower the rate. Under strong convexity and closed-sum hypotheses—typically met for reasonable marginals and costs—MPMS enjoys robust, linear convergence even in high dimensions (Eckstein et al., 27 Feb 2025).

The alternating projection structure, strong convexity, and explicit characterization of the dual functional guarantee both practical effectiveness and theoretical transparency. This structure unifies IPFP-type algorithms across optimal transport, adapted/martingale constraints, and multi-marginal extensions.