In this paper, we adopt conformal prediction, a distribution-free uncertainty quantification (UQ) framework, to obtain confidence prediction intervals with coverage guarantees for Deep Operator Network (DeepONet) regression. Initially, we enhance the uncertainty quantification frameworks (B-DeepONet and Prob-DeepONet) previously proposed by the authors by using split conformal prediction. By combining conformal prediction with our Prob- and B-DeepONets, we effectively quantify uncertainty by generating rigorous confidence intervals for DeepONet prediction. Additionally, we design a novel Quantile-DeepONet that allows for a more natural use of split conformal prediction. We refer to this distribution-free effective uncertainty quantification framework as split conformal Quantile-DeepONet regression. Finally, we demonstrate the effectiveness of the proposed methods using various ordinary, partial differential equation numerical examples, and multi-fidelity learning.

§ INTRODUCTION

Recent advancements in neural network technology have solidified its standing as a reliable and efficient tool for function approximation [1, 2]. This is particularly evident in scientific computing, where neural networks have demonstrated their power in approximating solutions to Partial Differential Equations (PDEs) or Ordinary Differential Equations (ODEs) [3, 4, 5, 6, 7]. A prevalent challenge in scientific computing involves approximating operators that map one function to another. Seminal research endeavors, such as those by Chen and Chen [8, 9], have explored the utilization of neural networks to approximate such operators, commonly referred to as neural operators.

In recent years, there has been a surge of interest in neural operators, with notable contributions from various scholars [10, 11, 12, 13, 14, 15, 16, 17]. Among these, two prominent approaches have emerged as particularly successful: Deep Operator Neural Networks (DeepONet) [10, 11, 18, 12, 19, 20] and Fourier Neural Operator (FNO) [21, 22, 14, 15, 16, 17, 23].

Compared to FNO, DeepONet represents a mesh-free neural operator, implying that the output functions do not require discretization. This characteristic enhances the flexibility of DeepONet in both training and testing phases [24, 25, 26, 11, 27, 28, 29, 22]. Furthermore, recent advancements in BelNet by [25, 26, 22] have extended DeepONet's capabilities, making it invariant to input function discretization. This means that input functions no longer need to conform to a shared discretization, thereby enhancing the versatility of DeepONet even further.

DeepONet and its extensions have been successfully applied to a diverse array of real-world applications. These include weather forecasting [23], sub-surface structure detection [13], electrical propagation on the left ventricle [30], geological carbon sequestration [31], disk-planet interactions in protoplanetary disks [32], hypersonic systems [33, 34], power systems [35], and optimization [36].

The need for Uncertainty Quantification (UQ) in scientific machine learning [37, 38] often arises from several factors. These include uncertainties in data due to measurement errors and numerical algorithm errors, uncertainties in model forms due to network architectures and physical models of varying fidelities, and uncertainties in parameters due to network training and system properties at different scales. As a result, researchers are tasked with providing confidence intervals for predicted outputs. A good interval should meet two fundamental criteria. Firstly, it should have a significant coverage rate, encompassing precise solutions or dynamics to the greatest extent possible. Secondly, the interval should be minimized in size.

Numerous works [21, 27, 37, 38, 39, 40, 41, 36, 42, 43] have established various uncertainty quantification (UQ) frameworks for DeepONet. These frameworks aim not only for accurate prediction but also for quantifying uncertainty. For example, a recent study [21] from the authors addresses this uncertainty issue in operator learning by framing the training process within a Bayesian framework. The study uses Langevin diffusion-based sampling methods [44, 45, 46, 47, 48] to generate ensembles, which aids in characterizing uncertainty in DeepONet predictions. To lessen the use of ensembles and increase efficiency, we introduced a probabilistic framework for UQ in DeepONets in [27], which provides an input-dependent standard deviation as a heuristic measure of uncertainty. However, despite the promising results from these methods, constructing rigorous confidence intervals for DeepONet predictions remains an unsolved challenge, which we aim to solve using conformal prediction.

Conformal prediction [49, 50, 51] is an alternative method for constructing confidence intervals that provide nonasymptotic, distribution-free coverage guarantees. Unlike the Bayesian Framework, this approach does not require prior distribution knowledge and avoids relying on strong assumptions of large-sample asymptotic approximations, which may be difficult to justify in practical scenarios. These characteristics make conformal prediction a valuable tool for addressing scientific machine learning problems, especially those arising from engineering applications with real observed data and a lack of prior knowledge.

There are two types of conformal prediction methods: full conformal prediction [51, 52] and split conformal prediction [51, 53]. In this paper, we focus on split conformal prediction, which operates with finite samples and enables building adaptive and validated confidence intervals. These properties have allowed split conformal prediction to become a popular approach for machine learning and uncertainty quantification tasks, including Language Modeling [54], Graph Neural Networks [55], time series [56], and quantile regression [57]. However, despite its potential, split conformal prediction remains underutilized in the scientific machine learning community, particularly in the context of Deep Operator Networks.

The goal of this paper is to address the challenge of uncertainty quantification (UQ) in operator learning. It aims to showcase the implementation of split conformal prediction for constructing confidence intervals and improving existing UQ algorithms tailored for DeepONet regression, with coverage guarantees.

To achieve this goal, we make the following contributions:

* We design a split conformal prediction algorithm (Section <ref>) that enhances the previously proposed uncertainty quantification frameworks, Bayesian DeepONet (B-DeepONet) and Probabilistic DeepONet (Prob-DeepONet). This novel algorithm allows B-DeepONet and Prob-DeepONet to generate confidence intervals with coverage guarantees.

* We then propose a new extension of DeepONet called Quantile-DeepONet (Section <ref>), which we use in a novel Conformal Quantile-DeepONet regression algorithm to construct rigorous confidence intervals.

* Finally, we demonstrate the effectiveness of our proposed conformalized-DeepONets by conducting multiple numerical experiments on ordinary differential equations, partial differential equations, and multi-fidelity settings.

The remaining sections of this paper are organized as follows. Section <ref> provides a review of the concepts of DeepONet, B-DeepONet, and Prob-DeepONet. In Section <ref>, we describe how to use split conformal prediction on B-DeepONet and Prob-DeepONet to build confidence intervals with coverage guarantees. Section <ref> presents the design of the Quantile DeepONet, which serves as the main component of the conformal Quantile-DeepONet regression algorithm for constructing rigorous confidence intervals. We present numerical experiments that validate the performance of all conformalized DeepONets in Section <ref>. Finally, Section <ref> discusses our results and future work, and Section <ref> concludes the paper.

§ BACKGROUND INFORMATION

This section reviews the Deep Operator Network (DeepONet) framework [10]. It also provides a summary of two uncertainty quantification (UQ) frameworks for DeepONet: Bayesian DeepONet (B-DeepONet [21]) and probabilistic DeepONet (Prob-DeepONet [27]).

§.§ Deep Operator Network (DeepONet)

The Deep Operator Network (DeepONet)  [10, 11, 12, 30] is an operator learning framework that aims to approximate the nonlinear operator $G$. The operator $G$ maps a function space $U$ with domain $K_1$ to a function space $V$ with domain $K_2$. DeepONet has been theoretically supported by universal approximation theorems [8, 26, 18, 58]. These theorems enable the approximation of the operator $G$ for a given $u \in U$ and $x \in K_2$ using the following linear trainable representation:

\begin{align*}
G(u)(x) \approx G_\theta(\hat{u})(x) = \sum_{k = 1}^{K} b_k(\hat{u})\tau_k(x),
\end{align*}

Here, $\theta$ is the vector of trainable parameters, while $b_k$ and $\tau_k$ are $K$ trainable coefficients and basis functions, respectively.

DeepONet consists of two sub-networks: the branch network and the trunk network (see Figure <ref>). The branch network maps the input function $u \in U$, which is discretized using $m$ sensors and denoted as $\hat{u}$, to a vector of $K$ trainable coefficients $b(\hat{u}) \in \mathbb{R}^K$. It is worth noting that various variants of DeepONet have been proposed to handle input functions with different discretizations. For instance, the new enhanced structures proposed in [25, 26] and the training algorithms proposed in [22].

Stacked version DeepONet $G_{\theta}$. $\bigotimes$ denotes the inner product in $\mathbb{R}^K$.

On the other hand, the trunk network maps a location $x \in K_2$, which belongs to the output function domain $K_2 \subset \mathbb{R}^d$, to a vector of $K$ trainable basis functions. This feature eliminates the need for discretization of the output function $G(u) \in V$ in DeepONet and its variants, making them mesh-free PDE approximate solvers.
To train DeepONet, one could optimize,

\begin{align*}
\mathcal{L}(\theta) = \frac{1}{N}\sum_{i = 1}^N\sum_{j = 1}^{N_j} \left |G^{(i)} - G_\theta \left(\hat{u}^{(i)}\right)\left(x^{(i)_j}\right) \right |^2.
\end{align*}

Without loss of generality and for better illustration of the algorithms, we assume $N_j = 1$, or $x^{(i)}\in K_2$ is the only evaluation point for one function $u^{(i)}$. That is, we have the dataset,
using the dataset of $N$ triplets: $\mathcal{D} = \left\{\hat{u}^{(i)}, x^{(i)}, G^{(i)}\right\}_{i=1}^N$, where $G^{(i)} = G\left(\hat{u}^{(i)}\right)\left(x^{(i)}\right)$ is the target operator value.

Despite the remarkable accuracy of DeepONet on multiple applications (e.g., [27, 23, 13, 30]), it only produces pointwise predictions. However, in cases where small training datasets or noisy inputs are involved, these pointwise predictions may be unreliable. Therefore, it is necessary to have some measure of the DeepONet's uncertainty. This requirement led the authors of this paper to develop two uncertainty quantification frameworks for DeepONet: Bayesian and probabilistic DeepONets, which we review next.

§.§ Bayesian Deep Operator Network (B-DeepONet)

This section reviews the first framework for uncertainty quantification (UQ) in DeepONet, known as Bayesian DeepONet (B-DeepONet) [21], developed by the authors of this paper. In particular, B-DeepONet is based on the Langevin theory that allows us to generate samples from an approximate posterior that converges to the Gibbs distribution proportional to $\exp(-U(\theta)/\tau)$, where $U$ is the non-convex energy function and $\tau >0 $ the temperature.

Measuring the uncertainty linked with constrained training datasets, noisy inputs, and over-parametrization of neural networks poses a considerable challenge. This challenge is even more pronounced in deep operator regression, as it involves mapping between functional spaces. To address this challenge, we proposed the Bayesian DeepONet (B-DeepONet) in [21]. B-DeepONet allows us to construct heuristic confidence intervals for the operator that maps an input function $u \in V$ to an output function $G(u)(x) \in \mathbb{R}$ evaluated at a point $x \in K_2 \subset \mathbb{R}^d$.

In B-DeepONet, our goal was to develop a method to sample from the predictive distribution $p(G|(u,x),\mathcal{D})$ of the target operator $G(u)(x)$, given a discretized input $\hat{u} \in \mathbb{R}^m$ at any $x \in K_2$ using the training dataset $\mathcal{D}$. To reach this goal, we first assumed a factorized Gaussian likelihood function for the dataset [44, 59]:

p(G|(\hat{u},x),\theta) = \mathcal{N}(G|G_\theta(\hat{u})(x), \text{diag}(\Sigma^2)) = \prod_{i=1}^N \mathcal{N}(G^{(i)}|G_\theta(\hat{u}^{(i)})(x^{(i)}), \sigma),

where the mean $G_\theta(\hat{u})(x)$ is the output of the underlying DeepONet with a vector of trainable parameters $\theta$ and $\text{diag}(\Sigma^2)$ a diagonal covariance matrix with $\Sigma^2 = (\sigma^2,\ldots,\sigma^2)$ on the diagonal [37], which we assume is given.

It is important to note that the target operator value $G(u)(x) \in \mathbb{R}$ for a given input function $u \in U$ at any location $x \in K_2$ and the training dataset $\mathcal{D}$ is the random variable $(G|(u,x),\mathcal{D})$. To calculate the density of this random variable, the DeepONet parameters $\theta$ need to be integrated out, as shown below:

p(G|(u,x), \mathcal{D}) = \int p(G|(u,x),\theta)p(\theta | \mathcal{D})d\theta.

In the above, $p(\theta | \mathcal{D})$ is the posterior distribution of the DeepONet trainable parameters. Using Bayes' rule, we can write the posterior of DeepONet parameters as

p(\theta|\mathcal{D}) \propto p(\mathcal{D}|\theta) p(\theta),

where $p(\theta)$ is the prior distribution of the parameters, and $p(\mathcal{D}|\theta)$ is the data likelihood, that is, $p(\mathcal{D}|\theta) = \prod_{i=1}^N p(G^{(i)}|(\hat{u}^{(i)}, x^{(i)}), \theta)$. We calculate this using the DeepONet forward pass and the i.i.d. training dataset $\mathcal{D}$.

Computing the posterior distribution in closed form using Bayes' rule is intractable [37] for deep operator networks. Thus, to generate samples from this distribution, we employ a Langevin diffusion-based sampling technique [47, 60]. This technique has been shown to produce samples that closely match the target distribution [61, 62]. Specifically, the Langevin diffusion process used to generate samples is as follows:

\begin{align*}
d\theta_t = -\nabla U(\theta_t) dt+\sqrt{2\tau}dW_t,
\end{align*}

where $\theta$ is the vector of parameters, $U$ is the energy function (i.e., the log data likelihood and prior), $\tau$ the temperature, and $W_t$ is the Brownian motion. By sampling from the posterior distribution using a discretized Langevin process, we obtained an ensemble of $M$ $\theta$ samples, denoted as $\{\theta_k\}_{k=1}^M$, as described next.

To obtain the ensemble of parameters, denoted as $\{\theta_k\}_{k=1}^M$, B-DeepONet uses a variation of the diffusion process mentioned earlier. Particularly, the authors consider the parallel tempering techniques.  These methods involve two chains and facilitate particle swapping between them, aiding in escaping local traps and addressing highly non-convex problems, thereby enhancing its capability to tackle complex UQ challenges. However, employing two chains doubles the cost of energy function evaluation. Despite efforts by the authors to mitigate these costs, achieving only up to a $25\%$ reduction in the doubled cost, this prompts the exploration of alternative UQ methods for managing uncertainty in neural operators.

To construct the confidence intervals, in our paper [21], we used the obtained $M$-ensemble $\{\theta_k\}_{k=1}^M$ of sampled parameters from the posterior distribution to fit a parametric predictive distribution, a common practice in the literature [37]. It is important to note that this assumption limits the applicability of these confidence intervals. Specifically, we used the Gaussian distribution $\mathcal{N}(\mu(\hat{u})(x), \sigma(\hat{u})(x))$ for an arbitrary $x$ location. The parameters of this distribution were obtained as follows:

\begin{align*}
\mu(\hat{u})(x) &= \frac{1}{M} \sum_{k=1}^M G_{\theta_k}(\hat{u})(x) \\
\sigma^2(\hat{u})(x) &= \frac{1}{M}\sum_{k=1}^M \left(G_{\theta_k}(\hat{u})(x) - \mu(\hat{u})(x)  \right)^2.
\end{align*}

Due to the assumption of fitting the parametric predictive distribution, B-DeepONet can only provide a heuristic estimate of uncertainty. In our paper [21], we demonstrated that constructing confidence intervals with B-DeepONet requires extensive hyperparameter optimization to achieve a high confidence level (e.g., 95%) for a specific test dataset. However, there is no guarantee that these heuristic confidence intervals, developed using B-DeepONet, will be reliable for other test examples. Additionally, constructing these heuristic confidence intervals requires the forward pass of $M$ DeepONets, making B-DeepONet a computationally expensive method.

§.§ Probabilistic Deep Operator Network (Prob-DeepONet)

One of the main drawbacks of quantifying uncertainty using B-DeepONet is its cost. B-DeepONet requires performing a forward pass across all the DeepONets within the collected ensemble of size $M$. To address this drawback, the authors of this paper developed the Probabilistic DeepONet (Prob-DeepONet) [27], which we will review in this section.

Prob-DeepONet is a probabilistic model whose output approximates the parameters of a normal distribution:

$$G|X=(\hat{u}, x) \sim \mathcal{N}\left(\mu_{\theta_{\mu}}(\hat{u})(x), \sigma_{\theta_{\sigma}}(\hat{u})(x) \right).$$

In the above, $\mu_{\theta_{\mu}}$ represents the mean of the probabilistic model and aims to approximate $G(u)(x) \in \mathbb{R}$, which is the target operator for a given input $u \in U$ at $x \in K_2$. This approximation is achieved using the following linear trainable representation:

$$\mu_{\theta_\mu}(\hat{u})(x) = \sum_{i=1}^K b^{\mu}_k(\hat{u})\cdot \tau_k^{\mu}(x).$$

Similarly, $\sigma_{\theta_{\sigma}}$ represents the standard deviation of the probabilistic model. It aims to estimate the uncertainty for $G(u)(x)$ through the following linear trainable representation:

$$\log(\sigma_{\theta_\sigma}(\hat{u})(x)) = \sum_{i=1}^K b^{\sigma}_k(\hat{u})\cdot \tau_k^{\sigma}(x).$$

Note that in the above, for numerical stability, we train the model to produce $\log \sigma$ and then recover the standard deviation $\sigma$ using the exponential function.

Similar to the classical DeepONet, Prob-DeepONet consists of two sub-networks: the branch network and the trunk network. In the branch network, we design two components: one for the mean and the other for the standard deviation. The branch component for the mean includes shared trainable layers and a few independent layers. It maps the discretized input $\hat{u} \in \mathbb{R}^m$ to a vector of $K$ trainable coefficients $b^{\mu}(\hat{u}) \in \mathbb{R}^K$ for the mean. Similarly, the branch component for the log standard deviation shares several layers with the mean component. It also has independent layers to process the map from the discretized input $\hat{u} \in \mathbb{R}^m$ to the vector of trainable coefficients $b^{\sigma}(\hat{u}) \in \mathbb{R}^K$.

On the other hand, the trunk network also has two components: one for the mean and another for the log standard deviation. Each component consists of a set of shared layers and a few independent layers. The trunk component for the mean maps the location $x \in K_2$ within the output function domain $K_2 \subset \mathbb{R}^d$ to the vector of $K$ trainable basis functions $\tau^{\mu}(x) \in \mathbb{R}^K$. Similarly, the trunk component for the log standard deviation maps the location $x \in K_2$ to the vector of trainable basis functions $\tau^{\sigma}(x) \in \mathbb{R}^K$.

To train the Prob-DeepONet parameters $\theta = \{\theta_{\mu},\theta_{\sigma}\}$, we minimize the negative log likelihood:

$$\mathcal{L}_{\text{nll}}(\theta) = -\frac{1}{2N} \left(\sum_{i=1}^N \left( \frac{ \left(G^{(i)} - \mu_{\theta_{\mu}}(\hat{u}^{(i)})(x^{(i)}) \right)^2}{\sigma_{\theta_{\sigma}}(\hat{u}^{(i)})(x^{(i)})^2} + 2 \log\sigma_{\theta_{\sigma}}(\hat{u}^{(i)})(x^{(i)})\right) + N\log 2 \pi\right),$$

using the dataset of $N$ triplets: $\mathcal{D} = \left\{\hat{u}^{(i)}, x^{(i)},G^{(i)} \right\}_{i=1}^N$ .

Similar to B-DeepONet, Prob-DeepONet only provides a heuristic estimate of uncertainty through the standard deviation. As demonstrated in our paper, constructing confidence intervals using Prob-DeepONet requires heavy hyper-parameter optimization to achieve a high confidence level (e.g., 95%) for a fixed test dataset. However, it is not guaranteed that these heuristic confidence intervals, developed using Prob-DeepONet and the estimated standard deviation, are reliable for other test examples. Therefore, the first objective of this paper is to enhance our previously proposed methods to transform these heuristic confidence intervals into rigorous prediction intervals. To accomplish this, we will employ split conformal prediction [53, 49], which will be described in the following section.

§ CONFORMAL PREDICTION FOR DEEP OPERATOR NETWORK REGRESSION

In this section, we will first demonstrate how to use split conformal prediction [53, 49] to create reliable confidence intervals using the heuristic uncertainty estimates $\sigma(\hat{u})(x)$ obtained from a trained Prob-DeepONet and the $M$-ensemble generated by B-DeepONet. Next, we will propose a novel Quantile-DeepONet extension that is better suited for developing a conformalized Quantile-DeepONet regression algorithm.

§.§ Split Conformal Prediction for DeepONet

Within this section, we employ split conformal prediction [53, 49], recognized as the predominant conformal prediction approach, facilitating Prob- and B-DeepONet in delivering robust confidence intervals. Formally, given any DeepONet test triplet $(\hat{u}_{\text{test}}, x_{\text{test}}, G_{\text{test}})$, our goal is to construct a marginal distribution-free confidence interval $\mathcal{C}(\hat{u}_{\text{test}}, x_{\text{test}})$ that likely contains the operator target $G_{\text{test}}$. That is, given a miscoverage rate $\alpha \in (0,1)$, we aim to achieve:

\begin{align} \label{eq:coverage-prob}
\mathbb{P}\{G_{\text{test}} \in \mathcal{C}(\hat{u}_{\text{test}}, x_{\text{test}})\} \ge 1-\alpha.
\end{align}

To achieve this goal, split conformal prediction (see Algorithm <ref>) proceeds as follows.
Split Conformal Prediction to Enhance B- and Prob-DeepONet

Input: Calibration data: $\{\hat{u}^{(i)}, x^{(i)},G^{(i)}\}_{i=1}^{n}$ where $G^{(i)} = G(u^{(i)})(x^{(i)})$ is the exact output function values, miscoverage level: $\alpha \in (0,1)$, and mean and standard deviation models: $\mu$ and $\sigma$, obtained from the $M$-ensemble of B-DeepONet or a trained Prob-DeepONet.
Use $\mu$ and $\sigma$ models to compute the calibration scores:

$$s_1 = s\left(\hat{u}^{(1)}, x^{(1)}, G^{(1)}\right), \dots, s_n = s\left(\hat{u}^{(n)}, x^{(n)}, G^{(n)}\right),$$

$$s(\hat{u},x,G) = \frac{|G - \mu(\hat{u} )(x)|}{\sigma( \hat{u} )(x)}.$$

Compute $\hat{q}$ as the $\frac{\lceil (n+1) (1-\alpha) \rceil}{n}-$th quantile of the calibration scores $s_1, \dots, s_n$.
Return: Rigorous confidence intervals (i.e., that satisfy (<ref>)) for any new test example $(\hat{u}_\text{test},x_{\text{test}})$:

$$\mathcal{C}(\hat{u}_\text{test},x_{\text{test}}) = \left[\mu(u_\text{test})(x_\text{test}) - \hat{q}\sigma(u_\text{test})(x_\text{test}),\mu(u_\text{test})(x_\text{test}) + \hat{q}\sigma(u_\text{test})(x_\text{test})\right].$$

Trained DeepONet model. First, we have access to trained B- or Prob-DeepONet models. These models predict the target operator $G_{\text{test}}$ by using the mean output $\mu(\hat{u}_{\text{test}})(x_{\text{test}})$ and provide an estimate of the uncertainty through the standard deviation output $\sigma(\hat{u}_{\text{test}})(x_{\text{test}})$. It is important to note that, as described in our previous works [21, 27], the standard deviation $\sigma$ produces a larger value when there is increased uncertainty about the test input.

Score function. Then, given $\mu$  and $\sigma$, we define a score function for any DeepONet data triplet $(\hat{u}, x, G)$:

\begin{align} \label{eq:score-function}
s(\hat{u},x,G) = \frac{|G - \mu(\hat{u})(x)|}{\sigma(\hat{u})(x)}.
\end{align}

Note that this score function acts as a correction factor for our DeepONet measure of uncertainty, i.e., $s(\hat{u},x,G) \cdot \sigma(\hat{u})(x) = |G - \mu(\hat{u})(x)|$.

Calibration quantile. Using a given calibration dataset of size $n$, denoted as $\{\hat{u}^{(i)}, x^{(i)}, G^{(i)}\}_{i=1}^n
$, we compute the calibration scores $s_1 = s(\hat{u}^{(1)}, x^{(1)}, G^{(1)}), \dots, s_n = s(\hat{u}^{(n)}, x^{(n)}, G^{(n)})$. Next, we calculate $\hat{q}$ as the $\frac{\lceil (n+1) (1-\alpha) \rceil}{n}$ quantile of the calibration scores $s_1, \dots, s_n$, where $\lceil \cdot \rceil$ denotes the ceil function. It is important to note that the calibration dataset must satisfy the exchangeability property [49], which is true in most DeepONet settings where the training, test, and calibration data are independently and identically distributed (i.i.d.).

Confidence intervals with coverage guarantees. We use the quantile $\hat{q}$ to construct rigorous confidence intervals for any new DeepONet test example $(\hat{u}_\text{test},x_{\text{test}})$:

\begin{align} \label{eq:confidence-intervals}
\mathcal{C}(\hat{u}_{\text{test}},y_{\text{test}}) = \{G:s(\hat{u}_{\text{test}},y_{\text{test}},G) \le \hat{q} \} \equiv \{G:|G - \mu(\hat{u}_\text{test})(x_\text{test})| \le \hat{q}\sigma(\hat{u}_\text{test})(x_\text{test}) \}.
\end{align}

These confidence intervals hold for any distribution of the data and are rigorous in the sense that they satisfy property (<ref>), which was demonstrated in the conformal coverage guarantee theorem <ref> replicated next for completeness.

Conformal calibration coverage guarantee theorem [50, 49]. Suppose $\{X_i,Y_i\}_{i=1}^n$ and $(X_\text{test}, Y_\text{test})$ are i.i.d. Then define $\hat{q}$ as

\hat{q} = \inf \left\{ q: \frac{|i:s(X_i,Y_i) \le q|}{n} \ge \frac{\lceil (n+1) (1-\alpha) \rceil}{n} \right\}

and the resulting prediction sets as

\mathcal{C}(X_\text{test}) = \{y:s(X_\text{test},y) \le \hat{q}\}.

$$\mathbb{P}\{ Y_\text{test} \in \mathcal{C}(X_\text{test})\} \ge 1 - \alpha.$$

We conclude this section by noting that the above theorem holds for any score function $s$.

§.§ Conformal Quantile DeepONet Regression

In the previous section, we demonstrated how to construct finite-sample, distribution-free, and adaptive confidence intervals of the form $\mu(\hat{u})(x) \pm\hat{q}\sigma(\hat{u})(x)$ for operator targets $G$ using the standard deviation $\sigma$ as a measure of uncertainty for DeepONets. We will show in Section <ref> that this approach effectively quantifies uncertainty in DeepONets. However, as stated in [49], there is no evidence to suggest that the heuristic measure of uncertainty $\sigma$ for B- and Prob-DeepONet is strongly correlated with the quantiles of the operator target distribution. Therefore, it is more natural to expect that a DeepONet method that estimates the quantiles for the target operator $G$ distribution will enhance the use of split conformal predictions. Let us now introduce this novel Quantile-DeepONet.

Quantile-DeepONet. Our goal is to design a Quantile-DeepONet that can approximate the conditional quantiles $t_{\alpha/2}$ and $t_{1-\alpha/2}$, given (i) a miscoverage level $\alpha \in (0,1)$, (ii) a discretized input function $\hat{u} \in \mathbb{R}^m$, and (iii) a location within the output function domain $x \in K_2$. To provide a complete understanding, let us recall the mathematical definition of a conditional quantile function. The $\gamma$th conditional quantile is defined as follows:

t_\gamma(x) = \inf \{y\in \mathbb{R}:F(y|X=x) \ge \gamma\},

F(y|X=x) = \mathbb{P}\{Y \leq y|X=x\},

is the conditional distribution of $Y$ given $X=x$. In other words, $t_\gamma$ provides us with information about the percentile of the distribution of $Y$ when conditioned on $X=x$.

Similar to Prob-DeepONet, Quantile-DeepONet approximates quantiles using the following linear trainable representations:

\begin{align*}
t_{\gamma}(\hat{u})(x) \approx \hat{t}_{\theta_\gamma}(\hat{u})(x) = \sum_{k = 1}^K b_k^{\gamma}(\hat{u})\cdot \tau_k^{\gamma}(x),
\end{align*}

for $\gamma \in \{\alpha/2,1-\alpha/2\}$. Therefore, Quantile-DeepONet consists of two sub-networks: the branch network and the trunk network.

The branch network comprises two components: one for the $(\alpha/2)$th quantile and another for the $(1-\alpha/2)$th quantile. Each component consists of a set of shared trainable layers and a few independent trainable layers. Additionally, each component maps the discretized input $\hat{u} \in \mathbb{R}^m$ to a vector of $K$ trainable coefficients $b^\gamma(\hat{u}) \in \mathbb{R}^K$, where $\gamma \in \{\alpha/2, 1-\alpha/2\}$. Similarly, the trunk network consists of two components, one for each quantile. These two trunk components also have shared and independent trainable layers. They map the location $x \in K_2$ to a vector of $K$ basis functions $\tau^\gamma(x) \in \mathbb{R}^m,$ where $\gamma \in \{\alpha/2,1-\alpha/2\}$.

We can train the Quantile-DeepONet parameters $\theta = \{\theta_{\alpha/2},\theta_{1-\alpha/2}\}$ by minimizing the following loss function for $\gamma \in \{\alpha/2, 1-\alpha/2\}$:

\mathcal{L}(\theta) = \frac{1}{N} \sum_{i=1}^N \rho_\gamma\left( G^{(i)},\hat{t}_{\theta_{\gamma}}\left(\hat{u}^{(i)}\right)\left(x^{(i)}\right)\right),

on the dataset of $N$ triplets: $\mathcal{D} = \left\{\hat{u}^{(i)}, x^{(i)},G^{(i)} \right\}_{i=1}^N$. In the above expression, the loss function $\rho_\gamma$ refers to the pinball loss [63, 64, 57]:

\rho_\gamma(y,\hat{y}) = \begin{cases} \gamma(y-\hat{y}) \qquad \qquad \text{if } y-\hat{y}>0, \\  (1-\gamma)(\hat{y}-y) \quad~~ \text{otherwise}.\end{cases}

Once trained, a Quantile-DeepONet can provide an estimate of the confidence interval for any new test sample:

\hat{C}(\hat{u}_\text{test},x_\text{test}) = \left[\hat{t}_{\theta^*_{\alpha/2}}\left(\hat{u}_\text{test}\right)\left(x_\text{test}\right), \hat{t}_{\theta^*_{1-\alpha/2}}\left(\hat{u}_\text{test}\right)\left(x_\text{test}\right) \right],

where $\theta^*_{\alpha/2}$ and $\theta^*_{1-\alpha/2}$ denote the optimizers of the corresponding losses.
However, the estimate provided above fails to meet the coverage probability requirement stated in (<ref>). To satisfy this requirement, we need to incorporate conformal prediction. It is worth noting that the combination of quantile and conformal prediction is referred to as Conformal Quantile Regression (CQR) in the literature [57], and we will now describe its application to Quantile-DeepONets.

Conformal Quantile-DeepONet Regression. The details of the split conformal Quantile-DeepONet regression procedure are provided in Algorithm <ref>. This procedure begins by defining a suitable score function.
Split Conformal Quantile DeepONet Regression

Input: Calibration data: $\{\hat{u}^{(i)}, x^{(i)},G^{(i)}\}_{i=1}^{n}$, miscoverage level: $\alpha \in (0,1)$, and trained quantile-DeepONet outputs $\hat{t}_{\alpha/2}$ and $\hat{t}_{1-\alpha/2}$.
Use the trained quantile-DeepONet to compute the calibration scores:

$$s_1 = s\left(\hat{u}^{(1)}, x^{(1)}, G^{(1)}\right), \dots, s_{n} = s\left(\hat{u}^{(n)}, x^{(n)}, G^{({n})}\right),$$

$$s(u,x,G) = \max \left\{\hat{t}_{\theta^*_{\alpha/2}}(u)(x) - G, G- \hat{t}_{\theta^*_{1-\alpha/2}}(u)(x)\right\}.$$

Compute $\hat{q}$ as the $\frac{\lceil (n+1) (1-\alpha) \rceil}{n}$ quantile of the calibration scores $s_1, \dots, s_{n}$.
Return: Rigorous confidence intervals for any new test example $(\hat{u}_\text{test},x_{\text{test}})$:

$$\mathcal{C}(\hat{u}_\text{test},x_{\text{test}}) = \left[\hat{t}_{\theta^*_{\alpha/2}}(\hat{u}_\text{test})(x_\text{test})-\hat{q}, \hat{t}_{\theta^*_{1-\alpha/2}}(\hat{u}_\text{test})(x_\text{test}) + \hat{q}\right].$$

Score function. Given the trained Quantile-DeepONet, we compute the score for any DeepONet data triplet $(\hat{u}, x, G)$ as follows:

s(u,x,G) = \max \left\{\hat{t}_{\theta^*_{\alpha/2}}(u)(x) - G, G- \hat{t}_{\theta^*_{1-\alpha/2}}(u)(x)\right\}.

Note that if the operator target $G$ satisfies $G<\hat{t}_{\theta^*_{\alpha/2}}(u)(x)$ or $G> \hat{t}_{\theta^*_{1-\alpha/2}}(u)(x)$, then the score represents the amount of error incurred by Quantile-DeepONet.

Calibration quantile. As before, using a given calibration i.i.d dataset of size $n$, denoted as $\{\hat{u}^{(i)}, x^{(i)}, G^{(i)}\}_{i=1}^N$, where $G^{(i)} = G(\hat{u}^{(i)})(x^{(i)})$, we compute the calibration scores $s_1 = s(\hat{u}^{(1)}, x^{(1)}, G^{(1)}), \dots, s_n = s(\hat{u}^{(n)}, x^{(n)}, G^{(n)})$. Then, we calculate $\hat{q}$ as the $\frac{\lceil (n+1) (1-\alpha) \rceil}{n}$  quantile of the calibration scores $s_1, \dots, s_n$.

Confidence intervals with coverage guarantees. We use the quantile $\hat{q}$ to construct rigorous confidence intervals for any new DeepONet test example $(\hat{u}_\text{test},x_{\text{test}})$:

\begin{align} \label{eq:confidence-intervals-CQR}
\mathcal{C}(\hat{u}_\text{test},x_{\text{test}}) = \left[\hat{t}_{\theta^*_{\alpha/2}}(\hat{u}_\text{test})(x_\text{test})-\hat{q}, \hat{t}_{\theta^*_{1-\alpha/2}}(\hat{u}_\text{test})(x_\text{test}) + \hat{q}\right].
\end{align}

These confidence intervals are adaptive, distribution-free, and satisfy the coverage probability property (<ref>). The validity of this property was demonstrated in [57] by extending theorem <ref> to the case of CQR.

§ NUMERICAL EXPERIMENTS

In this section, we assess the ability of conformalized-DeepONets (conformal prediction combined with B-DeepONet, Prob-DeepONet, or Quantile-DeepONet) to generate confidence intervals with coverage guarantees. We illustrate this through three experiments: (i) the nonlinear pendulum (Section <ref>), (ii) the diffusion-reaction system <ref>, and (iii) the viscous Burgers' equation <ref>.

Datasets. To train our DeepONet models, we used a dataset consisting of $N$ samples, each in the form of a triplet $\left(\hat{u}^{(i)},x^{(i)},G^{(i)}\right)$. The specific size of the training dataset, $N$, for each experiment is provided in Table <ref>. For the calibration dataset, we used $n = \frac{N}{10}$ samples, with details outlined in Table <ref>. Finally, for testing purposes, we evaluated all our DeepONet models using $100$ trajectories per experiment. Each trajectory included the target operator values $G(u)(X_\text{test})$ evaluated on a mesh $X_\text{test} \subset K_2$ composed of $100$ uniformly distributed location points within the output function domain $K_2$.
It should be noted that the testing mesh size is much more dense than the training mesh, i.e., our experiments also demonstrate the extrapolation of the proposed methods.

Experiment      # $N$ training      # $n$ calibration

The nonlinear pendulum equation    5000     500

The diffusion-reaction system     10000     1000

The viscous Burgers' equation       30000     3000

Summary of the number of training ($N$) and calibration ($n$) data samples for each experiment. Each sample is a triplet consisting of $\left(\hat{u}^{(i)},x^{(i)},G^{(i)}\right)$.

Neural Networks. We used the feed-forward neural network (FNN) architecture to build the branch and trunk sub-networks. The corresponding width and depth for each experiment, including the shared and independent layers for Prob- and Quantile-DeepONet, are listed in Table <ref>. We used the ReLU function as the activation function for each input and hidden layer.

Experiment      Trunk depth     Trunk width     Branch depth     Branch width

The nonlinear pendulum    3(1)     100     3(1)     100

The diffusion-reaction    4(1)     100     4(1)     100

The viscous Burgers'    5(1)     128     5(1)     128

Summary of the width and depth of the trunk and branch feed-forward neural networks used in each experiment. For the depth of Prob- and Quantile-DeepONets, we use the notation shared(independent) layers.

Optimizer and Sampler. For Prob- and Quantile-DeepONet, we used the Adam optimizer [65] with default hyperparameters and a learning rate $\eta = 10^{-3}$. We adjusted the learning rate using the reduced on plateau scheduler. For B-DeepONet, we used the stochastic gradient replica-exchange sampling algorithm [66, 21, 67] with standard hyperparameters.

Baselines. We compared the proposed Conformalized-DeepONets (which combine conformal prediction with B-DeepONet, Prob-DeepONet, or Quantile-DeepONet) against the following baselines: B-DeepONet [21], Prob-DeepONet [27], and Quantile-DeepONet (Section <ref>). Specifically, for each conformalized and baseline model, we calculated the coverage over the test dataset. This coverage is determined by

C_j = \frac{1}{n_\text{eval}} \sum_{i=1}^{n_\text{eval}} \mathbbm{1} \left \{ G^{(i,j)} \in \mathcal{C}_j \left(\hat{u}^{(j)}, x^{(i,j)} \right) \right\}, \text{ for } j=1,\dots,n_\text{traj}.

In other words, we calculate the coverage of $n_\text{traj}$ test trajectories and each $C_j$ is the ratio of location points whose reference solutions are in the established confidence interval. These are evaluated at $n_{\text{eval}}=100$ uniformly spaced points within the output function domain.

§.§ Experiment: The Nonlinear Pendulum

In this experiment, we consider the nonlinear pendulum subject to external forcing, represented by:

\begin{align*}
\frac{ds_1}{dt} &=s_2, \\
\frac{ds_2}{dt} &=-k \sin s_1 + u(t),
\end{align*}

with the initial condition $(s_1(0),s_2(0))^\top = (0,0)^\top$ and time domain $[0,1.0]$ (s). In the above, $(s_1(t),s_2(t))^\top$ is the state, $k=1.0$ is a constant determined by the acceleration due to gravity and the length of the pendulum, and $u(t)$ is the time-dependent external forcing with time domain $[0,1]$ (s).

Our aim is to approximate the solution operator $G:u(t) \mapsto s_1(t)$, which maps the external forcing input to the state $s_1$. We draw samples of the external forcing inputs for training, calibration, and testing from the mean-zero Gaussian Random Field (GRF) $u \sim \mathcal{G}(0,k_\ell(t_1, t_2))$. This GRF has a radial-basis function covariance kernel $k_\ell(t_1,t_2)$ with a length-scale of $\ell=0.1$. Finally, we discretize each sampled input $u$ using $m=100$ sensors.

Checking Coverage. We evaluated the ability of all conformalized-DeepONet models and baseline models to generate confidence intervals (with misscoverage rate $\alpha=0.05$) with coverage guarantees for trajectories in the test dataset. Figure <ref> illustrates the confidence intervals for a test trajectory randomly selected from the test dataset. Furthermore, Table <ref> presents the average coverage achieved by both conformalized-DeepONets and baseline models for all trajectories in the test dataset.

[ xlabel=$t$, ylabel=$G(u)(t) \equiv s_1(t)$,xmin=0, xmax=1,ymin=-0.07, ymax=0.01, legend pos=south west, line width=1pt,enlargelimits=false]

[solid, color=blue]
Reference Solution
[name path=prob_lower,  fill=none, draw=none,forget plot]
[name path=prob_upper, fill=none, draw=none,forget plot]

[name path=cop_lower, fill=none, draw=none,forget plot]
[name path=cop_upper, fill=none, draw=none,forget plot]

[fill=mycolor1, fill opacity=0.5] fill between[of=cop_upper and cop_lower, ];

[fill=red, fill opacity=0.5] fill between[of=prob_upper and prob_lower, ];

Conformal Prediction and Prob-DeepONet

[ xlabel=$t$, ylabel=$G(u)(t) \equiv s_1(t)$,xmin=0, xmax=1,ymin=-0.07, ymax=0.01, legend pos=south west, line width=1pt,enlargelimits=false]
[solid, color=blue]
Reference Solution
[name path=quantile_lower,  fill=none, draw=none,forget plot]
[name path=quantile_upper, fill=none, draw=none,forget plot]

[name path=coq_lower, fill=none, draw=none,forget plot]
[name path=coq_upper, fill=none, draw=none,forget plot]

[fill=mycolor1, fill opacity=0.5] fill between[of=coq_upper and coq_lower, ];

[fill=red, fill opacity=0.5] fill between[of=quantile_upper and quantile_lower, ];

Conformal Prediction and Quantile-DeepONet

Confidence intervals for a random test trajectory of the nonlinear pendulum experiment given a miscoverage rate $\alpha=0.05$. (a) Confidence intervals for the Conformalized Prob-DeepONet and the baseline Prob-DeepONet. (b) Confidence intervals for the Conformalized Quantile-DeepONet and the baseline Quantile-DeepONet. Please refer to Table <ref> for the summary of the average coverage rate for all testing trajectories with different approaches.

DeepONet model    Average coverage % ($\alpha=0.05$)

Conformalized Prob-DeepONet     94.69%

Conformalized Quantile-DeepONet     95.86%

Prob-DeepONet     3.80%

Quantile-DeepONet     55.31%

The average coverage percentage for the predicted confidence intervals (with a miscoverage rate of $\alpha=0.05$) calculated based on the 100 test trajectories in the nonlinear pendulum experiment.

Figure <ref> and Table <ref> offer significant insights. Firstly, the conformalized DeepONets (both Probabilistic and Quantile) deliver the desired coverage on average, with a miscoverage rate of $\alpha=0.05$. Moreover, the Quantile-DeepONet significantly outperforms the Probabilistic DeepONet. Lastly, the Probabilistic DeepONet displays the poorest average coverage. These findings corroborate the conclusions from our previous paper, which highlighted the necessity of extensive hyperparameter optimization to improve the quality of confidence intervals for Prob-DeepONet.

Finally, Figure <ref> shows the distribution of coverages (with a miscoverage rate of $\alpha=0.05$) for all confidence intervals obtained by applying the conformalized Quantile-DeepONet to all $100$ test trajectories of the nonlinear pendulum example. As anticipated, the coverages concentrate around $1-\alpha$, that is, with 95% coverage. Next, we present two ablation studies showing the adaptivity and effect of the calibration dataset size.

[width=*0.7, xlabel=Predicted confidence intervals' coverage (%), ylabel=Conuts, ymax=400,
[black, dashed] coordinates (0.95,0) (0.95,350);

] table [y index=0] figsData/statsData/coverages.txt;
[black, dashed] coordinates (0.95,0) (0.95,450);

The distribution of predicted confidence intervals' coverages (in percentages), obtained by using the conformalized Quantile-DeepONet on the 100 test trajectories of the nonlinear pendulum experiment, with a specified miscoverage rate of $\alpha=0.05$.

Ablation Study I: Adaptivity of Confidence Intervals. We tested the proposed conformalized Prob-DeepONet for its ability to provide adaptive confidence intervals – that is, confidence intervals of varying lengths for each $x \in K_2$, while still ensuring coverage guarantees.
The adaptive coverage will reduce the over-confidence issue of many algorithms.
To accomplish this, we applied the conformalized Prob-DeepONet to all $100$ test trajectories of the pendulum experiment, with the goal of analyzing the distribution of the confidence intervals' lengths. Figure <ref> illustrates the distribution of the lengths of confidence intervals. This distribution demonstrates that the confidence intervals are adaptive, meaning their length varies for each $x \in K_2$. If this wasn't the situation, the distribution would be focused around a specific length. Clearly, our experiment indicates that this is not the scenario.

[width=*0.7, xlabel=Confidence intervals' lengths ($\alpha=0.05$), ylabel=Counts,
xtick=0, 0.02, 0.04, 0.06, 0.08, 0.1, 0.12, 0.14,

xticklabels=0, $0.02$, $0.04$, $0.06$, $0.08$, $0.1$, $0.12$, $0.14$,

] table [y index=0] figsData/statsData/set_sizes.txt;

The distribution of the lengths of the confidence intervals (given a miscoverage rate $\alpha=0.05$) predicted using the proposed conformalized Prob-DeepONet for the $100$ test trajectories of the nonlinear pendulum experiment.

Ablation Study II: Testing the Effect of the Calibration Dataset Size $n$. We tested the effect of varying the calibration dataset size $n$ on the distribution of coverages for the nonlinear pendulum example and a miscoverage rate of $\alpha=0.05$.

To this end, we varied $n$ over the set $\{500, 1000, 5000, 10000\}$. For each $n$, we computed the distribution of coverages. This was done for $R$ rounds ($R=200$), where we generated a calibration dataset of size $n$ and a validation dataset of size $n_\text{val}$. Each dataset contains DeepONet triplets of the form $(\hat{u}, x, G)$. We then calculated the coverages as follows.

C_k = \frac{1}{n_\text{val}} \sum_{i=1}^{n_\text{val}} \mathbbm{1} \left\{ G^{(i,k)} \in \mathcal{C}_k\left(\hat{u}^{(i,k)}, x^{(i,k)}\right)\right\}, \quad \text{for } k=1,\dots,R.

Finally, we obtained the distribution of the obtained coverages $\{C_k\}_{k=1}^R,$ which we plot for each $n$.

Note that to prevent the generation of $R$ distinct datasets of size $n+n_\text{val}$, we used the approach proposed in [49]. This method creates just one dataset of size $n+n_\text{val}$ for all $R$ rounds. Then, in each round, we shuffled this dataset and used the initial $n$ samples for calibration and the remaining $n_\text{val}$ samples for validation.

Figure <ref> illustrates the obtained distribution of coverages for $n \in \{500, 1000, 5000, 10000\}$. Observe that for $n>500$, the distribution, as expected, centers around 95% ($\alpha=0.05$). Notably, for $n=1000$, we achieve a satisfactory distribution, and the gain from increasing the calibration dataset size is minimal.

Finally, for $n=500$, the distribution peaks at a value larger than 95%. This could be due to the small number of calibration samples, which may not allow for a proper characterization of the coverage. However, even in this extreme case, conformalized-DeepONet still produces satisfactory confidence intervals around $1-\alpha$.

[width=*0.7, xlabel=Average coverage $C_k$, ylabel=Counts, ymin=0, ymax=70,legend pos=north west]
[dashed,thick, color=black] coordinates (0.95,0) (0.95,70);
[solid, color=blue,line width=1pt]
[solid, color=orange,line width=1pt]
[solid, color=red,line width=1pt]

[solid, color=mycolor1,line width=1pt]

The empirical distribution of coverages $C_k$ for different values of the calibration dataset size $n \in \{500, 1000, 5000, 10000\}$. We observe that for $n>500$, the distribution, as expected, centers around $95$% ($\alpha=0.05$). Notably, for $n=1000$ and this nonlinear pendulum dataset, we achieve a satisfactory distribution, and the gain from increasing $n$ is minimal.

§.§ Experiment: The Diffusion-Reaction System

In the second experiment, we examine the following diffusion-reaction system with a source term, $u(x)$, as described below:

\begin{align*}
& \frac{\partial s}{\partial t} = D \frac{\partial^2 s}{\partial x^2} + k s^2 + u(x), \qquad x \in [0,1],~t\in[0,1],
\end{align*}

with zero boundary/initial conditions. In the above, $D=0.01$ is the diffusion coefficient, and $k=0.01$ is the reaction rate.

Our goal is to approximate the solution operator $G:u(x) \mapsto s(x,1.0)$, that is, the mapping from the source term to the solution at the terminal time. As in the nonlinear pendulum experiment, we sample the input function (the source term $u(x)$) from the non-zero Gaussian Random Field $u \sim \mathcal{G}(0, k_\ell(x_1, x_2))$ with a radial basis function kernel and a length-scale of $\ell=0.1$. Finally, we discretize the sampled source terms using $m=100$ sensors.

Checking Coverage. We assessed the ability of all conformalized and baseline DeepONet models to generate appropriate confidence intervals for a miscoverage rate of $\alpha=0.05$. Figure <ref> displays the confidence intervals predicted by the proposed conformalized-DeepONet and baseline models for a randomly selected trajectory from the test dataset. Furthermore, Table <ref> presents the average coverage of the confidence intervals predicted by both the conformalized and baseline DeepONet models for all $n_\text{traj}=100$ test trajectories of the diffusion-reaction system.

Figure <ref> shows that the confidence intervals predicted by the conformalized-DeepONet appear to be adequate. Table <ref> confirms this, indicating that the average coverage for the conformalized-DeepONets is as expected ($\approx 95$%), considering a miscoverage rate of $\alpha=0.05$. The baseline Quantile-DeepONet also has a good coverage level compared to the Prob-DeepONet. This superior coverage might enable the conformalized Quantile-DeepONet to perform slightly better than the conformalized Prob-DeepONet, which must compensate for its baseline's poor coverage.

[ xlabel=$x$, ylabel=$G(u)(x) \equiv s(x,1.0)$,xmin=0, xmax=1,ymin=-0.7, ymax=3.1, legend pos=north west, line width=1pt,enlargelimits=false]

[solid, color=blue]
Reference Solution
[name path=prob_lower,  fill=none, draw=none,forget plot]
[name path=prob_upper, fill=none, draw=none,forget plot]

[name path=cop_lower, fill=none, draw=none,forget plot]
[name path=cop_upper, fill=none, draw=none,forget plot]

[fill=mycolor1, fill opacity=0.5] fill between[of=cop_upper and cop_lower, ];

[fill=red, fill opacity=0.5] fill between[of=prob_upper and prob_lower, ];

Conformal Prediction and Prob-DeepONet

[ xlabel=$x$, ylabel=$G(u)(x) \equiv s(x,1.0)$,,xmin=0, xmax=1,ymin=-0.7, ymax=3.1, legend pos=north west, line width=1pt,enlargelimits=false]
[solid, color=blue]
Reference Solution
[name path=quantile_lower,  fill=none, draw=none,forget plot]
[name path=quantile_upper, fill=none, draw=none,forget plot]

[name path=coq_lower, fill=none, draw=none,forget plot]
[name path=coq_upper, fill=none, draw=none,forget plot]

[fill=mycolor1, fill opacity=0.5] fill between[of=coq_upper and coq_lower, ];

[fill=red, fill opacity=0.5] fill between[of=quantile_upper and quantile_lower, ];

Conformal Prediction and Quantile-DeepONet

Confidence intervals for a random test trajectory of the diffusion-reaction experiment given a miscoverage rate $\alpha=0.05$. (a) Confidence intervals for the Conformalized Prob-DeepONet and the baseline Prob-DeepONet. (b) Confidence intervals for the Conformalized Quantile-DeepONet and the baseline Quantile-DeepONet. The average coverage rate for all samples for all approaches are presented in Table <ref>.

DeepONet model    Average coverage % ($\alpha=0.05$)

Conformalized Prob-DeepONet     96.26%

Conformalized Quantile-DeepONet     95.64%

Prob-DeepONet     8.99%

Quantile-DeepONet     83.09%

The average coverage percentage for the predicted confidence intervals (with a miscoverage rate of $\alpha=0.05$) calculated based on the 100 test trajectories in the diffusion-reaction experiment.

§.§ Experiment: The Viscous Burgers' Equation

In this section, we explore a common example found in most operator learning literature: the viscous Burgers' equation:

\begin{align*}
&\frac{\partial u_s}{\partial t} + \frac{1}{2}\frac{\partial (u^2_s)}{\partial x} = \alpha \frac{\partial^2 u_s}{\partial x^2},\hspace{0.5em} x\in[0, 2\pi], \hspace{0.2em} t\in[0, 0.3]\\
&u_s(x, 0) = u^0_s(x),\\
&u_s(0, t) = u_s(2\pi, t),
\end{align*}

where $u^0_s(x)$ is the initial condition that depends on the parameter $s$. In the above, the viscosity is set to $\alpha = 0.05$.

Our objective is to estimate the solution operator $G:u_{s}^0(x) \mapsto u_s(x,0.3)$, which is the mapping from the initial condition to the solution at the terminal time. We generated initial conditions using a sum of two Gaussian distributions with a uniformly randomized weight from $[0, 5]$. The means and standard deviations of the two distributions were sampled uniformly from $[0, 2\pi]$ and $[0.1, 1]$, respectively.

Checking Coverage. We evaluated the proposed conformalized-DeepONet models' capacity to produce rigorous confidence intervals. These intervals have a specified miscoverage rate of $\alpha=0.05$ for test trajectories generated from the viscous Burgers' model. We also compared these proposed models with the baseline models. Figure <ref> presents a comparison between the conformalized and baseline models. This comparison uses the predicted confidence interval for a randomly chosen test trajectory from the viscous Burgers' example. The figure indicates that both the conformalized Prob-DeepONet and Quantile-DeepONet effectively capture the reference solution within the predicted confidence interval. In constrast, the baseline models do not capture several portions of the reference solution.

[ xlabel=$x$, ylabel=$G(u)(x) \equiv u_s(x,0.3)$, ,xmin=0, xmax=6.3,ymin=-0.02, ymax=0.05,legend pos=south west, line width=1pt,enlargelimits=false]

[solid, color=blue]
Reference Solution
[name path=prob_lower,  fill=none, draw=none,forget plot]
[name path=prob_upper, fill=none, draw=none,forget plot]

[name path=cop_lower, fill=none, draw=none,forget plot]
[name path=cop_upper, fill=none, draw=none,forget plot]

[fill=mycolor1, fill opacity=0.5] fill between[of=cop_upper and cop_lower, ];

[fill=red, fill opacity=0.5] fill between[of=prob_upper and prob_lower, ];

Conformal Prediction and Prob-DeepONet

[ xlabel=$x$, ylabel=$G(u)(x) \equiv u_s(x,0.3)$, xmin=0, xmax=6.3,ymin=-0.02, ymax=0.05,legend pos=south west, line width=1pt,enlargelimits=false]

[solid, color=blue]
Reference Solution
[name path=quantile_lower,  fill=none, draw=none,forget plot]                 tablefigsData/Burgers/Burgers_lower_traj_quantile56.txt;
[name path=quantile_upper, fill=none, draw=none,forget plot]

[name path=coq_lower, fill=none, draw=none,forget plot]
[name path=coq_upper, fill=none, draw=none,forget plot]

[fill=mycolor1, fill opacity=0.5] fill between[of=coq_upper and coq_lower, ];

[fill=red, fill opacity=0.5] fill between[of=quantile_upper and quantile_lower, ];

Conformal Prediction and Quantile-DeepONet

Confidence intervals for a random test trajectory of the viscous Burgers' experiment given a miscoverage rate $\alpha=0.05$. (a) Confidence intervals for the Conformalized Prob-DeepONet and the baseline Prob-DeepONet. (b) Confidence intervals for the Conformalized Quantile-DeepONet and the baseline Quantile-DeepONet. The average coverage for all testing trajectories and approaches are summarized and presented in Table <ref>.

Moreover, Table <ref> shows the average coverage of confidence intervals (with a targeted coverage rate of $\alpha=0.05$) generated from both conformalized and baseline models for all $n_\text{traj}=100$ test trajectories produced using the viscous Burgers' model. The results once again demonstrate that both the conformalized Prob-DeepONet and the conformalized Quantile-DeepONet provide (on average) the desired coverage guarantee. In contrast, both baseline models generally offer poor coverage guarantees. It's noteworthy that the conformalized Prob-DeepONet still produces satisfactory coverage results, even though these are based on the inferior results from the baseline Prob-DeepONet model.

DeepONet model    Average coverage % ($\alpha=0.05$)

Conformalized Prob-DeepONet     94.88%

Conformalized Quantile-DeepONet     95.96%

Prob-DeepONet     7.84%

Quantile-DeepONet     72.91%

The average coverage percentage for the predicted confidence intervals (with a miscoverage rate of $\alpha=0.05$) calculated based on the 100 test trajectories in the viscous Burgers' experiment.

§ DISCUSSION

On Our Results. Our results indicate that split conformal prediction [51, 49, 57] is an easy-to-use methodology that offers confidence intervals with coverage guarantees. We consider it easy-to-use as it does not require any distributional assumptions and only needs a finite number of DeepONet samples. This allowed us to build an adaptive, reliable, and efficient Uncertainty Quantification (UQ) and regression framework for Deep Operator Networks.

Our results also indicate that split conformal prediction is model agnostic. Note that we  applied it on top of two UQ DeepONet frameworks (B-DeepONet and Prob-DeepONet), proposed in our previous works, and a new UQ framework (Quantile-DeepONet) proposed here. Hence, we expect split conformal prediction to be compatible with many existing DeepONet extensions, such as Multi-fidelity DeepONet [12, 68, 69], BelNet [25], Fed-DeepONet [70], or D2NO [22]. To demonstrate this, we will present a simple experiment that highlights the use of split conformal prediction with Multi-fidelity DeepONets.

Conformal Prediction and Multi-Fidelity DeepONet.
In this experiment, we used a multi-fidelity conformalized Prob-DeepONet model to approximate the high-fidelity solution operator for the given 1D jump function:

\begin{equation}
\begin{aligned}
y_L(u)(x) &=
\begin{cases}
0.5(6x - 2)^2 \sin(u) + 10(x - 0.5) - 5 &  x \leq 0.5 \\
0.5(6x - 2)^2 \sin(u) + 10(x - 0.5) - 2 &  x > 0.5
\end{cases} \\
y_H(u)(x) &= 2y_L(u)(x) - 20x + 20 \\
u(x) &= ax - 4
\end{aligned}
\end{equation}

where $x \in [0, 1]$ and a $\in [10, 14]$. In the above, $y_L$ is the low-fidelity solution, $y_H$ the high-fidelity solution, and $u$ the input.

We built our Multi-fidelity Prob-DeepONet baseline model as follows. Initially, we used a trained classical DeepONet to approximate the low-fidelity solution operator $G:u(x) \mapsto \hat{y}{LF}(u)(x)$. Next, we trained a Prob-DeepONet to approximate the difference between the high-fidelity and low-fidelity solution operators $(y_{HF} - y_{LF})~|~X=(u,x) \sim \mathcal{N}(\mu(u)(x),\sigma(u)(x))$. Finally, we recovered the high-fidelity solution operator using the probabilistic model $\mathcal{N}(\mu(u)(x) + \hat{y}(u)(x), \sigma(u)(x))$.

We trained the baseline model using a dataset of $N_{LF}=3800$ low-fidelity DeepONet triplets and $N_{HF}=760$ high-fidelity DeepONet triplets. For calibration, we used the half the of number of high-fidelity data samples ($n=\frac{N_{HF}}{2} = 380$).

Checking Coverage. We tested the proposed Multi-Fidelity Conformalized Prob-DeepONet to create confidence intervals with coverage guarantees for a target miscoverage rate of $\alpha=0.05$. Figure <ref> shows the confidence intervals for a randomly chosen high-fidelity test trajectory. The Multi-fidelity Conformalized model provides a suitable confidence interval, showing increased uncertainty at the high-fidelity solution's jumps and turns. Additionally, Table <ref> displays the average coverage of confidence intervals (targeting a miscoverage rate of $\alpha=0.05$), generated by both the Multi-fidelity conformalized and baseline Prob-DeepONet models for all $n_\text{traj}=100$ high-fidelity test trajectories. The results suggest satisfactory coverage, which can be enhanced by using more calibration examples.

[width=*0.7, xlabel=$x$, ylabel=$y_H(u)(x)$,ymin=-4.5, ymax=14.0,legend pos=north west]

[dashed, color=blue,mark= *, mark size=1.2pt]
High-fidelity Reference Solution
[name path=prob_lower,  fill=none, draw=none,forget plot]
[name path=prob_upper, fill=none, draw=none,forget plot]

[name path=cop_lower, fill=none, draw=none,forget plot]
[name path=cop_upper, fill=none, draw=none,forget plot]

[fill=mycolor1, fill opacity=0.5] fill between[of=cop_upper and cop_lower, ];

[fill=red, fill opacity=0.5] fill between[of=prob_upper and prob_lower, ];

Confidence intervals predicted by the Multi-fidelity conformalized Prob-DeepONet and the corresponding baseline model for a random high-fidelity test trajectory of the 1D jump function, given a miscoverage rate of $\alpha=0.05$. We summarize the average coverage rate for all approaches in Table <ref>.

DeepONet model    Average coverage % ($\alpha=0.05$)

Multi-fidelity Conformalized Prob-DeepONet     97.02%

Multi-fidelity Prob-DeepONet     82.07%

The average coverage percentage for the predicted confidence intervals (with a given miscoverage rate of $\alpha=0.05$) calculated based on the $100$ high-fidelity test trajectories of the 1D jump function experiment.

On Our Future Work. The findings discussed in this paper are undoubtedly promising and applicable in scenarios where a reliable surrogate model is required. Thus, in our future work, we plan to investigate the application of conformalized-DeepONets in real-world settings. This includes power systems [27], rib-optimization in fluid systems [36], PDE-based topology optimization [12], and forecasting extreme events. Additionally, we plan to investigate the use of conformal prediction in DeepONet settings where the properties of exchangeability and independent and identically distributed (i.i.d.) do not apply. For instance, in DeepONet extrapolation [20], DeepONet for non-autonomous systems [19], or Federated DeepONet [70].

§ CONCLUSION

In the paper, we presented methods to apply split conformal prediction for providing confidence intervals for Deep Operator Network (DeepONet) prediction. These intervals come with coverage guarantees and do not require distributional assumptions. Specifically, we used split conformal prediction to enhance our previously proposed Probabilistic and Bayesian DeepONets, enabling them to predict rigorous confidence intervals at a predetermined miscoverage rate. Moreover, we designed a novel extension of DeepONet, known as Quantile-DeepONet, which estimated conditional quantiles for DeepONet predictions and provided a more natural setting for applying conformal prediction. By combining split conformal prediction with Quantile-DeepONet, we successfully developed an effective and distribution-free methodology for constructing confidence intervals with guaranteed coverage.

§ ACKOWNLEDGEMENT

G. Lin acknowledges the support of the National Science Foundation (DMS-2053746, DMS-2134209, ECCS-2328241, and OAC-2311848), and U.S. Department of Energy (DOE) Office of Science Advanced Scientific Computing Research program DE-SC0023161, the Uncertainty Quantification for Multifidelity Operator Learning (MOLUcQ) project (Project No. 81739) and DOE–Fusion Energy Science, under grant number: DE-SC0024583.
L. Lu was supported by the U.S. Department of Energy [DE-SC0022953].

[1]
Andrew R Barron.
Universal approximation bounds for superpositions of a sigmoidal function.
IEEE Transactions on Information theory, 39(3):930–945, 1993.

[2]
Ali Rahimi and Benjamin Recht.
Weighted sums of random kitchen sinks: Replacing minimization with randomization in learning.
Advances in neural information processing systems, 21, 2008.

[3]
Maziar Raissi, Paris Perdikaris, and George E Karniadakis.
Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.
Journal of Computational physics, 378:686–707, 2019.

[4]
Lu Lu, Xuhui Meng, Zhiping Mao, and George Em Karniadakis.
Deepxde: A deep learning library for solving differential equations.
SIAM review, 63(1):208–228, 2021.

[5]
George Em Karniadakis, Ioannis G Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang.
Physics-informed machine learning.
Nature Reviews Physics, 3(6):422–440, 2021.

[6]
Wing Tat Leung, Guang Lin, and Zecheng Zhang.
Nh-pinn: Neural homogenization-based physics-informed neural network for multiscale problems.
Journal of Computational Physics, page 111539, 2022.

[7]
Pratik Rathore, Weimu Lei, Zachary Frangella, Lu Lu, and Madeleine Udell.
Challenges in training pinns: A loss landscape perspective.
arXiv preprint arXiv:2402.01868, 2024.

[8]
Tianping Chen and Hong Chen.
Universal approximation to nonlinear operators by neural networks with arbitrary activation functions and its application to dynamical systems.
IEEE Transactions on Neural Networks, 6(4):911–917, 1995.

[9]
Tianping Chen and Hong Chen.
Approximations of continuous functionals by neural networks with application to dynamic systems.
IEEE Transactions on Neural networks, 4(6):910–918, 1993.

[10]
Lu Lu, Pengzhan Jin, Guofei Pang, Zhongqiang Zhang, and George Em Karniadakis.
Learning nonlinear operators via deeponet based on the universal approximation theorem of operators.
Nature Machine Intelligence, 3(3):218–229, 2021.

[11]
Lu Lu, Xuhui Meng, Shengze Cai, Zhiping Mao, Somdatta Goswami, Zhongqiang Zhang, and George Em Karniadakis.
A comprehensive and fair comparison of two neural operators (with practical extensions) based on fair data.
Computer Methods in Applied Mechanics and Engineering, 393:114778, 2022.

[12]
Lu Lu, Raphaël Pestourie, Steven G Johnson, and Giuseppe Romano.
Multifidelity deep neural operators for efficient learning of partial differential equations with application to fast inverse design of nanoscale heat transport.
Physical Review Research, 4(2):023210, 2022.

[13]
Min Zhu, Shihang Feng, Youzuo Lin, and Lu Lu.
Fourier-deeponet: Fourier-enhanced deep operator networks for full waveform inversion with improved accuracy, generalizability, and robustness.
arXiv preprint arXiv:2305.17289, 2023.

[14]
Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar.
Fourier neural operator for parametric partial differential equations.
arXiv preprint arXiv:2010.08895, 2020.

[15]
Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar.
Neural operator: Graph kernel network for partial differential equations.
arXiv preprint arXiv:2003.03485, 2020.

[16]
Gege Wen, Zongyi Li, Kamyar Azizzadenesheli, Anima Anandkumar, and Sally M Benson.
U-fno—an enhanced fourier neural operator-based deep-learning model for multiphase flow.
Advances in Water Resources, 163:104180, 2022.

[17]
Zongyi Li, Daniel Zhengyu Huang, Burigede Liu, and Anima Anandkumar.
Fourier neural operator with learned deformations for pdes on general geometries.
arXiv preprint arXiv:2207.05209, 2022.

[18]
Pengzhan Jin, Shuai Meng, and Lu Lu.
Mionet: Learning multiple-input operators via tensor product.
SIAM Journal on Scientific Computing, 44(6):A3490–A3514, 2022.

[19]
Guang Lin, Christian Moya, and Zecheng Zhang.
Learning the dynamical response of nonlinear non-autonomous dynamical systems with deep operator neural networks.
Engineering Applications of Artificial Intelligence, 125:106689, 2023.

[20]
Min Zhu, Handi Zhang, Anran Jiao, George Em Karniadakis, and Lu Lu.
Reliable extrapolation of deep neural operators informed by physics or sparse observations.
Computer Methods in Applied Mechanics and Engineering, 412:116064, 2023.

[21]
Guang Lin, Christian Moya, and Zecheng Zhang.
B-deeponet: An enhanced bayesian deeponet for solving noisy parametric pdes using accelerated replica exchange sgld.
Journal of Computational Physics, 473:111713, 2023.

[22]
Zecheng Zhang, Christian Moya, Lu Lu, Guang Lin, and Hayden Schaeffer.
D2no: Efficient handling of heterogeneous input function spaces with distributed deep neural operators.
arXiv preprint arXiv:2310.18888, 2023.

[23]
Jaideep Pathak, Shashank Subramanian, Peter Harrington, Sanjeev Raja, Ashesh Chattopadhyay, Morteza Mardani, Thorsten Kurth, David Hall, Zongyi Li, Kamyar Azizzadenesheli, et al.
Fourcastnet: A global data-driven high-resolution weather model using adaptive fourier neural operators.
arXiv preprint arXiv:2202.11214, 2022.

[24]
Chensen Lin, Zhen Li, Lu Lu, Shengze Cai, Martin Maxey, and George Em Karniadakis.
Operator learning for predicting multiscale bubble growth dynamics.
The Journal of Chemical Physics, 154(10), 2021.

[25]
Zecheng Zhang, Wing Tat Leung, and Hayden Schaeffer.
Belnet: Basis enhanced learning, a mesh-free neural operator.
arXiv preprint arXiv:2212.07336, 2022.

[26]
Zecheng Zhang, Wing Tat Leung, and Hayden Schaeffer.
A discretization-invariant extension and analysis of some deep operator networks.
arXiv preprint arXiv:2307.09738, 2023.

[27]
Christian Moya, Shiqi Zhang, Guang Lin, and Meng Yue.
Deeponet-grid-uq: A trustworthy deep operator framework for predicting the power grid’s post-fault trajectories.
Neurocomputing, 535:166–182, 2023.

[28]
Joel Hayford, Jacob Goldman-Wetzler, Eric Wang, and Lu Lu.
Speeding up and reducing memory usage for scientific machine learning via mixed precision.
arXiv preprint arXiv:2401.16645, 2024.

[29]
Yixuan Sun, Christian Moya, Guang Lin, and Meng Yue.
Deepgraphonet: A deep graph operator network to learn and zero-shot transfer the dynamic response of networked systems.
IEEE Systems Journal, 2023.

[30]
Minglang Yin, Nicolas Charon, Ryan Brody, Lu Lu, Natalia Trayanova, and Mauro Maggioni.
Dimon: Learning solution operators of partial differential equations on a diffeomorphic family of domains.
arXiv preprint arXiv:2402.07250, 2024.

[31]
Zhongyi Jiang, Min Zhu, Dongzhuo Li, Qiuzi Li, Yanhua O Yuan, and Lu Lu.
Fourier-mionet: Fourier-enhanced multiple-input neural operators for multiphase modeling of geological carbon sequestration.
arXiv preprint arXiv:2303.04778, 2023.

[32]
Shunyuan Mao, Ruobing Dong, Lu Lu, Kwang Moo Yi, Sifan Wang, and Paris Perdikaris.
Ppdonet: Deep operator networks for fast prediction of steady-state solutions in disk–planet systems.
The Astrophysical Journal Letters, 950(2):L12, 2023.

[33]
Zhiping Mao, Lu Lu, Olaf Marxen, Tamer A Zaki, and George Em Karniadakis.
Deepm&mnet for hypersonics: Predicting the coupled flow and finite-rate chemistry behind a normal shock using neural-network approximation of operators.
Journal of computational physics, 447:110698, 2021.

[34]
Patricio Clark Di Leoni, Lu Lu, Charles Meneveau, George Em Karniadakis, and Tamer A Zaki.
Neural operator prediction of linear instability waves in high-speed boundary layers.
Journal of Computational Physics, 474:111793, 2023.

[35]
Christian Moya, Guang Lin, Tianqiao Zhao, and Meng Yue.
On approximating the dynamic response of synchronous generators via operator learning: A step towards building deep operator-based power grid simulators.
arXiv preprint arXiv:2301.12538, 2023.

[36]
Izzet Sahin, Christian Moya, Amirhossein Mollaali, Guang Lin, and Guillermo Paniagua.
Deep operator learning-based surrogate models with uncertainty quantification for optimizing internal cooling channel rib profiles.
International Journal of Heat and Mass Transfer, 219:124813, 2024.

[37]
Apostolos F Psaros, Xuhui Meng, Zongren Zou, Ling Guo, and George Em Karniadakis.
Uncertainty quantification in scientific machine learning: Methods, metrics, and comparisons.
Journal of Computational Physics, 477:111902, 2023.

[38]
Zongren Zou, Xuhui Meng, Apostolos F Psaros, and George E Karniadakis.
Neuraluq: A comprehensive library for uncertainty quantification in neural differential equations and operators.
SIAM Review, 66(1):161–190, 2024.

[39]
Yibo Yang, Georgios Kissas, and Paris Perdikaris.
Scalable uncertainty quantification for deep operator networks using randomized priors.
Computer Methods in Applied Mechanics and Engineering, 399:115399, 2022.

[40]
Shailesh Garg and Souvik Chakraborty.
Vb-deeponet: A bayesian operator learning framework for uncertainty quantification.
Engineering Applications of Artificial Intelligence, 118:105685, 2023.

[41]
Ling Guo, Hao Wu, Wenwen Zhou, Yan Wang, and Tao Zhou.
Ib-uq: Information bottleneck based uncertainty quantification for neural function regression and neural operator learning.
arXiv preprint arXiv:2302.03271, 2023.

[42]
Zecheng Zhang, Christian Moya, Wing Tat Leung, Guang Lin, and Hayden Schaeffer.
Bayesian deep operator learning for homogenized1 to fine-scale maps for multiscale pde.

[43]
Christian Moya Calderon and Guang Lin.
Bayesian, multifidelity operator learning for complex engineering systems-a position paper.
Journal of Computing and Information Science in Engineering, pages 1–9, 2023.

[44]
Max Welling and Yee W Teh.
Bayesian learning via stochastic gradient langevin dynamics.
In Proceedings of the 28th international conference on machine learning (ICML-11), pages 681–688. Citeseer, 2011.

[45]
Arnak Dalalyan.
Further and stronger analogy between sampling and optimization: Langevin monte carlo and gradient descent.
In Conference on Learning Theory, pages 678–689. PMLR, 2017.

[46]
Arnak S Dalalyan and Avetik Karagulyan.
User-friendly guarantees for the langevin monte carlo with inaccurate gradient.
Stochastic Processes and their Applications, 129(12):5278–5311, 2019.

[47]
Maxim Raginsky, Alexander Rakhlin, and Matus Telgarsky.
Non-convex learning via stochastic gradient langevin dynamics: a nonasymptotic analysis.
In Conference on Learning Theory, pages 1674–1703. PMLR, 2017.

[48]
Haoyang Zheng, Wei Deng, Christian Moya, and Guang Lin.
Accelerating approximate thompson sampling with underdamped langevin monte carlo.
arXiv preprint arXiv:2401.11665, 2024.

[49]
Anastasios N Angelopoulos and Stephen Bates.
A gentle introduction to conformal prediction and distribution-free uncertainty quantification.
arXiv preprint arXiv:2107.07511, 2021.

[50]
Volodya Vovk, Alexander Gammerman, and Craig Saunders.
Machine-learning applications of algorithmic randomness.

[51]
Vladimir Vovk, Alexander Gammerman, and Glenn Shafer.
Algorithmic learning in a random world, volume 29.
Springer, 2005.

[52]
Glenn Shafer and Vladimir Vovk.
A tutorial on conformal prediction.
Journal of Machine Learning Research, 9(3), 2008.

[53]
Harris Papadopoulos, Kostas Proedrou, Volodya Vovk, and Alex Gammerman.
Inductive confidence machines for regression.
In Machine Learning: ECML 2002: 13th European Conference on Machine Learning Helsinki, Finland, August 19–23, 2002 Proceedings 13, pages 345–356. Springer, 2002.

[54]
Victor Quach, Adam Fisch, Tal Schuster, Adam Yala, Jae Ho Sohn, Tommi S Jaakkola, and Regina Barzilay.
Conformal language modeling.
arXiv preprint arXiv:2306.10193, 2023.

[55]
Kexin Huang, Ying Jin, Emmanuel Candes, and Jure Leskovec.
Uncertainty quantification over graph with conformalized graph neural networks.
Advances in Neural Information Processing Systems, 36, 2024.

[56]
Anastasios Angelopoulos, Emmanuel Candes, and Ryan J Tibshirani.
Conformal pid control for time series prediction.
Advances in Neural Information Processing Systems, 36, 2024.

[57]
Yaniv Romano, Evan Patterson, and Emmanuel Candes.
Conformalized quantile regression.
Advances in neural information processing systems, 32, 2019.

[58]
Beichuan Deng, Yeonjong Shin, Lu Lu, Zhongqiang Zhang, and George Em Karniadakis.
Approximation rates of deeponets for learning operators arising from advection–diffusion equations.
Neural Networks, 153:411–426, 2022.

[59]
Guanxun Li, Guang Lin, Zecheng Zhang, and Quan Zhou.
Fast replica exchange stochastic gradient langevin dynamics.
arXiv preprint arXiv:2301.01898, 2023.

[60]
Tianqi Chen, Emily Fox, and Carlos Guestrin.
Stochastic gradient hamiltonian monte carlo.
In International conference on machine learning, pages 1683–1691. PMLR, 2014.

[61]
RN Bhattacharya.
Criteria for recurrence and existence of invariant measures for multidimensional diffusions.
The Annals of Probability, pages 541–553, 1978.

[62]
Gareth O Roberts and Richard L Tweedie.
Exponential convergence of langevin distributions and their discrete approximations.
Bernoulli, pages 341–363, 1996.

[63]
Roger Koenker and Gilbert Bassett Jr.
Regression quantiles.
Econometrica: journal of the Econometric Society, pages 33–50, 1978.

[64]
Ingo Steinwart and Andreas Christmann.
Estimating conditional quantiles with the help of the pinball loss.
Bernoulli, 17(1):211–225, 2011.

[65]
Diederik P Kingma and Jimmy Ba.
Adam: A method for stochastic optimization.
arXiv preprint arXiv:1412.6980, 2014.

[66]
Guang Lin, Yating Wang, and Zecheng Zhang.
Multi-variance replica exchange sgmcmc for inverse and forward problems via bayesian pinn.
Journal of Computational Physics, 460:111173, 2022.

[67]
Wei Deng, Qi Feng, Liyao Gao, Faming Liang, and Guang Lin.
Non-convex learning via replica exchange stochastic gradient mcmc.
In International Conference on Machine Learning, pages 2474–2483. PMLR, 2020.

[68]
Amanda A Howard, Mauro Perego, George E Karniadakis, and Panos Stinis.
Multifidelity deep operator networks.
arXiv preprint arXiv:2204.09157, 2022.

[69]
Amirhossein Mollaali, Izzet Sahin, Iqrar Raza, Christian Moya, Guillermo Paniagua, and Guang Lin.
A physics-guided bi-fidelity fourier-featured operator learning framework for predicting time evolution of drag and lift coefficients.
Fluids, 8(12):323, 2023.

[70]
Christian Moya and Guang Lin.
Fed-deeponet: Stochastic gradient-based federated training of deep operator networks.
Algorithms, 15(9):325, 2022.