Artificial neural network solver for Fokker-Planck and Koopman eigenfunctions

TL;DR

This paper introduces a neural network-based method to compute eigenfunctions of Fokker-Planck and Koopman operators for stochastic differential equations, addressing high-dimensional PDE challenges with improved accuracy and less data.

Contribution

The paper presents a novel neural network solver that incorporates differential operators into the loss function for efficient eigenfunction computation of SDEs.

Findings

Effective in 2D, 3D, and 4D examples

Reduces need for large training datasets

Improves accuracy over traditional methods

Abstract

For a stochastic differential equation (SDE) that is an It\^{o} diffusion or Langevin equation, the Fokker-Planck operator governs the evolution of the probability density, while its adjoint, the infinitesimal generator of the stochastic Koopman operator, governs the evolution of system observables, in the mean. The eigenfunctions of these operators provide a powerful framework to analyze SDEs, and have shown to be particularly useful for systems of stochastic oscillators. However, computing these eigenfunctions typically requires solving high-dimensional PDEs on unbounded domains, which is numerically challenging. Building on previous work, we propose a data-driven artificial neural network solver for Koopman and Fokker-Planck eigenfunctions. Our approach incorporates the differential operator into the loss function, improving accuracy and reducing dependence on large amounts of accurate training data. We demonstrate our approach on several numerical examples in two, three, and four dimensions.

Equations

$$d\mathbf{X}=\mathbf{f}(\mathbf{X})dt+\mathbf{g}(\mathbf{X})d\mathbf{W}(t)$$

$$\rho (\mathbf{y},t|\mathbf{x},s)=\underset{|d\mathbf{y}|\to 0}{\lim }\frac{\mathbb{P}\left(\mathbf{X}(t)\in [\mathbf{y},\mathbf{y}+d\mathbf{y})|\mathbf{X}(s)=\mathbf{x}\right)}{|d\mathbf{y}|}$$

$$\frac{\partial }{\partial t}\rho =\mathcal{L}[\rho ]=-\sum _i{\partial }_{x_i}[\mathbf{f}\rho ]+\frac{1}{2}\sum _{ij}{\partial }_{x_i{x}_j}\{[\mathbf{g}{\mathbf{g}}^T{]}_{ij}\rho \}$$

$$-\frac{\partial }{\partial s}\rho ={\mathcal{L}}^{†}[\rho ]=\sum _i{\mathbf{f}}_i{\partial }_{y_i}[\rho ]+\frac{1}{2}\sum _{ij}[\mathbf{g}{\mathbf{g}}^T{]}_{ij}{\partial }_{y_i{y}_j}\{\rho \}$$

$$\begin{array}{rl}\mathcal{L}[P_{\lambda }(\mathbf{x})]& =\lambda P_{\lambda }(\mathbf{x})\\ {\mathcal{L}}^{†}[Q_{\lambda }^{\ast }(\mathbf{y})]& =\lambda Q_{\lambda }^{\ast }(\mathbf{y})\end{array}$$

$$\rho (\mathbf{y},t|\mathbf{x},s)=P_0+\sum _{\lambda \ne 0}{e}^{\lambda (t-s)}{P}_{\lambda }(\mathbf{y})Q_{\lambda }^{\ast }(\mathbf{x})$$

$$\rho (\mathbf{y},t|\mathbf{x},s)-P_0(\mathbf{y})\approx e^{{\lambda }_1(t-s)}{P}_{{\lambda }_1}(\mathbf{y})Q_{{\lambda }_1}^{\ast }(\mathbf{x})+\text{c.c.}$$

$$f(t;\mathbf{b})=b_1\exp (b_{5}t)\sin \left(\frac{2\pi t}{b_2}+\frac{2\pi }{b_3}\right)+b_4$$

$$F(\mathbf{t};\mathbf{b})=\sum _{j=t_s}^{t_f}({\tilde{\rho }}_{f,i}(t_j)-f(t_j;\mathbf{b}){)}^2$$

$$\rho (\mathbf{y},t|\mathbf{x},0)-P_0(\mathbf{y})\approx e^{{\lambda }_{1}t}{P}_{{\lambda }_1}(\mathbf{y})Q_{{\lambda }_1}^{\ast }(\mathbf{x})+\text{c.c.}=2\Re \left[e^{{\lambda }_{1}t}{P}_{{\lambda }_1}(\mathbf{y})Q_{{\lambda }_1}^{\ast }(\mathbf{x})\right]$$

$${\tilde{\rho }}_{f,i}(t_j)-P_0({\mathbf{x}}_i)=2\Re [\exp ({\tilde{\lambda }}_1{t}_j)P_{{\lambda }_1}({\mathbf{x}}_i)Q_{{\lambda }_1}^{\ast }({\mathbf{x}}_0)]$$

$${\tilde{P}}_{{\lambda }_1}({\mathbf{x}}_i)=2Q_{{\lambda }_1}^{\ast }({\mathbf{x}}_0)P_{{\lambda }_1}({\mathbf{x}}_i)$$

$${\tilde{\rho }}_{f,i}(t_j)-P_0({\mathbf{x}}_i)=\Re [\exp ({\tilde{\lambda }}_1{t}_j){\tilde{P}}_{{\lambda }_1}({\mathbf{x}}_i)]$$

$$\begin{split}\tilde{\rho}_{f,i}(t_{j})-P_{0}(\mathbf{x}_{i})-e^{\tilde{\mu}_{1}t_{j}}\Bigg{[}&\cos(\tilde{\omega}_{1}t_{j})P_{\lambda_{1}}^{R}(\mathbf{x}_{i})-\sin(\tilde{\omega}_{1}t_{j})P_{\lambda_{1}}^{I}(\mathbf{x}_{i})\Bigg{]}=0\end{split}$$

$${\mathbf{b}}_f=\left[\begin{array}{c}\exp (-{\tilde{\mu }}_1{t}_{T_s})[{\tilde{\rho }}_{f,i}(t_{T_s})-P_0({\mathbf{x}}_i)]\\ \exp (-{\tilde{\mu }}_1{t}_{T_s+1})[{\tilde{\rho }}_{f,i}(t_{T_s+1})-P_0({\mathbf{x}}_i)]\\ ⋮\\ \exp (-{\tilde{\mu }}_1{t}_{T_f})[{\tilde{\rho }}_{f,i}(t_{T_f})-P_0({\mathbf{x}}_i)]\end{array}\right],{\mathbf{u}}_f=\left[\begin{array}{c}{P}_{{\lambda }_1}^R({\mathbf{x}}_i)\\ P_{{\lambda }_1}^I({\mathbf{x}}_i)\end{array}\right]$$

$${\mathbf{A}}_f=\left[\begin{array}{cc}\cos ({\tilde{\omega }}_1{t}_{T_s})& -\sin ({\tilde{\omega }}_1{t}_{T_s})\\ \cos ({\tilde{\omega }}_1{t}_{T_s+1})& -\sin ({\tilde{\omega }}_1{t}_{T_s+1})\\ ⋮& \\ \cos ({\tilde{\omega }}_1{t}_{T_f})& -\sin ({\tilde{\omega }}_1{t}_{T_f})\end{array}\right]$$

$${\tilde{\rho }}_{b,i}(t_j)-{\int }_{\mathcal{B}}d\mathbf{y}{P}_0(\mathbf{y})=2\Re \left[\exp ({\tilde{\lambda }}_1{t}_j)\left({\int }_{\mathcal{B}}d\mathbf{y}{P}_{{\lambda }_1}(\mathbf{y})\right)Q_{{\lambda }_1}^{\ast }({\mathbf{x}}_i)\right]$$

$${\tilde{Q}}_{{\lambda }_1}^{\ast }({\mathbf{x}}_i)=2\left({\int }_{\mathcal{B}}d\mathbf{z}{P}_{{\lambda }_1}(\mathbf{z})\right)Q_{{\lambda }_1}^{\ast }({\mathbf{x}}_i)$$

$${\tilde{\rho }}_{b,i}(t_j)-{\int }_{\mathcal{B}}d\mathbf{z}{P}_0(\mathbf{z})=\Re [\exp ({\tilde{\lambda }}_1{t}_j){\tilde{Q}}_{{\lambda }_1}^{\ast }({\mathbf{x}}_i)]$$

$$\begin{split}\tilde{\rho}_{b,i}(t_{j})-\int_{\mathcal{B}}d\mathbf{z}\;P_{0}(\mathbf{z})-e^{\tilde{\mu}_{1}t_{j}}\Bigg{[}&\cos(\tilde{\omega}_{1}t_{j})Q_{\lambda_{1}}^{R}(\mathbf{x}_{i})-\sin(\tilde{\omega}_{1}t_{j})Q_{\lambda_{1}}^{I}(\mathbf{x}_{i})\Bigg{]}=0\end{split}$$

$${\mathbf{b}}_b=\left[\begin{array}{c}\exp (-{\tilde{\mu }}_1{t}_{T_s})[{\tilde{\rho }}_{b,i}(t_{T_s})-{\int }_{\mathcal{B}}d\mathbf{z}{P}_0(\mathbf{z})]\\ \exp (-{\tilde{\mu }}_1{t}_{T_s+1})[{\tilde{\rho }}_{b,i}(t_{T_s+1})-{\int }_{\mathcal{B}}d\mathbf{z}{P}_0(\mathbf{z})]\\ ⋮\\ \exp (-{\tilde{\mu }}_1{t}_{T_f})[{\tilde{\rho }}_{b,i}(t_{T_f})-{\int }_{\mathcal{B}}d\mathbf{z}{P}_0(\mathbf{z})]\end{array}\right],{\mathbf{u}}_b=\left[\begin{array}{c}{Q}_{{\lambda }_1}^R({\mathbf{x}}_i)\\ Q_{{\lambda }_1}^I({\mathbf{x}}_i)\end{array}\right]$$

$${\mathbf{A}}_b=\left[\begin{array}{cc}\cos ({\tilde{\omega }}_1{t}_{T_s})& -\sin ({\tilde{\omega }}_1{t}_{T_s})\\ \cos ({\tilde{\omega }}_1{t}_{T_s+1})& -\sin ({\tilde{\omega }}_1{t}_{T_s+1})\\ ⋮& \\ \cos ({\tilde{\omega }}_1{t}_{T_f})& -\sin ({\tilde{\omega }}_1{t}_{T_f})\end{array}\right]$$

$$\begin{split}\tilde{\rho}_{f,i}(t_{j})-P_{0}(\mathbf{x}_{i})&-\sum_{k=1}^{M}e^{\tilde{\mu}_{k}t_{j}}\Bigg{[}\cos(\tilde{\omega}_{k}t_{j})P_{\lambda_{k}}^{R}(\mathbf{x}_{i})-\sin(\tilde{\omega}_{k}t_{j})P_{\lambda_{k}}^{I}(\mathbf{x}_{i})\Bigg{]}=0\end{split}$$

$$\begin{split}\tilde{\rho}_{b,i}(t_{j})-\int_{\mathcal{B}}d\mathbf{z}\;P_{0}(\mathbf{z})&-\sum_{k=1}^{M}e^{\tilde{\mu}_{k}t_{j}}\Bigg{[}\cos(\tilde{\omega}_{k}t_{j})Q_{\lambda_{k}}^{R}(\mathbf{x}_{i})-\sin(\tilde{\omega}_{k}t_{j})Q_{\lambda_{k}}^{I}(\mathbf{x}_{i})\Bigg{]}=0\end{split}$$

$${\mathcal{F}}_{\mathcal{L}}({\lambda }_n,U_{{\lambda }_n})=\Vert (\mathcal{L}-{\lambda }_n)U_{{\lambda }_n}{\Vert }_{L^2(\mathcal{R})}^2+\Vert U_{{\lambda }_n}-{\tilde{U}}_{{\lambda }_n}{\Vert }_{L^2(\mathcal{R})}^2$$

$$\begin{array}{rl}{J}_{\mathcal{L}}(\chi )& =\frac{1}{N_y}\sum _{j=1}^{N_y}{\left(\mathcal{L}[{\mathcal{N}}^R({\mathbf{y}}_j;\chi )]-{\tilde{\mu }}_n{\mathcal{N}}^R({\mathbf{y}}_j;\chi )+{\tilde{\omega }}_n{\mathcal{N}}^I({\mathbf{y}}_j;\chi )\right)}^2\\ & +\frac{1}{N_y}\sum _{j=1}^{N_y}{\left(\mathcal{L}[{\mathcal{N}}^I({\mathbf{y}}_j;\chi )]-{\tilde{\mu }}_n{\mathcal{N}}^I({\mathbf{y}}_j;\chi )-{\tilde{\omega }}_n{\mathcal{N}}^R({\mathbf{y}}_j;\chi )\right)}^2\\ & +\frac{1}{N_x}\sum _{i=1}^{N_x}{\left({\mathcal{N}}^R({\mathbf{y}}_j;\chi )-\Re [{\tilde{U}}_{{\lambda }_n}({\mathbf{x}}_i)]\right)}^2+\frac{1}{N_x}\sum _{i=1}^{N_x}{\left({\mathcal{N}}^I({\mathbf{y}}_j;\chi )-\Im [{\tilde{U}}_{{\lambda }_n}({\mathbf{x}}_i)]\right)}^2\end{array}$$

$$\mathcal{L}[{\mathcal{N}}^R(\mathbf{x};{\chi }^{\ast })+i{\mathcal{N}}^I(\mathbf{x};{\chi }^{\ast })]={\lambda }_n[{\mathcal{N}}^R(\mathbf{x};{\chi }^{\ast })+i{\mathcal{N}}^I(\mathbf{x};{\chi }^{\ast })]$$

$$\underset{\beta \in {\mathbb{R}}^{2M}}{\min }\Vert \mathcal{X}\beta -\mathcal{Y}{\Vert }^2$$

$$\begin{array}{rl}\mathcal{X}& =\left[\begin{array}{cccccc}{e}^{{\mu }_1{t}_1}\cos ({\omega }_1{t}_1)& e^{{\mu }_1{t}_1}\sin ({\omega }_1{t}_1)& \cdots & \cdots & e^{{\mu }_M{t}_1}\cos ({\omega }_M{t}_1)& e^{{\mu }_M{t}_1}\sin ({\omega }_M{t}_1)\\ e^{{\mu }_1{t}_2}\cos ({\omega }_1{t}_2)& e^{{\mu }_1{t}_2}\sin ({\omega }_1{t}_2)& \cdots & \cdots & e^{{\mu }_M{t}_2}\cos ({\omega }_M{t}_2)& e^{{\mu }_M{t}_2}\sin ({\omega }_M{t}_2)\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ e^{{\mu }_1{t}_N}\cos ({\omega }_1{t}_N)& e^{{\mu }_1{t}_N}\sin ({\omega }_1{t}_N)& \cdots & \cdots & e^{{\mu }_M{t}_N}\cos ({\omega }_M{t}_N)& e^{{\mu }_M{t}_N}\sin ({\omega }_M{t}_N)\end{array}\right]\\ \mathcal{Y}& =[{\tilde{\rho }}_{f,i}(t_1)-P_0(\mathbf{x}),\cdots ,{\tilde{\rho }}_{f,i}(t_N)-P_0(\mathbf{x}){]}^T\end{array}$$

$${\beta }_0:=[P_{{\lambda }_1}^R({\mathbf{x}}_{\mathbf{i}}),P_{{\lambda }_1}^I({\mathbf{x}}_{\mathbf{i}}),\cdots ,P_{{\lambda }_M}^R({\mathbf{x}}_{\mathbf{i}}),P_{{\lambda }_M}^I({\mathbf{x}}_{\mathbf{i}}){]}^T$$

Full text

{MSCcodes}

60H10, 37M10, 65M75, 92C20, 47N40

1 Introduction

Many physical and biological systems have some degree of noise or uncertainty, and are widely modeled as systems of stochastic differential equations. Traditionally, such systems have been characterized by studying the time-dependent and stationary probability densities associated with the underlying stochastic dynamics [13, 14, 43]. The Fokker-Planck operator describes the evolution of the probability density analytically.

However, recent attention has focused on the infinitesimal generator of the stochastic Koopman operator (SKO)111For brevity, we will write “SKO” to mean the infinitesimal generator of the stochastic Koopman operator. as a complementary descriptor of stochastic processes. The SKO is adjoint to the Fokker-Planck operator, and describes the evolution of the expectation of system observables. Because the action of the SKO is linear on the space of $\mathcal{C}^{2}$-continuous system observables, the eigenspectrum of the SKO characterizes the original nonlinear, stochastic dynamics [9]. SKO approaches have gained significant traction in fluid dynamics [36, 44], control problems [26, 39], neuroscience [3, 34], and stochastic oscillators [29, 40, 41, 51].

The space of system observables is infinite-dimensional, so obtaining finite dimensional approximations of the SKO is important for practical computations [6]. While traditional PDE solvers (finite differences, finite elements) are effective for low-dimensional problems, they lose effectiveness for high-dimensional problems $(n\geq 3)$. Non-trivial boundary conditions on unbounded domains can further complicate numerical implementation, as can the potential presence of continuous spectra [35]. While much effort has been devoted to overcome these challenges, solving high-dimensional Fokker-Planck and Koopman eigenvalue problems is still a challenging task.

Machine learning is a powerful approach for solving high-dimensional PDEs and related eigenvalue problems. Some works have employed diffusion Monte Carlo methods [16, 46], deep Ritz methods [22, 33], physics-informed neural network (PINN) [42, 49], and BSDE methods [15, 50, 53], among other approaches [1, 2, 17, 54]. Recent work proposed a data-driven scheme that replaces boundary conditions by Monte Carlo simulation data when solving the Fokker-Planck equation [11, 31]. This method was extended into an artificial neural network (ANN) approach, and was found to be effective for solving high-dimensional problems [57]. The idea behind the method is to train an ANN to approximate a probability density function that minimizes both (i) the distance to the Monte Carlo reference data, and (ii) the norm of the residual of the forward operator $\mathcal{L}$. Importantly, it was found that relatively inaccurate Monte Carlo approximations at only a small number of reference points were sufficient to guide the ANN to an accurate solution of the invariant distribution. We also refer to [7] for a PINN-based Fokker-Planck solver that uses a different regularization term based on Monte Carlo simulations.

In this manuscript, we will extend previous results [31, 57] to develop a novel ANN approach to solve for the low-lying eigenfunctions and eigenvalues of the SKO and Fokker-Planck operators associated with stochastic oscillators. Leveraging results from [51], we use Monte Carlo approximations of the time-dependent probability densities to estimate the low-lying eigenvalues of a given system. Then, we obtain coarse approximations of the low-lying forward and Koopman eigenfunctions as solutions of a linear least squares problem. This Monte Carlo data for the eigenmodes serves as a reference in a loss function (as a replacement for the boundary conditions, which are not explicitly considered), which is minimized with an ANN, analogous to the approach considered in [57].

Our ANN approach is advantageous for several reasons. Existing ANN solvers for eigenfunction mainly use normalization to avoid converging to the trivial solution [16, 23]. In comparison, our approach directly estimates the value of an eigenfunction at numerous collocation points in the domain, leading to both fast convergence and improved solution quality. Other methods, such as extended dynamic mode decomposition (EDMD) and variants, have recently gained traction as data-driven methods to approximate the low-lying spectrum of the SKO [4, 5, 8, 37, 27, 55]. However, such methods rely on choosing a dictionary of functions; a poor choice of dictionary results in poor performance of the EDMD algorithm, and dictionaries which are too large often result in computational bottlenecks and ill-conditioned systems. In contrast, our novel ANN approach does not require a-priori system knowledge, nor a collection of dictionary functions. Instead, it requires only time-series data, making it more broadly applicable to a wide range of problems. Our ANN approach is also capable of simultaneously approximating several low-lying eigenfunctions of the Fokker-Planck and Koopman operators. Importantly, the user is free to define reference points in the domain at which solution values are required. This mesh-free approach is advantageous for solving high-dimensional problems because it requires relatively few reference points.

The remainder of this work is organized in the following manner. We first describe necessary assumptions and mathematical preliminaries in §2. Then, we detail our ANN approach and numerical algorithms in §3, discuss error associated with the eigenfunction approximation in §4, and apply our technique to several examples in two, three, and four dimensions in §5.

2 Mathematical preliminaries

We consider stochastic oscillators of the form

[TABLE]

with $\mathbf{X}(t)$ a stochastic process taking values in a domain $\Omega\subseteq\mathbb{R}^{n}$, $\mathbf{f}:\Omega\to\mathbb{R}^{n}$ a $\mathcal{C}^{2}$-continuous vector field, $\mathbf{g}$ an $n\times m$ matrix with $\mathcal{C}^{2}$-continuous entries such that $\mathbf{g}\mathbf{g}^{T}$ is nonsingular for all $\mathbf{x}\in\Omega$, and $d\mathbf{W}(t)$ an $m\times 1$ vector of increments of $m$ independent Wiener process. We interpret (1) in the sense of Itô.

We assume throughout that the time-dependent transition density of the Markov process

[TABLE]

exists for all $t>s$ and for all $\mathbf{x}$ and $\mathbf{y}$ in $\Omega$. The evolution of the density associated with (1) satisfies the Fokker-Planck (forward) equation, which defines the forward operator $\mathcal{L}$:

[TABLE]

and also the backward equation, which defines the backward operator $\mathcal{L}^{\dagger}$:

[TABLE]

for $1\leq i,j\leq n$. The operator $\mathcal{L}^{\dagger}$, aka the generator of the Markov process $\mathbf{X}$, is the infinitesimal generator of the stochastic Koopman operator [9]. We henceforth refer to $\mathcal{L}^{\dagger}$ as “the SKO” for brevity. We are concerned with the numerical solution of the forward and backward eigenvalue problems

[TABLE]

In particular, we focus on solving (5) for the lowest-lying eigenpairs.

We assume that (1) is robustly oscillatory: (i) all eigenvalues, apart from the simple trivial eigenvalue $\lambda_{0}\equiv 0$, have negative real parts, (ii) the first nontrivial eigenvalue, $\lambda_{1}=\mu\pm i\omega$, is a unique complex-conjugate pair, with quality factor $|\omega/\mu|\gg 1$, and (iii) all other eigenvalues $\lambda^{\prime}$ have more negative real parts, $\Re[{\lambda^{\prime}}]<\upsilon\,\Re[{\lambda_{1}}]$, where $0<\upsilon\in\mathbb{R}$. We further assume that (1) has a discrete set of eigenvalues with biorthonormal eigenfunctions, and has a unique stationary distribution. Under these assumptions, the forward transition density may be expressed as a linear combination of the forward and backward eigenmodes

[TABLE]

for times $t>s$. In cases where there is a large spectral gap $\left|\Re[\lambda_{1}]-\sup\Re[\lambda^{\prime}]\right|$, previous works have truncated the infinite sum at only the first term [40, 41, 51]

[TABLE]

While our focus is on the $Q$-function, the unique complex lowest-lying eigenpair, in practice we achieve a more faithful approximation of system dynamics by including the first several eigenpairs, i.e., by truncating (6) after several terms. In the following, we will use (7) to generate data-driven approximations of the forward and backward eigenfunctions.222Methods to extract approximations of the $Q$-function from time-series data were first introduced in [51] by Thomas and Lindner. By assuming knowledge of high-order eigenvalues, i.e., including more terms using (6), we will obtain more accurate solutions for the $Q$-function.

Finally, we recall that the forward eigenfunctions should satisfy “vanishing boundary conditions” at infinity, and that the backward eigenfunctions should have vanishing derivatives at infinity, when $\Omega$ is unbounded. In practice, it is not possible to implement these boundary conditions; instead, we restrict attention to a compact, rectangular computational domain, $\mathcal{R}\subset\Omega$.

3 Numerical methods and ANN implementation

We begin with an overview of our novel ANN approach to solve the eigenvalue problems (5). First, realizations of the stochastic process (1) are used to generate Monte Carlo approximations of the time-dependent probability density. The relaxation of the time-dependent probability density to an invariant distribution allows for an accurate estimation of the low-lying eigenvalue, $\lambda_{1}$, of the SKO and forward operator. Then, we formulate a linear least squares problem using the eigenfunction expansion ansatz (7). The solution of the least squares problem gives a coarse approximation of the forward and backward eigenmodes. Finally, the least squares data is fed into an ANN, which interpolates and smooths the data and increases the accuracy of the approximation.

3.1 Step 1: collocation points

Here, the goal is to generate two sets of collocation points: $\mathcal{X}=\{\mathbf{x}_{i},\;i=1,2,\dots,N_{x}\}$ and $\mathcal{Y}=\{\mathbf{y}_{j},\;j=1,2,\dots,N_{y}\}$. We call $\mathcal{X}$ the “training set” and $\mathcal{Y}$ the “reference set”. Ultimately, a neural network will learn a specific eigenfunction on the training set to replace boundary conditions, and will learn the differential operator, $\mathcal{L}$ or $\mathcal{L}^{\dagger}$, on the reference set. In other words, training on $\mathcal{Y}$ ensures that the output of the ANN satisfies (approximately) either the forward or backward equation, while training on $\mathcal{X}$ guides the ANN to a specific normalization of an eigenfunction and avoids the trivial solution.

To generate the training set, $\mathcal{X}$, we discretize the rectangular computational domain $\mathcal{R}$ into boxes of equal dimension, but we store in memory only $N_{x}$ of the boxes corresponding to the chosen collocation points. Each box is represented by the coordinates corresponding to the center of the respective box. This step is critical in higher dimensions, where storing a full mesh is not feasible. The points, $\mathbf{x}_{i}$, at the centers of the boxes will be associated with numerical estimates of the probability density. In contrast, we do not bin the reference set in boxes because these points are used solely for minimizing the residual of the operator loss.

The support of the forward eigenmodes is often restricted to the invariant distribution, which may be concentrated near low dimensional manifolds or in specific regions of the domain. While to the best of our knowledge this is not supported by known results, we expect concentration results of invariant probability measures [20, 21, 32] and quasi-stationary distributions [48] to be applicable after some modifications. In contrast, the support of the backward eigenmodes typically spans all of $\Omega$. The generation of the collocation points should take this difference into account.

To generate collocation points, $\mathcal{X}$, we first run (1) for a “burn-in” time, $T_{\text{burn}}$, so that the location of the terminal point of the process is effectively drawn from a distribution that is near the stationary distribution. Then, we choose a “ratio”, $\alpha\in[0,1]$, continue the realization of (1) starting from the terminal point of the burn-in simulation, and take only a fraction $\alpha$ (one should generally take $\alpha\in[0.5,0.9]$) of the collocation points from this trajectory. We establish a “sample” time, $t_{\text{gap}}$, and sample only at integer multiples of $t_{\text{gap}}$ so that points are not too close to each other. The other $1-\alpha$ of the collocation points are chosen from a uniform distribution on the computational domain $\mathcal{R}$. We associate each sampled point $\mathbf{x}_{i}$ with the box in $\mathcal{R}$ in which it lands provided it has not yet been taken by a collocation point. The process stops when $N_{x}$ collocation points are collected. We reiterate that only points $\mathbf{x}_{i}$ are binned in the boxes, and that one should take $\alpha=0$ when generating collocation points for the backward eigenmodes. This procedure ensures that many of the collocation points are in locations where the forward eigenmodes take on non-trivial values, and ultimately reduces the number of collocation points required for accurate estimation.

A similar procedure can be used to generate $\mathcal{Y}$, although these points are not binned in boxes because they are used to evaluate the differential operator and are not associated with numerical estimates of a probability density. Algorithms SM1.1 and SM1.2 in the supplement provide more details for the generation of $\mathcal{Y}$ and $\mathcal{X}$, respectively [57].

3.2 Step 2: density estimation

The goal of this step is to estimate time-dependent transition densities $\rho(\mathbf{y},t\,|\,\mathbf{x},s)$ associated with equation (1).

To numerically approximate a time-dependent density, we fix a sampling time interval, $t_{\text{gap}}$, and the number of sample times, $N_{t}$. Realizations of the process (1) are sampled at integer multiples of $t_{\text{gap}}$ so that samples are not too close together.333The parameter $t_{\text{gap}}$ should be much smaller than the mean period of oscillation of a given system. These samples are recorded in a matrix $\mathcal{D}$, which has $N_{t}$ rows (one for each time slice) and $N_{x}$ columns (one for each collocation point, $\mathbf{x}_{i}\in\mathcal{X}$).

The procedures for estimating the densities associated with forward and backward eigenfunctions differ. To distinguish between them, we denote by $\mathcal{D}_{F}$ and $\mathcal{D}_{B}$ the matrices containing the estimated time-dependent densities for the forward and backward eigenfunctions, respectively.

We begin by estimating the forward transition densities $\rho(\mathbf{x}_{i},t_{k}\mid\mathbf{x}_{0},0)$ for a fixed initial condition $\mathbf{x}_{0}$, where $\mathbf{x}_{i}\in\mathcal{X}$ and $t_{k}\in{t_{1},\dots,t_{N_{t}}}$. Using Monte Carlo simulation, we generate $K$ realizations of equation (1) starting from $\mathbf{x}_{0}$. At each sampling time $t_{k}$, we count how many trajectories fall into the bin centered at each collocation point $\mathbf{x}_{i}$. If the bin around $\mathbf{x}_{i}$ contains $K_{i}$ samples, the estimated transition probability is given by $K_{i}/(K\delta_{i})$, where $\delta_{i}$ is the volume of the bin around $\mathbf{x}_{i}$. These estimated values are then stored in the matrix $\mathcal{D}_{F}$.

Letting $\tilde{\rho}_{f,i}(t_{k})\equiv\tilde{\rho}(\mathbf{x}_{i},t_{k}\mid\mathbf{x}_{0},0)$ denote the estimated transition density, it follows that the $i$th column of $\mathcal{D}_{F}$ provides an approximation to the time series ${\tilde{\rho}_{f,i}(t_{k})}_{k=1}^{N_{t}}$. Since the initial condition $\mathbf{x}_{0}$ is fixed, the coefficient $Q_{\lambda_{1}}^{*}(\mathbf{x}_{0})$ in the eigenfunction expansion ansatz (7) is constant. If desired, one can also sample $\mathbf{x}_{0}$ from a prescribed distribution to reduce the influence of higher-order modes. Further implementation details can be found in Algorithm SM1.3 in the Supplement.

We now turn to the computation of the backward eigenmodes. The objective in this step is to approximate the forward transition probability $\rho(\mathcal{B},t\mid\mathbf{x}_{i},0)$ for each $\mathbf{x}_{i}\in\mathcal{X}$ at a sequence of time points $\mathbf{t}=\{t_{k}\}_{k=1}^{N_{t}}$. Here, $\mathcal{B}$ denotes a reference box, $\mathcal{B}\subset\mathcal{R}$, which is chosen to contain a nontrivial portion of the invariant density.444The reference box is introduced as a practical approximation. Ideally, one would track the probability mass arriving in an individual box, but the likelihood of a trajectory landing in a single small bin is typically very low, necessitating a prohibitively large number of samples. To improve sampling efficiency, we instead consider a larger region $\mathcal{B}$. We denote the estimated transition probability into $\mathcal{B}$ at time $t_{k}$, given initial condition $\mathbf{x}_{i}$, by $\tilde{\rho}_{b,i}(t_{k})\equiv\tilde{\rho}(\mathcal{B},t_{k}\mid\mathbf{x}_{i},0)$.

To compute $\tilde{\rho}(\mathcal{B},t_{k}\mid\mathbf{x}_{i},0)$, we simulate $K$ realizations of equation (1), all initialized at $\mathbf{x}_{i}$. For each time $t_{k}$, we count the number of trajectories, say $K_{j}$, that lie within the set $\mathcal{B}$. The estimated transition probability is then given by $\tilde{\rho}_{b,i}(t_{k})=K_{j}/K$, and the result is stored in the matrix $\mathcal{D}_{B}$. In our numerical examples, approximately $N_{\text{sample}}=10^{4}$ trajectories per $\mathbf{x}_{i}$ are sufficient to achieve accurate estimates. Additional implementation details are provided in Algorithm SM1.4 of the Supplement. It is important to note that in estimating $\tilde{\rho}_{b,i}(t_{k})$, the terminal region $\mathcal{B}$ is fixed. As a result, in the eigenfunction expansion ansatz (7), the integrals $\int_{\mathcal{B}}d\mathbf{y}\;P_{\lambda_{1}}(\mathbf{y})$ and $\int_{\mathcal{B}}d\mathbf{y}\;P_{0}(\mathbf{y})$ remain constant across all initial conditions $\mathbf{x}_{i}$. This property is exploited in Step 4 (see §3.4) to recover the backward eigenfunctions over the domain $\mathcal{X}$.

For low-dimensional problems, standard grid-based methods may be used (e.g., see [31]) to sample the probability density at the points $\mathbf{x}_{i}$. For higher dimensional problems, a grid-based approach is too slow and memory intensive. Instead, we opt for the memory-efficient Monte Carlo sampling algorithm introduced in [57] (see Appendix B in that work for a complete description).

3.3 Step 3: Eigenvalue Approximation

In this section, the goal is to estimate the low-lying eigenvalue, $\lambda_{1}=\mu_{1}+i\omega_{1}$, corresponding to the $Q$-function and related forward eigenfunction. The eigenvalue is then used to approximate the corresponding eigenfunction in Step 3.4.

The eigenvalue approximation relies on one of the numerically estimated densities from Step 2. Recall that the estimated densities arising from Algorithm SM1.3 are denoted by $\tilde{\rho}_{f,i}(t_{k})$, and the estimated densities arising from Algorithm SM1.4 are denoted by $\tilde{\rho}_{b,i}(t_{k})$. As before, $\mathbf{t}=\{t_{k}\}_{k=1}^{N_{t}}$ is understood to represent the $N_{t}$-dimensional vector of the time slices. For concreteness, we present results in terms of the $i$th density, $\tilde{\rho}_{f,i}(t_{k})$.555One should choose $i$ such that $\mathbf{x}_{i}$ is near a nontrivial portion of the invariant density. It is also practical to choose a reference box centered on a region with substantial invariant density and sum together each $\tilde{\rho}_{f,i}(t_{k})$ corresponding to locations in the reference box to reduce noise. The exact choice of box will influence $b_{1}$, $b_{3}$, and $b_{5}$ in (8), but will not influence the estimation of $\lambda_{1}$ (in the limit of large data), because $\lambda_{1}$ is a global property of the system.

Because the underlying system is a stochastic oscillator, we expect $\tilde{\rho}_{f,i}(t_{k})$ to have the form of a decaying oscillation. Our goal is to estimate the decay rate (real part of $\lambda_{1}$) and frequency (imaginary part of $\lambda_{1}$) of this oscillation. To that end, we define a function

[TABLE]

with parameters $\mathbf{b}=[b_{1},b_{2},b_{3},b_{4},b_{5}]$ to be determined. This function is chosen to model both contraction and oscillatory behavior, where:

• •

$b_{1}$ and $b_{4}$ adjust the amplitude and vertical shift,

• •

$b_{5}$ captures system damping, and approximates the real part of the eigenvalue,

• •

$b_{2}$ and $b_{3}$ describe the periodicity and phase shift of the oscillation, where $2\pi/b_{2}$ approximates the imaginary part of the eigenvalue.

• •

$b_{1}$, $b_{2}$, $b_{4}$ and $-b_{5}$ are strictly positive.

Then, we define an objective function

[TABLE]

which measures the least squares error between the data, $\tilde{\rho}_{f,i}(t_{k})$, and the model. By solving (9) in the (nonlinear) least squares sense666We use MATLAB’s “fminsearch” command to perform the fitting., we obtain estimates for the parameters $\mathbf{b}$. The parameter $b_{5}$ approximates the real part of $\lambda_{1}$, while $2\pi/b_{2}$ approximates the imaginary part of $\lambda_{1}$. We perform this fitting for different collocation points, and let $\tilde{b}_{5}$ and $\tilde{b_{2}}$ represent the average output of the fitting procedure. We let $\tilde{\lambda}_{1}=\tilde{\mu}_{1}+i\tilde{\omega}_{1}=\tilde{b}_{5}+(2\pi/\tilde{b}_{2})i$ denote the least squares estimate of the true eigenvalue, $\lambda_{1}$.

Notice that the sum is not necessarily taken over all the sample points $t_{k}$. Instead, it is taken over an intermediate time interval $t_{1}\leq t_{s}<t_{f}\leq t_{N_{t}}$. The formulation (9) assumes that significant contributions to the decaying oscillation arise only from the slowest decaying eigenvalue, $\lambda_{1}$. However, in practice, even if the spectral gap (recall the robustly oscillatory criteria in §2) is large, we find that the effects of higher-order eigenmodes cannot be ignored. As will be discussed in §4, large $t_{s}$ reduces the effect of higher-order eigenmodes but increases the Monte Carlo error. Therefore, the intermediate time interval $[t_{s},t_{f}]$ should be chosen so that the contributions from the higher-order modes decay but the decayed oscillation has not been overwhelmed by noise.

In some cases, it is useful to obtain estimates for the eigenvalues of higher-order eigenmodes. If the stochastic oscillator reduces to a stable limit cycle oscillator in the vanishing noise limit, the low-lying eigenvalue spectrum may be approximately parabolic with $\lambda_{n}\approx\mu_{1}n^{2}+i\omega_{1}n$ for integer $n\in\mathbb{Z}$. This structure is consistent with the existence of a change of coordinates under which system dynamics may be approximated as a diffusion on the ring: $d\phi=\omega_{1}dt+\sqrt{-2\mu_{1}}\,dW(t)$, whose eigenvalues lie exactly on this parabola [18, 51]. In such cases the approximated eigenvalue $\tilde{\lambda}_{1}$ may be used to estimate eigenvalues corresponding to higher-order modes.

3.4 Step 4: Eigenfunction approximation

We now discuss the approximation of the $Q$-function and the corresponding forward eigenfunction. In what follows, we assume knowledge of the stationary distribution, $P_{0}$, as the stationary distribution can be accurately computed by using methods in [11] or [57]. We also require the eigenfunction expansion ansatz (7), which we reproduce here for convenience (setting $s=0<t$)

[TABLE]

Forward eigenfunction

To numerically approximate the forward eigenfunction, we consider the ansatz

[TABLE]

where we replace the continuous terminal variable $\mathbf{y}$ with the collocation points: $\mathbf{y}\to\mathbf{x}_{i}\in\mathcal{X}$, replace the continuous initial variable $\mathbf{x}$ with points sampled from the starting box: $\mathbf{x}\to\mathbf{x}_{0}$, and replace continuous time with discrete time slices: $t\to t_{j}$. Here, $\{t_{j}\}_{j=T_{s}}^{T_{f}}\subset\{t_{k}\}_{k=1}^{N_{t}}$ is a subset of all the time-slices with $t_{1}\leq t_{T_{s}}<t_{T_{f}}\leq t_{N_{t}}$. We also replace the unknown true transition density with the known numerically approximated density, $\tilde{\rho}_{f,i}(t_{j})$, and the unknown true eigenvalue with the approximated eigenvalue $\tilde{\lambda}_{1}$. Notice that, as in Step 3 concerning the approximation of the eigenvalue $\lambda_{1}$, we consider only a subset of the possible time-slices due to contamination from higher-order modes.

Note that in (11), only the quantities $P_{\lambda_{1}}(\mathbf{x}_{i})$ and $Q_{\lambda_{1}}^{*}(\mathbf{x}_{0})$ are unknown. Our goal is to approximate $P_{\lambda_{1}}(\mathbf{x}_{i})$ for each $\mathbf{x}_{i}\in\mathcal{X}$. To that end, recall that $\tilde{\rho}_{f,i}(t_{k})\equiv\tilde{\rho}(\mathbf{x}_{i},t_{k}\>|\>\mathbf{x}_{0},0)$ where $\mathbf{x}_{0}$ is a fixed initial condition. Consequently, the quantity $Q_{\lambda_{1}}^{*}(\mathbf{x}_{0})$ remains constant across all indices $i$, and does not depend on the terminal point $\mathbf{x}_{i}$. We absorb this constant factor into the definition of a rescaled eigenfunction

[TABLE]

With this notation, we rewrite (11) as

[TABLE]

We find it convenient to introduce the notation $\tilde{P}_{\lambda_{1}}=P^{R}_{\lambda_{1}}+iP^{I}_{\lambda_{1}}$ to express (13) as

[TABLE]

which is a linear system $\mathbf{A}_{f}\mathbf{u}_{f}=\mathbf{b}_{f}$, with

[TABLE]

and

[TABLE]

Equation (14) is a linear least squares problem, and should be solved separately for each $\mathbf{x}_{i}\in\mathcal{X}$.777We use MATLAB’s “backslash” to perform this computation. The solution of (14) gives an approximation of the rescaled forward eigenfunction, $\tilde{P}_{\lambda_{1}}(\mathbf{x}_{i})$ for $\mathbf{x}_{i}\in\mathcal{X}$.

We note that the rescaled function $\tilde{P}_{\lambda_{1}}$ does not, in general, satisfy a biorthogonality condition as it absorbs the unknown scalar factor $2Q_{\lambda_{1}}^{*}(\mathbf{x}_{0})$.888If $Q_{\lambda_{1}}^{*}(\mathbf{x}_{0})=0$, then our estimation procedure fails. However, the $Q$-function for a robustly oscillatory stochastic oscillator typically has vanishing magnitude only at isolated points, namely at “phaseless points” that tend to lie close to deterministic phaseless points (fixed points of the deterministic system). One approach to ensure that $Q_{\lambda_{1}}^{*}(\mathbf{x}_{0})\neq 0$ is to choose $\mathbf{x}_{0}$ near a high density region of the invariant distribution. As a result, $\tilde{P}_{\lambda_{1}}$ approximates the true forward eigenfunction $P_{\lambda_{1}}$ only up to a complex constant. If normalization is required, e.g., to enforce biorthogonality, this scalar can be estimated by computing the inner product $\langle\tilde{P}_{\lambda_{1}},Q_{\lambda_{1}}^{*}\rangle$, provided that an approximation of $Q_{\lambda_{1}}^{*}$ is available. The function $\tilde{P}_{\lambda_{1}}$ can then be rescaled accordingly to recover a biorthonormal representative of the forward eigenfunction.

Backward eigenfunction

The computational procedure to approximate the backward eigenfunction is similar to that of the forward case. We consider the ansatz

[TABLE]

where we integrate the continuous terminal variable $\mathbf{y}$ over the reference box $\mathcal{B}$, replace the continuous initial variable $\mathbf{x}$ with the starting boxes: $\mathbf{x}\to\mathbf{x}_{i}\in\mathcal{X}$, and replace continuous time with discrete time-slices: $t\to t_{j}$ with $t_{j}$ as above. We also replace the true unknown density with the known numerically approximated density, $\tilde{\rho}_{b,i}(t_{j})$, and the unknown true eigenvalue with the approximated eigenvalue $\tilde{\lambda}_{1}$.

Note that in (17), only the quantities $\int_{\mathcal{B}}d\mathbf{y}\;P_{\lambda_{1}}(\mathbf{y})$ and $Q_{\lambda_{1}}^{*}(\mathbf{x}_{i})$ are unknown. Our goal is to approximate $Q_{\lambda_{1}}^{*}(\mathbf{x}_{i})$ for each $\mathbf{x}_{i}\in\mathcal{X}$. To that end, recall that $\tilde{\rho}_{b,i}(t_{k})\equiv\tilde{\rho}(\mathcal{B},t_{k}\>|\>\mathbf{x}_{i},0)$ where $\mathcal{B}$ is a fixed reference box for each initial condition $\mathbf{x}_{i}$, and is used as an approximation of a terminal box near the invariant distribution. Consequently, the quantity $\int_{\mathcal{B}}d\mathbf{y}\;P_{\lambda_{1}}(\mathbf{y})$ remains constant across all indices $i$, and does not depend on the starting point $\mathbf{x}_{i}$. We absorb this constant factor into the definition of a rescaled eigenfunction

[TABLE]

With this notation, we rewrite (17) as

[TABLE]

We find it convenient to introduce the notation $\tilde{Q}_{\lambda_{1}}^{*}=Q^{R}_{\lambda_{1}}+iQ^{I}_{\lambda_{1}}$ to express (19) as

[TABLE]

which is a linear system $\mathbf{A}_{b}\mathbf{u}_{b}=\mathbf{b}_{b}$, with

[TABLE]

and

[TABLE]

We use different notation for $\mathbf{A}_{f}$ and $\mathbf{A}_{b}$ to highlight that different intermediate time intervals may be used for the approximation of the forward and backward eigenfunctions, respectively. Equation (20) is a linear least squares problem, and should be solved separately for each $\mathbf{x}_{i}\in\mathcal{X}$. The solution of (20) gives an approximation of the rescaled backward eigenfunction, $\tilde{Q}_{\lambda_{1}}^{*}(\mathbf{x}_{i})$ for $\mathbf{x}_{i}\in\mathcal{X}$.

We note that, as in the case of the forward eigenfunction, the rescaled function $\tilde{Q}_{\lambda_{1}}^{*}$ does not, in general, satisfy a biorthogonality condition as it absorbs the unknown scalar factor $2\int_{\mathcal{B}}d\mathbf{z}\;P_{\lambda_{1}}(\mathbf{z})$.999Note that this constant factor is non-vanishing provided that the reference box overlaps a non-vanishing region of the invariant distribution. As a result, $\tilde{Q}_{\lambda_{1}}^{*}$ approximates the true backward eigenfunction $Q_{\lambda_{1}}^{*}$ only up to a complex constant. If needed, appropriate normalization can be enforced in a manner similar to that described for the case of the forward eigenfunction.

Computation of higher order eigenmodes

In some cases where the spectral gap is small, including more than one eigenmode in the least squares formulation can improve accuracy and allow for the approximation of higher-order eigenmodes, even if only the first eigenfunction is used in the artificial neural network training. We refer §4 for more analysis on how the linear regression error depends on higher-order modes. To account for higher-order eigenmodes, we expand (10) to include additional terms. We now list the augmented forward and backward least squares problems that incorporate $M$ eigenmodes:

[TABLE]

and

[TABLE]

respectively. We numerically investigate the impact of including higher order modes in §5. Note that our theory holds if a higher-order eigenmode and corresponding eigenvalue are purely real, as the imaginary parts are simply zero in the above expressions.

3.5 Step 5: ANN implementation

The approximations of the eigenfunctions from step 4 are only available at collocation points $\mathbf{x}_{i}\in\mathcal{X}$. Moreover, because these approximations rely on Monte Carlo simulations of the underlying probability density, the approximations are often noisy, especially in high dimensions. Here, our goal is to use an ANN to generate eigenfunction approximations at any point of interest in $\mathcal{R}$, and to reduce noise.

Building on the approach pioneered by [57], we consider an objective function corresponding to operator $\mathscr{L}$ with eigenpair $(\lambda_{n},U_{\lambda_{n}})$:

[TABLE]

Here, the generic operator $\mathscr{L}$ may play the role of either the forward operator, $\mathcal{L}$, or the backward operator, $\mathcal{L}^{\dagger}$. Accordingly, the eigenfunction $U_{\lambda_{n}}$ represents either a forward eigenfunction, $P_{\lambda_{n}}$, or a backward eigenfunction, $Q^{*}_{\lambda_{n}}$. The notation $\tilde{U}_{\lambda_{n}}$ is understood to represent a least-squares approximation of either the forward or backward eigenfunction from Step 4.

The first term in (25) ensures that $U_{\lambda_{n}}$ is a solution of the differential operator, i.e., one of the infinitely many eigenfunctions which are solutions to the eigenvalue problem without boundary conditions. The last term circumvents the need for boundary conditions by enforcing a particular normalization via reference data.

Our goal is to use an ANN to minimize the first term of the objective function (25) over the reference points, $\mathcal{Y}$, and to minimize the last term of the objective function over the training set, $\mathcal{X}$. To that end, consider an ANN with parameter set $\chi$ (further details are provided in the supplement) with two outputs, $\mathcal{N}^{R}(\mathbf{x};\chi)$ and $\mathcal{N}^{I}(\mathbf{x};\chi)$, corresponding to the real and imaginary parts of an eigenfunction, respectively. To approximate the eigenpair $(\lambda_{n},U_{\lambda_{n}})$, we introduce the loss function $J_{\mathscr{L}}(\chi)\equiv J_{\mathscr{L}}(\tilde{\lambda}_{n},\tilde{U}_{\lambda_{n}}(\mathcal{X}),\mathcal{Y};\chi)$:101010We break the loss function into real and imaginary parts because TensorFlow 2 cannot work directly with complex numbers.

[TABLE]

Minimization of the loss function $J_{\mathscr{L}}(\chi)$ with respect to $\chi$ generates parameters, $\chi^{*}$, for which the network solution approximately satisfies

[TABLE]

As noted in [57], our setup is similar to the so-called Physics-Informed Neural Network (PINN) [10], in which the first term over $\mathcal{Y}$ in the loss function (broken into real and imaginary parts) learns the “physics” of the problem via the differential operator, and the last term over $\mathcal{X}$ (broken into real and imaginary parts) is similar to boundary conditions.

The neural network training follows the “Alternating Adam” approach used in [57]. In each iteration, we shuffle two small batches from $\mathcal{X}$ and $\mathcal{Y}$ respectively, and update the parameters $\chi$ using the Adam optimizer. This is because the loss terms with the differential operators could have very different magnitudes than the loss terms corresponding to the reference data. The Adam optimizer is also almost invariant to rescaling [25].

4 Error decomposition of eigenfunction approximation

In this section we will study the error decomposition when estimating the forward eigenfunctions $P_{\lambda}(\mathbf{x})$ by solving a least squares problem, where $\mathbf{x}$ is a given collocation point. The case of the backward eigenfunctions is analogous. We start from equations (23) and (24) by assuming that the first $M$ eigenvalues denoted by $\lambda_{1}=\mu_{1}+i\,\omega_{1},\cdots,\lambda_{M}=\mu_{M}+i\,\omega_{N}$ (along with their complex conjugates), are known. This assumption helps us to focus on the effect of high-order modes, which is the main source of the error term in the eigenfunction. The first $M$ known eigenvalues can either come from the nonlinear regression method introduced in §3.3 or other approaches like the coupling method [12], EDMD [4, 5, 8, 27, 37, 55], or the Matrix Pencil Method [19, 30, 45].

Let $\{t_{j}\}_{j=1}^{N}$ be the sampling times with $t_{1}\equiv t_{T_{s}}$ and $t_{N}\equiv t_{T_{f}}$ as in §3.4. With the first $M$ eigenvalues known, the eigenfunction estimation problem becomes a linear least square problem for each $\mathbf{x}_{i}\in\mathcal{X}$

[TABLE]

for

[TABLE]

Define $\beta_{0}$ as the estimated value of the ground truth

[TABLE]

Assume further that $\mathscr{Y}$ has a decomposition

[TABLE]

where $\mathscr{Y}_{h}$ is the high-mode term satisfying

[TABLE]

Here, $C_{1}$ is an $\mathcal{O}(1)$ constant, $\phi$ is a phase-shift introduced so that $\mathscr{Y}_{m}$ remains real-valued, $\hat{\lambda}^{\pm}=\hat{\mu}\pm i\,\hat{\omega}$ is a known high-order eigenvalue satisfying $|\hat{\mu}|>|\mu_{i}|$ and $|\hat{\omega}|>|\omega_{i}|$ for each $i=1,\cdots,M$, and $\mathscr{Y}_{m}$ is an error term from the Monte Carlo simulation.

Let $\delta_{i}$ be the volume of the bin around $\mathbf{x}_{i}$.111111We assume uniform volume about each sample point $\mathbf{x}_{i}$: $\delta\equiv\delta_{i}$ for each $\mathbf{x}_{i}\in\mathcal{X}$. Assume $K$ samples are drawn from $X_{t_{j}}$ for each $j=1,\cdots,N$. Then the Monte Carlo estimate of $\rho_{f,i}(t_{j})$ is

[TABLE]

where $B(K,\rho_{f,i}(t_{j})\delta_{i})$ is a binomial random variable. If $K$ is sufficiently large, $\tilde{\rho}_{f,i}(t_{j})$ may be approximated by a normal random variable with mean $\rho_{f,i}(t_{j})$ and variance $\rho_{f,i}(t_{j})(1-\rho_{f,i}(t_{j}))/(K\delta_{i})$. Therefore, the $j$-th entry of $\mathscr{Y}_{m}$ is approximated by

[TABLE]

where each $\varepsilon_{j}$ is a standard normal random variable. Since each $\rho_{f,i}$ are $\mathcal{O}(1)$ quantities, we can write $\mathscr{Y}_{m}=D{\mathbf{\varepsilon}}$ for an $\mathcal{O}(1)$ diagonal matrix $D$ whose $j$-th entry is $\frac{\rho_{f,i}(t_{j})(1-\rho_{f,i}(t_{j}))}{K\delta_{i}}$.

Lemma 4.1**.**

Assume the real parts of $M$ known eigenvalues satisfy $0>\mu_{1}>\cdots\mu_{M}$ and $\hat{\lambda}^{\pm}=\hat{\mu}\pm i\,\hat{\omega}$ is a known high-order eigenvalue satisfying $|\hat{\mu}|>|\mu_{i}|$ and $|\hat{\omega}|>|\omega_{i}|$ for each $i=1,\cdots,M$. Let the time slices $\{t_{j}\}_{j=1}^{N}$ be an equipartition of $[T_{s},T_{f}]$ for $N$ sufficiently large. If Monte Carlo error terms $\{\varepsilon_{i}\}_{i=1}^{N}$ are i.i.d. standard normal random variables, then we have

[TABLE]

A detailed proof of Lemma 4.1 is provided in the supplement.

Lemma 4.1 shows that it is beneficial to include multiple eigenvalues into the linear regression, even if we use the neural network to solve for only the first eigenfunction. In addition, the choice of $T_{s}$ and $N$ should be properly balanced. Larger $T_{s}$ can reduce the error caused by higher modes, but increase the error caused by the Monte Carlo noise. The Monte Carlo error can be reduced by using more time slices. However, one should note that Lemma 4.1 assumes independent error terms at distinct $t_{j}$. In practice, the dependency among $\{\varepsilon_{j}\}_{j=1}^{N}$ cannot be ignored if the sampling gap $|t_{j+1}-t_{j}|$ is too small.

5 Examples

In this section, we apply our machine learning approach to approximate the SKO eigenfunctions and eigenvalues of several model systems. More examples and error analysis are provided in the supplement. When possible, we compare results with solutions generated by finite difference methods (see the supplement for more details).

5.1 2D Noisy Stuart-Landau Oscillator

We first study a noisy 2D Stuart-Landau oscillator

[TABLE]

where we take $\omega=2$ and $D=0.09473$. For the ANN approach, we consider a numerical domain $\mathcal{R}=[-2,2]\times[-2,2]$ discretized into $N^{2}$ boxes, with $N=200$. The leading SKO eigenvalue (computed via finite differences) is $\lambda_{1}=-0.1+2i$. Note that the stationary distribution, i.e., the eigenfunction $P_{0}$, is known exactly

[TABLE]

We first consider the approximation of the forward eigenmode $P_{\lambda_{1}}$ using our ANN approach. We run 400,000 trajectories and sample the forward density at $N_{t}=5000$ time-slices with timestep $dt=0.001$ and $t_{\text{gap}}=10$. Each realization has an initial condition drawn uniformly from the box containing the point $(1,0)$. The reference box for the eigenvalue approximation is $\mathcal{B}=[-2,0]\times[-2,0]$.

In the top row of Figure 1, we plot the decay of the forward time-dependent density in the reference box. The decaying density is fit to time-slices $t_{k}$ with $k\in[1000,5000]$. Our fitting returns an estimated eigenvalue $\tilde{\lambda}_{1}=-0.0996+2.0001i$. Both the real and imaginary parts have less than 1% relative error compared to the FD solution. We display the output of the ANN in the top row of Figure 2. We used $|\mathcal{Y}|=$500,000 reference points, all 40,000 training points $\mathcal{X}$, and trained the ANN over 30 epochs.

To compute the $Q$-function of the 2D Stuart-Landau model (36), we run 20,000 trajectories with initial conditions drawn uniformly on each of the $N^{2}=$ 40,000 boxes and sample the backward density at $N_{t}=200$ time-slices with timestep $dt=0.005$ and $t_{gap}=20$. The reference box is $\mathcal{B}=[-2,0]\times[-2,0]$.

In the bottow row of Figure 1, we plot the decay of (one of) the backward time-dependent densities in the reference box. The decaying density is fit to time-slices $t_{k}$ for $k\in[100,200]$. Our fitting returns an estimate eigenvalue $\tilde{\lambda}_{1}=-0.0992+1.9997i$. Once again, both the real and imaginary parts have less than 1% relative error compared to the FD solution. We also show the output of the ANN in the bottom row of Figure 2. We used $|\mathcal{Y}|=$500,000 reference points, all 40,000 training points $\mathcal{X}$, and trained the ANN over 30 epochs. Details on the error associated with the approximation of the forward and backward eigenfunctions may be found in the supplement.

5.2 4D Noisy Coupled Morris-Lecar Neurons

Here, we study a system of two noisy Morris-Lecar neurons coupled by a gap junction

[TABLE]

where functions and parameter values for the noiseless system are provided in the supplement, and $\iota=20$. The noise parameters are $D_{v_{1}}=15$, $D_{v_{2}}=50$, $\epsilon_{1}=0.3$, and $\epsilon_{2}=1$. We increase the noise of the second oscillator so that its isolated eigenvalues have more negative real parts, i.e., so that the dominant SKO eigenvalues of (38) correspond only to the first oscillator. We set $\kappa=0$ so that the oscillators are uncoupled, i.e., so that we can project a 4D solution onto a 2D plane and compare with known solutions generated by finite differences. For reference, the first two eigenvalues (computed with finite-differences) are $\lambda_{1}=-0.0879+1.9518i$ and $\lambda_{2}=-0.2521+3.9602i$.

Here, we focus only on the computation of the $Q$-function. To generate densities, we discretize the domain $\mathcal{R}=[-90,90]\times[0,1]\times[-90,90]\times[0,1]$ into $N^{4}$ boxes, with $N=100$. We store only 10,000 of the boxes in memory, chosen randomly and uniformly throughout $\mathcal{R}$. We take 30,000 initial conditions in each box, and keep track of the corresponding densities at $N_{t}=350$ time-slices. We set $t_{\text{gap}}=20$ and $dt=0.005$. In the top row of Figure 3, we depict the decaying probability density corresponding to one of the boxes. Our fitting procedure estimates that $\tilde{\lambda}_{1}=-0.0871+1.9573i$ taking time-slices $t_{k}$ with $k\in[125,350]$, which is in excellent agreement with the true value.

To generate approximations of the $Q$-function, we solve the least squares problem (20) and include the first two non-trivial eigenmodes in the expansion. We solve only with time-slices $t_{k}$ for $k\in[125,350]$. Once this data is generated, we train an ANN on the 10,000 collocation points $\mathcal{X}$, and $|\mathcal{Y}|=$1,000,000 reference points. The ANN was trained over 200 epochs, and iterations were divided into 128 batches, which were optimized separately to improve speed. In this case, the hidden layers of the ANN consist of only 12 nodes. We rescale voltage $\overline{v_{j}}=V_{j}/100$ for $j=1,2$ to be on the same order of magnitude as the gating variables $N_{j}$ when training the ANN. The results are shown in the bottom row of Figure 3.

5.3 Lorenz System with Additive Noise

As our final example, we consider a 3D Lorenz system in the chaotic regime that is subject to additive noise

[TABLE]

where $\sigma=10$, $\beta=8/3$, $\rho=28$, and $D=5$.

We are not aware of any solution methods for the SKO eigenfunctions of (39), apart from our machine learning approach.121212See these articles for a Koopman perspective on the deterministic Lorenz system [5, 24, 28, 39, 47]. The supplement contains information regarding the approximation of the stationary distribution of (39) when $D\neq 0$.

We now consider the approximation of the forward eigenmode, $P_{\lambda_{1}}$, using our ANN approach. We run 10,000,000 trajectories and sample the forward density in 20,000 boxes at $N_{t}=2000$ time-slices with timestep $dt=0.001$ and $t_{gap}=10$. Each realization is initialized with initial condition drawn uniformly from the box that contains the point $(10,10,10)$. The reference box is $\mathcal{B}=[-30,0]\times[-30,0]\times[-30,20]$. In the top row of Figure 4, we show the decay of the forward time-dependent density in the reference box. The decaying density is fit only for time-slices $t_{k}$ for $k\in[500,1700]$. Our fitting returns an estimated eigenvalue $\tilde{\lambda}_{1}=-0.4107+7.6684i$. We find the least squares approximation for the forward eigenfunction via (14), including only one eigenmode in the computation. We chose the same time interval as for the eigenvalue approximation, $k\in[500,1700]$. We also train an ANN, and show its output in the top row of Figure 5. The ANN was trained over 30 epochs. We used $|\mathcal{Y}|=$500,000 reference points, and all 20,000 training points $\mathcal{X}$.

Finally, we consider the $Q$-function of the noisy Lorenz system (39). We store 20,000 boxes in the domain $\mathcal{R}$. We take 20,000 initial conditions in each box, and keep track of the corresponding densities at $N_{t}=900$ time-slices. We set $t_{\text{gap}}=100$ and $dt=0.0001$. In the bottom row of Figure 4, we depict the decaying probability density corresponding to one of the boxes. Our fitting procedure estimates that $\tilde{\lambda}_{1}\approx-0.4145+7.6709i$ taking time-slices $t_{k}$ with $k\in[300,600]$.131313This estimate of $\lambda_{1}$ agrees to two decimal places with the approximation of the forward eigenmodes. We find the least squares approximation for the $Q$-function via (20), including only one eigenmode in the computation. We chose the same time interval as for the eigenvalue approximation, $k\in[300,600]$. We train an ANN over 40 epochs, with $|\mathcal{Y}|=$500,000 reference points, and all 20,000 training points $\mathcal{X}$. Iterations were broken up into 128 batches, which were optimized separately to improve speed. We interpolate both the least squares and ANN solutions over the “butterfly attractor”, i.e., the high-density region of the stationary distribution. The results are shown in the middle row of Figure 5. The bottom row of Figure 5 shows the stochastic asymptotic phase $\Psi(\mathbf{x})=\arg{Q_{\lambda_{1}}^{*}(\mathbf{x})}$.

6 Discussion

In this paper, we developed a data-driven machine learning framework for computing SKO and Fokker-Planck eigenfunctions. Our approach builds on prior work on data-driven Fokker-Planck solvers [31, 57], which leverage Monte Carlo simulation data to guide either finite difference schemes or neural network training. Unlike the Fokker-Planck equation, whose solution can be directly approximated using Monte Carlo data, the eigenfunctions cannot be sampled in this way. To address this, we extract the forward (resp. backward) eigenfunctions via nonlinear regression on time series data of the probability density function (resp. the expectation of test functions). A central challenge in neural network-based eigenfunction solvers is avoiding convergence to the trivial solution. While earlier methods often employ quadrature-based normalization techniques [16, 23], our regression-based approach directly estimates eigenfunction values at numerous collocation points, thereby accelerating convergence and improving solution quality. The method is flexible, requiring little problem-specific adjustments, and is both mesh-free and naturally suited to high-dimensional settings.

The framework presented here also extends existing theory pertaining to stochastic oscillators. Previous works [40, 41, 51] numerically studied the eigenfunctions of planar stochastic oscillators using finite differences, yet were not able to fully analyze systems in dimensions $n\geq 3$. For the first time, the theory developed in this paper allow for direct numerical investigation of higher-dimensional stochastic oscillators and systems of coupled stochastic oscillators, which has previously not been numerically feasible. A significant contribution of this paper is to show that upon the introduction of stochastic forcing, the 3D Lorenz system becomes a stochastic oscillator (in the sense introduced in [51]). The probability density exhibits decaying oscillations, and the slowest decaying SKO eigenmode is complex-valued. The complex argument of this eigenmode gives a well-defined stochastic asymptotic phase for the system. In contrast, the deterministic Lorenz system exhibits chaotic dynamics. It does not admit a periodic limit-cycle solution, and therefore does not have a well-defined phase reduction.141414However, an “optimal” phase description is possible for deterministic chaotic oscillators, including the Lorenz system [47].

In the deterministic case, systems with a spiral sink do not admit periodic-limit cycle solutions. Nevertheless, when perturbed by noise, the system behaves as a stochastic oscillator with a well-defined asymptotic stochastic phase [52]. Similarly, deterministic systems which admit heteroclinic cycles do not have periodic limit-cycle solutions and do not admit a well-defined phase reduction. But when perturbed by noise, such a system does exhibit a well-defined $Q$-phase reduction [51]. As far as we are aware, no previous work has shown that a deterministic chaotic system, when perturbed by noise, exhibits a well-defined $Q$-phase reduction. Here, for the first time, we demonstrate numerically that the Lorenz system (39) does have a well-defined $Q$-phase reduction when perturbed by additive noise. This observation supports the description of the $Q$-function as a universal description of stochastic oscillators in that it is agnostic to the underlying mechanism of oscillation [40].

Our machine learning tool allows for the study of high-dimensional systems of coupled stochastic oscillators. Several directions can be explored to further improve training quality. First, our sampling technique could be improved by, for example, using a forward-reverse sampler developed in [38], which can provide better estimate of the probability density function in high-dimensional setting. Second, our approach could benefit from integrating other data-driven techniques related to the Koopman operator. While extended dynamic mode decomposition (EDMD) has been used to compute eigenfunctions of the Fokker-Planck operator [56, 58], its performance depends strongly on the choice of basis functions, or “dictionary”. One promising direction is to incorporate EDMD into neural network training by using the approximated Koopman operator to provide additional regularization.

On the application side, our future work may focus on extending the current framework to compute a family of $Q$-functions and corresponding eigenvalues as a function of coupling strength. Another compelling direction is to adapt the ANN-based framework for systems with partial observations, enabling the computation of eigenfunctions for partially observed oscillators. Computing SKO eigenfunctions is also of significant interest to experimentalists, especially in scenarios where the underlying model is unknown. The framework developed here could be further modified to infer underlying dynamics directly from time-series data, making it applicable to a wide range of real-world problems.

Acknowledgments

This work was supported in part by NSF grant DMS-2052109 to PJT and DMS-2108628 to YL. This material is also based upon work supported in part by the National Science Foundation under Grant No. DMS-1929284 while the authors were in residence at the Institute for Computational and Experimental Research in Mathematics in Providence, RI, during the “Math + Neuroscience: Strengthening the Interplay Between Theory and Mathematics” program. This work was also supported in part by the Oberlin College Department of Mathematics.