Polynomial Chaos Expansion for Operator Learning

TL;DR

This paper introduces polynomial chaos expansion as a novel operator learning method in scientific machine learning, capable of approximating PDE solution operators with high accuracy and providing uncertainty quantification without extra cost.

Contribution

It establishes a mathematical framework for PCE in operator learning, enabling data-driven and physics-informed approximations, and demonstrates its effectiveness on various PDE problems.

Findings

High numerical accuracy in operator approximation

Efficient uncertainty quantification via PCE coefficients

Strong performance on diverse PDE problems

Abstract

Operator learning (OL) has emerged as a powerful tool in scientific machine learning (SciML) for approximating mappings between infinite-dimensional functional spaces. One of its main applications is learning the solution operator of partial differential equations (PDEs). While much of the progress in this area has been driven by deep neural network-based approaches such as Deep Operator Networks (DeepONet) and Fourier Neural Operator (FNO), recent work has begun to explore traditional machine learning methods for OL. In this work, we introduce polynomial chaos expansion (PCE) as an OL method. PCE has been widely used for uncertainty quantification (UQ) and has recently gained attention in the context of SciML. For OL, we establish a mathematical framework that enables PCE to approximate operators in both purely data-driven and physics-informed settings. The proposed framework reduces the task of learning the operator to solving a system of equations for the PCE coefficients. Moreover, the framework provides UQ by simply post-processing the PCE coefficients, without any additional computational cost. We apply the proposed method to a diverse set of PDE problems to demonstrate its capabilities. Numerical results demonstrate the strong performance of the proposed method in both OL and UQ tasks, achieving excellent numerical accuracy and computational efficiency.

Tables

Table 1. Table 1: Performance of PCE for operator learning and UQ across different example problems.

Problem	Operator Learning		Uncertainty Quantification
	MSE	Training	MAE		Comp. time (s)
		time (s)	$\mu$	$\sigma$	MCS	PC 2
Anti-derivative	4.55e-9	0.3	1.50e-3	2.10e-3	2.02	1
Advection-Diffusion	7.01e-8	0.4	1.58e-4	4.05e-5	1620	1.6
1D Burgers	3.38e-7	0.5	8.88e-5	1.11e-4	1737	1.7
2D Heat	6.21e-4	0.8	6.48e-16	6.19e-4	97200	56

Table 2. Table 2: Performance of PC 2 for operator learning and UQ across different example problems.

Problem	Operator Learning		Uncertainty Quantification
	MSE	Training	MAE		Comp. time (s)
		time (s)	$\mu$	$\sigma$	MCS	PC 2
Anti-derivative	1.51e-8	0.5	2.21e-3	2.09e-3	2.02	1.6
Advection-Diffusion	4.79e-5	30	4.68e-3	1.31e-3	1620	32
1D Burgers	4.41e-5	76	2.08e-3	1.1e-3	1737	84
2D Heat	6.24e-4	14	1.09e-6	5.3e-4	97200	116

Equations

$$\text{pazocal}G:\text{pazocal}U\to \text{pazocal}S,u(x)\mapsto s(y),$$

$${\theta }^{\ast }=\arg \underset{\theta }{\min }\frac{1}{N}\sum _{i=1}^N\text{pazocal}L\left(s^{(i)},\text{pazocal}{G}_{\theta }\left(u^{(i)}\right)\right),$$

$$\text{pazocal}{N}_{\nu }[s(\mathbf{x},t)]=f(\mathbf{x},t),\mathbf{x}\in \text{pazocal}D,t\in [0,T],$$

$$\text{pazocal}B[s(\mathbf{x},t)]=h(\mathbf{x},t),\mathbf{x}\in \partial \text{pazocal}D,t\in [0,T],$$

$$s(\mathbf{x},0)=g(\mathbf{x}),\mathbf{x}\in \text{pazocal}D,$$

$$\text{pazocal}G(u(\mathbf{x},t,\mathbf{\xi }\left(\omega \right)))(\mathbf{x},t)\approx \hat{s}(\mathbf{x},t,\mathbf{\xi }\left(\omega \right))=\sum _{\mathbf{\alpha }\in \text{pazocal}A}{s}_{\mathbf{\alpha }}\left(\mathbf{x},t\right){\Psi }_{\mathbf{\alpha }}\left(\mathbf{\xi }\left(\omega \right)\right),$$

$${\int }_{\text{pazocal}{D}_{\mathbf{\xi }}}{\Psi }_{\alpha }(\mathbf{\xi }){\Psi }_{\beta }(\mathbf{\xi })p_{\xi }(\mathbf{\xi })d\mathbf{\xi }=0(\mathbf{\alpha }\ne \mathbf{\beta }),$$

$$\text{pazocal}A=\left\{\mathbf{\alpha }\in {\mathbb{N}}^r:{\Vert \mathbf{\alpha }\Vert }_1\le p\right\}.$$

$$\mathrm{card}\text{pazocal}A=\frac{\left(r+p\right)!}{r!p!}\equiv P.$$

$$s_{\mathbf{\alpha }}\left(\mathbf{x},t\right)=\sum _{\mathbf{\beta }\in \text{pazocal}B}{c}_{\mathbf{\alpha }\mathbf{\beta }}{\Phi }_{\mathbf{\beta }}\left(\mathbf{x},t\right),$$

$$Q=\frac{\left(d+1+q\right)!}{(d+1)!q!}.$$

$$\hat{s}(\mathbf{x},t,\mathbf{\xi })=\sum _{\mathbf{\alpha }\in \text{pazocal}A}\sum _{\mathbf{\beta }\in \text{pazocal}B}{c}_{\mathbf{\alpha }\mathbf{\beta }}{\Phi }_{\mathbf{\beta }}\left(\mathbf{x},t\right){\Psi }_{\mathbf{\alpha }}(\mathbf{\xi }).$$

$$\hat{\mathbf{S}}=\mathbf{\Phi }\mathbf{C}\mathbf{\Psi },$$

$$\bm{\Phi}=\left[\begin{array}[]{ccc}\Phi_{1}\left((\mathbf{x},t)^{(1)}\right)&\ldots&\Phi_{Q}\left((\mathbf{x},t)^{(1)}\right)\\ \vdots&\ddots&\vdots\\ \Phi_{1}\left((\mathbf{x},t)^{(n)}\right)&\ldots&\Phi_{Q}\left((\mathbf{x},t)^{(n)}\right)\end{array}\right],$$

$$\bm{\Psi}=\left[\begin{array}[]{ccc}\Psi_{1}\left(\bm{\xi}^{(1)}\right)&\ldots&\Psi_{1}\left(\bm{\xi}^{(N)}\right)\\ \vdots&\ddots&\vdots\\ \Psi_{P}\left(\bm{\xi}^{(1)}\right)&\ldots&\Psi_{P}\left(\bm{\xi}^{(N)}\right)\end{array}\right].$$

$$\mathbf{C}^{*}=\dfrac{1}{n\cdot N}\underset{\tilde{\mathbf{C}}}{{\arg\min}}\Big{\|}\mathbf{S}-\hat{\mathbf{S}}\Big{\|}^{2},$$

$${\mathbf{C}}^{\ast }={\left({\mathbf{\Phi }}^{\top }\mathbf{\Phi }\right)}^{-1}\left({\mathbf{\Phi }}^{\top }\mathbf{S}{\mathbf{\Psi }}^{\top }\right){\left(\mathbf{\Psi }{\mathbf{\Psi }}^{\top }\right)}^{-1}.$$

$$\text{pazocal}{N}_{\nu }[\hat{s}(\mathbf{x},t,\mathbf{\xi })]=\sum _{\mathbf{\alpha }\in \text{pazocal}A}\sum _{\mathbf{\beta }\in \text{pazocal}B}{c}_{\mathbf{\alpha }\mathbf{\beta }}\text{pazocal}{N}_{\nu }\left[{\Phi }_{\mathbf{\beta }}\left(\mathbf{x},t\right)\right]{\Psi }_{\mathbf{\alpha }}(\mathbf{\xi })\approx f(\mathbf{x},t,\mathbf{\xi }),$$

$$\text{pazocal}B[\hat{s}(\mathbf{x},t,\mathbf{\xi })]=\sum _{\mathbf{\alpha }\in \text{pazocal}A}\sum _{\mathbf{\beta }\in \text{pazocal}B}{c}_{\mathbf{\alpha }\mathbf{\beta }}\text{pazocal}B\left[{\Phi }_{\mathbf{\beta }}\left(\mathbf{x},t\right)\right]{\Psi }_{\mathbf{\alpha }}(\mathbf{\xi })\approx h(\mathbf{x},t,\mathbf{\xi }),$$

$$\hat{s}(\mathbf{x},0,\mathbf{\xi })=\sum _{\mathbf{\alpha }\in \text{pazocal}A}\sum _{\mathbf{\beta }\in \text{pazocal}B}{c}_{\mathbf{\alpha }\mathbf{\beta }}{\Phi }_{\mathbf{\beta }}\left(\mathbf{x},0\right){\Psi }_{\mathbf{\alpha }}(\mathbf{\xi })\approx g(\mathbf{x},\mathbf{\xi }).$$

$${\mathbf{R}}_{\mathrm{PDE}}$$

$${\mathbf{R}}_{\mathrm{BC}}$$

$${\mathbf{R}}_{\mathrm{IC}}$$

$$L(\mathbf{C})$$

$$=\frac{1}{n_{\mathrm{PDE}}\cdot N}{\Vert {\mathbf{R}}_{\mathrm{PDE}}\Vert }^2+\frac{1}{n_{\mathrm{BC}}\cdot N}{\Vert {\mathbf{R}}_{\mathrm{BC}}\Vert }^2+\frac{1}{n_{\mathrm{IC}}\cdot N}{\Vert {\mathbf{R}}_{\mathrm{IC}}\Vert }^2.$$

$$\nabla L(\mathbf{C})=$$

$$+$$

$$\left[\frac{2}{n_{\mathrm{PDE}}}{\mathbf{\Phi }}_{\mathrm{PDE}}^{\top }{\mathbf{\Phi }}_{\mathrm{PDE}}+\frac{2}{n_{\mathrm{IC}}}{\mathbf{\Phi }}_{\mathrm{IC}}^{\top }{\mathbf{\Phi }}_{\mathrm{IC}}+\frac{2}{n_{\mathrm{BC}}}{\mathbf{\Phi }}_{\mathrm{BC}}^{\top }{\mathbf{\Phi }}_{\mathrm{BC}}\right]\mathbf{C}\mathbf{\Psi }{\mathbf{\Psi }}^{\top }$$

$$=\frac{2}{n_{\mathrm{PDE}}}{\mathbf{\Phi }}_{\mathrm{PDE}}^{\top }{\mathbf{F}}_{\mathrm{PDE}}{\mathbf{\Psi }}^{\top }+\frac{2}{n_{\mathrm{IC}}}{\mathbf{\Phi }}_{\mathrm{IC}}^{\top }{\mathbf{F}}_{\mathrm{IC}}{\mathbf{\Psi }}^{\top }+\frac{2}{n_{\mathrm{BC}}}{\mathbf{\Phi }}_{\mathrm{BC}}^{\top }{\mathbf{F}}_{\mathrm{BC}}{\mathbf{\Psi }}^{\top }$$

$$\begin{array}{rl}{\mathbf{C}}^{\ast }& ={\left[\frac{2}{n_{\mathrm{PDE}}}{\mathbf{\Phi }}_{\mathrm{PDE}}^{\top }{\mathbf{\Phi }}_{\mathrm{PDE}}+\frac{2}{n_{\mathrm{IC}}}{\mathbf{\Phi }}_{\mathrm{IC}}^{\top }{\mathbf{\Phi }}_{\mathrm{IC}}+\frac{2}{n_{\mathrm{BC}}}{\mathbf{\Phi }}_{\mathrm{BC}}^{\top }{\mathbf{\Phi }}_{\mathrm{BC}}\right]}^{-1}\\ & \left[\frac{2}{n_{\mathrm{PDE}}}{\mathbf{\Phi }}_{\mathrm{PDE}}^{\top }{\mathbf{F}}_{\mathrm{PDE}}{\mathbf{\Psi }}^{\top }+\frac{2}{n_{\mathrm{IC}}}{\mathbf{\Phi }}_{\mathrm{IC}}^{\top }{\mathbf{F}}_{\mathrm{IC}}{\mathbf{\Psi }}^{\top }+\frac{2}{n_{\mathrm{BC}}}{\mathbf{\Phi }}_{\mathrm{BC}}^{\top }{\mathbf{F}}_{\mathrm{BC}}{\mathbf{\Psi }}^{\top }\right]{\left[\mathbf{\Psi }{\mathbf{\Psi }}^{\top }\right]}^{-1}\end{array}$$

$${\mathbf{C}}^{(k+1)}={\mathbf{C}}^{(k)}-{\left[\text{pazocal}HL\left({\mathbf{C}}^{(k)}\right)\right]}^{-1}\nabla L\left({\mathbf{C}}^{(k)}\right),$$

$$\mu (\mathbf{x},t)=\mathbb{E}[\hat{\mathbf{S}}]=\mathbf{\Phi }(\mathbf{x},t){\mathbf{c}}_0\in {\mathbb{R}}^n$$

$$\mathbf{\Sigma }=\mathbf{\Phi }(\mathbf{x},t)\left({\mathbf{C}}^{\ast }\mathrm{diag}(0,1,\dots ,1){{\mathbf{C}}^{\ast }}^{\top }\right){\mathbf{\Phi }}^{\top }(\mathbf{x},t)\in {\mathbb{R}}^{n\times n}$$

$$\frac{d}{dx}s(x)$$

$$\Rightarrow \text{pazocal}G(u(x))(x)$$

Full text

Keywords Operator leaning  $\cdot$ DeepONet  $\cdot$ Polynomial chaos expansion  $\cdot$ Scientific machine learning  $\cdot$ Uncertainty quantification  $\cdot$ Physics-informed

1. Introduction

Operator learning (OL) has emerged as an influential field within scientific machine learning (SciML), focusing on learning mappings between infinite-dimensional functional spaces. In particular, OL has demonstrated significant potential in approximating the solution operator of partial differential equations (PDEs) by effectively mapping initial conditions, boundary conditions, source terms, and/or parameters to the corresponding PDE solution. The majority of existing work in OL has focused on neural network (NN) based approaches, with frameworks such as Deep Operator Networks (DeepONet) [1], Fourier Neural Operators (FNO) [2], Wavelet Neural Operators (WNO) [3], Laplace Neural Operators (LNO) [4], Convolutional Neural Operators (CNO) [5], and Graph Neural Operators (GNO) [6]. These models and their many variants have demonstrated significant potential across a wide range of scientific applications.

Recently, there has been growing research interest in exploring traditional machine learning (ML) methods within the OL framework. Particularly, Gaussian Processes (GPs) and other kernel-based approaches have demonstrated strong performance, especially in data-scarce settings where NN-based approaches may struggle. These methods offer interpretable solutions with built-in uncertainty quantification, as well as theoretical and computational guarantees [7, 8, 9]. Furthermore, hybrid approaches that integrate traditional ML methods into neural operator networks have also shown improved performance. For example, Lowery et al. [10] proposed the Kernel Neural Operator (KNO), which incorporates trainable kernel functions into neural operator architectures and achieves competitive accuracy with fewer parameters compared to FNO and GNO. Similarly, a hybrid GP/NN framework [8], which combines the expressiveness of NNs with the uncertainty quantification and interpretability of GPs, has shown enhanced performance. These traditional and hybrid frameworks offer a promising direction for OL, particularly in terms of computational efficiency, convergence guarantees, reduced data and resource requirements, and uncertainty quantification.

In this work, we introduce polynomial chaos expansion (PCE) for OL. PCE is a widely used method in the field of uncertainty quantification (UQ), where it approximates the solution of a stochastic system using a spectral expansion of orthogonal multivariate polynomials with deterministic coefficients. The choice of orthogonal polynomials depends on the probability distributions of the input stochastic variables, typically following the Wiener–Askey scheme [11]. Alternatively, for arbitrary input distributions, the polynomials can be constructed using approaches such as those proposed in [12]. Recently, Bahmani et al. [13] proposed a neural chaos method, where the orthogonal polynomial basis functions are replaced with NN parameterized basis functions in a purely data-driven manner. Once the basis is selected, the deterministic coefficients in PCE can be computed using either intrusive approaches, such as the stochastic Galerkin (SG) method [14], or non-intrusive approaches, such as the stochastic collocation method [15]. Non-intrusive methods are commonly categorized into three main types: interpolation [16, 17], regression [18], and pseudo-projection methods [19]. The resulting PCE coefficients can then be post-processed to efficiently estimate various response statistics, such as moments and sensitivities with respect to the input stochastic variables [20, 21, 22].

PCE has also been recently employed as an ML method in a purely data-driven setting, demonstrating competitive performance compared to NNs [23]. In SciML, we have recently proposed physics-constrained polynomial chaos expansion (

PC

2. [24, 25] – which incorporates physical constraints into the PCE framework by solving a constrained least squares regression problem to determine the PCE coefficients – and manifold PCE (m-PCE) [26, 27], which develops data-driven PCE on low-dimensional manifolds for high-dimensional UQ and has demonstrated competitive performance to DeepONet for complex PDEs [28]. While

PC

2 and m-PCE are not explicitly presented as OL methods, they lay the groundwork for the current work, which formally establishes PCE as an OL method.

This work presents the mathematical formulation for approximating the solution operator of a given PDE using PCE. The proposed framework reduces the task of OL to solving a system of equations for the PCE coefficients. Furthermore, the proposed method is capable of learning operators in both purely data-driven and physics-informed settings. Our framework is mathematically rigorous, easy to implement, and computationally efficient. Moreover, it provides UQ at no additional computational cost. We demonstrate the effectiveness of the proposed method in both OL and UQ tasks by applying it to a diverse set of PDE problems.

2. Methodology

Consider a general operator learning task, where the goal is to learn an operator $\pazocal{G}$ mapping input functions $u$ defined over a domain $\pazocal{D}$ to output functions $s$ defined over a domain $\pazocal{R}$:

[TABLE]

where $\pazocal{U}$ and $\pazocal{S}$ are function spaces (e.g., Banach or Hilbert spaces), $x\in\pazocal{D}\subseteq\mathbb{R}^{d_{x}}$ and $y\in\pazocal{R}\subseteq\mathbb{R}^{d_{y}}$. Here, $d_{x}$ and $d_{y}$ denote the dimensionality of the input and output domains, respectively.

The learning goal is to find an approximate operator $\pazocal{G}_{\theta}$, parameterized by $\theta$, that minimizes the discrepancy between predicted and true outputs over the given dataset $\left\{\left(u^{(i)},s^{(i)}\right)\right\}_{i=1}^{N},\text{ with }s^{(i)}=\pazocal{G}\left(u^{(i)}\right)$:

[TABLE]

where $\pazocal{L}$ is a suitable loss function (e.g., mean squared error) or any other metric suitable for the operator-learning task at hand.

In scientific and engineering applications, a common scenario involves learning an operator that is implicitly defined by a physical system governed by partial differential equations (PDEs). Consider a general PDE defined on a spatial domain $\pazocal{D}\subseteq\mathbb{R}^{d}$ and time domain $t\in$ $[0,T]$ as

[TABLE]

where $s(\mathbf{x},t)$ is the unknown PDE solution, $\pazocal{N}_{\nu}$ represents a general differential operator parameterized by $\nu$, $\pazocal{B}$ corresponds to the boundary operator, and $\partial\pazocal{D}$ denotes the boundary of the spatial domain. The functions $f(\mathbf{x},t),\ h(\mathbf{x},t)$, and $g(\mathbf{x})$ represent the source term, boundary condition, and initial condition, respectively. In the context of operator learning, any combination of these inputs, i.e., $f$, $h$, $g$, and $\nu$, may collectively form the input function $u$. The goal is to learn the mapping from $u$ to the corresponding PDE solution $s$.

For the task of learning the operator, these input functions are usually sampled from random fields $u(\mathbf{x},t,\omega)$, $\omega\in\Omega$, and $\Omega$ denotes the underlying sample space of the probability space ($\Omega,\mathscr{F},\mathscr{P})$ where $\mathscr{F}$ is the sigma-algebra of events and $\mathscr{P}$ is the associated probability measure. Once the operator is learned, it can provide the PDE solution corresponding to any given realization of the input random fields. This problem setup naturally fits within the framework of PCE, where PCE is used to learn the operator.

We can express the solution operator using a PCE as

[TABLE]

where $\bm{\xi}=\left(\xi_{1},\xi_{2},\ldots,\xi_{r}\right)\in\pazocal{D}_{\bm{\xi}}\subseteq\mathbb{R}^{r}$ is a vector of independent random variables with joint PDF $p_{\bm{\xi}}(\bm{\xi})=$ $\prod_{i=1}^{r}p_{\xi_{i}}\left(\xi_{i}\right)$. The multivariate polynomials $\left\{\Psi_{\alpha}(\bm{\xi})\right\}_{\alpha\in\pazocal{A}}$ are constructed by taking tensor-products of univariate orthogonal polynomials $\psi_{\alpha_{i}}\left(\xi_{i}\right)$. They satisfy the orthogonality relation

[TABLE]

and thus form an orthogonal basis of the weighted Hilbert space $L^{2}\left(\pazocal{D}_{\bm{\xi}},p_{\bm{\xi}}(\bm{\xi})d\bm{\xi}\right)$. These basis functions are typically selected according to the Wiener-Askey scheme [11], but can also be constructed for arbitrary distributions of $\bm{\xi}$ [12]. The coefficients $s_{\bm{\alpha}}\left(\mathbf{x},t\right)$ are deterministic functions of the spatio-temporal variables and must be determined. The index set $\pazocal{A}$ is a finite set of multi-indices with cardinality $P$.

A standard truncation scheme includes all multivariate polynomials in the $r$ input variables whose total degree does not exceed $p$, forming the index set:

[TABLE]

The cardinality of the truncated index set $\pazocal A$ is given by

[TABLE]

Other truncation schemes can be employed to reduce the cardinality of the index set, for example, by limiting the number of interaction terms through approaches such as hyperbolic truncation [29].

Next, we represent the deterministic coefficient function, $s_{\bm{\alpha}}\left(\mathbf{x},t\right)$, as a spectral expansion of orthogonal polynomials given as

[TABLE]

where $\Phi_{\bm{\beta}}\left(\mathbf{x},t\right)$ are multivariate orthogonal polynomial basis functions, constructed as tensor products of one-dimensional orthogonal polynomials with respect to some weight function, $w(x)$, defined over a specified interval for the spatial and temporal dimensions (e.g, Legendre polynomials defined on $[-1,1]$ with $w(x)=1$), and $c_{\bm{\alpha}\bm{\beta}}$ are the corresponding scalar coefficients.

We adopt the same truncation scheme as used for Eq. (4), i.e., $\pazocal B^{d+1,q}=\left\{\bm{\beta}\in\mathbb{N}^{d+1}:{\|\bm{\beta}\|}_{1}\leq q\right\},$ where $q$ denotes the total degree of the expansion, $d$ is the number of spatial dimensions, which can be at most three, along with one temporal dimension. The cardinality of set $\pazocal{B}$ is given by

[TABLE]

Substituting Eq. (7) in Eq. (4), we get the expansion of the solution operator as

[TABLE]

Rewriting Eq. (9) in matrix form, we get

[TABLE]

where

[TABLE]

[TABLE]

Here, $\bm{\Phi}\in\mathbb{R}^{n\times Q}$, $\bm{\Psi}\in\mathbb{R}^{P\times N}$, $n$ denotes the number of spatio-temporal coordinates at which the operator is evaluated, and $N$ is the number of samples of the input functions. The unknowns are the coefficients $\mathbf{C}\in\mathbb{R}^{Q\times P}$, which represent $\theta$ in Eq. (2).

The OL task reduces to learning the coefficient matrix of the PCE. Given the labeled training data, the coefficient can be obtained by minimizing the residual as

[TABLE]

where $\mathbf{S}$ is the true solution matrix and $\hat{\mathbf{S}}$ is the approximate PCE-based solution matrix, as defined in Eq. (10).

The closed-form solution for Eq. (13) reads as

[TABLE]

For a unique solution to Eq. (14), it is required that $n\geq Q$ and $N\geq P$. Once $\mathbf{C}^{*}$ is obtained, substituting it back into Eq. (10) yields the approximate solution for a given test input function at any desired spatial-temporal coordinates.

Next, we learn the operator in a physics-informed manner, i.e., without requiring labeled training data. We refer to this approach as the physics-constrained polynomial chaos expansion (

PC

2). It is worth noting that the present

PC

2 formulation differs from our previous work on this topic [24, 25], and is significantly more efficient and better tailored for operator learning. In this work, the stochastic and spatio-temporal basis functions are separated, leading to a far more compact expansion and efficient solution.

Given the general PDE described in Eq. (3), we replace the unknown PDE solution with the PCE approximation from Eq. (9). This gives

[TABLE]

The differential operators apply only to the physical basis functions, $\Phi_{\bm{\beta}}\left(\mathbf{x},t\right)$, as the constraints are defined with respect to the spatio-temporal variables rather than the stochastic variables. The constraints are enforced at randomly selected spatio-temporal virtual points for $N$ given realizations of the input random vector. Let the number of virtual points corresponding to the PDE, BC and IC be $n_{\mathrm{PDE}}$, $n_{\mathrm{BC}}$ and $n_{\mathrm{IC}}$, respectively. From Eq. (10) and Eq. (15), we can write the residuals as

[TABLE]

where matrices $\bm{\Phi}_{\mathrm{PDE}}\in\mathbb{R}^{n_{\mathrm{PDE}}\times Q}$, $\bm{\Phi}_{\mathrm{BC}}\in\mathbb{R}^{n_{\mathrm{BC}}\times Q}$ and $\bm{\Phi}_{\mathrm{IC}}\in\mathbb{R}^{n_{\mathrm{IC}}\times Q}$ are evaluated at virtual points corresponding to the PDE, BC and IC, respectively, after the application of the respective differential operators. The matrix $\bm{\Psi}\in\mathbb{R}^{P\times N}$ is evaluated for the given $N$ input realizations of the random vector. The matrices $\mathbf{F}_{\mathrm{PDE}}\in\mathbb{R}^{n_{\mathrm{PDE}}\times N}$, $\mathbf{F}_{\mathrm{BC}}\in\mathbb{R}^{n_{\mathrm{BC}}\times N}$ and $\mathbf{F}_{\mathrm{IC}}\in\mathbb{R}^{n_{\mathrm{IC}}\times N}$ contain the values of input functions corresponding to the PDE, BC, and IC, respectively, evaluated at their associated virtual points for each realization of the input random vector.

Now, using Eq. (16), we can write the loss function for

PC

2 as

[TABLE]

The coefficients can be learned by minimizing the loss function in Eq. (17), for which we set $\nabla L(\mathbf{C})=0$. For linear PDEs, this results in a system of linear equations, while for nonlinear PDEs, it leads to a system of nonlinear equations. Furthermore, a data loss term can be added to Eq. (17) if training data is available. Thus, the task of operator learning reduces to solving a system of equations.

For linear PDEs, we get

[TABLE]

This gives

[TABLE]

The optimal coefficients are obtained as

[TABLE]

For nonlinear PDEs, we solve a nonlinear system of equations using Newton-Raphson as

[TABLE]

where $\nabla L\left(\mathbf{C}^{(k)}\right)\in\mathbb{R}^{Q\times P}$ is the gradient of $L$ evaluated at the current iterate $\mathbf{C}^{(k)}$, and $\pazocal{H}L\left(\mathbf{C}^{(k)}\right)\in\mathbb{R}^{Q\times P}$ is the Hessian matrix of $L$ at $\mathbf{C}^{(k)}$. Closed-form expressions for the gradient and Hessian can be derived due to the polynomial structure of the expansion. We could also compute the optimal coefficients by minimizing the loss function using numerical optimization algorithms, such as gradient descent, ADAM [30], or L-BFGS.

Next, UQ can be efficiently performed by computing the first two moments at any desired spatio-temporal coordinate directly from the coefficient matrix $\mathbf{C}^{*}$, as follows:

[TABLE]

where $\mathbf{c}_{0}$ is the first column of $\mathbf{C^{*}}$ matrix.

The covariance matrix is given as:

[TABLE]

Additional output statistics, such as Sobol sensitivity indices [20, 31], can also be analytically derived from the coefficient matrix $\mathbf{C^{*}}$.

3. Numerical results

In this section, we demonstrate the capabilities of the PCE and

PC

2 in OL and UQ by considering four diverse ODE and PDE problems. For all example problems, the input functions are sampled from a Gaussian random field (GRF), which is discretized using the truncated Karhunen–Loève (KL) expansion [14]. We generate the training and testing data using the FEM solutions corresponding to different realizations of the input function. The training data is used solely for training PCE and is not required for

PC

2. The performance of both PCE and

PC

2 is evaluated on the same testing dataset. To validate the UQ results, we use Monte Carlo simulations (MCS) with 10,000 samples to compute reference estimates of the mean and standard deviation for all example problems. For OL, we assess the numerical accuracy of PCE and

PC

2 by computing the mean square error (MSE) with respect to FEM solutions on the testing samples. For UQ, we evaluate performance using the mean absolute error (MAE) in the predicted mean and standard deviation over the given domain, compared to reference results from MCS.

Table 1 summarizes the performance of PCE on all example problems for both OL and UQ. It can be observed that PCE excels in both OL and UQ tasks, providing exceptional numerical accuracy and computational efficiency. For OL, the MSE across all problems is on the order of $10^{-9}$ to $10^{-4}$, with training times of fractions of a second, which demonstrates excellent accuracy and efficiency of the PCE method. The UQ results show significantly lower MAE for both the mean and standard deviation across all examples, with a fraction of the computational cost compared to MCS, indicating superior performance.

From these results, it is evident that PCE excels in both OL and UQ tasks. However, as discussed in Section 2, the PCE solution (Eq. (14)) is unique when $N\geq P$. Moreover, $P$ increases significantly with the PCE order ($p$) and the stochastic dimensionality ($r$) (Eq. (6)), necessitating a larger amount of training data, which may become computationally infeasible to obtain. For example, in the case of the 2D heat problem, we get $P=2045$, and therefore we use $N=2500$ labeled training samples to ensure a unique solution. Acquiring this amount of labeled training data was computationally expensive ($\sim$7 h). Hence, it is essential to use

PC

2, which eliminates the need for labeled training data.

Table 2 shows the

PC

2 performance for both the OL and UQ tasks. As observed, the

PC

2 achieves excellent numerical accuracy and computational efficiency for both OL and UQ across all example problems.

PC

2 offers competitive performance relative to standard PCE without requiring any labeled training data, making it an ideal choice for OL and UQ, especially for high-dimensional problems. Thus, for all the example problems, we generate visualization plots using

PC

2. The details of the problem setup for each example are provided in the subsequent sections.

The codes were written in Python using the packages UQpy [32], JAX [33], and FEniCS [34]. All computations related to PCE and

PC

2 were performed on a single NVIDIA A100 GPU on the Rockfish HPC facility at Johns Hopkins University. FEM simulations were carried out on a Standard Compute node at Rockfish, equipped with two Intel Xeon Gold 6248R “Cascade Lake” processors (24 cores per CPU, 3.0 GHz base clock) and 192 GB DDR4 RAM.

3.1 Anti-derivative example

This example is adapted from the original DeepONet paper [1]. Here, we consider the problem of learning an anti-derivative operator. Let $u:[0,1]\rightarrow\mathbb{R}$ be an input function. The operator $\pazocal{G}$ maps $u$ to its antiderivative:

[TABLE]

where $u(x)\sim\operatorname{GP}\left(0,\operatorname{Cov}\left(x_{1},x_{2}\right)\right)$. The covariance function is defined by the radial basis functions (RBF) kernel, $k\left(x,x^{\prime}\right)=\sigma^{2}\exp\left(-\dfrac{\left(x-x^{\prime}\right)^{2}}{2\ell^{2}}\right)$, with $\sigma=1$ and $\ell=0.2$. The random field is discretized using six KL terms ($r=6$), which capture $99\%$ of its total variance. We set the polynomial orders to $p=3$ and $q=10$ for both PCE and

PC

2. For PCE, $N=100$ input samples with corresponding output solutions are used as labeled training data, while $N=1000$ samples are reserved for testing. As shown in Table 2, PCE provides excellent performance in terms of both accuracy and computational efficiency for both OL and UQ.

For

PC

2, we generate 1,000 input samples $u(x)$ for training without access to labeled outputs, i.e., relying solely on the PDE constraints. The performance is tested on the same testing dataset as PCE. As previously noted, only the

PC

2 results are used for visualization purposes. The OL plots are shown in Figure 1. For plotting, we select two representative input samples from the test set: one corresponding to the worst prediction error, and the other to the 75th percentile error. As shown in Figure 1, the

PC

2 provides nearly identical results to the ground truth, with MSE on the order of $10^{-8}$ for both cases. The training time for

PC

2 is 0.5 s, highlighting its excellent computational performance in the OL task. Figure 2 shows the UQ plots for this problem, including the predicted mean along with the $\pm 2\sigma$ confidence bounds for both MCS and

PC

2, as well as the absolute errors in the mean and standard deviation. From the plots, it is evident that

PC

2 provides excellent estimation of the mean and standard deviation as compared to MCS. The total computation time is 1.6 s, which includes preprocessing, training, and postprocessing. From this example, it is quite clear that

PC

2 performs both the OL and UQ tasks with high accuracy and exceptional computational efficiency.

3.2 Advection-Diffusion equation

In this example, we consider the 1D advection–diffusion equation with spatially varying advection velocity.

[TABLE]

where $D$ is the diffusion coefficient and the parameter $v(x)$ is the advection velocity.

Here, we consider $D=0.1$ and sample $v(x)$ from a Gaussian random field (GRF) with mean velocity, $\mu=1$, and covariance function defined using RBF kernel with $\sigma=0.05$ and $\ell=0.2$. We use six KL terms ($r=6$) to discretize the random field, capturing $99\%$ of its total variance. The given equation is advection-dominated, with a mean Péclet number, $\mathrm{Pe}=\frac{\mu L}{D}$ = 10, which is often challenging to model using physics-informed machine learning approaches [35]. We first learn the operator that maps the advection velocity $v(x)$ to the solution $u(x)$ using PCE and

PC

2 and then perform UQ. We use polynomial orders as $p=3$ and $q=14$. Similar to the previous example, for PCE, we generate $N=100$ labeled training samples and use $N=1000$ samples for testing. As observed in Table 1, PCE achieves exceptional numerical accuracy with MSE on the order of $10^{-8}$ with just a fraction of a second of training time, highlighting its effectiveness in OL. We get similar performance in UQ, with MAE on the order of $10^{-4}$ for mean and $10^{-5}$ for standard deviation, while requiring only a fraction of the computational time compared to MCS.

For

PC

2, we generate 1,000 input samples of $v(x)$ for training without using the labeled outputs. We use the same testing dataset as PCE for evaluating its performance. Figure 3 shows the OL plots with test samples corresponding to the worst error, and the 75th percentile error. As observed, the

PC

2 solution closely matches the FEM solution even in the worst case scenario, with an MSE on the order of $10^{-4}$. For the 75th percentile case, the MSE is on the order of $10^{-5}$, indicating excellent numerical accuracy in the operator learning (OL) task. In both cases, the maximum error occurs near the boundary. As shown in Table 2, the training time for

PC

2 in this example is 30 s, demonstrating superior computational efficiency.

Figure 4 shows the UQ plots for this problem, including the mean and standard deviation over the test domain for both MCS and

PC

2, along with the corresponding pointwise absolute errors. It can be observed from Figure 4 that

PC

2 provides accurate estimates of mean and standard deviation across the entire domain, compared to MCS with low absolute errors. The maximum absolute errors in the mean and standard deviation occur at the boundary and remain relatively low, indicating high numerical accuracy over the given domain. The computation time for this problem is 32 s, compared to 1620 s for MCS, as shown in Table 2, highlighting the superior computational efficiency of the proposed method. The results clearly demonstrate that

PC

2 performs exceptionally well in both OL and UQ tasks.

3.3 1D Burgers’ equation example

This example highlights the ability of

PC

2 to learn nonlinear operators. Consider the following 1D viscous Burgers’ equation with spatially varying source term as

[TABLE]

where $\nu$ is the fluid viscosity and $f(x)$ is the source term.

In this example, we set $\nu=0.001$, and the source term $f(x)$ is drawn from a zero-mean GRF with an RBF kernel, with $\sigma=0.1$ and $\ell=0.2$. Six KL terms ($r=6$) are used to discretize the random field, accounting for $99\%$ of its variance. We set polynomial orders as $p=3$ and $q=23$. Following the approach in the previous example, we generate $N=100$ labeled training samples and $N=1000$ samples for testing. As shown in Table 1, PCE demonstrates exceptional performance in OL, achieving the MSE on the order of $10^{-7}$ with a training time of just a fraction of a second. Similarly, we get excellent performance in UQ, with MAE on the order of $10^{-5}$ for the mean and $10^{-4}$ for the standard deviation, at a fraction of the computational cost compared to MCS.

For

PC

2, we use 1,000 samples of $f(x)$ for training without using the labeled outputs, and the performance is tested on the same testing dataset as PCE. Figure 5 shows the OL results corresponding to the worst case and the 75th percentile case. As shown in Figure 5,

PC

2 yields highly accurate results, with a low worst-case MSE of $2.82\times 10^{-4}$. This highlights the strong capability of the proposed method in learning a nonlinear operator. Furthermore, the maximum error in both cases is observed near the right boundary at $t=0.3$, corresponding to the steepening phase of the solution where sharp gradients emerge.

Figure 6 shows the UQ plots, where the estimated mean and standard deviation obtained using

PC

2 are compared with their reference counterparts from MCS. It can be observed from Figure 6 that

PC

2 provides accurate estimates of the mean and standard deviation over the given domain, as evidenced by the absolute error plots, with a maximum error on the order of $10^{-2}$. To achieve this level of accuracy,

PC

2 requires only 80 s of computation time, compared to 1737 s for MCS, as shown in Table 2, highlighting its superior computational efficiency. These results demonstrate the superior performance of PCE in UQ.

3.4 2D Heat equation

Consider the following 2D Heat equation with a spatially varying source term and initial condition as

[TABLE]

where $\alpha$ is the thermal diffusivity of the medium and $f(x,y)$ is the source term.

Here, we consider $\alpha=0.01$ and the source term $f(x,y)$ is sampled from a zero-mean 2D GRF with an RBF covariance kernel, using $\sigma=1$ and $\ell=0.2$. The random field is discretized using a Karhunen–Loève (KL) expansion with 20 KL terms, capturing $99\%$ of the total variance. The amplitude ($A$) of the initial condition is modeled as a standard normal random variable. So, in total, we have 21 stochastic variables ($r=21$). We use a polynomial order of $p=4$ for the stochastic basis with a hyperbolic truncation factor of 0.9 [29], and set the spatial-temporal expansion order to $q=16$. For training, we use $N=2500$ samples, and for testing, $N=100$ samples are employed. As observed in Table 1, PCE provides high numerical accuracy and computational efficiency for OL. For UQ, we compute the MAE in mean and standard deviation with respect to MCS at $t=1$. The MAE is on the order of $10^{-16}$ for the mean and $10^{-4}$ for the standard deviation. The computation time is just 56 s compared to $\sim$27 h for MCS, demonstrating exceptional computational efficiency.

For

PC

2, we use 5,000 input samples without labeled solutions for training, and test the performance on the testing dataset as PCE. Figure 7 presents the OL results for a test sample corresponding to the 75th percentile error. We plotted the results for different time stamps, i.e., $t=0.25,0.5,\text{and}\ 1$. As observed, the MSE for this case is $7.83\times 10^{-4}$, which suggests high numerical accuracy. The error tends to increase with time, as evident from the absolute error plots with a maximum error of $1.09\times 10^{-1}$ for $t=1$. As given in Table 2, the total MSE is $6.2\times 10^{-4}$, and

PC

2 achieves this error with only 14 s of training time, which is exceptional considering high stochastic dimensions ($r=21$) with three physical variables. This highlights the scalability of the proposed method. Furthermore,

PC

2 achieves accuracy identical to standard PCE, but without the need for training data, which is computationally expensive to generate for this example. This highlights the superior performance of

PC

2, particularly for high-dimensional problems where acquiring labeled training data is computationally expensive.

Figure 8 shows the UQ plots for the given problem. We obtain the estimated mean and standard deviation at $t=1$ and compare them with their reference counterparts from MCS. It can be seen from the plots that

PC

2 provides excellent estimates of the mean and standard deviation, with maximum absolute errors of $3.59\times 10^{-5}$ for the mean and $5.7\times 10^{-3}$ for the standard deviation. The computation time for

PC

2 is just 116 s, compared to $\sim$27 h for MCS, as shown in Table 2, highlighting the significant computational efficiency..

4. Conclusion

In this work, we have proposed a polynomial chaos expansion (PCE) framework for operator learning in a data-driven and physics-informed setting. In the proposed framework, the task of operator learning reduces to solving a standard system of linear/nonlinear equations. Furthermore, the framework provides uncertainty quantification at no extra computational cost. We have demonstrated the effectiveness of the proposed method in OL and UQ-related tasks by solving linear/nonlinear PDE problems with various input functions, including source terms, parameters, and initial conditions. Numerical results indicate that the proposed method provides high accuracy at minimal computational cost and resource requirements, highlighting its efficiency, practicality, and scalability.

5. Acknowledgements

This work was supported by the Defense Threat Reduction Agency, Award HDTRA12020001. LN gratefully acknowledges the support of the Czech Science Foundation under project No. 24-11845S. The international collaboration was partly supported by the Ministry of Education, Youth and Sports of the Czech Republic under project No. LUAUS24260.