Neural Spline Operators for Risk Quantification in Stochastic Systems

TL;DR

This paper introduces Neural Spline Operators, a physics-informed neural operator framework that efficiently learns risk probabilities in complex stochastic systems with functional dynamics, outperforming existing methods in accuracy and speed.

Contribution

The paper presents Neural Spline Operators (NeSO), a novel PINO framework using B-splines for risk quantification in systems with functional dynamics, with theoretical guarantees and practical advantages.

Findings

NeSO achieves better initial and boundary condition enforcement.

NeSO demonstrates significant online speed-up over existing methods.

Theoretical proof of universal approximation capability of NeSO.

Abstract

Accurately quantifying long-term risk probabilities in diverse stochastic systems is essential for safety-critical control. However, existing sampling-based and partial differential equation (PDE)-based methods often struggle to handle complex varying dynamics. Physics-informed neural networks learn surrogate mappings for risk probabilities from varying system parameters of fixed and finite dimensions, yet can not account for functional variations in system dynamics. To address these challenges, we introduce physics-informed neural operator (PINO) methods to risk quantification problems, to learn mappings from varying \textit{functional} system dynamics to corresponding risk probabilities. Specifically, we propose Neural Spline Operators (NeSO), a PINO framework that leverages B-spline representations to improve training efficiency and achieve better initial and boundary condition enforcements, which are crucial for accurate risk quantification. We provide theoretical analysis demonstrating the universal approximation capability of NeSO. We also present two case studies, one with varying functional dynamics and another with high-dimensional multi-agent dynamics, to demonstrate the efficacy of NeSO and its significant online speed-up over existing methods. The proposed framework and the accompanying universal approximation theorem are expected to be beneficial for other control or PDE-related problems beyond risk quantification.

Tables

Table 1. TABLE I : Test performance for recovery probability estimation over 10 test cases. Errors are reported as mean ± \pm standard deviation.

Metric	NeSO (Proposed)	FNO
MSE ($\times {10}^{-3}$)	$0.285\pm 0.190$	$5.298\pm 9.296$
MAE ($\times {10}^{-2}$)	$0.892\pm 0.380$	$2.827\pm 2.968$
Rel. Error ($\times {10}^{-2}$)	$5.697\pm 1.888$	$14.994\pm 10.093$
Train Time (s)	$756$	$1086$

Table 2. TABLE II : Comparison of online computation times.

Method	MC	PDE	NeSO (proposed)
Time (s)	4618.48	2079.71	56.47

Equations

$$d{x}_t=f(x_t,t)dt+\sigma d{W}_t,$$

$${\mathcal{C}}_{\alpha }=\left\{x\in {\mathbb{R}}^n:{\varphi }_{\alpha }(x)\ge 0\right\},$$

$$F(x,t)=\mathbb{P}(x_{\tau }\in {\mathcal{C}}_{\alpha },\forall \tau \in [0,t]\mid x_0=x).$$

$$\frac{dF(x,t)}{dt}=f(x,t)\frac{dF(x,t)}{dx}+\frac{1}{2}{\sigma }^{\top }\sigma \mathrm{tr}\left[\frac{d^{2}F(x,t)}{d{x}^2}\right],$$

$$F(x,0)=1,\forall x\in {\mathcal{C}}_{\alpha },F(x,t)=0,\forall x\in \partial {\mathcal{C}}_{\alpha },\forall t>0.$$

$$G_{\theta }(h):\mathcal{H}({\mathbb{R}}^m)\to \mathcal{H}({\mathbb{R}}^l),$$

$$\int k(x,y)f(x)dx+b(y),$$

$$B_{i,d}(x)=\frac{x-{\hat{x}}_i}{{\hat{x}}_{i+d}-{\hat{x}}_i}{B}_{i,d-1}(x)+\frac{{\hat{x}}_{i+d+1}-x}{{\hat{x}}_{i+d+1}-{\hat{x}}_{i+1}}{B}_{i+1,d-1}(x),$$

$$B_{i,0}(x)=\{\begin{array}{ll}1,& {\hat{x}}_i\le x<{\hat{x}}_{i+1},\\ 0,& \text{otherwise}.\end{array}$$

$$c:=[c_1,c_2,\dots ,c_{\ell }],$$

$$B_d(x):=[B_{1,d}(x),B_{2,d}(x),\dots ,B_{\ell ,d}(x){]}^{\top }.$$

$$\hat{F}(x,t)=c{B}_d(x).$$

$$\hat{F}(x_1,$$

$$\sum _{i_t=1}^{{\ell }_t}$$

$$W_{\hat{F}}(x,t,f)=$$

$$-\frac{1}{2}{\sigma }^{\top }\sigma \mathrm{tr}\left[\frac{d^2\hat{F}(x,t)}{d{x}^2}\right].$$

$${\mathcal{L}}_p=\frac{1}{|\mathcal{P}|}\sum _{x\in \mathcal{P}}|W_{\hat{F}}(x,t,f){|}^2,$$

$${\mathcal{L}}_d=\frac{1}{|\mathcal{D}|}\sum _{(x,t)\in \mathcal{D}}|F(x,t)-\hat{F}(x,t){|}^2,$$

$$\mathcal{L}=w_p{\mathcal{L}}_p+w_d{\mathcal{L}}_d,$$

$$\frac{d^p}{d{x}^p}{B}_{i,d}(x)=$$

$$\frac{(d-1)!}{(d-p-1)!}\sum _{k=0}^p(-1{)}^k(\genfrac{}{}{0ex}{}{p}{k})\frac{B_{i+k,d-p}(x)}{{\prod }_{j=0}^{p-1}\left({\hat{x}}_{i+d-j-1}-{\hat{x}}_{i+k}\right)}.$$

$$\displaystyle\big{|}F(x,t)-G_{\theta}(f,\alpha)(x,t)\big{|}\leq\epsilon,$$

$$G_{\theta }(f,\alpha )(x,t)=G_{\theta }(f,\alpha )\cdot B_{\ell ,d}(x,t),$$

$$f(x,t)=A_1\sin \left(2\pi \frac{{\omega }_1}{10}t+{\psi }_1\right)+A_2\sin \left(2\pi \frac{{\omega }_2}{10}t+{\psi }_2\right),$$

$$F(x,t):=\mathbb{P}\left(\exists \tau \in [0,t],\text{ s.t. }{x}_{\tau }\in {\mathcal{C}}_{\alpha }\mid x_0=x\right).$$

$$F(x,0)=0,\forall x\in {\mathcal{C}}_{\alpha }^c,F(x,t)=1,\forall x\in \partial {\mathcal{C}}_{\alpha },\forall t>0,$$

$$F(x,t)={\int }_0^t\frac{(4-x)}{\sqrt{2\pi {\tau }^3}}\exp \left(-\frac{((4-x)-S(\tau ){)}^2}{2\tau }\right)d\tau ,$$

$$S(t)={\int }_0^{t}f(x_{\tau },\tau )d\tau .$$

$$d{x}_t^{(k)}=\left(A{x}_t^{(k)}+B{u}_t^{(k)}\right)dt+{\sigma }_k{w}_t^{(k)},$$

$$y_t^{(k)}=C{x}_t^{(k)},\text{where }C=\left[\begin{array}{cc}1& 0\end{array}\right].$$

$$u_t^{(k)}=-\sum _{i\in \mathcal{I}(k)}{l}_{k,i}{y}_t^{(i)},$$

$$L=\left[\begin{array}{ccccccc}5& -1& -1& -1& -1& -1& 0\\ -1& 3& 0& -1& 0& 0& -1\\ -1& 0& 2& 0& -1& 0& 0\\ -1& -1& 0& 4& -1& -1& 0\\ -1& 0& -1& -1& 4& -1& 0\\ -1& 0& 0& -1& -1& 3& 0\\ 0& -1& 0& 0& 0& 0& 1\end{array}\right].$$

Full text

I Introduction

Accurately quantifying long-term risk probabilities online across a wide range of stochastic systems is crucial to designing safe control policies [1]. Estimating these probabilities is non-trivial due to the rare nature of risk events and the long time horizons that need to be accounted for, resulting in heavy computational cost for sampling-based methods [2]. Partial differential equation (PDE)-based approaches alleviate this computational burden by directly characterizing the evolution of risk probabilities of different initial states and time horizons through PDEs, but they often struggle to scale efficiently online with varying forms of system dynamics [3]. Physics-informed neural network (PINN) methods can quantify risk probabilities online for multiple systems by accepting parameters of system dynamics as input [4, 5]. However, these approaches are inherently limited as they cannot accommodate inputs with varying parameter dimensions or system dynamics of different functional forms. In this paper, we aim to efficiently and accurately estimate long-term risk probabilities (or equivalently, safety probabilities) for diverse stochastic systems characterized by varying dynamic functions and safe region descriptions. We specifically address the limitations of existing PINN-based methods and introduce an approach capable of simultaneously handling systems with different functional dynamics.

Long-term safety probabilities in stochastic systems can often be described by solutions of certain convection-diffusion equations [3]. Traditional numerical methods for solving these PDEs include grid-based and sampling-based approaches, both of which face significant challenges in accounting for varying functional forms of system dynamics. Recent PINN methods overcome some limitations by learning parameterized PDEs corresponding to safety probabilities of parameterized system dynamics. These networks are trained to approximate PDE solutions directly, enabling the prediction of safety probabilities across the entire state space and over various long-term time horizons. However, these PINNs remain restricted to fixed finite dimensional parameterized forms of system dynamics and cannot directly handle dynamics described by arbitrary functions.

To overcome these challenges, we introduce physics-informed operator learning methods [6, 7] into risk quantification, enabling direct learning of mappings from general functional system dynamics to risk probability values. While operator learning and neural operators significantly broadens applicability in many problems [8], direct application to risk quantification is insufficient due to difficulties in enforcing PDE constraints effectively, which is crucial for accuracy in safety-critical scenarios. Thus, we propose Neural Spline Operators (NeSO), a novel physics-informed neural operator (PINO) learning framework leveraging B-spline representations [9] to efficiently enforce PDE constraints during training. The spline-based representation inherently offers enhanced smoothness, inital and boundary condition enforcement, and in particular analytical derivative calculations, which is essential but previously computationally expensive for physics-informed learning [10, 7]. We also theoretically demonstrate the universal approximation capability of our proposed NeSO framework. The primary contributions of this paper are summarized as follows:

• •

Introduction of physics-informed operator learning methods to risk quantification to enable simultaneous handling of varying functional system dynamics.

• •

Development of the Neural Spline Operator (NeSO), employing a spline-based representation for efficient physics-informed operator learning, accompanied by theoretical guarantees on universal approximation.

• •

Demonstration of effectiveness and significant computational speed-ups across various systems, including those with diverse functional dynamics and high-dimensional multi-agent settings.

II Problem Formulation

We consider the following stochastic dynamical systems

[TABLE]

where $x\in\mathcal{X}\in\mathbb{R}^{n}$ is the state, $W_{t}\in\mathbb{R}^{n}$ is the standard Wiener process with $W_{0}=\mathbf{0}$, and $\sigma\in\mathbb{R}^{n}$ is the magnitude of the noise, and $f:\mathbb{R}^{n+1}\rightarrow\mathbb{R}^{n}$ is the (closed-loop) dynamics of the system. We define the safe set as

[TABLE]

where $\phi_{\alpha}:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is a function parameterized by $\alpha$. Specifically, in this paper we consider safe sets (or sub safe sets for multi-agent systems introduced in section VI-B) with $\phi_{\alpha}=-\sum_{k=1}^{n}\mathds{1}(x_{k}\notin[\alpha_{k}^{-},\alpha_{k}^{+}])$ where $\alpha_{k}^{\pm}\in\mathbb{R}$ for all $k=1,\cdots,n$ and $\mathds{1}(\cdot)$ is the indicator function, *i.e., *the safe region is characterized by the $n$-dimensional bounding box with edge values $\alpha_{k}^{\pm}$. This characterization of safe set is used in many scenarios such as lane keeping in autonomous driving [11], dynamic walking of humanoid-robot [12], etc. The risk quantification problem aims to estimate the following long-term safety probability under dynamics (1).

[TABLE]

Eq. (3) characterizes the probability of the system (1) remaining safe under time horizon $t$ starting from initial state $x$. Note that the safety probability $F$ in (3) is also dependent on the system dynamics (1) and safe set function (2), but we drop the notations for conciseness. This long-term safety probability is critical for designing safe control for stochastic systems [1], and the goal is to accurately estimate such value for all $x\in\mathcal{X}$ and different time horizon $t>0$, with varying dynamics function $f$ and safe set parameter $\alpha$.

III Related Work

III-A Risk Quantification

A standard method to estimate (3) is to use sampling-based methods such as Monte Carlo (MC), where the system dynamics (1) is simulated starting from state $x$ for time horizon $t$ with multiple independent trajectories, and the safety probability is empirically estimated via the ratio of safe trajectories. Multiple techniques such as importance sampling [13, 14], subset simulation [15], multi-level Monte Carlo [16] are further developed to reduce the sample complexity of MC, but such methods are intrinsically designed to quantify risk of single dynamics with fixed safe set. In this work, we propose methods to learn mappings of varying system dynamics and safe sets to the safety probabilities.

Another method to estimate the long-term safety probability is to solve a corresponding partial differential equation. Specifically, from [3] we know that the safety probability (3) is the solution of the following convection diffusion equation

[TABLE]

with the following initial and boundary conditions (ICBCs)

[TABLE]

When the system dynamics is fixed and simple (e.g., $f$ is a linear function), the PDE associated with the safety probability either has analytical solution or can be solved effectively with numerical solvers (*e.g., *finite element method-based solvers [17] and finite volume method-based solvers [18]). However, when the system dynamics is changing and is a general function, traditional solvers usually require special treatment to maintain stability [19, 20] and cannot solve the safety probability PDEs in efficient and generic manners. Physics-informed learning methods leverage both data samples and PDE models, and are capable of estimating the safety probability mapping from fixed or parameterized dynamics function [4]. However, how to effectively estimate the safety probability with varying and general functional system dynamics remains an open challenge. We propose the neural spline operator method to address this challenge, by solving (4) with varying functional $f$ efficiently.

III-B Neural Operators and B-Splines

Neural operators are generalizations of neural networks by learning mappings between infinite-dimensional function spaces rather than finite-dimensional Euclidean spaces [6]. Specifically, a neural operator $G_{\theta}$, parameterized by $\theta$, maps an input function $h$ from one function space to another, formally written as:

[TABLE]

where $\mathcal{H}(\mathbb{R}^{k})$ denotes an appropriate function space defined on a given domain $\mathbb{R}^{k}$, for $k\in\mathbb{N}^{+}$. The key distinction between neural operators and classical neural networks lies in their internal operations. The standard linear transformation $Wx+b$ (weighted sum of inputs $x$, with weights $W$ and biases $b$) in multilayer perceptrons in neural networks is replaced by a continuous integral operator of the form:

[TABLE]

where $f(x)$ is the input function, $k(x,y)$ is an integral kernel, and $b(y)$ is a bias function. This integral operator formulation allows neural operators to accept input and produce output as functions at arbitrary resolutions. These properties make neural operators very competitive solutions to applications where input is of multiple resolutions and cannot be parameterized, such as weather forecast [21], turbulence modeling [22] and in particular learning surrogate maps for the solution operators of PDEs under the physics-informed neural operator learning setting [7, 23]. In this work, we leverage physics-informed neural operators (PINOs) to learn solutions of PDEs corresponding to safety probabilities, and address risk quantification problems for systems with varying functional dynamics. To the best of our knowledge, this is the first work to explore physics-informed operator learning in the context of risk quantification. In addition, although there are multiple existing architectures for physics-informed neural operator learning such as Fourier neural operators (FNO) [24], graph neural operators (GNO) [25], spectral neural operators (SNO) [26] and DeepONets [27], expensive derivative calculations through automatic differentiation or finite difference methods are often required [7], and the trained models may not always comply with initial and boundary conditions (ICBCs) that the PDE needs to satisfy [28]. We develop a specific neural operator architecture based on B-splines that provides analytical derivative calculations and direct enforcement of ICBC constraints, to improve training efficiency for PINOs in risk quantification problems.

B-splines are piece-wise polynomial functions derived from slight adjustments of Bezier curves, aimed at obtaining polynomial curves that tie together smoothly [9]. B-splines are used in combination with finite element methods [29] and through the variational dual formulation [30] to solve PDEs. B-splines are used as basis functions for surface reconstruction [31], nonlinear system modeling [32], image segmentation [33], and controller design for dynamical systems [34]. B-splines with weights learned from neural networks (NNs) are used to approximate fixed or parameterized ODEs [35] and PDEs [36, 37], and specifically PDEs for safety probabilities [5], but these methods cannot account for functional inputs, thus not varying functional system dynamics in risk quantification. We leverage the advantages of B-splines and incorporate it into a physics-informed neural operator learning framework, to learn mappings of safety probabilities from functional dynamics and varying safe sets, for efficient risk quantification of diverse stochastic systems.

IV Proposed Method

In this section, we present our proposed neural spline operator (NeSO) method, which combines both data and PDE model for learning on functional spaces, to efficiently address the risk quantification problem with varying functional dynamics. Specifically, we aim to design neural operators that learn a mapping from the dynamics function $f:\mathbb{R}^{n+1}\rightarrow\mathbb{R}^{n}$ and the safe set parameter $\alpha\in\mathbb{R}^{\xi}$, to the safety probability $F:\mathbb{R}^{n+1}\rightarrow\mathbb{R}$.

IV-A B-Splines and Basis Functions

In this section, we introduce B-splines as basis functions to approximate the safety probabilities. We start by considering the one-dimensional cases. For one-dimensional state $x\in\mathbb{R}$ (assuming fixed time $t$), the B-spline basis functions are given by the Cox-de Boor recursion formula:

[TABLE]

and

[TABLE]

Here, $B_{i,d}(x)$ denotes the value of the $i$-th B-spline basis of order $d$ evaluated at $x$, and $\hat{x}_{i}\in(\hat{x}_{i})_{i=1}^{\ell+d+1}$ is a non-decreasing vector of knot points. Since a B-spline is a piece-wise polynomial function, the knot points determine in which polynomial the parameter $x$ belongs. While there are multiple ways of choosing knot points, we use $(\hat{x}_{i})_{i=1}^{\ell+d+1}$ with $\hat{x}_{1}=\hat{x}_{2}=\cdots=\hat{x}_{d+1}$ and $\hat{x}_{\ell+1}=\hat{x}_{\ell+2}=\cdots=\hat{x}_{\ell+d+1}$, and for the remaining knot points we select equispaced values. Note that for the proposed framework we define the B-spline basis functions in a normalized domain $[0,1]$, and the basis functions for any safe set domain of interest $[\alpha_{k}^{-},\alpha_{k}^{+}]$ can be obtained through rescaling. We then define the control points

[TABLE]

and the B-spline basis functions vector

[TABLE]

Then, we can approximate the safety probability $F(x,t)$ at a fixed time $t$ with

[TABLE]

Note that with our choice of knot points, we ensure the initial and final values of $\hat{F}(x)$ coincide with the initial and final control points $c_{1}$ and $c_{\ell}$. This property will be used later to directly impose initial conditions and Dirichlet boundary conditions with the neural spline operator.

More generally, for a $n$-dimensional state $x=[x_{1},\cdots,x_{n}]\in\mathbb{R}^{n}$ and time $t\in\mathbb{R}^{+}$, we can generate B-spline basis functions based on the Cox-de Boor recursion formula (8) and (9) along each dimension $x_{k}$ and $t$ with order $d_{k}$ for $x_{1},\cdots,x_{n}$ and $d_{t}$ for $t$. We use $B_{i_{k},d_{k}}(x_{k})$ to denote the B-spline basis of order $d_{k}$ for the $i_{k}$-th function of $x_{k}$, and $B_{i_{t},d_{t}}(t)$ for $i_{t}$-th function of $t$ with order $d_{t}$.

The $n+1$-dimensional control point tensor will be given by $C=[c_{i_{1},\cdots,i_{n},i_{t}}]_{\ell_{1}\times\cdots\times\ell_{n}\times l_{t}}$, where $i_{k}$ is the $k$-th state index of the control point, and $\ell_{k}$ is the number of control points along the $k$-th dimension. Similarly $i_{t}$ and $\ell_{t}$ are the index and the number of control points along $t$. The safety probability on the $n+1$-dimensional space can then be approximated with B-splines and control points via

[TABLE]

We will use $C\cdot B_{\ell,d}(x,t)$ to denote the weighted sum of the right hand side of (13) when the context is clear for the rest of the paper.

IV-B Neural Spline Operators

In this section, we introduce our proposed neural spline operator (NeSO). The overall diagram of the framework is shown in Fig. 1. The framework composites a coefficient neural functional that learns the control point tensor $C$ given the dynamics function $f$ and parameters $\alpha$, and the B-spline layers that produce basis functions $B_{\ell,d}$.

During the forward pass, the coefficient neural functional with parameter $\theta$, denoted by $G_{\theta}(f,\alpha)$, learns the control point tensor $C$ from the input dynamics function $f$ and parameters $\alpha$. The learned control points are then multiplied with the B-spline basis functions $B_{\ell,d}$ via (13) to get the approximation $\hat{F}(x,t)=G_{\theta}(f,\alpha)\cdot B_{\ell,d}(x,t)$, and we use $\hat{F}(x,t)=G_{\theta}(f,\alpha)(x,t)$ for simplicity. Note that inside the coefficient neural functional $C=G_{\theta}(f,\alpha)$, any initial condition or Dirichlet boundary condition of a PDE can be directly specified using the control point values. For example, in a 2D case where the initial condition is given by $F(x,0)=0,\forall x$, we can set the first column of the control point tensor $c_{i_{1},1}=0$ for all $i_{1}=1,\cdots,\ell_{1}$ and this will ensure that the initial condition is met.

For the backward pass, two losses in PINO training are imposed to efficiently and effectively train the neural spline operator. We define the PDE residual

[TABLE]

We first impose a physics model loss

[TABLE]

where $\mathcal{P}$ is the set of points sampled to evaluate the governing PDE. When data is available, we can additionally impose a data loss

[TABLE]

to capture the mean square error of the approximation, where $F$ is the data point for the safety probability, $\mathcal{D}$ is the data set, and $\hat{F}$ is the prediction from the neural spline operator. The total loss is given by

[TABLE]

where $w_{p}$ and $w_{d}$ are the weights for physics and data losses. Specific choices of $w_{p}$ and $w_{d}$ will be given in the case study section, and the general rule is to increase $w_{p}$ when the PDE model is accurate and increase $w_{d}$ when the data are reliable, and choose the weights so that the magnitude of the gradients matches the choice of the learning rate for training.

The coefficient neural functional composites a neural operator layer and a linear neural network layer, to achieve learning from functional spaces to control point tensors. The neural operator layer can take many architectures such as Fourier neural operators (FNO) [24], graph neural operator (GNO) [25], U-shaped neural operator [38, 39] and low-rank neural operator (LNO) [6]. We treat the safe set parameters $\alpha$ as constant functions to adapt to existing neural operator architectures [40, 41]. In the case studies in section VI, we specifically consider FNO as the neural operator layer, since it is widely accepted in the operator learning community. Due to space limits, we refer the readers to [24] for the implementation details of FNO.

Note that there are several advantages of the proposed neural spline operator.

The neural operator formulation enables mapping from dynamics functionals to safety probabilities, which covers a wide variety of system dynamics for risk quantification. 2. 2.

The B-spline functions admit analytical derivatives [42], enables efficient training with physics-informed losses (15). Specifically, the $p$-th derivative of the $d$-th ordered B-spline is given by

[TABLE] 3. 3.

Control points of B-splines allow direct assignment of ICBCs without any learning involved [5], which improves accuracy and efficiency in learning PDEs for risk quantification.

V Theoretical Analysis

In this section, we provide theoretical analysis on the universal approximation properties of the proposed neural spline operators.

We assume that the state space $\mathcal{X}$ is a compact set in $\mathbb{R}^{n}$, and the time horizon $t\in\mathcal{T}$ where $\mathcal{T}$ is a compact set in $\mathbb{R}$. Let $\Omega_{1}:=\mathcal{X}\subset\mathbb{R}^{n}$ be the space for the arguments of the input functions, and $\Omega_{2}:=\mathcal{X}\times\mathcal{T}\subset\mathbb{R}^{n}\times\mathbb{R}$ be the space for the arguments of the output functions. We assume the dynamics function $f\in\mathcal{K}_{1}\subset C(\Omega_{1};\mathbb{R}^{n})$, where $C(\Omega_{1};\mathbb{R}^{n})$ is the space of continuous functions that map $\Omega_{1}$ to $\mathbb{R}^{n}$. Since the parameter $\alpha$ is input to the neural operator as a function, with slight abuse of notation, in this section, we denote $\alpha$ as a constant function on $\mathcal{X}$ and we assume $\alpha\in\mathcal{K}_{2}\subset C(\Omega_{1};\mathbb{R}^{\xi})$. Similarly, we assume the safety probability $F\in\mathcal{K}_{3}\subset C(\Omega_{2};\mathbb{R})$. We assume that $\mathcal{K}_{1}$, $\mathcal{K}_{2}$ and $\mathcal{K}_{3}$ are continuous Sobolev spaces [43], which are subsets of Banach spaces with norm $\|\cdot\|_{C}=\max_{\Omega}\|\cdot\|$.111With slight abuse of notation, here $C$ denotes the space of continuous functions instead of control points. Continuity of the Sobolev spaces can be guaranteed if $kp>n$, where $k$ is the order of the derivative, $p$ is integrability component, and $n$ is the dimension [44, Section 5.6.3]. We further assume the mapping between $(f,\alpha)\rightarrow F$ is continuous.222Such continuity holds with regularity conditions on $f$ and $\mathcal{C}$ [45].

Theorem 1**.**

Given dynamical system (1) and safe set (2), for any $\epsilon>0$, there exists the number of control points $\ell_{1},\cdots,\ell_{n},\ell_{t}\in\mathbb{N}^{+}$, order of B-spline basis $d_{1},\cdots,d_{n},d_{t}\in\mathbb{N}^{+}$, and a parameter set $\theta$ for the coefficient neural functional $G$, such that for any function $f\in\mathcal{F}$ and safe set parameter $\alpha\in\mathbb{R}^{\xi}$, and any $x\in\mathcal{X}$ and $t\in\mathcal{T}$, and the corresponding safety probability $F(x,t)$ defined in (3), the following holds

[TABLE]

where

[TABLE]

is the dot product of the control points and B-spline basis functions evaluated at $(x,t)$.

Proof.

See Appendix -A. ∎

VI Case Studies

In this section, we present two case studies for the proposed neural spline operator methods to estimate safety probabilities, specifically on a 1-dimensional system with varying functional dynamics in section VI-A, and on a 14-dimensional multi-agent system with varying dynamics and safe sets in section VI-B. Code is available at https://github.com/jacobwang925/NeSO.

VI-A Functional Dynamics

We consider system (1) with $x\in\mathbb{R}$ and safe set (2) with $\alpha=4$ and $\phi_{\alpha}(x)=x-4$ fixed. The dynamics functional of the system is given by

[TABLE]

where $A_{1},A_{2}\sim\mathcal{U}(-1,1)$ are random amplitudes, $\omega_{1},\omega_{2}\sim\mathcal{U}(0.5,2)$ are randomly chosen frequencies, and $\psi_{1},\psi_{2}\sim\mathcal{U}(0,2\pi)$ are random phase shifts. Here, $\mathcal{U}(a,b)$ denotes uniform distribution in $[a,b]$. With probability 0.5, we set $A_{2}=0$, reducing the function to a single sine wave. The goal is to estimate the following long-term recovery probability, which is a slight variant of (3).

[TABLE]

This characterizes the probability of the system recovering to the safe region for at least once under time horizon $t$ starting from an initial state $x$. Such probability can be estimated either via MC, or through solving (4), with the following modified initial and boundary condition [3]

[TABLE]

where $\mathcal{C}_{\alpha}^{c}$ is the complement of $\mathcal{C}_{\alpha}$, *i.e., *unsafe region. In this case study, our aim is to estimate (22) in the state-temporal space $(x,t)\in[-10,4]\times[0,10]$, with varying system dynamics $f$.

The ground truth recovery probability can be calculated via the cumulative integral of the hitting time density:

[TABLE]

where the integrated effect of the dynamics function is

[TABLE]

For the choice of neural operators, we compare the proposed NeSO method with Fourier neural operator (FNO) [24] as the backbone of the neural coefficient functional, against an architecture only of FNO, to learn the mappings from the dynamics function $f$ directly to the recovery probability $F$. In both methods, we use 8 Fourier modes in spatial and temporal dimensions with width of 32, and 3 spectral convolution layers in the FNO module. We train both NeSO and FNO with 10 independent random functions and corresponding ground truth solutions for data loss function, and impose physics loss based on (4), with weights $w_{d}=3$ and $w_{p}=1$. Additional loss on the residual of initial condition and boundary condition (ICBC) (23) is imposed for FNO with weight $w_{\text{ICBC}}=10$. We train both models with Adam optimizer with an initial learning rate of 0.001 for 500 epochs. Both methods converge during training, and Fig. 2 visualizes the prediction results on a random testing function $f$. It can be seen that both methods generalize to unseen dynamics to provide accurate estimation, due to the functional learning properties of neural operators. Table I summarizes the training time and prediction errors on 10 independent random testing functions for both methods. Mean square errors (MSE), mean absolute errors (MAE) and average relative errors are reported. It can be seen that with the B-spline representation, the proposed NeSO outperforms FNO in both training efficiency and prediction accuracy, due to the analytical derivative calculations during training and the direct assignment of ICBCs.

VI-B High-Dimensional Multi-Agent Dynamics

We consider a multi-agent system consisting of $N=7$ agents. Each agent $\Sigma_{k}$ is a mass-spring-damper system of the following dynamics,

[TABLE]

where $x^{(k)}:=[p^{(k)},v^{(k)}]\in\mathbb{R}^{2}$ is the state, $u^{(k)}\in\mathbb{R}$ is the control input, and $w^{(k)}\in\mathbb{R}^{2}$ is the standard Wiener process to capture the noise, and $\sigma_{k}=0.2\;I_{2}\in\mathbb{R}^{2}$ is the noise magnitude. Here, $I_{N}$ is the $N$-dimensional identity matrix and we will use this notation for the rest of this section. The system matrices are defined as $A=\begin{bmatrix}0&1\\ -\beta_{1}&-\beta_{2}\end{bmatrix}$ and $B=\begin{bmatrix}0\\ 1\end{bmatrix}$, where $\beta_{1},\beta_{2}>0$ are the physical parameters governing the spring and damping effects, respectively. The output of each agent is given by

[TABLE]

The control input $u^{(k)}_{t}$ of each dynamic agent is determined by the outputs $y^{(i)}_{t}$ of its neighboring agents, thereby enabling mutual interaction among the agents. Specifically, the control input is given by

[TABLE]

where $\mathcal{I}(k)\subseteq\{1,2,\ldots,N\}$ denotes the set of neighboring agents of $\Sigma_{k}$, and $l_{k,i}$ is the feedback control gain that characterizes the interaction force of agent $i$ on agent $k$.

The agents interact through an undirected graph, encoded in a Laplacian matrix $L\in\mathbb{R}^{N\times N}$ whose $(i,k)$-th entry is defined by $L_{i,k}=l_{i,k}$, representing the interaction network among the agents. Specifically, in this case study we consider

[TABLE]

The overall system dynamics for the multi-agent system are obtained by stacking the states of all agents into the collective state vector x_{t}=\bigl{[}x^{(1)\top}_{t},\dots,x^{(N)\top}_{t}\bigr{]}^{\top}\in\mathbb{R}^{2N}, yielding

[TABLE]

where $w_{t}\in\mathbb{R}^{2N}$ is the standard Wiener process and $\sigma=\text{diag}(\text{diag}(\sigma_{1}),\cdots,\text{diag}(\sigma_{N}))=0.2\,I_{2_{N}}$. The matrix $H$ governs the coupled system dynamics and is defined as

[TABLE]

where $\otimes$ denotes the Kronecker product, and $BC$ describes the coupling between agents through their outputs.

The safe set of the multi-agent system is defined as (2) with

[TABLE]

where $\mathbf{t}_{k}$ is the $k$-th eigenvector of the Laplacian matrix $L$, and $\alpha_{k}\in\mathbb{R}^{+}$ is the corresponding safety threshold for agent $k$. The goal is to estimate the safety probability (3) for the entire multi-agent system (30).

Let $\lambda_{k}$ be the corresponding eigenvalue of the Laplacian matrix $L$ for mode $k$, let

[TABLE]

and we define safe sets for sub-systems as

[TABLE]

The safety probability for sub-system $x^{(k)}$ is defined as

[TABLE]

and we know its value evolves according to the following PDE [46]

[TABLE]

with the following ICBC

[TABLE]

Here, $\gamma_{k}=\beta_{1}+\lambda_{k}$, and the spatial domain is constrained by a threshold $[-\alpha_{k},\alpha_{k}]$. It is shown in [46, Theorem 1] that the overall safety probability can be calculated via

[TABLE]

which is the multiplication of the safety probabilities for all sub-systems.

We consider fixed $\beta_{1}=1$, and varying $\beta_{2}\sim\mathcal{U}(0.5,2)$ and $\alpha_{k}\sim\mathcal{U}(1,2)$ for different agent $k$. We then generate training data on 20 random systems with ground truth safety probability for each sub-system (35) generated through solving (36) with numerical solvers. With that, we train a NeSO that maps $(\beta_{2},\lambda_{k},\alpha_{k})$ to $F_{k}$ for all $x^{(k)}\in[-\alpha_{k},\alpha_{k}]$ and $t\in[0,10]$. After training, we test on $\beta_{2}=1$ and $[\alpha_{1},\cdots,\alpha_{7}]=[2,2,1,1,1,1,1]$, which is not seen during training. Fig. 3 shows the safety probability prediction results for the full 14-dimensional system, via prediction from NeSO on (35) and multiplication in (38). The ground truth is obtained by running MC on the 14-dimensional multi-agent system (30) with 10000 trajectories and empirically estimate (3). It can be seen that NeSO reconstructs the safety probability of the high-dimensional multi-agent system with high accuracy, even though it did not train on data from this system. The time-averaged absolute difference between the MC and NeSO estimation is 0.0074, indicating the efficacy of the NeSO prediction. Table II shows the computation time for safety probability estimation of the 14-dimensional system through MC, solving PDEs (36) for sub-systems, and the proposed NeSO. It can be seen that NeSO achieves a 2 orders of magnitude reduction in online computation time, compared to both the MC method and the PDE method. The training time for NeSO is 3999$\mathrm{s}$ on an L4 GPU. Based on this, we know that the total time for evaluating safety probabilities of $n$ dynamics and safe set instances is $4618n$ seconds for MC, $2079n$ seconds for the PDE method, and $3999+57n$ seconds for NeSO (including training). It can be seen that even with $n=1$, the proposed NeSO method is faster than MC, and it will obtain better efficiency than the PDE method when the number of evaluated dynamics $n>2$.

VII Conclusions

In this paper, we consider the problem of risk quantification for stochastic systems with varying functional dynamics. We introduce physics-informed operator learning methods to this problem to directly find mappings from varying dynamic functionals and related parameters to the safety probabilities, while respecting all PDE constraints. Specifically, we propose neural spline operators, a framework for operator learning based on B-spline representations that improves training efficiency and estimation accuracy for risk quantification. We provide a theoretical analysis of the neural spline operator as universal approximators. We also show in experiments the efficacy of the proposed framework dealing with functional dynamics and high-dimensional multi-agent systems, where tremendous speedup is obtained compared to existing methods. The proposed framework and the accompanying universal approximation theorem are expected to benefit other related problems beyond risk quantification.

-A Proof of Theorem 1

Proof.

(Theorem 1)

We state the general proof of the universal approximation property of neural operators with basis functions, for learning continuous mappings between Sobolev function spaces. In the end, we show that this proof applies to our setting, where B-spline basis is considered and mappings to safety probabilities are studied.

First, we prove the coefficient neural functional is a universal approximator of the optimal basis function weights. We prove by construction. We start by showing the neural operator in the coefficient neural functional can approximate the target function, and show that a unique linear mapping can be learned by the linear layer in the coefficient neural functional to find the optimal weights.333Note that this is not the only possible setting for the coefficient neural functional to approximate the optimal basis function weights.

We denote the input Sobolev function space as

[TABLE]

where $k\in\mathbb{N}^{+}$ and $D^{r}h$ denotes the weak derivative of order $r$. Similarly, we denote the output function space as

[TABLE]

Let $W:\mathcal{H}\to\mathcal{Y}$ be a continuous non-linear mapping, and let $h=(f,\alpha)\in\mathcal{K}_{1}\times\mathcal{K}_{2}\subset\mathcal{H}$. The goal is to approximate $g=W(h)$ via a neural operator, where $g\in\mathcal{K}_{3}\subseteq\mathcal{Y}$ is the target function.

From the universal approximation theorem for neural operators (*e.g., *FNO [43], DeepONet [27]), there exists a continuous approximation of $W$ indicated as $\hat{W}$ defined over the same spaces, such that, for any $\epsilon_{1}>0$

[TABLE]

For a more compact notation, let us define $g_{h}\triangleq W(h)$ and define $\hat{g}_{h}\triangleq\hat{W}(h)$.

We then take $N<\infty$ orthogonal elements $\phi_{i}\in\mathcal{Y}$, and define a finite-dimensional subspace444The proof naturally generalizes to cases where $\phi_{i}\in\mathcal{Y}$ are not orthogonal, but $B_{{\phi}_{i,j}}=\langle\phi_{j},\phi_{i}\rangle$ is positive definite. In our case B-spline basis functions are orthogonal to each other.

[TABLE]

For any function $g_{h}\in\mathcal{K}_{3}\subseteq\mathcal{Y}$, we want to find a unique element $P_{V_{N}}g_{h}\in V_{N}$ such that

[TABLE]

where $P_{V_{N}}$ is a projection from $\mathcal{Y}$ to $V_{N}$. Since Sobolev spaces are Hilbert spaces, they are equipped with a norm $\|\cdot\|_{\mathcal{Y}}$ induced by the inner product. Considering

[TABLE]

we want

[TABLE]

This equation can be written as

[TABLE]

By defining the matrix $B_{{\phi}_{i,j}}=\langle\phi_{j},\phi_{i}\rangle$, which is invertible since $\phi_{j}$ are orthogonal thus linearly independent, we can find the coefficients of $P_{V_{N}}$ denoted by $C^{*}$ with

[TABLE]

where the elements of $b(g_{h})$ are $b_{i}(g_{h})=\langle g_{h},\phi_{i}\rangle$.

For the uniqueness of the projection, consider the presence of $v_{1}$ and $v_{2}$ in $V_{N}$ such that

[TABLE]

From the linearity of the inner product, we can rewrite the above equation as

[TABLE]

Since $v_{1}-v_{2}\in V_{N}$, we can take $v=v_{1}-v_{2}$, then we have

[TABLE]

This implies that $v_{1}=v_{2}$, which demonstrates the uniqueness of the projection.

For the vector $b(g_{h})$, we have

[TABLE]

By the Cauchy-Schwarz inequality

[TABLE]

we get

[TABLE]

Hence,

[TABLE]

Since $B_{\phi}\in\mathbb{R}^{N\times N}$ is symmetric and positive-definite, the operator norm of $B_{\phi}^{-1}$ is

[TABLE]

where $\lambda_{\min}(B_{\phi})>0$ denotes the smallest eigenvalue of $B_{\phi}$. For any $g_{h},\hat{g}_{h}\in\mathcal{K}_{3}\subseteq\mathcal{Y}$, we have:

[TABLE]

We then define the following Lipschitz constants

[TABLE]

and

[TABLE]

where $\lambda_{\max}(B_{\phi})$ is the largest eigenvalue of $B_{\phi}$. From the universal approximation of the neural operator layer (41), we have

[TABLE]

and we can choose $\epsilon_{1}=\frac{\epsilon}{4L_{1}L_{2}}$ where $\epsilon>0$. Hence,

[TABLE]

Since $C^{*}(\cdot)$ is a unique bounded linear mapping, it can be approximated by a neural network. By the universal approximation theorem [47], there exists a neural network $\hat{C}^{*}$ such that

[TABLE]

By using the triangular inequality, and choosing $\epsilon_{2}=\epsilon/(4L_{2})$, we can write

[TABLE]

Then, we show that given the coefficient neural functional is a universal approximator of the optimal weights, and the basis functions has sufficient representation, the overall framework is a universal approximator.

Let the basis $\{\phi_{i}\}_{i=1}^{N}$ and the number of basis functions $N$ are chosen so that the subspace $V_{N}$ is sufficiently rich to approximate any target function within a given tolerance $\epsilon_{3}>0$:

[TABLE]

Recalling that $g=W(h)$, the output of the operator $W(\cdot)$ with $h$ as input function, we can write

[TABLE]

By selecting $\epsilon_{3}=\epsilon/2$, we obtain

[TABLE]

In our case, we use B-splines as basis functions (13), therefore

[TABLE]

and $C^{*}(g_{h})$ represents the optimal control points tensor denoted by $C^{*}$. Since B-splines are universal approximators [5], there exists the number of control points $\ell_{1},\cdots,\ell_{n},\ell_{t}\in\mathbb{N}^{+}$, order of B-spline basis $d_{1},\cdots,d_{n},d_{t}\in\mathbb{N}^{+}$, such that

[TABLE]

Since the coefficient neural functional produces $\hat{C}^{*}(\hat{g}_{h})$ and is denoted by $G_{\theta}(h)=G_{\theta}(f,\alpha)$, we have

[TABLE]

and for our setting, $g_{h}=F$ is the safety probability.

∎