Solving Partial Differential Equations with Random Feature Models

To introduce random feature method, we first give a short introduction to kernel method, which is also known as kernel trick. It is one of the popular techniques for capturing nonlinear relations between features and targets. Let $\mathbf{x},\mathbf{x}^{\prime}\in X\subset\mathbb{R}^{d}$ be samples and $\phi:X\to\mathcal{H}$ be a feature map transforming samples to a high-dimensional (even infinite-dimensional) reproducing kernel Hilbert space $\mathcal{H}$ where the mapped data can be learned by a linear model. In practice, the explicit expression of feature map $\phi$ is not necessarily known to us. The inner produce between $\phi(\mathbf{x})$ and $\phi(\mathbf{x}^{\prime})$ endowed by $\mathcal{H}$ can be computed by using a kernel function $k(\cdot,\cdot):X\times X\to\mathbb{R}$, i.e.

$$k(\mathbf{x},\mathbf{x}^{\prime})=\langle\phi(\mathbf{x}),\phi(\mathbf{x}^{\prime})\rangle_{\mathcal{H}}.$$

Due to the ease of computing the inner product, kernel method is effective for nonlinear learning problems with a wide range of successful applications [9]. However, kernel method does not scale well to extremely large datasets. For example, given $m$ training samples, kernel regression requires $\mathcal{O}(m^{3})$ training time and $\mathcal{O}(m^{2})$ to store the kernel matrix, which is often computationally infeasible when $m$ is large. Random feature method is one of the most popular techniques to overcome the computational challenges of kernel method. The theoretical foundation of random (Fourier) feature builds on the following classical result from harmonic analysis.

Bochner’s theorem guarantees that the Fourier transform of kernel $k(\delta)$ is a proper probability distribution if the kernel is scaled properly. Denote the probability distribution by $\rho(\boldsymbol{\omega})$, we have

$$k(\mathbf{x},\mathbf{x}^{\prime})=\int_{\mathbb{R}^{d}}\exp(i\langle\boldsymbol{\omega}_{k},\mathbf{x}-\mathbf{x}^{\prime}\rangle)d\rho(\boldsymbol{\omega})=\int_{\mathbb{R}^{d}}\exp(i\langle\boldsymbol{\omega}_{k},\mathbf{x}\rangle)\overline{\exp(i\langle\boldsymbol{\omega}_{k},\mathbf{x}^{\prime}\rangle)}d\rho(\boldsymbol{\omega}).$$

(1)

Using Monte Carlo sampling technique, we randomly generate $N$ i.i.d samples $\{\boldsymbol{\omega}_{k}\}_{k\in[N]}$ from $\rho(\boldsymbol{\omega})$ and define a random Fourier feature map $\phi:\mathbb{R}^{d}\to\mathbb{C}^{N}$ as ²²2It depends on the i.i.d samples $\{\boldsymbol{\omega}_{k}\}_{k\in[N]}$ from $\rho(\boldsymbol{\omega})$. To simplify the notation, we omit the dependency on $\{\boldsymbol{\omega}_{k}\}_{k\in[N]}$ when we define $\phi$.

$$\phi(\mathbf{x})=\frac{1}{\sqrt{N}}\Big{[}\exp(i\langle\boldsymbol{\omega}_{1},\mathbf{x}\rangle),\dots,\exp(i\langle\boldsymbol{\omega}_{N},\mathbf{x}\rangle)\Big{]}^{T}\in\mathbb{C}^{N}.$$

(2)

Using random Fourier feature map, we can define a kernel function $\hat{k}(\mathbf{x},\mathbf{x}^{\prime}):X\times X\to\mathbb{R}$ as

$$\hat{k}(\mathbf{x},\mathbf{x}^{\prime}):=\langle\phi(\mathbf{x}),\phi(\mathbf{x}^{\prime})\rangle.$$

Kernel function $\hat{k}$ is finite dimensional since the random Fourier feature map $\phi$ transforms data to a finite dimensional space $\mathbb{C}^{N}$. We can also use random cosine features to approximate any shift-invariant kernel. Specifically, setting random feature $\boldsymbol{\omega}\in\mathbb{R}^{d}$ generated from $\rho(\boldsymbol{\omega})$ and $b\in\mathbb{R}$ sampled from uniform distribution on $[-\pi,\pi]$, the random cosine feature is defined as $\cos(\langle\boldsymbol{\omega},\mathbf{x}\rangle+b)$. Similarly, we can define a random cosine feature map $\phi:\mathbb{R}^{d}\to\mathbb{R}^{N}$ as

$$\phi(\mathbf{x})=\frac{1}{\sqrt{N}}\Big{[}\cos(\langle\boldsymbol{\omega}_{1},\mathbf{x}\rangle+b_{1}),\dots,\cos(\langle\boldsymbol{\omega}_{N},\mathbf{x}\rangle+b_{N})\Big{]}^{T}\in\mathbb{R}^{N},$$

(3)

using random i.i.d samples $\{\boldsymbol{\omega}_{k},b_{k}\}_{k\in[N]}$, and hence define a kernel function $\hat{k}:X\times X\to\mathbb{R}$ using random cosine feature map. We could approximate the original shift-invariant kernel function $k$ defined in (1) by the finite-dimensional kernel $\hat{k}$. Utilizing random features allows efficient learning with $\mathcal{O}(mN^{2})$ time and $\mathcal{O}(mN)$ storage capacity.

Now, we set ourselves to the regression setting where our aim is to learn a function $f:\mathbb{R}^{d}\to\mathbb{R}$ using training samples $\{(\mathbf{x}_{j},y_{j})\}_{j\in[m]}$. To implement the random feature model, we first draw $N$ i.i.d random features $\{\boldsymbol{\omega}_{k}\}_{k\in[N]}\subset\mathbb{R}^{d}$ from a probability distribution $\rho(\boldsymbol{\omega})$, and then construct an approximation for target function $f$ taking the form

$$f^{\sharp}(\mathbf{x})=\sum_{k=1}^{N}c_{k}^{\sharp}\phi(\mathbf{x},\boldsymbol{\omega}_{k}).$$

(4)

We assume that the $m$ sampling points $\mathbf{x}_{j}$’s are drawn from a certain distribution with the corresponding output values

$$y_{j}=f(\mathbf{x}_{j})+e_{j},\quad\mbox{ for all }j\in[m],$$

where $e_{j}$ is the measurement noise. Let $\mathbf{A}\in\mathbb{R}^{m\times N}$ be the random feature matrix defined component-wise by $\mathbf{A}_{j,k}=\phi(\mathbf{x}_{j},\boldsymbol{\omega}_{k})$ for $j\in[m]$ and $k\in[N]$. Training the random features model (4) is equivalent to finding the coefficient vector $\mathbf{c}^{\sharp}\in\mathbb{R}^{N}$ such that $\mathbf{A}\mathbf{c}^{\sharp}\approx\mathbf{y}$, where $\mathbf{c}^{\sharp}=[c_{1}^{\sharp},\dots,c_{N}^{\sharp}]^{\top}\in\mathbb{R}^{N}$ and $\mathbf{y}=[y_{1},\dots,y_{m}]\in\mathbb{R}^{m}$.

In the under-parameterized regime where we have more measurements than features ($m\geq N$), the coefficients are trained by solving the (regularized) least squares problem:

$$\mathbf{c}_{\lambda}^{\sharp}\in\mathop{\mathrm{argmin}}_{\mathbf{c}\in\mathbb{R}^{N}}\|\mathbf{A}\mathbf{c}-\mathbf{y}\|_{2}^{2}+\lambda\|\mathbf{c}\|_{2}^{2},$$

where $\lambda>0$ is the regularization parameter. It is also referred to ridge regression since the ridge regularization term $\lambda\|\mathbf{c}\|_{2}^{2}$ is added.

Recently, over-parametrized models have received great attention since those trained models not only fit the training samples exactly but also predict well on unseen test data [11, 12]. In the over-parametrized regime, we have more features than measurements ($N\geq m$), and we consider training the coefficient vector $\mathbf{c}^{\sharp}\in\mathbb{R}^{N}$ using the min-norm interpolation problem:

$$\mathbf{c}^{\sharp}\in\mathop{\mathrm{argmin}}_{\mathbf{c}\in\mathbb{R}^{N}}\|\mathbf{c}\|_{2}^{2}\quad\mbox{ subject to }\mathbf{A}\mathbf{c}=\mathbf{y}.$$

This problem is also referred to ridgeless regression problem since the solution $\mathbf{c}^{\sharp}$ can be viewed as the limit of $\mathbf{c}_{\lambda}^{\sharp}$ as $\lambda\to 0$.

The generalization analysis of random features models have been of recent interest [13, 14, 15, 16, 17, 18, 19]. In [13], the authors showed that the random feature model yields a test error of $\mathcal{O}(N^{-\frac{1}{2}}+m^{-\frac{1}{2}})$ when trained on Lipschitz loss functions. Therefore, the generalization error is $\mathcal{O}(N^{-\frac{1}{2}})$ for large $N$ if $m\asymp N$. However, the model is trained by solving a constrained optimization problem which is rarely used in practice. In [14], it was shown that for $f$ in an RKHS, using $N=\mathcal{O}(\sqrt{m}\log(m))$ features is sufficient to achieve a test error of $\mathcal{O}(m^{-\frac{1}{2}})$ with squared loss. In [15], the authors showed that a regularized model can achieve $N^{-1}+m^{-\frac{1}{2}}$ risk provided that the target function belonging to an RKHS. Nevertheless, results in [14, 15] depend on the assumptions of kernel and a certain decay rate of second moment operator, which may be difficult to verify in practice. Extending the results of random feature models from squared loss to 0-1 loss, the authors of [16] showed that the support vector machine with random features $N\ll m$ can achieve the learning rate faster than $\mathcal{O}(m^{-\frac{1}{2}})$ on a training set with $m$ samples. In [17], the authors computed the precise asymptotic bound of the test error, in the limit $N,m,d\to\infty$ with $N/d$ and $m/d$ fixed. In [18], the authors derived non-asymptotic bounds including both the under-parametrized setting using (regularized) least square problem and the over-parametrized setting using min-norm problem or sparse regression. Their results relied on the condition numbers of the random feature matrix and indicated double descent behavior in random feature models. However, the target function space is a subset of a RKHS, which limits the approximation ability of random feature model. In [19], the authors consider a RKHS as target function space and derived a similar non-asymptotic bound by utilizing different proof techniques.

Furthermore, the random feature model can be viewed as a two-layer (one hidden layer) neural network where the weights connecting input layer and hidden layer and biases are sampled randomly and independently from a known distribution. Only the weights of the output layer are trainable using training samples, and hence it leads to solving a convex optimization problem when training random feature models. There are other types of randomized neural network, including the random vector functional link (RVFL) network [20, 21], and the extreme learning machine (ELM) [22, 23, 24], among others. Sharing the same structure as the random feature model, RVFL network is also a shallow neural network where the input-to-hidden weights and biases are randomly selected. However, the motivations of two models are different. RVFL networks were designed to address the difficulties associated with training deep neural networks and weights are usually sampled from uniform distribution. Random feature models were originally used to approximate large-scale kernel machines, and hence random weights depend on the kernel function, see Table 1 for some examples of commonly used kernels and the corresponding densities. Extreme learning machine is one type of deep neural network (more than two hidden layers) in which all the hidden-layer weights are randomly selected and then fixed. Only the output-layer coefficients are trained. Theoretical guarantees on RVFL and ELM suggest that they are universal approximators, see [25].

In this section, we present our framework for solving PDEs using random feature models. Let us consider the PDE problem of the general form

		$\displaystyle\mathcal{P}[u](\mathbf{x})=0,$		$\displaystyle\mathbf{x}\in\Omega$		(5)
		$\displaystyle\mathcal{B}[u](\mathbf{x})=0,$		$\displaystyle\mathbf{x}\in\partial\Omega,$

where $\Omega\subset\mathbb{R}^{d}$ is the domain with the boundary $\partial\Omega$, $\mathcal{P}$ is the interior differential operator and $\mathcal{B}$ is the boundary differential operator. For the sake of brevity, we assume that the PDE is well-defined pointwise and has a unique strong solution throughout this paper. We propose to solve the PDE (5) by using random feature model. More precisely, let $\{\mathbf{x}_{j}\}_{j\in[M]}$ be a collection of collocation points such that $\{\mathbf{x}_{j}\}_{j\in[M_{\Omega}]}$ is a collection of points in the interior of $\Omega$ and $\{\mathbf{x}_{j}\}_{j=M_{\Omega}+1}^{M}$ is a set of points on the boundary $\partial\Omega$. Random features $\{\boldsymbol{\omega}_{k}\}_{k\in[N]}$ are randomly drawn from a known distribution $\rho(\boldsymbol{\omega})$. The random feature model takes the form

$$u^{\sharp}(\mathbf{x})=\sum_{k=1}^{N}c_{k}^{\sharp}\phi(\mathbf{x},\boldsymbol{\omega}_{k})$$

(6)

We train the random feature model by solving the following optimization problem:

		$\displaystyle\mathop{\mathrm{minimize}}_{\mathbf{c}\in\mathbb{R}^{N}}$		$\displaystyle\|\mathbf{c}\|_{2}^{2}$		(7)
		s.t.		$\displaystyle\mathcal{P}[u^{\sharp}](\mathbf{x}_{j})=0,\quad\mbox{ for }j=1,\dots,M_{\Omega}$
		$\displaystyle\mathcal{B}[u^{\sharp}](\mathbf{x}_{j})=0,\quad\mbox{ for }j=M_{\Omega}+1,\dots,M$

We wish to find an approximation of the true solution $u$ with the min-norm random feature model satisfying the PDE and boundary data at collocation points.

However, addressing (7) directly can be complicated if the PDE is nonlinear. Therefore, we may solve an unconstrained optimization problem with regularization instead,

$$\mathop{\mathrm{minimize}}_{\mathbf{c}\in\mathbb{R}^{N}}\|\mathbf{c}\|_{2}^{2}+\lambda_{1}\sum_{j=1}^{M_{\Omega}}\left(\mathcal{P}[u^{\sharp}](\mathbf{x}_{j})\right)^{2}+\lambda_{2}\sum_{j=M_{\Omega}+1}^{M}\left(\mathcal{B}[u^{\sharp}](\mathbf{x}_{j})\right)^{2},$$

(9)

where $\lambda_{1},\lambda_{2}>0$ are regularization parameters. When $\lambda_{1},\lambda_{2}\to 0$, the solution of (9) converges to the solution of (7). In practice, we use variants of stochastic gradient descent and the modern automatic differential libraries to solve problem (9). In scenarios where PDEs are challenging, a large number of collocation points are required to capture the solution details. Compared with the framework using the standard kernel matrix to construct the solution approximation in [3], our proposed random feature method approximates the kernel matrix, and hence reduces the computational cost and accelerate computation involving the standard kernel matrix (and its inverse) when dealing with massive collocation points.

In this section, we show the convergence analysis of our method. Our convergence analysis relies on the standard convergence analysis of kernel method in [6] and the kernel approximation by using random features. Recall the kernel-based method for solving PDEs in [5, 6], a reproducing kernel Hilbert space $\mathcal{H}$ is chosen and we aim to solve the following

		$\displaystyle\mathop{\mathrm{minimize}}_{u\in\mathcal{H}}$		$\displaystyle\|u\|_{\mathcal{H}}$		(10)
		s.t.		$\displaystyle\mathcal{P}[u^{\sharp}](\mathbf{x}_{j})=0,\quad\mbox{ for }j=1,\dots,M_{\Omega}$
		$\displaystyle\mathcal{B}[u^{\sharp}](\mathbf{x}_{j})=0,\quad\mbox{ for }j=M_{\Omega}+1,\dots,M$

We first state the main assumptions on the domain $\Omega$ and its boundary $\partial\Omega$, the PDE operators $\mathcal{P}$ and $\mathcal{B}$, and the reproducing kernel Hilbert space $\mathcal{H}$.

Item (C1) is a standard assumption when analyzing PDEs. Item (C2) assumes that the PDE to be Lipschitz well-posed with respect to the right hand side/source term. It relates to the analysis of nonlinear PDEs and is independent of our numerical scheme. Assumption (C3) dictates the choice of the RKHS $\mathcal{H}$ , and in turn the kernel, which should be carefully selected based on the regularity of the strong solution $u$. We are now ready to state the first theorem which concerns the convergence rate of kernel method

With the kernel minimizer $\hat{u}\in\mathcal{H}$ at hand, our next step is approximating $\hat{u}$ using random features. We first adopt an alternative representation of preselected RKHS $\mathcal{H}$. Denote the corresponding random Fourier feature map (or random cosine feature map) by $\phi:X\to\mathbb{C}^{N}(\mathbb{R}^{N})$, then we define the following function space

$$\mathcal{F}(\rho):=\left\{f(\mathbf{x})=\int_{\mathbb{R}^{d}}\alpha(\boldsymbol{\omega})\phi(\mathbf{x},\boldsymbol{\omega})d\rho(\boldsymbol{\omega}):\|f\|^{2}_{\rho}=\mathbb{E}_{\boldsymbol{\omega}}[\alpha(\boldsymbol{\omega})^{2}]<\infty\right\},$$

(11)

where $\rho(\cdot)$ is the Fourier transform density associated with kernel $k$. Notice that the completion of $\mathcal{F}(\rho)$ is a Hilbert space equipped with RKHS norm $\|f\|_{\rho}$. Recall the Proposition 4.1 in [26], it is indeed the reproducing kernel Hilbert space $\mathcal{H}$ with associated kernel function $k$. The endowed norms $\|\cdot\|_{\mathcal{H}}$ and $\|\cdot\|_{\rho}$ are equivalent.

In the next theorem, we address the approximation ability of finite sum random feature model taking the form (6).

The last theorem in this section presents the convergence rate of our proposed random feature model.

We test with the instance of nonlinear elliptic PDE

	$\displaystyle-\Delta u(\mathbf{x})+u(\mathbf{x})^{3}$	$\displaystyle=f(\mathbf{x}),\quad$		$\displaystyle\mathbf{x}\in\Omega$
	$\displaystyle u(\mathbf{x})$	$\displaystyle=g(\mathbf{x}),\quad$		$\displaystyle\mathbf{x}\in\partial\Omega,$

where $\Omega=[0,1]^{2}$. The solution $u(\mathbf{x})=\sin(\pi x_{1})\sin(\pi x_{2})+4\sin(4\pi x_{1})\sin(4\pi x_{2})$, and the right-hand side $f(\mathbf{x})$ is computed accordingly via the solution $u(\mathbf{x})$. The boundary condition is $g(\mathbf{x})=0$. This example was also considered in [5, 3]. We randomly sample 1000 random features ($N=1000$) from normal distribution $\mathcal{N}(0,\sigma^{2}\mathbf{I})$ with variance $\sigma^{2}=100$. The random feature model is trained on $M_{\Omega}=900$ interior collocation points and 124 collocation points on the boundary $\partial\Omega$. We take the uniform grids of size $100\times 100$ as our test points. In Figure 1, we show predicted solution using our proposed random feature model, true solution, and entry-wise absolute errors at test points. We observe that our proposed random feature model provides an accurate prediction of the true solution.

We also compare the training epochs, test errors and training times of our proposed method with PINN, see results in Table 2. The PINN model has 2 hidden layers and each layer has 64 neurons. The nonlinear activation function is tanh function. PINN is trained on the same training samples and is tested on the uniform grid as well. We observe that training PINN is complicated, which requires more epochs, and hence longer training time. If we set the same number of epochs, then the PINN gives a bad prediction compared with random feature method. Moreover, our proposed method has smaller test error. Overall, our proposed random feature model outperforms PINN in this example.

Finally, we numerically verify the convergence rate of nonlinear elliptic PDE. We first fix the number of random features to be $N=100$ and varies the number of collocation points. We sample $M_{\Omega}=400,900,1600$ points uniformly in the domain, and $M_{\partial\Omega}=84,124,164$ points uniformly on the boundary. We sample another set of 100 test points and evaluate the test errors. In the second experiment, we fix $M_{\Omega}=400$ interior collocation points and $M_{\partial\Omega}=84$ points on the boundary. We take different numbers of random features, i.e. $N=100,200,300$. We sample another set of 100 test points and evaluate the test errors. In Figure 2, we show the test errors as a function of the number of collocation points and a function of the number of random features, respectively. For each point in the figure, we repeat the experiments 10 times and take the average.

In this section, we test our proposed method with high-dimensional nonlinear Poisson PDEs. We consider the domain $\Omega=[-1,1]^{d}$ and the following problem defined on $\Omega$,

	$\displaystyle-\nabla\cdot(a(u)\nabla u)$	$\displaystyle=f(\mathbf{x}),\quad$		$\displaystyle\mathbf{x}\in\Omega$
	$\displaystyle u(\mathbf{x})$	$\displaystyle=g(\mathbf{x}),\quad$		$\displaystyle\mathbf{x}\in\partial\Omega,$

where $a(u)=u^{3}-u$. The solution is crafted as $u(\mathbf{x})=\exp(-\frac{1}{d}\sum_{i=1}^{d}x_{i})$. The function $g(\mathbf{x})=u(\mathbf{x})$ on the boundary, and the right-hand side $f(\mathbf{x})$ is computed using the true solution, which has the explicit expression

$$f(\mathbf{x})=\frac{1}{d}\left[-3\exp\left(-\frac{3}{d}\sum_{i=1}^{d}x_{i}\right)+2\exp\left(-\frac{2}{d}\sum_{i=1}^{d}x_{i}\right)\right].$$

We first test with the 2D nonliear Poisson PDE and show the predicted solution, true solution, and entry-wise absolute errors in Figure 3. The training samples of size 1024 contains 900 interior collocation points and 124 boundary points, which are uniformly generated over the domain $\Omega=[-1,1]^{2}$ and its boundary $\partial\Omega$, respectively. We randomly generate 500 test samples to evaluate the performance. We generate 1000 random features, following the standard Normal distribution, to construct the random feature model.

Furthermore, we compare between our proposed model and PINN model. The neural network has 2 hidden layers and 64 neurons at each layer. We select tanh function as the activation function. We test 2D nonlinear Poisson PDE as well as high-dimensional nonlinear Poisson PDEs ($d=4$ and $d=8$). We aim to compare the test errors and training times. We also report the number of random features for our model and the number of epochs for both models. All numerical results are summarized in 3. We observe that our proposed random feature model achieves similar performance or even beats the PINN models in terms of the test error.

We note that the accuracy of our model is related to the choice of variance of Gaussian random feature, but it is not sensitive to the variance. Usually, we can tune the hyperparameter using cross-validation. In Table 4, we report the variances and the corresponding test errors for $d=8$ dimension nonlinear Poisson PDE. In this example, it is better to use small variance, but the test error is not very sensitive to the choice of variance when it is smaller than some threshold.

In Figure 4, we report the empirical convergence rates for the nonlinear Poisson equations. Figure 4(a) shows the test error as a function of the number of collocation points. For each dimension, we fix $N=100$ random features and varies the number of collocation points. We uniformly sample $M_{\Omega}=100,400,900$ interior points, and $M_{\partial\Omega}=44,84,124$ boundary points. We sample a different set of 100 test points to evaluate the test errors. Figure 4(b) shows the test error as a function of the number of random features. For each dimension, we fix the collocation points of size 484 with 400 interior collocation points. We take $N=100,400$ random features. We average over 10 experiments to produce each point in Figure 4.

Next, we consider a 2D stationary Allen-Cahn equation with a source function and Dirichlet boundary conditions, i.e.

$$\Delta u+\gamma(u^{m}-u)=f(\mathbf{x}),\quad\mathbf{x}\in[0,1]^{2},$$

where $\gamma=1$ and $m=3$. The solution takes the form $u(\mathbf{x})=\sin(2\pi ax_{1})\cos(2\pi ax_{2})$, and the function $f(\mathbf{x})$ is computed using the solution $u(\mathbf{x})$. Positive parameter $a$ controls the frequency of the solutions. We test three cases $a=1$, $a=10$, and $a=20$.

We first compare the performance between random feature method and PINN. In each case, we randomly sample $M=1024$ points with $M_{\Omega}=900$ interior collocation points to train models. The test performances for both models are evaluated on 100 test samples, which are uniformly generated over the domain and boundary. We report the parameter selections and summarize the numerical results in Table 5. From the numerical results, we observe that the random feature models beat PINNs cross all cases, especially when the frequency parameter $a$ is large. It might be related to the known ”spectral bias” of neural networks [28, 29, 30]. Precisely, neural networks trained by gradient descent fit a low frequency function before a high frequency one. Therefore, it is difficult for PINNs to learn the high frequency PDE solutions. To alleviate ”spectral bias” and learn high frequency functions, previous work proposed to use Fourier features and both theoretically and empirically showed that a Fourier feature mapping can improve the performance [31]. Fourier features have been used to solve high frequency PDEs, see [32]. As our results suggest, higher frequency in the PDE solutions leads to a larger variance of Gaussian random feature. Moreover, it requires more random features and epochs to train the random feature model as the frequency increasing. In Figure 5, we show predicted solution, true solution and entry-wise errors at test points.

Figure 6 illustrates the numerical verifications of the convergence rate of Allen-Cahn equation. We first show the test error as a function of the number of collocation points. In this experiment, we fix the number of random features to be $N=100$. To produce the collocation points, we uniformly sample $M_{\Omega}=400,900,1600$ points in the domain, and $M_{\partial\Omega}=84,124,164$ points on the boundary. We sample another set of 100 test points to evaluate the test errors. In the second experiment, where we show the test error as a function of the number of random features, we uniformly sample and then fix $M_{\Omega}=400$ interior collocation points and $M_{\partial\Omega}=124$ points on the boundary. We sample another set of 100 test points and evaluate the test errors. We generate $N=100,200,400$ random features from Gaussian distribution. In Figure 6, we show the test errors as a function of the number of collocation points and a function of the number of random features, respectively. For each point in the figure, we repeat the experiments 10 times and take the average.

Finally, we test our method with advection diffusion equation. Consider the initial boundary value problem on the spatial-temporal domain $(x,t)\in[-1,1]\times[0,1]$, the PDE is

	$\displaystyle u_{t}-u_{xx}+u_{x}$	$\displaystyle=f(x,t),\quad$		$\displaystyle(x,t)\in[-1,1]\times[0,1]$
	$\displaystyle u(x,t)$	$\displaystyle=g(x,t),\quad$		$\displaystyle(x,t)\in\{-1,1\}\times[0,1],$
	$\displaystyle u(x,0)$	$\displaystyle=h(x),\quad$		$\displaystyle x\in[-1,1].$

The true solution is employed as $u(x,t)=\sin(x)\exp(-t)$. Functions $f(x,t)$, $g(x,t)$, and $h(x)$ are set according to the true solution. When we simulate this problem, we treat the time variable $t$ in the same way as the spatial variable $x$. We uniformly generate 1000 training samples over $[-1,1]\times[0,1]$. We enforce the boundary condition on 100 collocations points and the initial condition on 200 collocations points. Gaussian random features ($N=100$) are randomly sampled from standard Normal distribution (variance $\sigma^{2}=1$). The PINN model has 2 hidden layers with 64 neurons at each layer. We report number of epochs, the test errors and training times for both models in Table 6. In this example, our proposed method achieves similar test error as PINN, but the training is simpler in the sense that it requires less epochs and the training time of random feature model is around 50% of that of PINN model. In Figure 7, we compare the predicted solution and true solution at various times $t=0.1,0.5,0.9$ to further highlight the ability of our method in learning the true solution.

Finally, we perform a convergence study for advection-diffusion equation. In Figure 8(a), we show the test errors as a function of the number of collocation points. We use $M=100,200,400$ collocation points, which are uniformly generated from $[-1,1]\times[0,1]$. We use 100 points on the boundary and 200 points for initial values. In Figure 8(b), we use 200 interior points, 100 boundary points, and 200 points for initial condition, respectively. We vary the number of random features, i.e $N=100,200,400$. For each point in the figure, we take the average over 10 repetitions of the experiment.