1 Introduction↩︎

Approximating the solutions of partial differential equations (PDEs) with numerical algorithms is a fundamental problem in scientific and engineering applications. Traditional methods include finite difference, finite element, and spectral method which require access to the governing equation. In recent years, the use of machine learning for solving PDEs has attracted attention, with several directions for addressing the approximation problems. Physics-informed neural network (PINN) [1] trains a network to fit a single PDE task using a least squares fit on the data and on the PDE using a collocation method. However, any changes to the original problem leads to retraining the neural network, which is computationally expensive. Operator Learning trains a network to approximate the solution operator between the input and output function spaces directly. For example, we aim to approximate an operator between boundary condition and the PDE solution itself. The approximation of operators using neural network was first introduced in [2]. With the development of deep neural networks, recent results include various neural operators such as Deep Operator Nets [3], Fourier Neural Operators [4] and several others [5], [6]. Neural operators have been used in several scientific applications, for example, solving spatio-temporal dynamics [3], [7]–[9], dynamical system with control [2], reduced order modeling (ROM) [10], climate predictions [11] and uncertainty quantification [12]. While these neural networks based operator learning methods are often used in practice, the theoretical analysis relies on universal approximation property of deep neural networks, which shows the existence of a network of a requisite size achieving a certain error rate. However, the existence result does not guarantee that the network is obtained in practice, [2], [13], [14]. Some prior results go beyond universal approximation by providing theoretical guarantees for training algorithms and sample complexity, see [15]–[18] and references therein.

To address the limited theoretical results in operator learning, a kernel or Gaussian Process (GP) based framework was proposed in [19], [20] along with a priori error estimates and convergence guarantees. In several benchmark tests, it was shown that the kernel-based approach matches or outperforms neural operator methods in terms of test accuracy and computational costs. Kernel methods operate on the kernel matrix (Gram matrix) of the data with size \(m\times m\), where \(m\) represents the number of samples in the dataset. This leads to poor scaling with the size of training dataset, i.e., large training sets incur significant computational and storage costs.

The class of random features is one of the most popular techniques to speed up kernel methods in large-scale problems. Rather than forming and computing solutions using the kernel matrix directly, the random feature method (RFM) [21], [22] maps data into a relatively low-dimensional randomized feature space, which significantly reduces the computation needed for training. Additionally, the random feature model can be viewed as a randomized two-layer neural network whose weights connecting input layer and single hidden layer are randomly generated from a known distribution rather than trainable parameters. Several approaches utilized RFMs for approximating solutions to PDEs [23]–[25]. In [23], the authors utilized the random features as a randomized Galerkin method for solving the PDE. In [24] and a subsequent review paper [25], the authors took the operator learning perspective to approximate the solutions of the PDEs. Their analysis relies on the theory of vector-valued reproducing kernels. As a consequence of [26], a high probability non-asymptotic error bound was derived. However, this quantitative result is dependent on the assumption that the target operator belongs to a reproducing kernel Hilbert space corresponding to an operator-valued kernel, which may not hold in practice. Moreover, the random feature maps were carefully designed and adapted to different problems. However, training also requires samples in the frequency domain, which may not be available for the original problem.

In this paper, we propose a novel random features based operator learning method. We use the random Fourier features [21] and randomly generate features from a known probability distribution to approximate the solution operator. We provide the theoretical guarantees and experimental validation of our proposed method. Our contributions are summarized below.

• Random Feature Method for Operator Learning. We follow the existing operator learning framework in [18], [19], [27]. Rather than using neural networks or kernel methods, we introduce the random features method for the operator learning problem addressing one of the open problems in [19].

• Error Analysis. We derive an error bound for proposed approach in Section 4. Our error analysis relies on the generalization error of random feature model for function approximation (Section 3), which has been of recent interest as well, see [28]–[38]. We consider the overparametrized setting (the number of random features \(N\) is larger than the number of training samples \(m\)) and the min-norm interpolation problem (or the ridgeless limit as the regularization parameter \(\lambda\to0^+\)). Our analysis depends on the condition number of random feature matrix, which is similar as the theoretical results in [31]. To better match the expected smoothness of the solutions, we propose Cauchy random features which lead to mixed Sobolev spaces.

• Improved Numerical Performance. In Section 5, we compare our method with kernel based method and neural operator methods empirically. Numerical experiments on benchmark PDE problems show several advantages of our method: 1) our method is easy to implement, 2) expensive computational resources are not required, 3) the training involves solving several convex optimization problems (in parallel) which significantly reduces the training time, and 4) it is competitive in terms of test-accuracy.

The paper is organized as follows. In Section 2, we introduce the proposed random feature operator learning method. Then, we prove the generalization error bounds of the random feature method for function learning in Section 3. Using the obtained error bounds in section 3, we derive generalization error bounds for random feature operator learning in Section 4. Lastly, Section 5 presents the numerical results and comparisons.

Notations. We let \(\mathbb{R}\) be the set of all real number and \(\mathbb{C}\) be the set of all complex number where \(i=\sqrt{-1}\) denotes the imaginary unit. We define the set \([N]\) to be all natural numbers smaller than of equal to \(N\), i.e., \(\{1,2,\dots,N\}\). We use bold letters to denote vectors or matrices, and denote the identity matrix of size \(n\times n\) by \(\mathbf{I}_n\). For any two vectors \(\mathbf{x},\mathbf{y}\in\mathbb{C}^d\), the inner product is denoted by \(\langle \mathbf{x},\mathbf{y}\rangle = \sum_{j=1}^d x_j\Bar{y}_j\), where \(\mathbf{x}= [x_1,\dots,x_d]^\top\) and \(\mathbf{y}= [y_1,\dots,y_d]^\top\). For a vector \(\mathbf{x}\in\mathbb{C}^d\), we denote by \(\|\mathbf{x}\|_p\) the \(\ell_p\)-norm of \(\mathbf{x}\) and for a matrix \(\mathbf{A}\in\mathbb{C}^{m\times N}\) the (induced) \(p\)-norm is written as \(\|\mathbf{A}\|_p\). The conjugate transpose of a matrix \(\mathbf{A}\in\mathbb{C}^{m\times N}\) is denoted by \(\mathbf{A}^*\). We use \(\mathcal{U}\) to denote the input function set with domain \(D_\mathcal{U}\), and \(\mathcal{V}\) to denote the output function set with domain \(D_\mathcal{V}\). The operator mapping functions in \(\mathcal{U}\) to functions in \(\mathcal{V}\) is denoted by \(G\). The \(L^2\) norm of a function over domain \(D_\mathcal{U}\) is defined as \(\|u\|_{L^2(D_\mathcal{U})} = \sqrt{\int_{D_\mathcal{U}} |u(\mathbf{x})|^2 d\mathbf{x}}\). We write \(B_\delta(\mathcal{U})\) for the ball of radius \(\delta>0\) in a normed space \(\mathcal{U}\) equipped with the \(L^2\) norm.

2 Problem Statement↩︎

The goal is to learn an unknown Lipschitz operator \(G:\mathcal{U}\to\mathcal{V}\) between two function sets \(\mathcal{U}\) and \(\mathcal{V}\) from training samples \(\{(u_j, v_j)\}_{j\in[M]} \subset \mathcal{U}\times\mathcal{V}\) such that \(G(u_j)=v_j\) for all \(j\in[M]\). We consider Lipschitz operators in the following sense.

For our problem, we consider the setting where the input/output functions are only partially observed through a finite collection of function values at given collocation points. Denote collocation points in \(D_\mathcal{U}\) and \(D_\mathcal{V}\) by \(\{\mathbf{x}_j\}_{j\in[n]} \subset D_\mathcal{U}\) and by \(\{\mathbf{y}_j\}_{j\in[m]} \subset D_\mathcal{V}\), respectively. Then, we define sampling operators \(S_\mathcal{U}:\mathcal{U}\to\mathbb{R}^n\) and \(S_\mathcal{V}:\mathcal{V}\to\mathbb{R}^m\) as \[S_\mathcal{U}(u) = [u(\mathbf{x}_1),\dots,u(\mathbf{x}_n)]^\top\in\mathbb{R}^n,andS_\mathcal{V}(v) = [v(\mathbf{y}_1),\dots,v(\mathbf{y}_m)]^\top\in\mathbb{R}^m.\] Indeed sampling operators are linear, and we further assume that \(S_\mathcal{U}\) and \(S_\mathcal{V}\) are bounded operators throughout this paper. We could also generalize to the situation where functions are accessible through arbitrary linear measurements. For the sake of simplicity, we will focus on point evaluations in this paper. We formalize the problem statement as follows:

2.1 Operator learning framework↩︎

Our problem statement gives rise to a diagram for operator learning, which is depicted in Figure 1. Here, the map \(f:\mathbb{R}^n\to\mathbb{R}^m\) is defined explicitly as \[\label{operator95f} f = S_\mathcal{V}\circ G\circ R_\mathcal{U},\tag{1}\] where \(R_\mathcal{U}:\mathbb{R}^n\to\mathcal{U}\) is a recover map. Notice that the mapping \(f:\mathbb{R}^n\to\mathbb{R}^m\) is between finite-dimensional spaces, which is more amenable to numerical approximation. We make some assumptions on \(f\) throughout this paper.

After constructing recovery maps \(R_\mathcal{V}:\mathbb{R}^m\to\mathcal{V}\), and \(\hat{f}:\mathbb{R}^n\to\mathbb{R}^m\), we are able to define operator \(\Hat{G}\) to approximate the target operator \(G\) as \[\Hat{G} = R_\mathcal{V}\circ \Hat{f} \circ S_\mathcal{U}.\] In [19], the authors proposed the same framework to do operator learning where they used the optimal recovery maps in reproducing kernel Hilbert spaces (RKHS) as recovery maps \(R_\mathcal{V}\), \(R_\mathcal{U}\) and \(\hat{f}\). Precisely, the map \(\hat{f}:\mathbb{R}^n\to\mathbb{R}^m\), which is viewed as an approximation of \(f\), is the optimal recovery map in a vector-valued RKHS. The optimal recovery map \(R_\mathcal{V}:\mathbb{R}^m\to\mathcal{V}\) takes function values as inputs and returns a function in \(\mathcal{V}\). As the theory of optimal recovery suggests, the optimal recovery map returns the kernel interpolation [39]. Other function approximation/reconstruction techniques can be applied as well. For example, in [40], the authors used neural networks to approximate \(f:\mathbb{R}^n\to\mathbb{R}^m\).

2.2 Random Feature Operator Learning↩︎

The key idea of the random feature method is to draw \(N\) random frequencies \(\boldsymbol{\omega}_k\in\mathbb{R}^d\) from a density \(\rho(\boldsymbol{\omega})\), and then to construct an approximation as \[\label{RF95model} f(\cdot,\mathbf{c}^\sharp) = \sum_{k=1}^N c_k^\sharp \exp \left(i\langle \boldsymbol{\omega}_k,\cdot\rangle\right).\tag{2}\] The random feature model can be viewed as a one layer neural network where the weights in the first layer are drawn from a known distribution and only the weights in the second layer are trainable parameters. Instead of using the exponential function as an activation function, we could also use rectified linear unit (ReLU) activation function [41], or trigonometric functions [30]. Originally, the random feature method was proposed to speed up large-scale kernel machines, see [21], [42] for details.

We define the recovery map \(R_\mathcal{U}:\mathbf{u}\mapsto R_\mathcal{U}(\mathbf{u}):=f(\cdot, \mathbf{c}^\sharp(\mathbf{u}))\), where \[\label{opt95Ru} \mathbf{c}^\sharp(\mathbf{u}) \in \mathop{\mathrm{argmin}}_{\mathbf{c}\in\mathbb{R}^N} \|\mathbf{c}\|_2 \qquad subject toS_\mathcal{U}(f(\cdot,\mathbf{c})) = \mathbf{u}.\tag{3}\] The coefficient vector \(\mathbf{c}^\sharp\in\mathbb{R}^N\) is indeed trained by solving the min-norm interpolation problem. Recall that \(S_\mathcal{U}:\mathcal{U}\to\mathbb{R}^n\) is the point evaluation operator, then the constraint reads as \(R_\mathcal{U}(\mathbf{u})(\mathbf{x}_i) = u_i\) for \(i\in[n]\). Therefore, we can rewrite the min-norm interpolation problem and the optimal coefficient vector \(\mathbf{c}^\sharp\in\mathbb{R}^N\) is trained by solving \[\mathbf{c}^\sharp(\mathbf{u}) \in \mathop{\mathrm{argmin}}_{\mathbf{A}_x \mathbf{c}=\mathbf{u}} \|\mathbf{c}\|_2,\] where matrix \(\mathbf{A}_x\) is defined (component-wise) by \((\mathbf{A}_x)_{j,k} = \exp(i\langle \boldsymbol{\omega}_k, \mathbf{x}_j\rangle)\). The solution is given by \(\mathbf{c}^\sharp = \mathbf{A}_x^*(\mathbf{A}_x\mathbf{A}_x^*)^{-1}\mathbf{u}= \mathbf{A}_x^\dagger\mathbf{u}\), where \(\mathbf{A}_x^\dagger\) is the pseudo-inverse of matrix \(\mathbf{A}_x\). We will prove that the matrix \(\mathbf{A}_x\mathbf{A}_x^*\) is invertible with high probability in Section 3, and hence pseudo-inverse \(\mathbf{A}_x^\dagger\) is well-defined. Following the same arguments, we define the recovery map \(R_\mathcal{V}:\mathbb{R}^m\to\mathcal{V}\). The optimal coefficient vector \(\mathbf{c}^\sharp\in\mathbb{R}^N\) is trained by solving the following min-norm interpolation problem ¹ \[\mathbf{c}^\sharp \in \mathop{\mathrm{argmin}}_{\mathbf{A}_y \mathbf{c}=\mathbf{v}} \|\mathbf{c}\|_2,\] where the matrix \(\mathbf{A}_y\) is defined (component-wise) by \((\mathbf{A}_y)_{j,k} = \exp(i\langle \boldsymbol{\omega}_k,\mathbf{y}_j\rangle)\). The optimal coefficient vector is computed by \(\mathbf{c}^\sharp = \mathbf{A}_y^*(\mathbf{A}_y\mathbf{A}_y^*)^{-1}\mathbf{v}= \mathbf{A}_y^\dagger \mathbf{v}\), where \(\mathbf{A}_y^\dagger\) is the pseudo-inverse of matrix \(\mathbf{A}_y\).

Suppose that Assumption 2 holds, we propose to approximate each \(f_j\) separately using random feature model, and hence the vector-valued random feature map \(\Hat{f}:\mathbb{R}^n\to\mathbb{R}^m\) is defined as \[\label{operator95rf} \Hat{f}(\mathbf{u}) = \left[ \Hat{f}_1(\mathbf{u}), \dots, \Hat{f}_m(\mathbf{u})\right]^\top\in\mathbb{R}^m.\tag{4}\] Each \(\Hat{f}_j:\mathbb{R}^n\to\mathbb{R}\) is a random feature approximation of \(f_j:\mathbb{R}^n\to\mathbb{R}\) and takes the form \[\label{hat95f} \Hat{f}_j(\mathbf{u})= \sum_{k=1}^{N} c_k^{(j)} \exp(i\langle \boldsymbol{\omega}_k,\mathbf{u}\rangle).\tag{5}\]

To compute the coefficient vectors \(\mathbf{c}^{(j)}\in\mathbb{R}^N\) for all \(j\in[m]\), we solve the following min-norm interpolation problem \[\min_{\mathbf{c}^{(1)},\dots,\mathbf{c}^{(m)}\in\mathbb{R}^N} \sum_{j=1}^{m}\|\mathbf{c}^{(j)}\|_2 \qquad s.t.\Hat{f}(\mathbf{u}_\ell) = \mathbf{v}_\ell \quadfor all\ell\in[M]\] where \(\mathbf{u}_\ell\in\mathbb{R}^n\) and \(\mathbf{v}_\ell\in\mathbb{R}^m\) contain the function values at collocation points of the \(\ell\)-th data pair \((u_\ell,v_\ell)\in \mathcal{U}\times \mathcal{V}\), i.e \[\mathbf{u}_\ell = S_u(u_\ell) = [u_\ell(\mathbf{x}_1), \dots, u_\ell(\mathbf{x}_n)] \in\mathbb{R}^n,and\mathbf{v}_\ell = S_v(v_\ell) = [v_\ell(\mathbf{y}_1), \dots, v_\ell(\mathbf{y}_m)] \in\mathbb{R}^m.\] Due to Assumption 2, the training process can be done simultaneously in all components of \(f\). Therefore, we consider a sequence of parallel min-norm interpolation problems \[\label{parallel95c} \min_{\mathbf{c}^{(j)}\in\mathbb{R}^N} \|\mathbf{c}^{(j)}\|_2 \quad s.t.\mathbf{A}\mathbf{c}^{(j)} = \mathbf{v}^{(j)},\tag{6}\] where the matrix \(\mathbf{A}\in\mathbb{R}^{M\times N}\) is defined (component-wise) by \(\mathbf{A}_{\ell,k}=\exp(i\langle \boldsymbol{\omega}_k,\mathbf{u}_\ell\rangle)\) and \(\mathbf{v}^{(j)} = [v_1(\mathbf{y}_j),\dots,v_M(\mathbf{y}_j)]\in\mathbb{R}^M\). The solution to the min-norm interpolation problem is given by \(\mathbf{c}^{(j)} = \mathbf{A}^*(\mathbf{A}\mathbf{A}^*)^{-1}\mathbf{v}^{(j)} = \mathbf{A}^\dagger \mathbf{v}^{(j)}\) for each \(j\in[m]\), where \(\mathbf{A}^\dagger\) is the pseudo-inverse of \(\mathbf{A}\).

3 Random Feature Generalization Error↩︎

In this section, we provide bounds on the generalization error of approximating a function using random feature model.

3.1 Set-up and Notations↩︎

For a probability distribution \(\rho\) associated with random feature \(\boldsymbol{\omega}\) in \(\mathbb{R}^d\), we define function space \[\label{target95space} \mathcal{F}(\rho) := \left\{ f(\mathbf{x}) = \int_{\mathbb{R}^d} \hat{f}(\boldsymbol{\omega})\exp(i\langle \boldsymbol{\omega}, \mathbf{x}\rangle)d\boldsymbol{\omega}: \|f\|^2_\rho = \mathbb{E}_\omega \left|\frac{\hat{f}(\boldsymbol{\omega})}{\rho(\boldsymbol{\omega})} \right|^2 <\infty \right\},\tag{7}\] where \(\hat{f}(\boldsymbol{\omega})\) is the Fourier transform of function \(f\). Let \(D\subset\mathbb{R}^d\) be a compact domain with finite volume, i.e., \(\mathop{\mathrm{vol}}(D)<\infty\). Function \(f:D\to\mathbb{C}\) has finite \(\rho\)-norm if it belongs to \(\mathcal{F}(\rho)\). By Proposition 4.1 in [38], the function space \(\mathcal{F}(\rho)\) is a reproducing kernel Hilbert space (RKHS) with associated kernel function \[\label{kernel} k(\mathbf{x},\mathbf{y}) = \int_{\mathbb{R}^d} \exp(i\langle \boldsymbol{\omega}, \mathbf{x}\rangle)\exp(-i\langle \boldsymbol{\omega}, \mathbf{y}\rangle)d\rho(\boldsymbol{\omega}).\tag{8}\] We consider Cauchy random features which are drawn from the tensor-product Cauchy distribution with scaling parameter \(\gamma>0\), whose density function is \[\rho(\boldsymbol{\omega}) := \prod_{j=1}^d \frac{1}{\pi\gamma(1+\omega_j^2/\gamma^2) }.\] Here, \(\omega_j\) is the \(j\)-th entry of \(\boldsymbol{\omega}\in\mathbb{R}^d\). Using the tensor-product Cauchy distribution, the kernel function defined in 8 is indeed the Laplace kernel, i.e. \(k(\mathbf{x},\mathbf{y}) = \exp(-\gamma\|\mathbf{x}-\mathbf{y}\|_1)\). Hence, function space \(\mathcal{F}(\rho)\) defined in 7 is the reproducing kernel Hilbert space corresponding to Laplace kernel. The function space \(\mathcal{F}(\rho)\) also relates to the well-known Sobolev space of mixed smoothness, which is an appropriate function space for studying PDEs [43]. Precisely, let \(p>0\), we define the Sobolev spaces of mixed smoothness by [44], \[H^p_{\text{mix}}(\mathbb{R}^d):= \left\{f:D\rightarrow \mathbb{C}^d \,\, \Big\rvert \,\, \|f\|_{H^p_{\text{mix}}(\mathbb{R}^d)}< \infty \right\}\] where the associated norm is defined by: \[\|f\|^2_{H^p_{\text{mix}}(\mathbb{R}^d)}:=\int_{\mathbb{R}^d} \, |\hat{f}(\boldsymbol{\omega})|^2 \prod_{i=1}^d \left(1+|\omega_i|^2\right)^p d\boldsymbol{\omega}.\] If the scaling parameter of tensor-product Cauchy distribution \(\gamma=1\), then we have \(\|f\|_\rho = \|f\|_{H^1_{\text{mix}}}\). For arbitrary \(\gamma>0\), we could derive that \[\|f\|_\rho^2 \leq \|f\|^2_{H^1_{\text{mix}}(\mathbb{R}^d)} \max\left(\frac{\pi}{\gamma}, \pi\gamma\right)^d.\]

We consider the regression problem which is to find an approximation of \(f\in\mathcal{F}(\rho)\). Given a set of random weights \(\{ \boldsymbol{\omega}_k\}_{k\in[N]}\) which are i.i.d. samples from tensor-product Cauchy distribution \(\rho(\boldsymbol{\omega})\), we train a random feature model \[\label{rf95train} f^\sharp(\mathbf{x}) = \sum_{k=1}^N c_k^\sharp \exp(i\langle \boldsymbol{\omega}_k,\mathbf{x}\rangle)\tag{9}\] using samples \((\mathbf{x}_j, f(\mathbf{x}_j))\) for \(j\in[m]\). Let \(\mathbf{A}\in\mathbb{C}^{m\times N}\) be the random feature matrix with \(\mathbf{A}_{j,k}=\exp(i\langle \boldsymbol{\omega}_k, \mathbf{x}_j\rangle)\) for \(j\in[m]\) and \(k\in[N]\). The coefficient vector \(\mathbf{c}^\sharp\in\mathbb{R}^N\) is trained by solving the min-norm interpolation problem: \[\mathbf{c}^\sharp \in \mathop{\mathrm{argmin}}_{\mathbf{A}\mathbf{c}= \mathbf{y}} \|\mathbf{c}\|_2,\] where \(\mathbf{y}= [f(\mathbf{x}_1), \dots,f(\mathbf{x}_m)]\in\mathbb{R}^m\). Our main goal of this section is to bound the generalization error \(\|f-f^\sharp\|_{L^2(D)}\), which is defined as \[\|f-f^\sharp\|_{L^2(D)} = \sqrt{ \int_{D} |f(\mathbf{x})-f^\sharp(\mathbf{x})|^2 d\mathbf{x}}.\] We decompose the generalization error into two parts \[\|f-f^\sharp\|_{L^2(D)} \leq \underbrace{ \| f - f^* \|_{L^2(D)} }_{ \text{Approximation Error} } + \underbrace{ \|f^* - f^\sharp\|_{L^2(D)} }_{ \text{Estimation Error} }\] by triangular inequality. The approximation \(f^\sharp\) is the random feature model defined in 9 and \(f^*\) is the best random feature model that will be defined later. We first notice that we can write the target function \(f(\mathbf{x})\) as \[\label{target95f} f(\mathbf{x}) = \int_{\mathbb{R}^d} \hat{f}(\boldsymbol{\omega})\exp(i\langle \boldsymbol{\omega}, \mathbf{x}\rangle)d\boldsymbol{\omega}= \mathbb{E}_{\boldsymbol{\omega}\sim\rho} \left[ \frac{\hat{f}(\boldsymbol{\omega})}{\rho(\boldsymbol{\omega})}\exp(i\langle \boldsymbol{\omega}, \mathbf{x}\rangle) \right].\tag{10}\] To simplify the notation, we let \(\alpha(\boldsymbol{\omega}) = \frac{\hat{f}(\boldsymbol{\omega})}{\rho(\boldsymbol{\omega})}\) (defined on the \(\boldsymbol{\omega}\in \text{supp}(\rho)\) and zero otherwise) and then define \(\alpha_{\leq T}(\boldsymbol{\omega}) = \alpha(\boldsymbol{\omega})\mathbb{1}_{\left|\alpha(\boldsymbol{\omega})\right|\leq T}\). Therefore, we could write \(\alpha_{>T} = \alpha(\boldsymbol{\omega}) - \alpha_{\leq T}(\boldsymbol{\omega})\), i.e. \[\alpha_{\leq T}(\boldsymbol{\omega}) = \begin{cases} \alpha(\boldsymbol{\omega}) \quad &if\left|\alpha(\boldsymbol{\omega}) \right| \leq T \\ 0 \quad &otherwise\end{cases} \quadand\quad \alpha_{> T}(\boldsymbol{\omega}) = \begin{cases} \alpha(\boldsymbol{\omega}) \quad &if\left|\alpha(\boldsymbol{\omega}) \right| > T \\ 0 \quad &otherwise . \end{cases}\] We define \(f^*(\mathbf{x})\) as \[\label{rf95best} f^*(\mathbf{x}) = \frac{1}{N}\sum_{k=1}^N \alpha_{\leq T}(\boldsymbol{\omega}_k) \exp(i\langle \boldsymbol{\omega}_k,\mathbf{x}\rangle)\tag{11}\] where \(\boldsymbol{\omega}_k\) are i.i.d samples from distribution \(\rho(\boldsymbol{\omega})\). Note that the expectation is given by: \[\mathbb{E}f^*(\mathbf{x}) = \mathbb{E}_{\boldsymbol{\omega}} \left[ \alpha_{\leq T}(\boldsymbol{\omega})\exp(i\langle \boldsymbol{\omega},\mathbf{x}\rangle) \right].\]

3.2 Main Results↩︎

Before showing the main results of this section, we first introduce fill-in distance of training samples and state assumptions on the samples and random features. Specifically, the samples are well-separated and the random features are randomly generated from the tensor-product Cauchy distribution. The assumptions are more amenable to solving PDE problems where the Sobolev space is an appropriate space and collocation points are usually non-random.

We first state the main results of Section 3 where we provide a generalization error bound for random features model.

3.3 Approximation Error Bound↩︎

Theorem 2 provides an upper bound of the approximation error, which is the error between target function \(f\) and best random feature approximation \(f^*\).

The proof of Theorem 2 is based on the following two lemmas.

3.4 Estimation Error Bound↩︎

In this section, we aim to bound the estimation error \(\|f^*-f^\sharp\|_{L^2(D)}\). The error bound depends on the condition number of random feature matrix \(\mathbf{A}\), which can be bounded by using the concentration properties of random feature matrix. We refer readers to [31], [35] for more details about the concentration properties of random feature matrix. We first state the concentration results, and then derive an upper bound for the estimation error.

The direct consequence of this result is that each eigenvalue of \(\frac{1}{N}\mathbf{A}\mathbf{A}^*\) is close to 1. Specifically, if the conditions in Theorem 3 are satisfied, then \[\left| \lambda_k\left(\frac{1}{N}\mathbf{A}\mathbf{A}^*\right) -1 \right| \leq 2\eta\] with probability at least \(1-\delta\) for all \(k\in[m]\). Here, \(\lambda_k(\frac{1}{N}\mathbf{A}\mathbf{A}^*)\) is the \(k\)-th eigenvalue of matrix \(\frac{1}{N}\mathbf{A}\mathbf{A}^*\). Therefore, the matrix \(\mathbf{A}\mathbf{A}^*\) is invertible with high probability and hence the pseudo-inverse \(\mathbf{A}^\dagger\) is well-defined.

Utilizing the concentration property of random feature matrix \(\mathbf{A}\in\mathbb{C}^{m\times N}\), we derive an upper bound for the estimation error \(\|f^* - f^\sharp\|_{L_2(D)}\), which is summarized in Theorem 4.

Now, we are ready to prove Theorem 4 which states an upper bound for the estimation error. Recall the result of Lemma 3, it remains to bound \(\|\mathbf{A}\|_2^2\) and \(\|\mathbf{c}^*-\mathbf{c}^\sharp\|_2^2\).

4 Error Analysis for Random Feature Operator Learning↩︎

In this section, we present an error bound for our proposed random feature operator learning method. Let \(G\) be the target operator, the generalization error of estimator \(\Hat{G}\) is defined as \(\|G(u) - \Hat{G}(u)\|_{L^2(D_\mathcal{V})}\) for any \(u\in\mathcal{U}\), where \(\mathcal{U}\) is chosen to be the function space \(\mathcal{F}(\rho)\) defined in 7 . Using triangle inequality, we decompose the generalization error of our proposed estimator \(\Hat{G} = R_\mathcal{V}\circ\Hat{f}\circ S_\mathcal{U}\) into \[\begin{align} & \|G(u) - R_\mathcal{V}\circ\Hat{f}\circ S_\mathcal{U}(u)\|_{L_2(D_\mathcal{V})} \\ & \leq \underbrace{\|G(u) - R_\mathcal{V}\circ f\circ S_\mathcal{U}(u)\|_{L_2(D_\mathcal{V})}}_{\text{Approximation Error}} + \underbrace{ \|R_\mathcal{V}\circ f\circ S_\mathcal{U}(u) - R_\mathcal{V}\circ\Hat{f}\circ S_\mathcal{U}(u)\|_{L_2(D_\mathcal{V})} }_{\text{Estimation Error}}. \end{align}\]

To bound the estimation error, we need to introduce the following vector-valued random feature map \(\Bar{f}:\mathbb{R}^n\to\mathbb{R}^m\) where each \(\Bar{f}_j:\mathbb{R}^n\to\mathbb{R}\) is a random feature approximation taking the form \[\label{operator95temp} \Bar{f}_j(\mathbf{u}) = \sum_{k=1}^N \Bar{c}_k^{(j)} \exp(i\langle \boldsymbol{\omega}_k,\mathbf{u}\rangle),\tag{17}\] where the coefficient vector \(\Bar{\mathbf{c}}^{(j)}\in\mathbb{R}^N\) for all \(j\in[m]\) is trained by solving the min-norm interpolation problem \[\min_{\mathbf{c}^{(j)}\in\mathbb{R}^N} \|\mathbf{c}^{(j)}\|_2 \quads.t.\Bar{f}_j(\mathbf{u}_\ell) = f_j(\mathbf{u}_\ell) \quadfor all\ell\in[M].\] Notice that we can never compute \(\Bar{f}\) in practice since \(f:\mathbb{R}^n\to\mathbb{R}^m\) is not known to us.

Recall that \(R_v:\mathbb{R}^m\to\mathcal{V}\) is a random feature recovery map and denote \(f\circ S_\mathcal{U}(u)\), \(\Hat{f}\circ S_\mathcal{U}(u)\), and \(\Bar{f}\circ S_\mathcal{U}(u)\) by vectors \(\mathbf{z}\), \(\Hat{\mathbf{z}}\), and \(\Bar{\mathbf{z}}\in\mathbb{R}^m\), respectively. Then, we can write \[\begin{align} R_\mathcal{V}(\mathbf{z}) = \sum_{k=1}^N c_k \exp(i\langle \boldsymbol{\omega}_k,\mathbf{y}\rangle),R_\mathcal{V}(\Hat{\mathbf{z}}) = \sum_{k=1}^N \Hat{c_k} \exp(i\langle \boldsymbol{\omega}_k,\mathbf{y}\rangle),R_\mathcal{V}(\Bar{\mathbf{z}}) = \sum_{k=1}^N \Bar{c_k} \exp(i\langle \boldsymbol{\omega}_k,\mathbf{y}\rangle), \end{align}\] where \(\mathbf{c}= \mathbf{A}^\dagger_y \mathbf{z}\), \(\Hat{\mathbf{c}} = \mathbf{A}^\dagger_y\Hat{\mathbf{z}}\), \(\Bar{\mathbf{c}} = \mathbf{A}^\dagger_y \Bar{\mathbf{z}}\in\mathbb{R}^N\), and the matrix \(\mathbf{A}_y\) is defined (component-wise) by \((\mathbf{A}_y)_{j,k} = \exp(i\langle \boldsymbol{\omega}_k, \mathbf{y}_j\rangle)\) with \(\mathbf{A}^\dagger_y=\mathbf{A}_y^*(\mathbf{A}_y\mathbf{A}_y^*)^{-1}\) being its pseudo-inverse.

Let \(\mathcal{U}\) and \(\mathcal{V}\) are set of functions defined on compact domain \(D_\mathcal{U}\) and \(D_\mathcal{V}\), respectively. Distinct collocation points are denoted by \(X=\{\mathbf{x}_j\}_{j\in[n]}\subset D_\mathcal{U}\) and \(Y=\{\mathbf{y}_j\}_{j\in[m]}\subset D_\mathcal{V}\). The finite training data \(\{ (u_j,v_j) \}_{j\in[M]} \subset \mathcal{U}\times \mathcal{V}\) satisfies \(G(u_j) = v_j\) for all \(j\in[M]\) and is accessible via sampling operators \(S_\mathcal{U}\) and \(S_\mathcal{V}\), i.e., denote the point evaluations of function \(u_j\) and \(v_j\) by \[\mathbf{u}_j = S_\mathcal{U}(u_j) = [u_j(\mathbf{x}_1), \dots, u_j(\mathbf{x}_n)]^\top \in\mathbb{R}^n, \quad \mathbf{v}_j = S_\mathcal{V}(v_j) = [v_j(\mathbf{y}_1), \dots, v_j(\mathbf{y}_m)]^\top \in\mathbb{R}^m.\] Let \(U\subset \mathbb{R}^n\) be a compact set with Lipschitz boundary and \(\{ \mathbf{u}_j\}_{j\in[M]}\) be the set of all point evaluation vectors with fill-in distance \[h_\mathbf{u}:= \max_{\mathbf{u}'\in U}\min_{1\leq j\leq M} \|\mathbf{u}_j - \mathbf{u}'\|_2.\]

5 Numerical Experiments↩︎

Numerical experiments are used to benchmark and compare our proposed random feature operator learning method to kernel-based method [19] and neural operator methods, including DeepONet [3] and FNO [4]. We measure the model performance by using an empirical relative error on a test set \(U_{test}\): \[\mathrm{Error}(\Hat{G}) = \frac{1}{|U_{test}|} \sum_{u_k\in U_{test}} \frac{\|G(u_k)-\hat{G}(u_k)\|_{L^2(D_\mathcal{V})}}{\|G(u_k)\|_{L^2(D_\mathcal{V})}}\] where \(G:\mathcal{U}\to\mathcal{V}\) is the true operator. For any \(v\in\mathcal{V}\), we take \(\|v\|_{L^2(D_\mathcal{V})} = \sqrt{ \int_{D_\mathcal{V}} |v(x)|^2 dx}\), which in turn is estimated by the function values at collocation points. We outline the setup of each problem in the following subsections. In all problems, \(\mathcal{U}\) and \(\mathcal{V}\) are spaces of real-valued functions defined on domains \(D_\mathcal{U}, D_\mathcal{V}\subset\mathbb{R}^d\) for \(d=1,2\). All the experiments are conducted using Python and our source codes are available online³.

For the kernel method, we consider RBF kernel and Matérn kernel. Specifically, the RBF kernel function with parameter \(\gamma>0\) is \(K(\mathbf{x},\mathbf{x}') = \exp(-\gamma\|\mathbf{x}- \mathbf{x}'\|_2^2)\). The Matérn kernel with scaling parameter \(\sigma>0\) and smoothness parameter \(\nu>0\) is defined by \[K_{\nu, \sigma}(\mathbf{x},\mathbf{x}') = \frac{2^{1-\nu}}{\Gamma(\nu)} \left( \frac{\sqrt{2\nu}\|\mathbf{x}-\mathbf{x}'\|_2}{\sigma} \right)^\nu K_\nu\left( \frac{\sqrt{2\nu}\|\mathbf{x}-\mathbf{x}'\|_2}{\sigma} \right),\] where \(\Gamma\) is the gamma function, and \(K_\nu\) is the modified Bessel function of the second kind of order \(\nu\). The Gaussian random feature model with parameter \(\gamma>0\) means that the random features \(\boldsymbol{\omega}\) are drawn from normal distribution \(\mathcal{N}(0,2\gamma)\), which approximates RBF kernel with parameter \(\gamma\). The Cauchy random feature model with parameter \(\gamma>0\) means that the random features \(\boldsymbol{\omega}\) are drawn from tensor-product Cauchy distribution with scaling parameter \(\gamma>0\).

5.1 Advection Equations , , and↩︎

Consider the one-dimensional advection equation: \[\begin{align} \frac{\partial v}{\partial t} + \frac{\partial v}{\partial x} &= 0 \quad && x\in(0,1), t\in(0,1] \\ v(x,0) &= u(x) \quad && x\in(0,1), \end{align}\] we aim to learn the operator \(G\) mapping from the initial condition \(u(x)=v(x, 0)\) to the solution at time \(t=0.5\), which is denoted by \(v(\cdot, 0.5)\). The domains of function spaces \(\mathcal{U}\) and \(\mathcal{V}\) are \(D_\mathcal{U}= D_\mathcal{V}= (0,1)\). This problem was considered in [45], [46] with different distributions \(\mu\) for the initial condition. We consider, referred as Advection here, the initial condition which is a square wave centered at \(x = c\) of width \(b\) and height \(h\): \[u(x) = h1_{ \{c- 2b ,c + 2b\} }.\] The parameters \((c, b, h)\) are randomly generated following the uniform distribution on \([0.3, 0.7] \times [0.3, 06] \times [1, 2]\). The spatial grid is of resolution 40, and we use 1000 instances for training and 200 instances for testing. Figure 4 shows an example of training input, training output, true test example and its Cauchy random feature approximation, and the corresponding pointwise error.

Advection equation takes a more complicated initial condition, i.e. \[u(x) = h_1 1_{ \{c_1 - w, c_1 + w\} } + \sqrt{ \max(h_2^2 - a_2^2(x - c_2)^2, 0)}.\] The resolution is of 40 and we use 1000 samples to train the random feature models and another 1000 samples to test the performance. Similarly, we show an example of training input, training output, true test example, its Cauchy random feature approximation, and the corresponding pointwise error in Figure 5.

A different initial condition was considered in [46]. We call it Advection equation here. Suppose that \(\tilde{u}_0\) is generated from a Gaussian process \(GP(0,(-\Delta+32\mathbf{I})^{-2})\) where \(\Delta\) and \(\mathbf{I}\) represent the Laplacian and the identity, respectively. The initial condition is defined as \[u=-1+21_{ \{\tilde{u}_0 \geq 0\} }.\] The resolution is of 200 and we use 1000 samples for training and another 1000 instances for testing. In Figure 6, we show an example of training input and output, true test function, its Cauchy random feature approximation, and the corresponding pointwise error.

In Table 1, we compare the relative test errors and computational times of RBF kernel-based model, Gaussian random feature model, and Cauchy random feature model. The scaling parameter \(\gamma\) for each model and the number of features \(N\) for random feature models are included as well. We observe that our proposed random feature method is reliable in terms of accuracy and even outperforms vanilla kernel method. Moreover, the training time is reduced significantly. As discussed in [19], since the Advection tests have more jumps, the prediction task poses challenges for smooth models such as our model and the kernel models around the discontinuities.

5.2 Burgers’ Equation↩︎

Consider the one-dimensional Burgers’ equation: \[\begin{align} \frac{\partial w}{\partial t} + w\frac{\partial w}{\partial x} = \nu \frac{\partial^2 w}{\partial x^2},& \quad && (x,t)\in(0,1)\times (0,1] \\ w(x,0) = u(x),& \quad && x \in (0,1). \end{align}\] We set the viscosity parameter \(\nu=0.1\). Our goal is to learn the operator \(G:u(x)\mapsto w(x,1)\), which maps the initial condition \(u(x) = w(x,0)\) to the solution \(w(x,t)\) at time \(t=1\). The data are generated by sampling the initial condition \(u\) from a Gaussian Process \(\mu\sim GP(0, 625(-\Delta + 25I)^{-2}))\), where \(\Delta\) and \(I\) represent the Laplacian and the identity, respectively. As in [19], [45], we use a spatial resolution with 128 grid points to represent the input and output functions, and use 1848 instances for training and 200 instances for testing. Figure 7 shows an example of training input, training output, true test function and its approximation using Cauchy random feature model along with the pointwise error.

5.3 Darcy Problem↩︎

Consider the two-dimensional Darcy flow problem: \[\begin{align} -\mathop{\mathrm{div}} e^u \nabla v &= w, &&inD, \\ v &= 0, &&on\partial D \end{align}\] with \(D=(0,1)^2\) and zero Dirichlet boundary condition, we are interested in learning the operator \(G:u\mapsto v\) where \(u\) is the permeability field and \(v\) is the solution. Hence, the domains of function \(u\) and \(v\) are the same, i.e., \(D_\mathcal{U}= D_\mathcal{V}= (0,1)^2\). The source term \(w\) is assumed to be fixed. Here, we consider the piecewise constant permeability field. Specifically, the coefficient \(u\) is sampled from \(u \sim \log \circ h (\mu)\), where \(\mu\) is sampled from Gaussian Process \(\mu \sim GP(0, (-\Delta+9I)^{-2})\) and \(h\) is binary function mapping positive inputs to 12 and negative inputs to 3. Therefore, the permeability field \(e^u\) is piecewise constant, see an example of permeability field (training input) and the corresponding solution (training output) in Figure 8. The grid of resolution is \(29\times29\) and we use 1000 samples for training, and 200 samples for testing. For our random feature model, we randomly generate \(N=10^4\) random features from the tensor-product Cauchy distribution with scaling parameter \(\gamma=2\times10^{-4}\). Figure 8 shows one example of true test solution \(v\), its Cauchy random feature approximation, and the pointwise error.

5.4 Helmholtz’s Equation↩︎

We consider the Helmholtz PDE \[\label{Helmholtz} \begin{align} & \left( - \Delta - \frac{\omega^2}{u^2(x)} \right)v = 0, \quad x\in(0,1)^2 \\ & \frac{\partial v}{\partial n} = 0, \quad x\in\{0,1\}\times[0,1]\cup[0,1]\times\{0\} \quadand\quad \frac{\partial v}{\partial n} = v_N, \quad x\in[0,1]\times\{1\}, \end{align}\tag{20}\] where frequency \(\omega\) and wave speed field \(u: D_\mathcal{U}\to \mathbb{R}\) are assumed to be fixed and given. The excitation field \(v: D_\mathcal{V}\to \mathbb{R}\) solves equation 20 . The domains of function \(u\) and \(v\) are \(D_\mathcal{U}= D_\mathcal{V}= (0,1)^2\) in this example. In the numerical experiments, we take \(\omega=10^3\) and \(v_N=I_{\{0.35\leq x\leq 0.65\}}\). Our goal is to learn operator \(G:u\mapsto v\) from the wave speed field \(u\) to the excitation field \(v\). The wave speed field \(u\) takes the form \(u(x) = 20+\tanh(\tilde{u}(x))\), where \(\tilde{u}\) is drawn from Gaussian Process \(GP(0,(-\Delta+9I)^{-2})\). Here, denote \(\Delta\) by the Laplacian on \(D_\mathcal{U}\) subject to homogeneous Neumann boundary conditions on the space of spatial-mean zero functions. As described in [19], [46], samples are generated by solving the equation 20 using a Finite Element Method on a discretization grid of size \(26\times26\) of the unit square. We use 1000 functions for training and 1000 functions for testing. In Figure 9, we show an example of training input, training output, true test function, Cauchy random feature approximation, and the corresponding pointwise error.

5.5 Structural Mechanics↩︎

The system that governs the displacement vector \(w\) in an elastic solid undergoing infinitesimal deformations is defined as \[\begin{align} & \nabla\cdot \sigma = 0 \quad &&inD \\ & w=\Bar{w} \quad &&on\Gamma_w \\ & \sigma\cdot n = u \quad && on\Gamma_u \end{align}\] where \(\sigma\) is the (Cauchy) stress tensor and \(n\) is the outward unit normal vector. The computational domain is denoted by \(D=(0,1)^2\), and its boundary \(\partial D\) is split in \([0,1]\times 1 = \Gamma_u\) and its complement \(\Gamma_w\). The prescribed displacement \(\Bar{w}\) and the surface traction \(u\) are imposed on the domain boundaries \(\Gamma_w\) and \(\Gamma_u\), respectively. We aim to learn the operator that maps the one-dimensional load \(u\) on \(\Gamma_u\) to the two-dimensional von Mises stress field \(v\) on \(D_\mathcal{V}=(0,1)^2\). The load \(u\) is drawn from Gaussian Process \(GP(100, 400^2(-\Delta+\tau^2\mathrm{\mathop{Id}})^{-d})\) with \(\Delta\) being the Laplacian subject to homogeneous Neumann boundary conditions on the space of spatial-mean zero functions. The inverse length scale of the random field is taken to be \(\tau=3\) and \(d = 1\) determines its regularity (upto 1/2 a fractional derivative for samples from this measure). The load \(u\) is interpolated on a 41 grid and extruded in the \(y\) direction and the stress field is interpolated on a 21 \(\times\) 21 grid via radial basis function interpolation. We use 1000 functions for training and 1000 functions for testing. Figure 10 shows an example of training data, true test function and its approximation using Cauchy random features along with the corresponding pointwise error.

5.6 Navier-Stokes Equation↩︎

We consider the vorticity-stream \((\omega,\psi)\) formulation of the incompressible Navier-Stokes equation: \[\begin{align} & \frac{\partial w}{\partial t} + (c\cdot \nabla)w - v\Delta w = u \\ & w=-\Delta\psi \qquad \int_{D}\psi=0 \\ & c = (\frac{\partial\psi}{\partial x_2}, -\frac{\partial\psi}{\partial x_1}), \end{align}\] where computational domain is \(D=[0,2\pi]^2\) and periodic boundary conditions are considered. We are interested in learning the operator \(G:u\mapsto w(\cdot,T)\), where \(u\) is the forcing term \(u\) and \(w(\cdot,T)\) is the vorticity field at a given time \(t=T\). The forcing term \(u\) is sampled from Gaussian Process \(GP(0, (-\Delta+\tau^2\mathrm{\mathop{Id}})^{-d})\). As previous examples, \(\Delta\) denotes the Laplacian on \(D\) subject to periodic boundary conditions on the space of spatial-mean zero functions, \(\tau=3\) denotes the inverse length scale of the random field and \(d = 4\) determines its regularity; the choice of \(d\) then leads to up to 3 fractional derivatives for samples from this measure. We fix the initial condition \(w(\cdot,0)\) which is generated from the same distribution. Given the constant viscosity \(\nu=0.025\), the equation is solved on a \(16\times16\) grid with a pseudo-spectral method and Crank-Nicholson time integration. The size of the training dataset is 2000 and the test dataset is of size 2000. Figure 11 shows an example of training data, true test function, its approximation using Cauchy random features, and the corresponding pointwise error.

5.7 Results↩︎

In Table 2, we report the relative test errors of random feature methods and of kernel methods. Even though random feature methods can be viewed as an approximation of a kernel method, we observe that the random feature method achieves similar relative test errors or even improves the accuracy over the kernel method in some tests. Both the random feature method and kernel method are easy to implement in practice. However, the implementation of kernel methods requires high-performance numerical computing and machine learning python library such as JAX, see [19]. The implementation of our proposed random feature method can be done by using scikit-learn library and GPU computing is not required. Therefore, our proposed algorithm achieves similar accuracy even if the computational resources are limited. As concluded in [19], the kernel method either matches the performance of neural network methods or outperforms them in several benchmarks. Therefore, the random feature method can match the performance or outperform neural network methods as well.

We further compare the performances of the kernel method and the random feature method provided that each model is trained over a similar time period. Since the training time of the kernel method depends on the number of training samples, we further subsample a small dataset to train the kernel methodsİn Table [Table:RF95kernel2], we report the relative test errors and training times of kernel methods and of random feature methods. We observe that our proposed random feature method outperforms kernel method cross all benchmarks provided that the same amount of training time is allocated to each model.

We also compare our proposed random feature method with DeepONet and FNO directly in terms of the accuracy, see the summary of test relative errors in Table [Table:Test95error]. We report the training times of random feature methods as well. The relative errors of DeepONet and FNO are cited from [19], [45], [46]. We observe that the random feature methods can be trained fast even if the problem is complicated. For example, there are 20000 training samples in the Navier-Stokes problem and each training function is interpolated on a \(41\times 41\) grid. The random feature method matches the performance of DeepONet and FNO in the structural mechanics example, and outperforms DeepONet in the Helmholtz example, and outperforms both DeepONet and FNO in the Navier-Stokes example.

The choice of discretizations is important in operator learning. A true operator learning method should demonstrate the invariance of the expected relative test error to the mesh resolution used for training and testing. We empirically verify that our proposed method is stable to discretization changes, and hence is discretization-invariant, see also [47] for a discussion on discretization-invariant operator networks. We present the relative test errors at different resolutions for Darcy problem and Helmholtz problem in Tables 3 and 4, respectively. For each problem, we use 1000 functions for training and 1000 functions for testing.

In Figure 12, we empirically verify the error bound we obtained in Theorem 5 using some benchmarks. As the number of random features \(N\) increases, we observe that the test errors are all convergent as a rate of \(1/\sqrt{N}\). We also observe that the growth rate of training time is between \(\sqrt{N}\) and \(N\).

6 Conclusion↩︎

We propose a random feature method for learning nonlinear Lipschitz operators related to PDE. We provide a detailed error analysis for our proposed method with Cauchy random features. The theory suggests that the generalization error decays as a rate of \(1/\sqrt{N}\), where \(N\) is the number of random features. The results hold for both Gaussian and Cauchy distributions. Numerical experiments show that the random feature method not only enjoys the easy and fast implementation, but also matches or outperforms kernel method and neural network methods in benchmark 2D PDE examples. In addition, the theoretical and computational analysis shows the benefits of randomization. The use of Cauchy features supports the application to PDE. These methods may also be applicable to the recent multi-operator learning approaches, for example, the transformer based models [48]–[52] or DeepONet based approaches [53]. Future work could explore other heavy tail random features and noisy data.

Acknowledgments↩︎

HS was partially supported by NSF DMS 2331033. DN was partially supported by NSF DMS 2408912.

References↩︎