1 Introduction↩︎

The application of machine learning (ML) techniques to partial differential equations (PDEs) has seen remarkable progress in recent years, presenting novel strategies to address persistent difficulties in scientific computing [1]–[4]. Within this landscape, operator networks have gained attention as an effective and practical approach owing to their capacity to provide rapid solution predictions once training is completed. In contrast to traditional numerical methods that iteratively compute solutions for each new PDE data (e.g., boundary conditions, initial conditions, and external forcing terms), operator networks learn the underlying solution operator from PDE data to the corresponding solutions, thereby enabling fast solution predictions for varying parameters. This advantage positions operator networks as a compelling paradigm for studying parametric PDEs. Representative contributions in this area include the Deep Operator Network (DeepONet) [5] and the Fourier Neural Operator (FNO) [6]. While both DeepONet and FNO enable fast prediction of solutions under varying PDE data, some challenges remain when extending their use to real-world scenarios. A primary obstacle lies in the necessity of a large collection of pre-computed training datasets of solutions. This is typically generated through classical numerical methods for PDEs, which becomes particularly burdensome for nonlinear or highly complex systems. To address this limitation, hybrid frameworks such as the Physics-Informed DeepONet (PIDeepONet) [7] and the Physics-Informed Neural Operator (PINO) [8] have been proposed. These approaches integrate the benefits of Physics-Informed Neural Networks (PINNs) and operator learning by incorporating governing equations directly into the loss functions of neural operators. Nevertheless, such models still encounter difficulties, including reduced accuracy for intricate geometries, challenges in handling stiff problems, and considerable generalization errors caused by limited training data [9]–[11]. In addition, the use of neural networks to represent solution spaces complicates the enforcement of diverse boundary conditions, which in turn affects the reliability of the resulting solutions [12].

To overcome these challenges, an unsupervised operator network founded on classical finite element methods (FEMs), referred to as the Finite Element Operator Network (FEONet), has been proposed in [13]. Within the FEM framework, the numerical solution \(u_h(x)\) with a mesh size \(h>0\) is expressed as a linear combination of nodal basis functions \(\phi_k(x)\), which are piecewise polynomials defined over meshes. Formally, this can be written as \[\label{FEM95sol} u_h(x) = \sum_{k=1}^{N_h} \alpha_k \phi_k(x), \quad x \in \Omega,\tag{1}\] where \(N_h\) is the number of basis functions depending on \(h>0\). Extending this formulation, instead of solving a linear algebraic system to determine \(\{\alpha_k\}^{N_h}_{k=1}\) in 1 , FEONet predicts the coefficients \(\{\widehat{\alpha}_k\}_{k=1}^{N_h}\) using neural networks to construct the approximate solution for the given PDE. In this setting, compared to DeepONet [5], FEONet relies only on a branch network without a trunk network, which significantly reduces the number of parameters. The loss function of FEONet, motivated by classical FEM, is designed using the residual of the Galerkin approximation, which guarantees both accurate PDE solutions and strict enforcement of boundary conditions. Due to the inherent capability of FEM in incorporating boundary conditions, the solutions generated by FEONet precisely satisfy these constraints.

A distinguishing feature of FEONet is its ability to solve parametric PDEs without requiring paired input-output training data, marking a significant advance in computational efficiency and general applicability. To be more specific, FEONet predicts PDE solutions under varying inputs such as initial and boundary conditions, external forcings, and variable coefficients. As shown in [13], it demonstrates flexibility in managing diverse PDE instances across complex domains, while avoiding reliance on pre-generated data. One further advantage of FEONet lies in its versatility, namely, its ability to directly incorporate the techniques from classical numerical analysis. For instance, consider a singularly perturbed problem whose solution develops boundary or interior layers [14], [15]. Such problems are classical yet notoriously difficult in numerical analysis, and many specialized numerical methods have been proposed to treat them [16], [17]. One such approach is the enriched finite element method, which employs corrector basis functions derived from asymptotic analysis to capture the stiff behavior of the solution and yield improved FEM approximations [18]–[20]. This idea extends naturally to the FEONet framework: by incorporating a boundary-layer element into the finite element space via an appropriate corrector basis function, one can construct an enriched FEONet basis that effectively captures sharp variations. Further details and a rigorous convergence analysis of FEONet are given in [13], [21].

While FEONet demonstrates flexibility in delivering accurate and efficient solutions across a wide range of scenarios, there still remain some computational challenges. One of the main problems is the issue of computational cost. As previously noted, within the FEONet framework, the dimension of the neural network output coincides with the number of elements employed in the solution representation. For relatively simple problems, we can obtain accurate solutions via FEM with a modest number of bases; however, for complex problems, more refined computations often necessitate a substantially larger number of elements. In such cases, we observed that FEONet’s performance deteriorates as the number of bases grows. A primary cause of this limitation lies in the approximation and generalization issues that arise when the finite element coefficients are represented through neural network approximation. As proved in [22], as the mesh size \(h\) decreases (equivalently, as the number of elements increases), the error of FEONet initially decreases, but beyond a certain threshold it begins to increase again. This phenomenon directly corresponds to the key discussion in [22], where the main cause of this phenomenon was identified. A second major limitation arises from the scalability issue, which is commonly encountered in large-scale deep-learning tasks. While the errors are mitigated thanks to the theory developed in [22], the computational cost remains substantial for domains with a large number of elements, leading to additional challenges in terms of computation and memory. Problems requiring a large number of elements inevitably suffer from considerable computational overhead, manifested in substantially increased computational cost and impeded training efficiency. This constitutes a crucial challenge from a practical perspective, and it must be addressed if FEONet is to be deployed more broadly in real-world scenarios in a manner comparable to FEM.

This discussion raises a fundamental question about how to handle computational cost efficiently in operator learning methods, which typically require a large amount of computation. One important idea is to impose a suitable sparsity structure on the neural networks used. If one can design an appropriate sparsity pattern by taking into account the architecture, loss function, and training procedure of a given operator network, then it should be possible to achieve a significant improvement in computational efficiency while maintaining accuracy. We refer to such approaches collectively as Sparse Neural Operator netWorks (SNOW). For widely used architectures such as DeepONet and the FNO, developing SNOW-type methods to effectively control computational cost is, in our view, one of the important future directions in operator learning.

In this perspective, the main contributions of this paper are as follows. First, we propose a SNOW approach for FEONet that substantially improves computational efficiency, even when a large number of elements is required. As we show later, our method is inspired by the observation in FEM that only neighboring elements exert strong interactions, whereas the influence between distant elements is comparatively minor. Building upon this motivation, we introduce a new strategy using sparse neural networks that enables effective FEONet computations with a significantly reduced number of parameters. Next, we theoretically establish that the proposed sparse architecture has sufficient approximation capacity for the target operator and admits stable training. In particular, we prove the universal approximation property of the proposed sparse network and provide the stability analysis that guarantees the robust training and solution prediction. Lastly, we validate the efficiency of the proposed method through a series of numerical experiments. Across a variety of benchmark scenarios, we demonstrate that the proposed approach achieves a substantial reduction in the number of trainable parameters and memory usage while maintaining high accuracy. We also show that the proposed sparse structure outperforms randomly constructed sparse networks in terms of accuracy. In addition, in regimes where very fine meshes are required due to the sharp-transition or high-frequency nature of the solution, the original FEONet becomes practically untrainable, whereas the proposed sparse architecture provides a fast and accurate solution prediction.

The remainder of the paper is structured as follows. Section 2 reviews the preliminaries required for the development of our approach. Section 3 presents the proposed methodology in detail. Section 4 provides a theoretical analysis supporting the validity of the method, while Section 5 demonstrates the efficiency of the proposed method through extensive experiments. Finally, Section 6 offers concluding remarks and discusses future research directions.

2 Preliminaries↩︎

The objective of this section is to introduce FEONet, which forms the baseline of our proposed sparse methodology. Since FEONet is built upon the classical FEM, we provide a brief overview of the setting of FEM, and subsequently provide a detailed description of FEONet. As a model problem, we shall consider the general second-order linear elliptic PDE of the form \[\label{main95eq} \begin{align} -\,{\mathrm{div}}\,(\boldsymbol{a}(x)\nabla u)+\boldsymbol{b}(x)\cdot\nabla u+c(x)u&=f(x)\quad{\mathrm{in}}\,\,D, \\ u(x)&=0\quad\hskip+14pt {\mathrm{on}}\,\,\partial D. \end{align}\tag{2}\] Here, for simplicity, let us assume the following: \[\label{coef95cond951} \boldsymbol{a}\in L^{\infty}(D)^{d\times d},\quad \boldsymbol{b}\in W^{1,\infty}(D)^d,\quad c\in L^{\infty}(D),\quad f\in H^{-1}(D).\tag{3}\] For the diffusion tensor \(\boldsymbol{a}=(a_{ij})\), we further assume uniform ellipticity; that is, there exists a constant \(\tilde{a}>0\) such that \[\label{coef95cond952} \sum_{i,j=1}^d a_{ij}(x)\xi_i\xi_j \;\geq\; \tilde{a}\sum_{i=1}^d \xi_i^2,\quad \forall\, \xi=(\xi_1,\ldots,\xi_d)\in\mathbb{R}^d,\; x\in\overline{D}.\tag{4}\] To ensure the well-posedness of the problem 2 , we additionally assume that \[\label{coef95cond953} c(x)-\tfrac{1}{2}\,\mathrm{div}\,\boldsymbol{b}(x)\;\geq\;0,\quad x\in\overline{D}.\tag{5}\] The weak formulation of the problem is then given as follows: we seek \(u\in H^1_0(D)\) such that \[B[u,v] := \int_{D}\boldsymbol{a}(x)\nabla u\cdot\nabla v\,\mathrm{d}x + \int_{D}\boldsymbol{b}(x)\cdot\nabla u\, v\,\mathrm{d}x + \int_{D}c(x)uv\,\mathrm{d}x= \int_D f(x)v\,\mathrm{d}x =:\ell(v),\] for arbitrary \(v\in H^1_0(D)\). Under the assumptions 3 5 , there exist constants \(c_0, c_1, c_2>0\) such that \[\label{lax95cond} B[v,v]\geq c_0\|v\|^2_{H^1(D)}, \,\, |B[u,v]|\leq c_1\|u\|_{H^1(D)}\|v\|_{H^1(D)}, \,\, |\ell(v)|\leq c_2\|v\|_{H^1(D)}.\tag{6}\] Therefore, the existence and uniqueness of the weak solution follow directly from the classical Lax–Milgram theorem (see, e.g., [23]). For brevity, the above discussion has been restricted to the linear case. However, as will be explained in detail later, FEONet employs the residual of the variational formulation as the loss functional, and therefore, it can be directly applied to nonlinear problems as well. In fact, in the numerical experiments section, we will also solve some nonlinear equations using the proposed method.

2.1 Finite element methods↩︎

We first provide a brief overview of the classical FEM. As a preliminary step, we introduce the finite element space that will serve as the foundation throughout the paper. Let \(\mathcal{G}_h\) denote a shape-regular partition of the physical domain \(\overline{D}\), where \(h_E\) represents the diameter of an element \(E\in\mathcal{G}_h\), and \(h=\max_{E\in\mathcal{G}_h}h_E\). We assume that there exists a constant \(\gamma>0\), independent of \(h\), such that \[\max_{E\in\mathcal{G}_h}\frac{h_E}{\rho_E}\leq \gamma,\] where \(\rho_E\) is the supremum of the diameters of inscribed balls in \(E\in\mathcal{G}_h\). For a given partition \(\mathcal{G}_h\), the corresponding finite element space is defined as \[\mathbb{V}_h=\mathbb{V}(\mathcal{G}_h)\;\mathrel{\mathop{:}}=\;\bigl\{V\in C(\overline{D}):\,V|_E\in\hat{\mathbb{P}}_{\mathbb{V}},\; E\in \mathcal{G}_h,\; V|_{\partial D}=0\bigr\},\] where \(\hat{\mathbb{P}}_{\mathbb{V}}\subset W^{1,\infty}(\hat{E})\) is a finite-dimensional subspace. We further assume that \(\mathbb{V}_h\) admits a finite, locally supported basis. Specifically, for each \(h>0\), there exists \(N_h\in\mathbb{N}\) such that \(\mathbb{V}_h=\mathrm{span}\{\phi_1,\ldots,\phi_{N_h}\}\). For a basis function \(\phi_i\), there exists a patch \(\mathcal{P}_i\) such that \(\textrm{supp}(\phi_i) \subset \mathcal{P}_i\) and \(\phi_i(x) = 0\) for \(x\not\in \mathcal{P}_i\). If \(\mathcal{P}_i\) is the set of elements that contain the node \(i\), then \(\mathcal{P}_i = \textrm{supp}(\phi_i)\).

We now define the Galerkin approximation: we seek a discrete solution of the form \[\label{basis95exp} u_h=\sum_{k=1}^{N_h}\alpha_k\phi_k\in \mathbb{V}_h\tag{7}\] satisfying \[\label{var95for} B[u_h,v_h]=\ell(v_h)\quad \text{for all } v_h\in\mathbb{V}_h.\tag{8}\] Introducing the matrices \(A=(A_{ij})_{1\leq i,j\leq N_h}\in\mathbb{R}^{N_h\times N_h}\) and the vector \(F=(F_j)_{1\leq j\leq N_h}\in\mathbb{R}^{N_h}\) defined by \[\label{vec95form} A_{ij}=B[\phi_i,\phi_j]\quad{\rm{and}} \quad F_j=\ell(\phi_j),\tag{9}\] the discrete scheme 8 can be reformulated as the linear algebraic system \[\label{LAS} A\alpha=F,\quad \alpha=(\alpha_k)_{1\leq k\leq N_h}\in\mathbb{R}^{N_h}.\tag{10}\] It follows that \(B[\phi_i, \phi_j] = 0\) if \(\mathcal{P}_i \cap \mathcal{P}_j\) is of measure zero. As a result, the linear system 10 has a sparse structure. By solving this system, we determine the coefficient \(\alpha\), which in turn allows us to compute an approximate solution of the given PDE 2 via the basis expansion 7 .

2.2 Finite element operator networks↩︎

We now introduce the finite element operator network (FEONet), originally proposed in [13]. As discussed earlier, the input to FEONet can represent various types of PDE data, such as external forces, variable coefficients, or boundary conditions. For clarity, we present a prototype setting in which the input corresponds to an external forcing term, though the framework can be naturally extended to other types of inputs. Given a forcing function \(f\), rather than computing the coefficients \(\alpha\) through the linear system 10 , FEONet predicts them via deep neural networks. Specifically, the neural network takes as input the forcing term \(f\), which is parameterized by a random variable \(\omega\) defined on the probability space \((\Omega,\mathcal{T},\mathbb{P}_{\Omega})\). Typical examples include Gaussian random fields or random forcing functions of the form \[f(x,\omega)=\omega_1\sin(2\pi\omega_2 x)+\omega_3\cos(2\pi\omega_4 x),\] with \(\omega=(\omega_1,\omega_2,\omega_3,\omega_4)\) drawn i.i.d. uniformly with \(\omega_i\in[a_i,b_i]\) for \(i=1,2,3,4\). Once the input feature \(\omega\in\Omega\) passes through the neural network, it outputs coefficients \(\{\widehat{\alpha}_i\}_{i=1}^{N_h}\), and the solution prediction is reconstructed as \[\label{sol95recon} \widehat{u}_h(x,\omega)=\sum_{i=1}^{N_h}\widehat{\alpha}_i(\omega)\phi_i(x).\tag{11}\]

For training, the population loss is defined by the residual of the variational formulation 8 : \[\label{pop95loss} \mathcal{L}(\alpha)=\mathbb{E}_{\omega\sim\mathbb{P}_{\Omega}} \left[\sum_{i=1}^{N_h}\big|B[\widehat{u}_h(x,\omega),\phi_i(x)]-\ell(\phi_i(x))\big|^2\right]^{1/2}.\tag{12}\] In practice, the empirical loss function is employed, obtained via Monte Carlo sampling of 12 : \[\label{emp95loss} \mathcal{L}^M(\alpha)=\frac{|\Omega|}{M}\sum_{j=1}^M \left[\sum_{i=1}^{N_h}\big|B[\widehat{u}_h(x,\omega_j),\phi_i(x)]-\ell(\phi_i(x))\big|^2\right]^{1/2},\tag{13}\] where \(\{\omega_j\}_{j=1}^M\) are i.i.d. samples drawn from \(\mathbb{P}_{\Omega}\). At each training epoch, the network parameters are updated to minimize \(\mathcal{L}^M\), until the empirical loss becomes sufficiently small. The final prediction is then obtained via 11 , after sufficient training has been performed.

A crucial feature of FEONet is that training relies solely on random samples from \(\Omega\), without requiring any precomputed input-output pairs. Moreover, since the method is based on the basis expansion 11 , exact boundary conditions can be imposed in the same way as in classical FEM. A schematic overview of the FEONet structure is given in Figure 1. Extensive numerical experiments on benchmark problems are reported in [13], confirming the effectiveness of the approach in terms of accuracy, generalization capability, and computational efficiency, and the rigorous convergence analysis of FEONet was conducted in [22]. As noted in these papers, a major limitation of FEONet is that its computational cost grows rapidly as the number of finite elements increases, which makes large-scale applications challenging. In some cases, training may not progress at all. As discussed before, this observation provides a key motivation for our work. To address the computational cost, training efficiency, and memory bottlenecks, we propose in this paper a new sparse architecture tailored to FEONet. As a simple motivating example, Figure 2 summarizes learning curves for different values of \(N_h\), the number of degrees of freedom in the finite element discretization. When \(N_h\) is small, the original FEONet and our proposed method exhibit comparable behavior. However, for relatively large \(N_h\), training of the original FEONet often fails to make sufficient progress and is terminated early (e.g., via early stopping), whereas our sparse architecture consistently yields robust training performance across all tested values of \(N_h\). As will be further demonstrated in the experiments, many practical settings require fine meshes with a large number of elements; therefore, our method is expected to advance FEONet toward scalable, real-world applicability.

3 Methodology: Sparse FEONet↩︎

As in [13], FEONet can take fully-connected (FC) neural networks or convolutional neural networks (CNN) as model architectures, but in this paper, we shall consider the FC case. In an FC layer with \(N\) input nodes and \(M\) output nodes, the number of weights and biases equals \((N+1)M\). In this case, the dimension of the parameter space grows rapidly as the number of nodes and layers increases. Thus, from an optimization viewpoint, it is crucial to reduce the number of parameters. This motivates us to introduce a new approach that assigns sparse weights to each layer. In this approach, we represent each layer’s weight matrix in a sparse-matrix format.

We now describe how to design a sparse network for the FEONet. To clearly convey the idea, we shall consider the unit-square domain, and the case of homogeneous Dirichlet boundary condition as in 2 , which allows us to consider only the interior nodes. In Figure 3, we need to determine the coefficients \(\{\alpha_k\}\) in 7 for the 25 interior nodes. Hence, the input layer has 25 nodes, which matches the number of unknowns. For simplicity, we also set 25 nodes in each hidden layer and in the output layer. The parameter \(C_\ell\) determines the level of support expansion for the basis function associated with each node. The key idea is that, in the FEM, interactions occur only between neighboring elements. In particular, \(B[\phi_i, \phi_j] = 0\) whenever the support of \(\phi_i\) and \(\phi_j\) do not overlap. Accordingly, we design the neural network to be sparse by introducing connectivity only between nodes corresponding to elements that influence each other (i.e., adjacent elements). More precisely, for \(C_\ell = 1\), we connect a node only to the nodes contained in the support of its basis function. For \(C_\ell = 2\), new nodes are further connected to all nodes contained in the union of the supports of the nodes selected at \(C_\ell = 1\). Thus, the number of connections increases as \(C_\ell\) increases. Figure 3 shows the FC layer and the sparse layer when \(C_\ell = 1\). For example, the nodes contained in the support of the basis function at node 7 are nodes 2, 3, 6, 7, 8, 11, and 12. We therefore connect input node 7 to those output nodes. To quantify the degree of sparsity more precisely, we introduce the notion of sparsity measure. With this sparse connectivity above, the number of weights is 137 with the sparsity 0.7808, where \[S = 1 - \frac{\textrm{\# of (nonzero) weights for a sparse layer }}{\textrm{\# of weights for a FC layer}}.\] Table 1 shows the number of weights and the sparsity \(S\) for the FC layer and for sparse layers at several values of \(C_\ell\). For fixed \(N_h\), the sparsity decreases as \(C_\ell\) increases; when \(C_\ell\) is large enough, the sparse layer coincides with the FC layer, as in the case \(N_h = 25\). As \(N_h\) increases, the sparsity increases, and for sufficiently large \(N_h\), the sparsity becomes less sensitive to \(C_\ell\), which makes our proposed sparsity more practical.

To provide a more formal illustration of the proposed method, we consider the following representative example. For simplicity, let \(\Omega = (0,1)^2 \subset \mathbb{R}^2\) be the unit square. Note that the same construction extends to general dimensions and more general domains. We consider a uniform Cartesian grid with \(h = 1/n\) in both \(x\) and \(y\) axis directions, with integer indices \(i,\;j \in \{0, 1, \dots, n\}.\) On each grid square, we take the standard right isosceles triangle split (e.g., cut along a diagonal), yielding a triangulation of \(\Omega\) for a conforming piecewise linear finite element method. Since we consider homogeneous Dirichlet boundary conditions, only the values corresponding to interior nodes \((i,j)\) with \(1\le i, \;j \le n-1\) are unknowns. Here, we assign a single global node number by the row-major mapping \(k = (n-1)(j-1) + i\) (see, e.g., Figure 3). Let \(\eta_{i,j} = (ih, jh)\) denote the coordinate of node \((i, j)\) and \(\mathcal{V}_h = \{(i,j): 1\le i, \;j\le n-1\}\) be the set of interior nodes with the size \(N_h = |\mathcal{V}_h| = (n-1)^2\). Also, we denote the corresponding nodal basis centered at \(\eta_{i,j}\) by \(\phi_{i,j}\). For the global indices \(k\) and \(\ell\), we may write \(\phi_k = \phi_{p,q}\) and \(\phi_\ell = \phi_{r,s}\) with \(p,q,r,s \in \{1,\dots,n-1\}\). Under this notation, 9 can be rewritten as \(A_{k\ell} = A_{(p,q),(r,s)}\) and \(F_k = F_{(p,q)}\).

In our setting, every layer has width \(N_h\) and the neurons in every layer are indexed by the same interior-node indices \((i,j) \in \mathcal{V}_h\). We let \(z^{(0)} = F\in\mathbb{R}^{N_h}\) be an input of a neural network and, for \(\ell = 0,1,\dots,L-1\), we write \[\label{val95layer} z^{(\ell+1)}_{i,j} = \sigma\left(\sum_{(p,q)\in \mathcal{V}_h} W^{(\ell+1)}_{(i,j),(p,q)} z^{(\ell)}_{p,q}+b^{(\ell+1)}_{{(i,j)}} \right),\tag{14}\] where \(\sigma\) is an activation function, \(W^{(\ell)}\) is the weight matrix, and \(b^{(\ell)}\) is the bias vector for the \(\ell\)-th layer. The final output is \(z^{(L)} \in \mathbb{R}^{N_h}\), which is the prediction for the finite element coefficients. In 14 , the weights are not fully-connected. Let \(\mathcal{P}^{(1)}_{i,j} = \textrm{supp}(\phi_{i,j})\). For \(t \ge 1\), we define the level-\((t+1)\) patch recursively by \[\label{level95patch} \mathcal{P}^{(t+1)}_{i,j} := \bigcup_{(p,q) \in \mathcal{V}_h}\{\text{supp}(\phi_{p,q}) \;: \; |\textrm{supp}(\phi_{p,q}) \cap \mathcal{P}^{(t)}_{i,j}| \neq 0\}.\tag{15}\] We then define the level-\(t\) neighborhood by \[\mathcal{V}_{t}(i,j) = \{(p,q) \in \mathcal{V}_h : |\textrm{supp}(\phi_{p,q}) \cap \mathcal{P}^{(t)}_{i,j}| \neq 0 \}.\] Accordingly, for interior nodes \((i,j), (p,q) \in \mathcal{V}_h\), we allow a nonzero weight as follows: \[\label{sparse95weight} W^{(k)}_{(i,j),(p,q)} \neq 0 \quad \text{is allowed only if}\;(p,q) \in \mathcal{V}_{C_\ell}(i,j),\tag{16}\] for some given connectivity level \(C_\ell \in \mathbb{N}\), and we set the weight to zero otherwise to enforce sparsity. For a fixed \(C_\ell\) and a stack of \(L\) such layers, the dependence of \(z^{(L)}_{i,j}\) on the input \(z^{(0)}\) is supported on \(\mathcal{V}_{L C_\ell}(i,j)\). Hence, the effective receptive field grows linearly with the depth \(L\) (see, e.g., Figure 4). Details on the effect of connectivity are provided in the appendix. Remarkably, accurate approximations are achieved even when the effective receptive field does not fully cover the entire domain.

4 Theoretical Analysis↩︎

In this section, we provide the theoretical background for the proposed sparse FEONet approach. As discussed earlier, our method significantly reduces the number of parameters compared to the original FEONet. While this is clearly advantageous from a computational point of view, it raises a natural theoretical question on the approximation capability of the resulting sparse architecture. In addition, we need to compare our finite-element-guided sparsification strategy with more naive approaches that impose sparsity in a random or purely heuristic manner. In particular, we discuss in what sense the finite-element-guided design leads to theoretical advantages, and how these advantages are reflected in practice. Moreover, we investigate, from a theoretical perspective, whether the proposed sparsity pattern not only improves computational efficiency but also leads to highly stable and efficient training. The goal of this section is to address these theoretical questions for sparse FEONet and to establish rigorous guarantees for its approximation and stability properties. We conclude the section by presenting brief numerical experiments that support the theoretical findings.

4.1 Universal approximation theorem↩︎

This section is devoted to providing the theoretical justification for the proposed sparse neural networks. In particular, we address the question of whether the proposed sparse neural network constitutes a suitable approximation class. While the universal approximation property of dense neural networks is well established (see, e.g., [24]–[26]), our approach relies on a sparse architecture in which connectivity is restricted to nodes associated with neighboring elements. Thus, it is essential to demonstrate that neural networks with such a structure can still approximate the target function effectively. Accordingly, in this section, we shall first establish a universal approximation theorem for the proposed sparse network.

Note that for our proposed method, the number of nodes in each layer of our neural network coincides with the number of degrees of freedom \(N_h\in\mathbb{N}\). However, in such a setting, it is well known in the literature that even dense neural networks cannot, in general, be guaranteed to exhibit the universal approximation property. In fact, according to the results known to date (e.g., [27]–[29]), a class of neural networks attains the universal approximation property only when the number of nodes in its hidden layer is taken to be sufficiently large. Consequently, one has even less justification for expecting any universal approximation capability from a sparse network defined in the previous section. Note, however, that the mapping we seek to approximate by means of a neural network does not belong to an arbitrary class of continuous functions. Rather, it is the mapping that assigns, to the parameters defining a given PDE problem, the corresponding coefficients of the associated finite element discretization. As described in 10 , the finite element coefficient can be characterized by a linear algebraic system. Thus, in our framework, the essential task is to represent the linear map \(x\mapsto A^{-1}x\) where \(A\) denotes the finite element matrix, by a sparse neural network. As will be verified in the subsequent proof, this structural feature aligns with our setting in a particularly precise and favorable manner, and the desired conclusion indeed follows.

To simplify the theoretical setting of the network connectivity, we first describe our proposed method within a graph-theoretic framework. Let us begin with the following definitions.

For a given simple undirected graph \(G\), we define a sparse matrix associated with the connectivity of \(G\). Here, \(\mathrm{GL}_N(\mathbb{R})\) denotes the group of invertible matrices over \(\mathbb{R}\) with respect to the matrix multiplication.

The first key result in our analysis is the following.

To show this equivalence, let us first prove the implication in one direction, which is encapsulated in the following proposition.

In order to prove the converse, we first recall some standard notation for elementary matrices.

In classical linear algebra, the following fact is well-known, which describes a generating set for \(\mathrm{GL}_N(\mathbb{R})\) in terms of transvections and diagonal matrices. More precisely, the general linear group \(\mathrm{GL}_N(\mathbb{R})\) can be generated by all transvection matrices \(E_{ij}(t)\) with \(i\neq j\), \(t\in \mathbb{R}\) and all invertible diagonal matrices \({\rm{diag}}(\lambda_1,\dots,\lambda_N)\). Therefore, to complete the proof of Theorem 1, it remains to show that \(H_G(\mathbb{R})\) contains all transvection matrices and invertible diagonal matrices provided that \(G\) is connected. It is obvious that any invertible diagonal matrices \({\rm{diag}}(\lambda_1,\dots,\lambda_n)\) are \(G\)-sparse, and hence contained in \(H_G(\mathbb{R})\). Moreover, if \(\{i,j\}\in E\), it is easy to see directly from the definition that \(E_{ij}(t)\in H_G(\mathbb{R})\) for all \(t\in\mathbb{R}\). Therefore, what remains to prove is that \(E_{ij}(t)\in H_G(\mathbb{R})\) for all \(t\in\mathbb{R}\) even if \(\{i,j\}\notin E\). To do this, let us introduce a commutator identity: for distinct indices \(i, j, k\), and \(s,t\in\mathbb{R}\), there holds \[\label{comm95id} E_{ij}(st)=E_{ik}(s)\, E_{kj}(t)\, E_{ik}(-s)\, E_{kj}(-t),\tag{17}\] which follows immediately by directly computing the matrix multiplications using the identity \(e_{pq} e_{rs} = \delta_{qr} e_{ps}\). We now use the connectivity of \(G\) to show that \(H_G(\mathbb{R})\) contains every transvection matrices \(E_{ij}(t)\) even if \(\{i,j\}\notin E\).

Let us now apply the above theory to our sparse neural network setting to prove the universal approximation theorem. For simplicity, we present the argument in the two-dimensional setting described above. It is, however, straightforward to verify that our approach extends to arbitrary spatial dimensions. Note that the weight matrix \(W^{(\ell)}\) in the sparse layer 14 can be reformulated using the above notation. More precisely, we shall define a simple undirected graph \(G_h = (\mathcal{V}_h,E_h)\). For the set of edges \(E_h\), we set the vertices \(u=(i,j)\), \(v=(p,q)\in\mathcal{V}_h\) to be adjacent if \(v=(p,q)\in\mathcal{V}_{C_\ell}(i,j)\) for a prescribed constant \(C_\ell>0\). The important first step is the following lemma. It follows immediately because, for fixed \(C_\ell\) and mesh size \(h > 0\), the entire domain \((0,1)^2\) can be covered by finitely many graph balls \[\mathcal{B}_G(u; C_\ell) := \{v \in \mathcal{V}_h \;: \;d_G(u, v) \le C_\ell\}.\] By the patch construction 15 , \(\mathcal{B}_G(u; C_\ell)\) coincides with \(\mathcal{V}_{C_\ell}(u)\), which in turn corresponds to \(\mathcal{P}^{(C_\ell)}_u\).

An important observation is that our sparse weight matrix \(W^{(\ell)}\) in 14 can be characterized by the graph \(G_h\). To be more specific, if we construct the network sparsely according to the procedure proposed in this Section 3, then our weight matrices \(W^{(\ell)}\) are \(G_h\)-sparse for \(\ell=1,2,\ldots, L\). Now, we are ready to prove the universal approximation property of our sparse network. Since the graph \(G_h\) corresponding to the sparse layer defined above is connected and the finite element matrix \(A\) is invertible, by Theorem 1, we see that \(A^{-1}\) can be represented as a product of \(G_h\)-sparse matrices. Based on this fact, we shall prove that our target function \(x\mapsto A^{-1}x\) can be represented as a ReLU network with sparse layers. In other words, each \(G_h\)-sparse matrix can be realized as a single affine layer of a neural network with the prescribed sparsity pattern, which allows us to translate the above matrix factorization into a neural network representation.

The idea of the proof is to implement each factor of \(M\) by a small ReLU subnetwork that acts as the identity on the relevant compact subset of intermediate representations. The key observation is that ReLU coincides with the identity on the positive half-line.

Theorem 5 means that any given invertible linear mapping can be represented exactly by a ReLU sparse network. A natural subsequent question is whether an analogous property holds for more general activation functions. In this case, as in the classical universal approximation theorem, we can approximate any given invertible linear mapping to arbitrary accuracy, which is encapsulated in the following theorem.

As we can see from Theorem 5 and Theorem 6, under a given mesh size \(h\) and a fixed connectivity constant \(C_{\ell}\), the connectivity of the underlying graph plays a crucial role in guaranteeing the above results. In particular, when the graph is connected so that information from the input can propagate to the final layer, we are able to establish a universal approximation property. From this viewpoint, it is natural to wonder whether one might lose such desirable guarantees if, instead of designing the network as proposed here, one imposes sparsity in an ad hoc manner to construct a sparse neural network. To investigate this, we performed some numerical experiments. For \(N_h=225\) and \(N_h=961\), we constructed two sparse networks, one using our proposed method and the other using a randomly imposed sparsity pattern, and then carried out the FEONet experiments with both architectures. Figure 5 illustrates the structure of the weight matrices used in each case. The blue entries indicate nonzero values, whereas the remaining white entries correspond to zeros, yielding sparse matrices. On a uniform mesh, the proposed sparse weight matrices have a banded structure, which is shown in (b) and (d) of Figure 5. For a fair comparison, we construct a random sparse weight matrix with the same number of nonzero entries as the proposed one, and we exclude any case with an all-zero row or an all-zero column (see (a) and (c) of Figure 5).

As we can see from Table 2, the sparsity pattern introduced by the supports of the FEM basis functions yields substantially higher efficiency than a random sparse pattern. For each \(N_h\), we generated 10 different random sparse weight matrices and built ten corresponding networks. These were trained for 10,000 epochs under exactly the same settings as the sparse FEONet. The relative \(\ell_2\) Errors are calculated with the finite element solution at the same \(N_h\) as the reference solution. This confirms that the connectivity condition discussed above plays a central role, both in the theory and in practical numerical performance.

4.2 Stability↩︎

In this section, in order to further highlight the efficiency of the proposed method, we present a theoretical study of the network’s stability in training and inference. The analysis relies on layer-wise operator norms and on an activation with Lipschitz continuity, and it yields stability bounds with explicit dependence on depth and resolution. We compare dense (FC) and sparse connectivity and find that the sparse network preserves locality and yields stronger stability guarantees. We shall also present some brief numerical tests to support the theoretical result.

By the Marchenko–Pastur law [30], if each entry of \(W\in\mathbb{R}^{N\times N}\) is independently identically distributed random variables with mean 0 and variance \(\sigma^2 < \infty\), i.e., \(W_{ij} \sim \mathcal{N}(0,\sigma^2)\) for all \(i,j\in\{1,2,\cdots,N\}\), it is known, with high probability, that \[\label{Marchenko95estimate} \|W\|_2 \approx 2\sigma \sqrt{N}.\tag{22}\] On the other hand, the key observation is that for a sparse matrix \(W \in\mathbb{R}^{N\times N}\) constructed via 16 , we have \[\label{sparse95bound} \|W\|_2 \le (\|W\|_1 \|W\|_\infty)^{1/2} \le \omega \gamma = \mathcal{O}(1),\tag{23}\] where \(\gamma =\max_{i,j} |\mathcal{V}_{C_\ell}(i,j)|\) denotes the maximum number of nonzero entries per row. Note that \(\gamma\) depends on the constant \(C_\ell\) and the mesh structure, but it is independent of \(N\). With this in mind, we compare the stability of the dense (FC) network with that of our proposed sparse network, which is encapsulated in the following theorem.

In order to demonstrate Theorem 8, we test depth \(L = 6\) under two connectivities: FC layer and sparse layer that satisfies 16 . We consider both untrained i.i.d. Gaussian weights and trained weights, under identical architectures and training conditions. We generated 3,000 input samples and added random noise with a magnitude with \(1\%\) of the maximum norm of the input dataset. We then compared the model outputs for the original and perturbed inputs. We define the empirical sensitivity by \[\widehat{C}_S(f,\delta) = \frac{\|\mathcal{N}_L(f) - \mathcal{N}_L(\hat{f})\|_2}{\|\delta(f)\|_2}.\] For each resolution \(N_h\), we report the mean (and standard deviation) of \(\widehat{C}_S\) over the 3,000 samples and compare it with the theoretical upper bound \(C_S\).

By the Marchenko–Pastur estimate 22 , each \(C_{W^{(\ell)}}\) for the FC network doubles when \(N_h\) is quadrupled. So the right-hand side in 24 increases by a factor \(2^{L} = 2^6\). Figure 6 and Table 3 confirm that the upper bound \(C_S\) grows with slope \(\approx 6\) on a log scale, whereas the empirical sensitivity remains strictly below the upper bound for all \(N_h\). We note that the FC network could not be instantiated due to memory limitations at \(N_h = 255^2\). For the sparse network, \(C_\ell = \mathcal{O}(1)\) is independent of \(N_h\).

5 Numerical Experiments↩︎

In this section, we present the experimental results comparing our proposed sparse network with the FC neural network within the FEONet framework. We tested various PDEs with Dirichlet boundary conditions, from coarse to fine mesh resolutions.

For each experiment, we utilized 3,000 training samples and 3,000 test samples from randomly generated external forcing terms. As the mesh becomes finer, the number of trainable parameters increases significantly, which leads to memory bottlenecks and computational infeasibility in FC architectures. However, our model overcame these problems, improving efficiency and achieving higher accuracy even in a high-resolution regime. To evaluate the robustness of our sparse model across diverse geometries, we tested the model not only on structured triangular meshes in a square domain, but also on irregular triangular meshes in both square and circular-hole domains generated via FEniCS [31]. These tests demonstrate that the proposed model shows strong performance across diverse domains and mesh geometries. Both the FC model and the proposed sparse model were trained under identical settings on the same dataset for a fair comparison. We increased the connectivity step by step to find out the minimum number of connections needed to achieve accuracy comparable to the fully connected baseline model. Each model employed five hidden layers and was optimized using the Adam optimizer. We used the Swish activation function and applied a cosine decay scheduler to gradually reduce the learning rate for better training convergence. The training was performed in an unsupervised manner by minimizing the weak-form residual.

For the performance evaluation of the proposed method, we first measure the relative \(L^2\) error and the relative \(H^1\) semi-norm error against the FEM solution \(u_h^*\) on a fine mesh (\(n=1,024\)) which can be regarded as a true solution. To be more specific, we measure the following errors: \[\textrm{Rel. L^2 Err.} := \frac{\|u_{\textrm{pred}} - u_h^*\|_{L^2(\Omega)}}{\|u_h^*\|_{L^2(\Omega)}}, \quad \textrm{Rel. H^1 Semi Err.} := \frac{\|\nabla u_{\textrm{pred}} - \nabla u_h^*\|_{L^2(\Omega)}}{\|\nabla u_h^*\|_{L^2(\Omega)}}.\] In addition, to measure the error introduced by incorporating neural-network approximation into the original FEM, we also evaluated the accuracy of the FEM solution computed on the same mesh as the solution predicted by FEONet. These values are shown in parentheses in the tables below. Moreover, we measured the memory usage of trainable parameters to evaluate computational efficiency.

Our model was implemented in JAX(v0.4.7), using Flax and Optax libraries. All experiments were conducted on a workstation with a single NVIDIA RTX 3090 GPU(24GB VRAM), running CUDA 11.4 and CUDNN 8.2.4.

5.1 Advection-diffusion-reaction equation↩︎

For the basic performance evaluation, we first consider the 2D advection-diffusion-reaction equation defined as \[\begin{align} -0.1\;\Delta u(x,y) + a\cdot\nabla u(x,y) + 20\;u(x,y)&= f(x,y) && (x,y) \in \Omega, \\ u(x,y) &= 0 && (x,y) \in \partial \Omega, \end{align}\]

where \(\Omega=[-1, 1]^2\), \(a=(-1,0)^T\). Moreover, we set external forces as inputs of neural networks, which are given by \[\label{f95input} f(x,y)=m_0\sin(n_0x+n_1y)+m_1\cos(n_2x+n_3y),\tag{25}\] where \(m_0, \, m_1\) and \(n_0, \, n_1, \, n_2, \, n_3\) are random samples from \([0,1)\) and \([0,1)\times\pi\), respectively.

All models consist of five hidden layers. We consider a uniform Cartesian grid with \(h = 2/n\) in both \(x\) and \(y\) axis directions, and take the standard right isosceles triangle split yielding a triangulation of \(\Omega\) for a conforming piecewise linear finite element method. Corresponding to the mesh resolutions of \(n=16,\;32\), and \(64\), the connectivity is \(C_{\ell}=1\), \(C_{\ell}=2\), and \(C_{\ell}=3\) for each respective resolution. We determined these connectivity values empirically by testing increasing connectivity levels and selecting those that provide stable training and accurate predictions. To ensure a fair comparison, each model was trained for 10,000 epochs using the same optimizer and learning rate scheduler.

Table 4 shows the number of parameters, memory, weak form loss, relative \(L^2\) error, and relative \(H^1\) semi-norm error on the test set. The error and loss values are averaged over 3000 data samples per epoch. The sparse neural network performs comparably to, or even better than the dense network with up to about 99% fewer parameters. In the case of the coarse grid (\(n=16\)), the errors of both models are similar to those of the finite element method, which is the value in the parentheses. As the resolution increases, the sparse model shows better convergence. At \(n=64\), the dense model achieves 0.0074 in relative \(L^2\) error, while the sparse model achieves 0.0052, which is close to the accuracy of the finite element solution at that resolution. Figure 7 shows the input function, the finite element solution, the sparse FEONet prediction, and the absolute error between the finite element solution and the prediction in \(n = 64\).

5.2 Helmholtz equation↩︎

Some problems require high-resolution solutions with many basis functions due to stiffness or highly oscillatory behavior. In such regimes, the original FEONet often faces substantial computational difficulties, whereas our proposed method can predict solutions effectively. To illustrate this point, we consider the 2D Helmholtz equation following the formulation in [32], given by \[\begin{align} \Delta u(x,y) + k^2 u(x,y) &=q(x,y) && (x,y) \in \Omega, \\ u(x,y) &= g(x,y) && (x,y) \in \partial \Omega, \end{align}\]

where \(\Omega=[-1, 1]^2\). We construct data using the manufactured solution. More precisely, we prescribe the exact solution \[u(x,y)=\sin(a_1 \pi x)\sin(a_2 \pi y),\] and obtain the corresponding forcing term \(q(x,y)\) and boundary data \(g(x,y)\) given by \[\begin{align} q(x,y) &= - (a_1 \pi)^2 \sin(a_1 \pi x)\sin(a_2 \pi y) - (a_2 \pi)^2 \sin(a_1 \pi x)\sin(a_2 \pi y) + k^2 \sin(a_1 \pi x)\sin(a_2 \pi y), \\ g(x,y) &=u\big|_{\partial \Omega}, \end{align}\] with the random samples \(a_1, \, a_2 \in [1, 10),\) where we fixed \(k = 1\).

To make the problem more challenging, we sampled relatively large values of \(a_1\), and \(a_2\). This produces high oscillations, which require fine-mesh resolutions for accurate solution prediction. This setting allows us to evaluate the model’s capability in this high-resolution regime.

We conducted the experiments using the same setting as in the advection-diffusion-reaction equation, except that each model was trained for longer epochs because of the increased difficulty caused by the nonhomogeneous Dirichlet boundary conditions. We tested the model on mesh resolutions, \(n=16,\;32,\;64\), and \(128\). The sparse model uses the same local connectivity strategy as before with \(C_{\ell}=5\) for \(n=32\), \(C_{\ell}=7\) for \(n=64\), and \(C_\ell=10\) for \(n=128\). Because the solution exhibits high oscillations, the coarse grid cannot capture the overall behavior of the solution. Although finer mesh resolutions can resolve these oscillations, the dense model has significant optimization difficulties and fails to converge in this regime. Also, for \(n=128\), the number of parameters in the dense model results in the GPU memory overflow, which prevented the experiment from being conducted. However, the sparse model converges well on \(n=32, 64\), and even on the high-resolution case \(n=128\). On \(n=128\), our model achieves stable convergence with only about 1.91% of the parameters required by the dense model. Table 5 summarizes the results, and Figure 8 shows the plot of the prediction of the sparse model at \(n=128\).

5.3 Nonlinear Burgers’ equation↩︎

We next show that our method also performs well on nonlinear equations. To do this, we tested our model on the 1D nonlinear Burgers’ equation to assess the model’s capability in handling nonlinear terms. The equation can be written as \[\begin{align} -0.1u_{xx}+uu_x &= f(x) && x \in [-1,1], \\ u(x) &= 0 && x \in \{-1, 1\}, \end{align}\] where \(f(x)=m_0\sin(n_0x)+m_1\cos(n_1x)\) with the random samples \(m_0, \, m_1 \in [0,1)\) and \(n_0, \, n_1\in [0,1)\times\pi\). For the sparse model, we set the connectivity to \(C_\ell=8\) for \(n=64\), \(C_\ell=13\) for \(n=128\), and \(C_\ell=30\) for \(n=256\).

Table 6 presents the quantitative results, and Figure 9 shows the prediction of the sparse model at \(n=256\). The sparse model achieves convergence in the relative \(H^1\) semi-norm error that is comparable to that of the finite element solution across all resolutions. In contrast, the dense model shows a noticeable difference from the FEM as \(n\) increases. This behavior suggests that the optimization of the dense model becomes more difficult at higher resolutions, and that the training may not be sufficient within the fixed number of epochs. For the relative \(L^2\) error, the sparse model still achieves small errors, but a larger gap with respect to the finite element solution is observed. This difference is likely due to the weak-form-residual-based training, which directly controls the gradient of the solution and therefore favors accuracy in the \(H^1\) semi-norm, while the \(L^2\) error is only indirectly controlled. In addition, the nonlinear and advective nature of the Burgers’ equation can make the \(L^2\) error more sensitive to small shifts in the solution profile. Compared to the two-dimensional case, the ratio of parameters between the dense and sparse models is relatively large in one dimension. However, since each node has fewer neighboring degrees of freedom, the overall reduction in the number of parameters becomes less significant. These results indicate that the proposed sparse model remains effective for problems with nonlinear terms.

5.4 Irregular triangular meshes and complex domain↩︎

To evaluate the performance of the proposed sparse model for a general case, we test the model on irregular triangular meshes and examine whether the sparse connectivity strategy remains effective when the mesh structure is not aligned with a regular Cartesian grid. To do this, we solve the 2D Poisson equation with the homogeneous Dirichlet boundary condition and the forcing terms given by 25 . Using FEniCS, we generate irregular triangular meshes on a standard square domain and on a square domain with a circular hole. The two domains have 919 and 679 degrees of freedom, respectively. Except for the structure of the mesh and the domain geometry, all training settings are kept consistent with the previous experiments.

Table 7 shows the experimental results. In both cases, the sparse model converges accurately with fewer parameters compared to the dense model. These results indicate that the local connectivity-based sparse structure remains effective on unstructured meshes and complex geometries. Figure 10 visualizes the prediction of the sparse model on the domain with a circular hole.

5.5 Stokes equations↩︎

Finally, we extend our evaluation to the 2D Stokes equations: \[\begin{align} -\Delta \mathbf{u}(x,y) + \nabla p(x,y) &= \mathbf{f}(x,y) && (x,y) \in \Omega, \\ \nabla \cdot \mathbf{u}(x,y) &= 0 && (x,y) \in \Omega, \\ \mathbf{u}(x,y) &= \mathbf{0} && (x,y) \in \partial \Omega, \end{align}\] where \(\Omega=[-1, 1]^2\), \(\mathbf{u}=(u_1, u_2)\) represents the velocity field, and \(p\) denotes the kinematic pressure. To ensure uniqueness of the pressure, we also impose the zero-mean condition, \(\int_\Omega p\;{\rm{d}}x = 0\). The forcing term \(\mathbf{f}=(f_1, f_2)\) is of the form \[f_i(x,y)=m_{i,0}\sin(n_{i,0}x+n_{i,1}y)+m_{i,1}\cos(n_{i,2}x+n_{i,3}y),\] with the random samples \(m_{i,0}, \, m_{i,1} \in [1,2],\, n_{i,0}, n_{i,1}, n_{i,2}, n_{i,3} \in [0,1)\times\pi\) for \(i=1,2\).

We use \(P_2\) elements for the velocity components and \(P_1\) elements for the pressure field. This significantly increases the total number of degrees of freedom compared to the previous problems, making the learning process much more challenging. The discrete finite element system consists of three components related to the momentum equation \(-\Delta \mathbf{u} + \nabla p = \mathbf{f}\), the incompressibility condition \(\nabla\cdot\mathbf{u} = 0\), and the zero-mean pressure condition \(\int_\Omega p\;{\rm{d}}x = 0\). We define the training loss as a weighted sum of the \(l_2\) norms of the residuals of these three components. We use the constant weights to balance the scales of the residual terms.

We conducted the experiments on mesh resolutions of \(n=8, 16\) and \(32\). Due to the high computational and memory requirements for this problem, training a dense model becomes impractical. Therefore, we focus on evaluating whether the proposed sparse FEONet converges well in this problem. The results are summarized in Table 8. The sparse model shows high computational efficiency. At \(n=32\), it requires only \(9.18 \%\) of the parameters relative to a dense model. At \(n=16\), the relative \(H^1\) semi-norm error for the velocity field is \(0.0311\), which is close to the error between the reference solution at \(n=128\) and the finite element solution at \(n=16\). At the finer resolution, \(n=32\), the problem becomes more challenging due to the rapid increase in the number of degrees of freedom and physical constraints. To train the model for this problem, we reduced the learning rate and extended the training epochs. Although this improves convergence, the prediction errors increase moderately compared with those on the coarse mesh. Nevertheless, the sparse model still converges reliably and achieves acceptable accuracy under this challenging setting. The same model architecture is used, and the weight settings are the same across the mesh resolutions without additional hyperparameter tuning. Figure 11 visualizes the prediction of the sparse model at \(n=32\).

6 Concluding Remarks↩︎

In this work, we consider the FEONet, an unsupervised operator-learning framework for parametric PDEs based on the classical FEM. While FEONet demonstrates strong accuracy and robustness over a wide range of problems, its computational burden grows rapidly with mesh refinement, and its performance can degrade as the number of elements increases, which limits its applicability to large-scale settings. To overcome these challenges, we proposed a new sparse network architecture guided by the intrinsic locality and connectivity structure of finite elements. The proposed design significantly reduces the computational cost and improves memory efficiency, while preserving accuracy comparable to the original FEONet across extensive numerical experiments. Beyond empirical validation, we established theoretical guarantees: we proved that the sparse architecture can approximate the target operator effectively, and we provided a stability analysis that supports reliable training and prediction. Taken together, these results suggest that incorporating finite-element structure into operator networks is a principled and practical strategy for scaling FEONet to finer discretizations.

Looking ahead, an important future direction is to deepen the theoretical understanding of the proposed method by exploiting properties that are specific to the sparse FEONet and are not present in the original FEONet. In particular, by leveraging the finite-element-induced sparsity pattern and its associated locality structure, it should be possible to carry out a convergence analysis that more directly reflects the role of the new architecture, clarifying how the sparse connectivity affects approximation error and generalization error as the mesh is refined. Such a result would not only strengthen the theoretical foundation of the method but also provide principled guidelines for designing and tuning sparse architectures in large-scale regimes. Another promising avenue is to move beyond multi-layer perceptrons and develop appropriate sparse designs for more general neural network classes. Constructing FEM-relevant sparse structures for broader architectures may significantly widen the applicability of the approach. Establishing conditions under which such sparse generalizations preserve approximation power and stability, while maintaining computational advantages, would be an important and intriguing step toward making the proposed framework a more universal tool for operator learning in parametric PDEs.

Appendix↩︎

In this appendix section, we investigate the effect of connectivity \(C_\ell\) for both the advection-diffusion-reaction equation and the Helmholtz equation. In all experiments, we consider a uniform mesh on the unit square domain described in Section 3, and we report results for mesh resolutions \(n = 32\) and \(n = 64\). The number of hidden layers is fixed to five throughout all experiments. As illustrated in Figure 4, information from a single node can propagate across the entire domain through successive layers, even with local connectivity. Based on this observation, one can estimate the minimum value of \(C_\ell\) required for global information propagation. The resulting values are \(C_\ell=3\) for \(n=16\), \(C_\ell=5\) for \(n=32\), \(C_\ell=11\) for \(n=64\), and \(C_\ell=21\) for \(n=128\).

Tables 9–12 summarize the effect of the connectivity parameter \(C_\ell\) on the performance of the sparse model. Tables 9 and 10 show the results for the advection-diffusion-reaction equation at mesh resolutions \(n = 32\) and \(n = 64\), respectively, while Tables 11 and 12 present the corresponding results for the Helmholtz equation. In all cases, the relative errors are evaluated on the test set using the finite element solution \(u_h\) computed on the same mesh. The results show that increasing the connectivity generally improves the accuracy of the sparse model, although accurate approximations are already achieved at connectivity levels below the corresponding estimated minimum connectivity. In addition, the Helmholtz equation, which involves higher-frequency solution components, requires larger connectivity values to achieve comparable accuracy. These results indicate that the minimum connectivity required for accurate approximation increases with both problem complexity and mesh resolution.

References↩︎

---

1. Department of Mathematics, Inha University, 22212, Incheon, Republic of Korea. Email: scko@inha.ac.kr↩︎

2. Department of Mathematics, Ajou University, 16499, Suwon, Republic of Korea. Email: gkim0201@ajou.ac.kr↩︎

3. Department of Mathematics, Ajou University, 16499, Suwon, Republic of Korea. Email: dws@ajou.ac.kr↩︎