Table of Contents
Fetching ...

Newton-Krylov Methods for Computing Steady States of Particle Timesteppers via Optimal Transport

Hannes Vandecasteele, Nicholas Karris, Alexander Cloninger, Ioannis G. Kevrekidis

TL;DR

The paper addresses computing steady states of stochastic particle timesteppers, where direct fixed-point equations are ill-defined due to randomness. It introduces a Wasserstein-distance-based residual and reinterprets particle ensembles as probability measures, enabling matrix-free Newton–Krylov methods to converge to steady-state distributions despite high noise. To overcome noise-induced obstacles in higher dimensions, it develops smooth distribution representations—notably ICDF-to-ICDF in 1D and CDF- or sliced-Wasserstein-based timesteppers in higher dimensions—and demonstrates fast, robust convergence on representative 1D and 2D problems, including bacterial chemotaxis, economic agents, and a half-moon potential. The framework provides a unified variational approach bridging deterministic and stochastic timesteppers, with practical implications for variance reduction, continuation, and scalable computation of distributional steady states in complex systems.

Abstract

Timesteppers constitute a powerful tool in modern computational science and engineering. Although they are typically used to advance the system forward in time, they can also be viewed as nonlinear mappings that implicitly encode steady states and stability information. In this work, we present an extension of the matrix-free framework for calculating, via timesteppers, steady states of deterministic systems to stochastic particle simulations, where intrinsic randomness prevents direct steady state extraction. By formulating stochastic timesteppers in the language of optimal transport, we reinterpret them as operators acting on probability measures rather than on individual particle trajectories. This perspective enables the construction of smooth cumulative- and inverse-cumulative-distribution-function ((I)CDF) timesteppers that evolve distributions rather than particles. Combined with matrix-free Newton-Krylov solvers, these smooth timesteppers allow efficient computation of steady-state distributions even under high stochastic noise. We perform an error analysis quantifying how noise affects finite-difference Jacobian action approximations, and demonstrate that convergence can be obtained even in high noise regimes. Finally, we introduce higher-dimensional generalizations based on smooth CDF-related representations of particles and validate their performance on a non-trivial two-dimensional distribution. Together, these developments establish a unified variational framework for computing meaningful steady states of both deterministic and stochastic timesteppers.

Newton-Krylov Methods for Computing Steady States of Particle Timesteppers via Optimal Transport

TL;DR

The paper addresses computing steady states of stochastic particle timesteppers, where direct fixed-point equations are ill-defined due to randomness. It introduces a Wasserstein-distance-based residual and reinterprets particle ensembles as probability measures, enabling matrix-free Newton–Krylov methods to converge to steady-state distributions despite high noise. To overcome noise-induced obstacles in higher dimensions, it develops smooth distribution representations—notably ICDF-to-ICDF in 1D and CDF- or sliced-Wasserstein-based timesteppers in higher dimensions—and demonstrates fast, robust convergence on representative 1D and 2D problems, including bacterial chemotaxis, economic agents, and a half-moon potential. The framework provides a unified variational approach bridging deterministic and stochastic timesteppers, with practical implications for variance reduction, continuation, and scalable computation of distributional steady states in complex systems.

Abstract

Timesteppers constitute a powerful tool in modern computational science and engineering. Although they are typically used to advance the system forward in time, they can also be viewed as nonlinear mappings that implicitly encode steady states and stability information. In this work, we present an extension of the matrix-free framework for calculating, via timesteppers, steady states of deterministic systems to stochastic particle simulations, where intrinsic randomness prevents direct steady state extraction. By formulating stochastic timesteppers in the language of optimal transport, we reinterpret them as operators acting on probability measures rather than on individual particle trajectories. This perspective enables the construction of smooth cumulative- and inverse-cumulative-distribution-function ((I)CDF) timesteppers that evolve distributions rather than particles. Combined with matrix-free Newton-Krylov solvers, these smooth timesteppers allow efficient computation of steady-state distributions even under high stochastic noise. We perform an error analysis quantifying how noise affects finite-difference Jacobian action approximations, and demonstrate that convergence can be obtained even in high noise regimes. Finally, we introduce higher-dimensional generalizations based on smooth CDF-related representations of particles and validate their performance on a non-trivial two-dimensional distribution. Together, these developments establish a unified variational framework for computing meaningful steady states of both deterministic and stochastic timesteppers.
Paper Structure (29 sections, 3 theorems, 82 equations, 17 figures)

This paper contains 29 sections, 3 theorems, 82 equations, 17 figures.

Key Result

Theorem 2.1

Let $\sigma,\mu\in\mathcal{P}_2(\mathbb{R}^d)$. If $\sigma$ has a density with respect to the Lebesgue measure, then the optimal coupling $\gamma$ that satisfies eq:wass-dist-defn is unique, and there exists a $\sigma$-a.e. unique map $T:\mathbb{R}^d\to\mathbb{R}^d$ such that $\gamma = (\mathop{\mat Moreover, there exists a $\sigma$-a.e. unique (up to additive constant) convex function $\phi$ such

Figures (17)

  • Figure 1: Numerical results for the Wasserstein-Adam optimizer on the chemotaxis model. Left: Wasserstein loss (blue) and gradient-norm (orange) per epoch. Right: Histogram of the initial particles in blue, the Wasserstein-Adam optimized particles in orange, and the analytic steady-state density (see equation \ref{['eq:chemotaxis_ss']}) in red.
  • Figure 2: Numerical results for the Wasserstein-Adam optimizer on the economic agents model. Left: Wasserstein loss (blue) and gradient-norm (orange) per epoch. Right: Histogram of the initial particles in blue, the Wasserstein-Adam optimized particles in orange, and the analytic steady-state density (see Figure \ref{['fig:agent_pde_ss']}) in red.
  • Figure 3: A color plot of the half-moon potential, see equation \ref{['eq:halfmoonpotential']} for more details.
  • Figure 4: Numerical results for the Wasserstein-Adam optimizer on the half-moon potential. Left: Wasserstein loss (blue) and gradient-norm (orange) per epoch. Right: Colormap of the 2D histogram of the optimized particles in orange.
  • Figure 5: Numerical results for the Newton--Krylov optimizer for bacterial chemotaxis. Left: Norm of the objective function, $\left\lVert\tilde{F}(X_k)\right\rVert$ per iteration. Right: Initial particles (blue), optimized particles (orange) and analytic steady-state distribution (dashed red).
  • ...and 12 more figures

Theorems & Definitions (6)

  • Theorem 2.1: Brenier brenier1991polar
  • Theorem 2.2: Proposition 2.1 in peyre2019computational
  • Remark 1
  • Remark 2
  • Theorem A.1
  • proof