Learning solutions of parameterized stiff ODEs using Gaussian processes

TL;DR

The paper tackles the challenge of learning parameterized solutions to stiff ODEs with Gaussian processes, where nonstationarity hampers accuracy and data efficiency. It introduces a data-driven arc-length based curve reparameterization, coupled with a flexible, monotone σ function built from Hermite splines, to render the target solution approximately stationary without altering GP internals. The method extends to multi-parameter inputs and demonstrates substantial convergence gains (typically 2–6×) across benchmarks like the Van der Pol oscillator, Tunnel Diode Oscillator, and Brusselator, with negligible overhead. This approach enables efficient, uncertainty-aware surrogate modeling of parameterized stiff ODEs and offers promising directions for higher-dimensional reparameterizations and PDE applications.

Abstract

Stiff ordinary differential equations (ODEs) play an important role in many scientific and engineering applications. Often, the dependence of the solution of the ODE on additional parameters is of interest, e.g.\ when dealing with uncertainty quantification or design optimization. Directly studying this dependence can quickly become too computationally expensive, such that cheaper surrogate models approximating the solution are of interest. One popular class of surrogate models are Gaussian processes (GPs). They perform well when approximating stationary functions, functions which have a similar level of variation along any given parameter direction, however solutions to stiff ODEs are often characterized by a mixture of regions of rapid and slow variation along the time axis and when dealing with such nonstationary functions, GP performance frequently degrades drastically. We therefore aim to reparameterize stiff ODE solutions based on the available data, to make them appear more stationary and hence recover good GP performance. This approach comes with minimal computational overhead and requires no internal changes to the GP implementation, as it can be seen as a separate preprocessing step. We illustrate the achieved benefits using multiple examples.

Figures (10)

• Figure 1: Left: stiff ODE solution with regions of rapid variation that have complicated shapes. Right: one particular solution for $p_1 = 10^{-6}$ (indicated by the dashed line on the left) illustrating the extremely stiff behavior more closely.
• Figure 1: Exemplary sequential designs based on the $k_\infty$ Matérn (Gaussian) kernel for the VDPO with $p = 1$ (left) and $p = 10$ (right). The plots show the original solution $x(t)$, the prediction $\hat{x}(t)$, the posterior standard deviation $\hat{\sigma}(t)$ as well as the observations used for training $(t_i, x_i)$. The plots on the right illustrate the issues for stiff solutions.
• Figure 1: Reparameterized solutions $\tilde{x}(s)$ using the linear definition of $\sigma$, both for the non-stiff case (left) and the stiff case (right). The original solutions are indicated in light grey with the reparameterized solutions in solid green and the reparameterizations $\tilde{\tau}^{-1}(s)$ as dashed green lines. The reparameterizations are scaled by a factor to fit the plots.
• Figure 1: Convergence of sequential designs for the reparameterized solutions (top) and reparameterizations (bottom), based on different parameter values $p$ and different choices for $\sigma$ (indicated by $\mu$).
• Figure 2: Convergence of the sequential designs, both for the non-stiff case with $p = 1$ (left) and the stiff case with $p = 10$ (right).

Theorems & Definitions (10)

• Proposition 3.1
• Proposition 3.2
• Remark 3.3: Relation to warping
• Remark 3.4
• Proposition A.1
• Proof 1
• Proposition A.2
• Proof 2
• Proof 3
• Proof 4