1 Introduction↩︎

Parameter-dependent partial differential equations (PDEs) arise naturally across a broad spectrum of scientific and engineering applications, including forward simulations, control theory, uncertainty quantification, and inverse problems. There, the repeated numerical solution of the underlying PDEs often represents the dominant computational cost, particularly for high-dimensional parameter spaces. Consequently, the efficient approximation of the solution operator, which maps parameter inputs such as initial values or forcing terms to the PDE solution, has emerged as a central problem in numerical analysis and scientific computing.

Recent studies indicate that learned approaches, see [1]–[3] for an overview, can outperform classical techniques [4]. As the overall computational effort is often dominated by the number of required PDE solves for generating the training samples rather than by the training itself, methods should be sample efficient. To formalize operator learning in the PDE context, let \(\mathcal{F}\) denote a (potentially infinite-dimensional) parameter space, \(\mathcal{U}\) a target space of functions on some domain \(\Omega \subset \mathbb{R}^d\), and \(G \colon \mathcal{F} \to \mathcal{U}\) the unknown solution operator. As in reduced-order modeling [5] and consistent with the finite capacity of the hidden representations used by neural operators [6], we project the input parameters onto a finite set of basis coefficients \(\boldsymbol{a} \in \mathbb{R}^n\), yielding a finite-dimensional operator \(G_n \colon \mathbb{R}^n \to \mathcal{U}\). Given training data \(\{(\boldsymbol{a}_i, \boldsymbol{x}_i, u_{\boldsymbol{a}_i}(\boldsymbol{x}_i))\}_{i=1}^M \subset \mathbb{R}^n \times \Omega \times \mathbb{R}\), consisting of parameter vectors paired with solution evaluations (often with multiple spatial evaluations per parameter \(\boldsymbol{a}\)), the goal is to learn an accurate approximation of \(G_n\). A central difficulty is the curse of dimensionality: the sample and computational complexity may scale exponentially in the input dimension \(n\).

To construct a tractable approximation of \(G_n\), one must first choose a suitable parametric architecture. In approximation theory, a standard approach for handling the infinite-dimensional range space \(\mathcal{U}\) is to choose a basis \(\{\Phi_j\}_{j \in J} \subset \mathcal{U}\), such as the Fourier or Chebyshev basis, and to work with the truncated expansion \[\label{eq:OperatorApprox} G_n[\boldsymbol{a}](\boldsymbol{x}) \approx \sum_{j \in J} b_j(\boldsymbol{a}) \Phi_j(\boldsymbol{x}) \qquad \forall\, \boldsymbol{x} \in \Omega.\tag{1}\] This expansion is also used by the DeepONet approach [7]. However, we restrict 1 to a fixed, non-learnable basis. Then, parameterizing \(b \colon \mathbb{R}^n \to \mathbb{R}^{|J|}\) by a neural network \(b(\cdot\,;\boldsymbol{c})\) with parameters \(\boldsymbol{c} \in \mathbb{R}^P\), as in DeepONet’s branch network, leads to the training problem \[\label{eq:NN} \mathop{\mathrm{arg\,min}}_{\boldsymbol{c} \in \mathbb{R}^{P}} \sum_{i=1}^M \biggl\vert \sum_{j \in J} b_j(\boldsymbol{a}_i; \boldsymbol{c}) \Phi_j(\boldsymbol{x}) - u_{\boldsymbol{a}_i}(\boldsymbol{x}_i) \biggr\vert^2.\tag{2}\] By expanding a neural operator [6] output in the basis \(\{\Phi_j\}_{j \in J}\), we can formally interpret its learning within the same framework. While expressive, such nonlinear architectures typically require many training samples and offer limited theoretical guarantees. In contrast, we adopt a simple linear feature model of the form \(b_j(\boldsymbol{a}; \boldsymbol{c}) = \sum_{k \in K} c_{k,j} \Psi_k(\boldsymbol{a})\), where \(\{\Psi_k\}_{k \in K}\) is a fixed dictionary of features \(\Psi_k \colon \mathbb{R}^n \to \mathbb{R}\) and the coefficients \(c_{k,j} \in \mathbb{R}\) are the learnable parameters. Setting \(\boldsymbol{k} = (k,j) \in \Gamma \mathrel{\vcenter{:}}= K \times J\) and \(T_{\boldsymbol{k}}(\boldsymbol{a}, \boldsymbol{x}) \mathrel{\vcenter{:}}= \Psi_k(\boldsymbol{a}) \Phi_j(\boldsymbol{x})\), the resulting approximation 1 can be compactly written as \(G_n[\boldsymbol{a}](\boldsymbol{x}) \approx \sum_{\boldsymbol{k} \in \Gamma} c_{\boldsymbol{k}} T_{\boldsymbol{k}}(\boldsymbol{a}, \boldsymbol{x})\). Then, the generic operator learning task 2 reduces to the least-squares problem \[\label{eq:LeastSquares} \mathop{\mathrm{arg\,min}}_{\boldsymbol{c} \in \mathbb{R}^{\vert \Gamma \vert}} \sum_{i=1}^M \biggl\vert \sum_{\boldsymbol{k} \in \Gamma} c_{\boldsymbol{k}} T_{\boldsymbol{k}}(\boldsymbol{a}_i, \boldsymbol{x}_i) - u_{\boldsymbol{a}_i}(\boldsymbol{x}_i) \biggr\vert^2,\tag{3}\] which is theoretically well analyzed. However, since the number of unknowns \(|\Gamma|\) typically far exceeds the number of training samples \(M\), solving 3 directly leads to severe overfitting. Sparse approximation methods [8]–[11] address this by assuming that \(G_n\) can be represented accurately using only a small number of nonzero coefficients \(c_{\boldsymbol{k}}\). This assumption is not merely a regularization device: the analysis of specific PDEs shows that the analytical solutions indeed satisfy such sparsity constraints [12]. Together, these observations motivate the constrained training problem \[\label{eq:SparseRec} \mathop{\mathrm{arg\,min}}_{I \subset \Gamma:\, |I| \leq s}\; \min_{\boldsymbol{c} \in \mathbb{R}^{|I|}} \sum_{i=1}^M \biggl\vert \sum_{\boldsymbol{k} \in I} c_{\boldsymbol{k}} T_{\boldsymbol{k}}(\boldsymbol{a}_i, \boldsymbol{x}_i) - u_{\boldsymbol{a}_i}(\boldsymbol{x}_i) \biggr\vert^2.\tag{4}\] The problem 4 can be solved, e.g., with greedy algorithms such as orthogonal matching pursuit (OMP) [13]. Under suitable conditions on the feature dictionary and the sample distribution, robust recovery guarantees for 4 are available; see, e.g., [14]–[16].

While sparse recovery addresses the overfitting problem, the cardinality of the search space \(\Gamma\) can itself be prohibitively large. For a product-type basis for \(\Phi_j\) and \(\Psi_k\) with, e.g., \(64\) basis functions per direction, the number of candidate indices already grows to \(|\Gamma| \approx 64^{n+d}\), rendering an exhaustive search computationally infeasible. As a remedy, the dimension-incremental framework of [17]–[19] constructs the support set \(I\) dimension by dimension: beginning with a one-dimensional candidate set, it successively extends the current support rather than searching \(\Gamma\) directly. The extension step involves solving several sparse recovery problems of the form 4 , each with a tractable number of candidates. In this way, the method can efficiently exploit the underlying sparsity structure while considerably reducing both computational and sample complexity compared to an exhaustive search. A recent alternative for managing large (continuous) candidate sets is generative feature training [20], which adaptively resamples the candidates for 4 in place of a greedy approach.

1.0.0.1 Contributions

Our contributions can be summarized as follows.

• We combine the dimension-incremental support detection framework of [17], [18] with OMP-based sparse recovery [16] to compute sparse approximations of PDE solution operators from carefully selected samples. Compared to cubature-based recovery [19], this significantly reduces the sample size and consequently the number of PDE solves.

• By design, the cubature-based recovery in [12] requires an independent PDE solve for each sample. In contrast, sparse recovery methods can naturally incorporate multiple samples with different spatial locations obtained from a single solve, which improves the computational efficiency even further.

• We numerically validate the proposed methods for several PDEs, comparing approximation quality, runtime, and sampling complexity against cubature-based sparse approximation methods and Tensorized Fourier Neural Operators (TFNOs) as a baseline.

• Unlike purely data-driven black-box approaches, our approach provides interpretable information about the parameter dependence of the solution via the detected index sets and their interaction structure.

Figure 1 provides a schematic overview of the proposed approximation strategies together with the methods used for comparison. Moreover, Table 1 summarizes typical parameter sizes occurring in the numerical experiments of Section 3, as well as the corresponding numbers of PDE solves required by the different approaches. The Python implementation and numerical experiment scripts are publicly available at [21].

1.0.0.2 Outline

The remainder of this work is organized as follows. Section 2 discusses the parametric representation of PDE solution operators and their sparse approximation. Subsequently, we recall the dimension-incremental support detection framework and combine it with OMP-based coefficient recovery and an improved sampling strategy for PDE-based problems. Section 3 presents numerical experiments for the heat equation, Burgers’ equation, and a parametric diffusion equation, comparing approximation quality, runtime, and sampling complexity against cubature-based approaches and TFNOs. Finally, Section 3.4 discusses the numerical results and the trade-offs between approximation quality and computational cost.

1.0.0.3 Related work

The approximation of (infinite-dimensional) operators from training data has recently attracted significant attention, with approaches broadly falling into three categories: linear expansions, neural operators, and (random) kernel methods.

From the viewpoint of approximation theory, operator learning can be cast as finding the coefficients of a basis expansion [10], [22]. This perspective is closely related to sparse polynomial approximation [9], [23]: if the solution operator admits a sparse representation in a suitable basis, only a small number of coefficients need to be recovered [24]. To this end, a PDE-oriented recovery strategy was investigated in [12]. Alternatively, one can directly learn mappings between the basis expansion coefficients of input and output functions using neural networks [25], [26] or approximate the Green’s function instead [27].

Neural operators are typically trained in a supervised setting using (potentially physics-informed [28]) loss functions. Prominent architectures include DeepONets [7], [29] and Fourier Neural Operators (FNOs) [30], many of which are implemented as part of the NeuralOperator library [31]. Spectral neural operators and related architectures replace purely Fourier-based representations with polynomial or spectral basis expansions [32], [33], a perspective closely related to this work. Other recent directions include the incorporation of attention mechanisms [34], [35] and the aim for foundation models [36]. Comprehensive surveys covering practical and theoretical aspects of neural operators can be found in [1], [2], [37].

Operator-valued kernels provide a principled framework [38], though the direct implementation is computationally prohibitive in general. Random feature approaches [39] offer a scalable alternative: by fixing the random parameters during training, the optimization reduces to a linear problem. Interestingly, kernel-based methods can achieve accuracy competitive with neural operators [40]. A related Gaussian process framework was recently proposed in [41].

Finally, we remark that there is also growing interest in topics such as error estimates [6], [42], discretization consistency [43] and sample complexity [44] of neural operators.

2 Theory↩︎

The presentation of the general framework follows [12], where operator learning for partial differential equations (PDEs) is investigated from a sparse high-dimensional approximation perspective. For the purposes of this paper, we adapt their approach and incorporate dedicated sparse recovery algorithms. While we restrict ourselves to time-dependent PDEs as a running example below, we stress that the proposed method applies to a broad class of operators, which also includes stationary and parametric PDEs.

2.1 Time-dependent partial differential equations↩︎

Let \(\Omega \subset \mathbb{R}^d\) be a bounded spatial domain and \(T>0\). Further, assume that we are given function spaces \(\mathcal{U}\) and \(\mathcal{F}\) defined on the space-time domain \(\Omega \times (0,T)\), a function space \(\mathcal{H}\) on \(\partial \Omega \times (0,T)\), and a function space \(\mathcal{U}_0\) on \(\Omega\). Here, the pairs \((\mathcal{U}, \mathcal{U}_0)\) and \((\mathcal{U}, \mathcal{H})\) should be compatible in the sense that \(u(\cdot,0) \in \mathcal{U}_0\) and \(u|_{\partial\Omega \times (0,T)} \in \mathcal{H}\) for any \(u \in \mathcal{U}\). Within this setting, we consider differential equations of the form \[\begin{align} \label{eq:general} L u = f \qquad \qquad \text{s.t. } u|_{\partial\Omega \times (0,T)} = h \text{ and } u(\cdot,0) = u_0, \end{align}\tag{5}\] where \(L \colon \mathcal{U}\to\mathcal{F}\) is a (possibly nonlinear) differential operator, \(f\in\mathcal{F}\) is an inhomogeneity, and \(u_0 \in \mathcal{U}_0\) is an initial value. Throughout, we assume that 5 is well-posed in the sense that for every \((f,h,u_0) \in \mathcal{F} \times \mathcal{H} \times \mathcal{U}_0\), there exists a unique solution \(u \in \mathcal{U}\) depending continuously on the data \((f,h,u_0)\). Then, 5 induces a solution operator \[\label{eq:solutionOp} G\colon \mathcal{F} \times \mathcal{H} \times \mathcal{U}_0 \to \mathcal{U}, \qquad (f,h,u_0) \mapsto u.\tag{6}\] In practice, often only a subset of the data \((f,h,u_0)\) is varied, while the remaining components are fixed. In particular, a common situation is that the solution operator \(G\) reduces to a mapping from the initial data \(u_0\) to the solution \(u\), with \(f\) and \(h\) being fixed. Then, our goal is to approximate \(G\) by a surrogate model that enables fast evaluations for new input data \(u_0\) without the need for solving the PDE 5 .

2.2 Parametric operator representation and sparse approximation↩︎

In order to obtain a finite-dimensional representation of the domain of \(G\), we expand the initial condition \(u_0\) as \[\begin{align} u_0(\boldsymbol{x}) \approx \sum_{j=1}^n a_j A_j(\boldsymbol{x}), \qquad \boldsymbol{x} \in \Omega, \label{eq:u095approx} \end{align}\tag{7}\] with basis functions \(A_j \colon \Omega \to \mathbb{R}\), \(j=1,\ldots,n\), and coefficients \(\boldsymbol{a} = (a_1,\dots,a_n)\in[-1,1]^n\). Note that the chosen functions \(A_j\) should lead to a sufficiently accurate approximation in \(\mathcal{U}_0\). Based on 7 , we get the induced solution operator \(G_n \colon \mathbb{[}-1,1]^n \to \mathcal{U}\) with \(\boldsymbol{a} \mapsto u\).

In general, any appropriate discretization \(G_n\) of \(G\) can be identified with a function defined on an extended bounded product domain . In the case of 7 , we get \(\mathcal{D} = \Omega \times [0,T] \times [-1,1]^n\). For the next step, we assume that the induced function admits an expansion in a bounded orthonormal product basis (BOPB) \(\{\Phi_{\boldsymbol{k}} : \boldsymbol{k} \in \mathbb{N}^{d_\mathrm{ex}}\} \subset L_2(\mathcal{D},\mu)\), where \(\mu\) is a product measure. Recall that a BOPB consists of tensor-products \(\Phi_{\boldsymbol{k}}(\boldsymbol{y}) = \prod_{j=1}^{d_\mathrm{ex}} \phi_{j,k_j}(y_j)\). More details can be found, e.g., in [19].

For the operator \(G_n\) induced by 7 , the BOPB expansion reads as \[\begin{align} u(\boldsymbol{x}, \tau; \boldsymbol{a}) = \sum_{\boldsymbol{k} \in \mathbb{N}^{d_\mathrm{ex}}} c_{\boldsymbol{k}}\,\Phi_{\boldsymbol{k}}(\boldsymbol{x}, \tau, \boldsymbol{a}), \qquad (\boldsymbol{x}, \tau, \boldsymbol{a})\in \mathcal{D}, \label{eq:basexp95of95u} \end{align}\tag{8}\] with \(c_{\boldsymbol{k}}\in \mathbb{C}\) being the corresponding basis coefficients \[\begin{align} \label{eq:integral95ck} \hat{c}_{\boldsymbol{k}} = \int_\mathcal{D} g(\boldsymbol{x}, \tau, \boldsymbol{a}) \overline{\Phi_{\boldsymbol{k}}(\boldsymbol{x}, \tau, \boldsymbol{a})}\mathrm d (\boldsymbol{x}, \tau, \boldsymbol{a}). \end{align}\tag{9}\] Then, our aim is to accurately approximate \(u\) (if possible) by a sparse truncation \[\begin{align} S_I u(\boldsymbol{x}, \tau; \boldsymbol{a}) = \sum_{\boldsymbol{k}\in I} c_{\boldsymbol{k}}\,\Phi_{\boldsymbol{k}}(\boldsymbol{x}, \tau, \boldsymbol{a}), \label{eq:basexp95of95u95trunc1} \end{align}\tag{10}\] where \(I \subset \mathbb{N}_0^{d_\mathrm{ex}}\) is a finite index set. Since the exact \(c_{\boldsymbol{k}}\) are unknown, we replace them by approximations \(\hat{u}_{\boldsymbol{k}}\) (for which we outline a computation procedure below) and obtain \[\begin{align} S_I^\mathcal{A}u(\boldsymbol{x}, \tau, \boldsymbol{a}) = \sum_{\boldsymbol{k}\in I} \hat{u}_{\boldsymbol{k}}\,\Phi_{\boldsymbol{k}}(\boldsymbol{x}, \tau, \boldsymbol{a}). \label{eq:basexp95of95u95trunc} \end{align}\tag{11}\] In summary, 11 yields a continuous approximation of \(G\) by first projecting \(u_0\) onto its finite coefficient vector \(\boldsymbol{a}\) and then evaluating \(S_I^\mathcal{A}u(\cdot,\cdot,\boldsymbol{a})\). The resulting approximation error can be split into a truncation part and a coefficient recovery part. More precisely, using the uniform boundedness of the basis functions \(\Phi_{\boldsymbol{k}}\), we obtain \[\begin{align} \left\Vert u-S_I^\mathcal{A}u\right\Vert_\infty &\leq \left\Vert u-S_Iu\right\Vert_\infty + \left\Vert S_Iu-S_I^\mathcal{A}u\right\Vert_\infty \notag\\ &\leq \sum_{\boldsymbol{k}\notin I} \left\vert c_{\boldsymbol{k}}\right\vert \left\Vert \Phi_{\boldsymbol{k}}\right\Vert_\infty + \sum_{\boldsymbol{k}\in I} \left\vert c_{\boldsymbol{k}}-\hat{u}_{\boldsymbol{k}}\right\vert \left\Vert \Phi_{\boldsymbol{k}}\right\Vert_\infty \notag\\ &\leq \sup_{\boldsymbol{k} \in I} \left\Vert \Phi_{\boldsymbol{k}}\right\Vert_\infty \biggl( \sum_{\boldsymbol{k}\notin I} \left\vert c_{\boldsymbol{k}}\right\vert + \sum_{\boldsymbol{k}\in I} \left\vert c_{\boldsymbol{k}}-\hat{u}_{\boldsymbol{k}}\right\vert \biggr). \end{align}\] Hence, a good approximation requires both identifying a suitable index set \(I\) and an accurate method for estimating the coefficients \(c_{\boldsymbol{k}}\) from samples \(u(\boldsymbol{x}_j,\tau_j,\boldsymbol{a}_j)\). For our running example, we keep in mind that samples are obtained by numerically solving the PDE 5 for the parameters \(\boldsymbol{a}_j\), which makes sampling costly.

Given a finite index set \(I\), the corresponding coefficients \(c_{\boldsymbol{k}},\, \boldsymbol{k} \in I,\) can be estimated using a cubature formula based on 9 . Alternatively, they may be obtained by solving the linear system that arises from inserting the samples \(u(\boldsymbol{x}_j,\tau_j,\boldsymbol{a}_j)\) into 8 , for instance via a standard least squares approach. However, since the index set \(I\) is not known a priori, one typically has to work with the product space \(\Gamma \supset I\) instead. Without imposing strong a priori assumptions, such search spaces grow rapidly in high-dimensional settings and render the aforementioned approaches computationally infeasible. We therefore turn to the problem of determining a suitable index set \(I \subset \Gamma \subset \mathbb{N}_0^{d_\mathrm{ex}}\).

2.3 Dimension-incremental support detection and sparse recovery↩︎

Since the candidate set \(\Gamma\) has product structure, its cardinality scales as \(10^{d_\mathrm{ex}}\) even with only 10 basis functions per dimension. For \(d_\mathrm{ex} = 12\), this already yields \(10^{12}\) candidate index tuples, making an exhaustive search for a suitable \(I\) in general computationally infeasible.

As an alternative, we employ the dimension-incremental support detection strategy proposed in [19]. The main idea is to build \(I\) progressively, starting from one-dimensional projections and then combining the detected candidates across dimensions. Below we provide a brief summary; see Algorithm 2 and [19] for further details and a detection guarantee for a suitable index set \(I\). To simplify the notation, we switch from the problem-specific function \(u\) to a generic function \(g\), and correspondingly collect all variables into a single \(d_\mathrm{ex}\)-dimensional input vector.

Let \(g\colon \mathbb{R}^{d_\mathrm{ex}}\to\mathbb{C}\) be an (unknown) target function with (\(s\)-sparse) expansion \[\label{eq:expansion} g = \sum_{\boldsymbol{k}\in I} \hat{g}_{\boldsymbol{k}}\Phi_{\boldsymbol{k}},\tag{12}\] where \(I\) with \(\vert I \vert =s\) is contained in a prescribed finite candidate space \(\Gamma \subset \mathbb{N}_0^{d_\mathrm{ex}}\). Typically, \(\Gamma\) is simply a cube in \(\mathbb{N}_0^{d_\mathrm{ex}}\). For each coordinate \(t\in\{1,\dots,d_\mathrm{ex}\}\), we consider the one-dimensional projected candidate set \(\mathcal{P}_{\{t\}}(\Gamma) \mathrel{\vcenter{:}}= \{ l \in \mathbb{N}_0 \,|\,\exists \boldsymbol{k} \in \Gamma: k_t=l\}\) and for each \(l \in \mathcal{P}_{\{t\}}(\Gamma)\) the so-called projected coefficients \[\begin{align} \hat{g}_{\{t\},l}(\boldsymbol{\theta}) \mathrel{\vcenter{:}}= \int_{\mathcal{D}_t} g_{\{t\}}(\xi,\boldsymbol{\theta}) \overline{\phi_{t,l}(\xi)} \,\mathrm{d}\xi = \sum_{\boldsymbol{k}\in I : k_t = l} \hat{g}_{\boldsymbol{k}} \prod_{j \neq t} \phi_{j, k_j}(\theta_j), \label{eq:proj95coef951d95new} \end{align}\tag{13}\] where \(\boldsymbol{\theta} \in \mathbb{R}^{d_\mathrm{ex}-1}\) is an anchor and \(g_{\{t\}}(\xi,\boldsymbol{\theta})\) denotes the value of \(g\) with \(\xi\) in the \(t\)-th dimension and \(\boldsymbol{\theta}\) for the remaining coordinates. Due to the representation 12 , at most \(s\) coefficients \(\hat{g}_{\{t\},l}(\boldsymbol{\theta})\) should be nonzero. Thus, among all coefficients \(l \in \mathcal{P}_{\{t\}}(\Gamma)\), we retain the \(s\) largest ones (in terms of \(|\hat{g}_{\{t\},l}(\boldsymbol{\theta})|\)) exceeding a prescribed threshold, and discard all others as negligible. All detected candidates \(l\) are collected in the temporary index sets \(I_{\{t\}}\).

Next, for \(t\geq 2\), the candidate sets are combined incrementally. Given \(I_{\{1,\dots,t-1\}}\), one forms the sets \[\begin{align} \mathcal{P}_{\{1,\ldots,t\}}(\Gamma) &\mathrel{\vcenter{:}}= \{ \boldsymbol{l} \in \mathbb{N}_0^t \,|\,\exists \boldsymbol{k} \in \Gamma \subset \mathbb{N}_0^{d_\mathrm{ex}}: \boldsymbol{k}_{\{1,\ldots,t\}}= \boldsymbol{l}\}\\ K &\mathrel{\vcenter{:}}= \bigl(I_{\{1,\dots,t-1\}}\times I_{\{t\}}\bigr) \cap \mathcal{P}_{\{1,\dots,t\}}(\Gamma), \label{eq:candidate95set} \end{align}\tag{14}\] and investigates the corresponding projected coefficients on \(K\) as in the one-dimensional case. In particular, we again consider the projected coefficients \[\begin{align} \hat{g}_{\{1,\ldots,t\},\boldsymbol{l}}(\boldsymbol{\theta}) \mathrel{\vcenter{:}}=\!\!\! \int\limits_{\prod_{j=1}^t \mathcal{D}_j} g_{\{1,\ldots,t\}}(\boldsymbol{\xi},\boldsymbol{\theta}) \overline{\Phi_{\{1,\ldots,t\},\boldsymbol{l}}(\boldsymbol{\xi})} \,\mathrm{d}\boldsymbol{\xi} = \!\!\!\! \sum_{\boldsymbol{k}\in I : \boldsymbol{k} _{\{1,\ldots,t\}} \!\!= \boldsymbol{l}} \hat{g}_{\boldsymbol{k}} \!\!\!\prod_{j \notin \{1,\ldots,t\}} \Phi_{j, k_j}(\theta_j), \label{eq:proj95coef95hd95new} \end{align}\tag{15}\] estimate their absolute value, retain the \(s\) largest ones exceeding the threshold, and collect them in the index set \(I_{\{1,\dots,t\}}\). In this way, higher-order interactions are identified successively until the final index set \(I_{\mathrm{est}} \mathrel{\vcenter{:}}= I_{\{1,\dots,d_\mathrm{ex}\}}\) is obtained.

In [12], the projected coefficients defined in 13 and 15 are approximated via cubature rules, which typically require a large number of sampling points \(\smash{\boldsymbol{\xi}^{(j)}}\) at carefully chosen locations. Specifically, so-called rank-1 lattices are employed, which provide accurate (and in some cases even exact) cubature for functions with suitable smoothness or frequency structure. When using large lattices, this enables reliable and efficient recovery of the projected coefficients for all candidate indices, including the ones that are actually negligible.

In contrast to this exhaustive approximation, we propose to directly exploit the system’s inherent sparsity. More precisely, 12 implies that also the projected expansions \[\label{eq:SparseRecProb} g_{\{t\}}(\xi,\boldsymbol{\theta}) = \sum_{l \in \mathcal{P}_{\{t\}}(\Gamma)} \hat{g}_{\{t\},l}(\boldsymbol{\theta}) \Phi_{\{t\},l}(\xi),\tag{16}\] and similarly the ones with 15 , are \(s\) sparse. Classical sparse recovery algorithms, such as the orthogonal matching pursuit (OMP) [8] discussed in Section 2.4, typically require only around \(M=5s\) random sample locations \(\boldsymbol{\xi}^{(j)}\), thereby considerably improving the computational efficiency of Algorithm 2 compared to [12]. Since our approach treats \(g\) as a black box, any numerical scheme, such as finite elements, finite differences, or spectral methods, can be used to compute the function values \(g_{\{1,\ldots,t\}}(\boldsymbol{\xi}^{(j)}, \boldsymbol{\theta})\). In turn, the accuracy of these evaluations directly affects the quality of the computed approximation 11 . Finally, we emphasize that the sampled PDE input parameters \(\boldsymbol{a}_i\) must conform to the anchor decomposition: each \(\boldsymbol{a}_i\) is partitioned into active coordinates, sampled randomly, and fixed anchor coordinates.

2.4 Orthogonal Matching Pursuit↩︎

Now, we briefly discuss OMP, which is used to reconstruct the relevant projected coefficients in lines 7 and 19 of Algorithm 2. This algorithm is a greedy approach developed for sparse approximation and compressed sensing. For a broad overview, we refer to [8], [13], [45]. Assume, we are in the \(t\)-th dimension-incremental step of Algorithm 2 and consider the candidate set \(K\) as given in 14 and sampling values \(g_{\{1,\ldots,t\}}(\boldsymbol{\xi}) \mathrel{\vcenter{:}}= (g_{\{1,\ldots,t\}}(\boldsymbol{\xi}^{(j)}))_{j=1}^M \in \mathbb{R}^M\). Starting from an empty active set, OMP iteratively selects at each step the product of basis functions \(\phi_{\{1,\ldots,t\}, \boldsymbol{l}}\) for which the corresponding column of the basis matrix \[\begin{align} \label{eq:OMP95basis95mat} \boldsymbol{\Phi} \mathrel{\vcenter{:}}= (\phi_{\{1,\ldots,t\},\boldsymbol{l}}(\boldsymbol{\xi}^{(j)})_{j=1,\ldots,M,\, \boldsymbol{l} \in K} \in \mathbb{C}^{M \times \vert K \vert} \end{align}\tag{17}\] is most strongly correlated with the current residual. The active set is then updated with the index \(\boldsymbol{l}\), and the basis coefficients \({\hat{g}}_{\{1,\ldots,t\},\boldsymbol{l}}(\boldsymbol{\theta})\) restricted to this set are determined via a least squares solve. This procedure is repeated until the prescribed sparsity \(s\) is reached or a suitable stopping criterion, e.g., an absolute or relative error tolerance, is satisfied.

Under suitable assumptions on the sampling points , the resulting Algorithm 3 is guaranteed to recover the correct support of the basis coefficients \({\hat{g}}_{\{1,\ldots,t\},\boldsymbol{l}}(\boldsymbol{\theta})\), \(\boldsymbol{l} \in K,\) with high probability provided that the number of samples scales like \(M \gtrsim s \log |K|\). The complexity of a single iteration can be bounded by \(\mathcal{O}(M |K| + M s + s^2)\), cf. [46], which is dominated by the cost of the matrix-vector multiplication (\(M |K|\)). Performing \(s\) iterations leads to the overall computational complexity \(\mathcal{O}(s M |K|)\). Our implementation is based on the methodology presented in [16], with a Python implementation being publicly available¹.

2.5 Improving Sample Efficiency (OMP+)↩︎

In Algorithm 2, each evaluation of \(u\) at a node \((\boldsymbol{x}^{(j)}, \tau^{(j)}, \boldsymbol{a}^{(j)}) \in \mathcal{D}\) requires a separate numerical solve of 5 if the parameter instances \(\boldsymbol{a}^{(j)}\) differ. As described so far, Algorithm 2 incorporates random samples \(\{\boldsymbol{\xi}^{(j)}\}_{j=1}^M\) conditional to a random anchor, and thus in principle only a single value from each PDE solve, even though the solution is typically available at many grid points \((\boldsymbol{x}, \tau)\). Given the substantial cost of PDE solves, this represents a significant inefficiency. Fortunately, in contrast to the cubature-based approach in [12], OMP imposes no structural constraints on the sampling nodes \(\{\boldsymbol{\xi}^{(j)}\}_{j=1}^M \in \mathbb{R}^{t}\). Thus, we can systematically incorporate multiple samples from a computed solution, as detailed below.

Consider a \(t\)-th dimension-incremental step with \(t > d+1\) in Algorithm 2. Then, we have with \(\boldsymbol{a}^{(j)}_1 \in \mathbb{R}^{t - d - 1}\) and the anchor \(\boldsymbol{\theta} = \boldsymbol{a}_2 \in \mathbb{R}^{n + d + 1 - t}\). Instead of evaluating the PDE solution \(u(\cdot,\cdot,\boldsymbol{a}^{(j)})\) for , \(j = 1,\ldots,M\), only at the sampled \((\boldsymbol{x}^{(j)}, \tau^{(j)})\), we additionally evaluate them at \(b^{d+1}\) randomly chosen points \((\boldsymbol{x},\tau) \in \mathbb{R}^{d+1}\), which of course differ for each \(\boldsymbol{a}^{(j)}\). We recommend \(b=3\) as default. The number of evaluations per solve must be carefully chosen to balance the exploration of the parameter and space-time domains. A smaller amount of space-time evaluations requires a higher number of costly PDE solves to achieve a given target sample size. Conversely, too many evaluations per parameter result in insufficient coverage of the parameter domain (the extreme case being only a single parameter). Moreover, since the complexity of OMP scales linearly with \(M\), adding more space-time evaluations is only recommended when they contain non-redundant information. In our experiments, we refer to this modified version of OMP as OMP+.

3 Numerics↩︎

We benchmark Algorithm 2 for several PDEs, in both its standard version OMP and the modified variant OMP+ from Section 2.5. For \(\boldsymbol{a} \in [-1,1]^n\), we measure the approximation accuracy in terms of the absolute \(\ell_2\)-error \[\begin{align} \label{eq:abs95err} \text{err}_{\text{abs}}(\boldsymbol{a}) &\mathrel{\vcenter{:}}= \left\Vert S_{I_{\mathrm{est}}}^\mathcal{A}u(\boldsymbol{x},\tau,\boldsymbol{a}) - u(\boldsymbol{x},\tau,\boldsymbol{a})\right\Vert_{\ell_2} = \biggl(\sum_{j=1}^{G}\left\vert S_{I_{\mathrm{est}}}^\mathcal{A}u(\boldsymbol{x}^{(j)},\tau^{(j)},\boldsymbol{a}) - u(\boldsymbol{x}^{(j)},\tau^{(j)},\boldsymbol{a})\right\vert^2\biggr)^{\frac{1}{2}}, \end{align}\tag{18}\] which coincides with the quality criterion in Algorithm 2, and the (scale-invariant) relative \(\ell_2\)-error \[\begin{align} \label{eq:rel95err} \text{err}_{\text{rel}}(\boldsymbol{a}) &\mathrel{\vcenter{:}}= \frac{\left\Vert S_{I_{\mathrm{est}}}^\mathcal{A}u(\boldsymbol{x},\tau,\boldsymbol{a}) - u(\boldsymbol{x},\tau,\boldsymbol{a})\right\Vert_{\ell_2}}{\left\Vert u(\boldsymbol{x},\tau,\boldsymbol{a})\right\Vert_{\ell_2}} \notag\\ &= \frac{\left(\sum_{j=1}^{G}\left\vert S_{I_{\mathrm{est}}}^\mathcal{A}u(\boldsymbol{x}^{(j)},\tau^{(j)},\boldsymbol{a}) - u(\boldsymbol{x}^{(j)},\tau^{(j)},\boldsymbol{a})\right\vert^2\right)^{\frac{1}{2}}}{\left(\sum_{j=1}^{G}\left\vert u(\boldsymbol{x}^{(j)},\tau^{(j)},\boldsymbol{a})\right\vert^2\right)^{\frac{1}{2}}}. \end{align}\tag{19}\] In 18 and 19 , the \((\boldsymbol{x}^{(j)},\tau^{(j)})\), \(j=1,\ldots,G\), are equidistant grid points in the space-time domain \(\Omega \times [0,T]\) (or spatial domain \(\Omega\) for stationary problems). We statistically compare these quantities for several randomly drawn coefficients \(\boldsymbol{a}\) in terms of range (min and max), median, first quartile, and third quartile.

Throughout, we employ the tensorized Chebyshev polynomial basis, which is given by \[\begin{align} T_{\boldsymbol{k}}(\boldsymbol{z}) \mathrel{\vcenter{:}}= \prod_{j=1}^{d_\mathrm{ex}} T_{k_j}(z_j) \quad\text{with}\quad T_{k_j}(z_j) \mathrel{\vcenter{:}}= \begin{cases} 1 & k_j=0 \\ \sqrt2 \cos(k_j \arccos(z_j)) & k_j \not=0 \end{cases}, \end{align}\] where \(d_\mathrm{ex} = d+1+n\) denotes the dimension of the space-time-parameter domain as in Section 2.2. This choice is motivated by the fact that several classes of parameter-dependent PDEs admit sparse or approximately sparse representations in terms of Chebyshev polynomials [12]. If not stated otherwise, Algorithm 2 uses the following parameters:

• the search space \(\Gamma\) is the (non-negative) full grid \([0,N]^{d+1+n}\) in \(d+1+n\) dimensions (or \(d+n\) for stationary examples) with extension \(N\);

• the detection threshold \(\delta_+ = 10^{-12}\);

• \(r=5\) detection iterations for OMP and \(r=2\) for OMP+, see also Remark 3;

• an individual sparsity level \(s\) for each experiment.

See [19] for more detailed information on these parameters and how they affect the behavior of the support detection. For the one-dimensional detection steps, we use \(s\) samples. This suffices because typically \(N \ll s\), which ensures that the number of candidates in these steps (namely \(N+1\)) remains much smaller than \(s\). For the coupled detection steps, the standard OMP approach utilizes an oversampling factor of \(c=5\), leading to \(5s\) samples and PDE solves. In contrast, the OMP+ approach utilizes \(3^{d+1}\) (or \(3^d\) for stationary problems) additional space-time evaluations as discussed in Section 2.5, reducing the requirement to only \(s\) PDE solves. In both cases, the resulting sparse recovery problems 16 are solved using the OMP Algorithm 3 detailed in Section 2.4.

We compare our methods against the dimension-incremental method using Chebyshev multiple rank-1 lattices (cMR1L) [12] that inspired our work. Further, we consider Tensorized Fourier Neural Operators (TFNO), implemented in PyTorch using the NeuralOperator library [31]. Given training samples , the underlying Fourier Neural Operator (FNO) framework learns the solution operator \(G \colon u_0 \mapsto u\) from data by minimizing the supervised Sobolev loss \[\begin{align} \label{eq:TFNO95min} \min_{\theta} \;\frac{1}{M} \sum_{j=1}^M \|G_\theta(u_0^{(j)}) - u^{(j)}\|_{H^1}^2. \end{align}\tag{20}\] The neural operator \(G_\theta\) is realized with a pointwise lifting layer, followed by a composition of multiple mappings (frequently called Fourier layers) of the form \[\begin{align} v_{k+1}(\boldsymbol{x}) = \sigma \left( W_k v_k(\boldsymbol{x}) + \mathcal{F}^{-1} \big( R_k \cdot \mathcal{F}(v_k) \big)(\boldsymbol{x}) \right), \end{align}\] where \(R_k\) are learned Fourier multipliers on a truncated set of modes and \(W_k\) are pointwise linear mappings, and a final pointwise projection layer. For more details on FNOs and their tensorized variants, we refer to [30], [49]. All TFNO models use \(4\) Fourier layers with \(64\) hidden channels and are trained using the AdamW optimizer with initial learning rate \(10^{-3}\) and weight decay \(10^{-4}\). Training is performed for 1000 epochs with a batch size of 64 using random parameter samples \(\boldsymbol{a}_i\) with corresponding solutions \(u_{\boldsymbol{a}_i}\) evaluated at resolution \(100\). Further, we incorporate a cosine annealing schedule that decays the learning rate to \(10^{-8}\) and is stepped after each epoch.

All experiments are performed in Python^TM and can be found together with the algorithm in [21]. Unless stated otherwise, all computations are carried out on two AMD EPYC 9534 64-Core processors. To accelerate the sampling process, parallelization with up to \(96\) workers is employed. The TFNO models are trained separately on a NVIDIA GeForce RTX 4070 Ti Super GPU. Individual runtimes are discussed in the corresponding sections.

3.1 Heat equation↩︎

Our first numerical example is the heat equation with homogeneous boundary conditions on the domain \(\Omega = (0,1)\) with final time \(T=1\), namely \[\label{eq:heat951d} \begin{alignedat}{3} \partial_\tau u - \alpha^2 \partial_{xx} u &= 0, & \qquad x \in (0,1),\, \tau&\in (0,1) \\ u(x,0) &= f(x), & \qquad x &\in (0,1)\\ u(0,\tau) = u(1,\tau) &= 0, & \qquad \tau &\in (0,1). \end{alignedat}\tag{21}\] For the experiments, we choose \(\alpha = 0.25\). As in [12], we use initial conditions given by a truncated sine expansion, namely \[\begin{align} \label{eq:heat95initial} u(x,0) = f(x) = \sum_{\ell=1}^9 a_\ell \sin (\ell \pi x), & & x \in [0,1] \end{align}\tag{22}\] with coefficients \(a_\ell \in [-1,1], \ell \in \mathbb{N}\), which leads to the analytical solution \[\begin{align} \label{eq:heat95exact} u(x,\tau) = \sum_{\ell=1}^9 a_\ell \sin \left(\ell \pi x\right) \exp \left(-\ell^2\pi^2\alpha^2 \tau \right), & & x \in [0,1], \tau\in [0,1]. \end{align}\tag{23}\]

Clearly, the truncation together with the condition \(a_\ell \in [-1,1]\) restricts the set of admissible initial conditions \(f\). Functions far outside this class require a modified parameterization, for instance by increasing the number of retained modes the range of the coefficients. As explained in Remark 2, we shift both PDE variables \(x\) and \(t\) via the invertible transformation \(\mathcal{T}\colon [-1,1]\to[0,1]\) with \(\mathcal{T}(z) = \frac{1}{2}(z+1)\) to be able to use the Chebyshev basis.

For the sampling, we solve the PDE 21 using SciPy’s solve_ivp with the implicit Radau IIA Runge–Kutta method of order five, a relative tolerance of \(10^{-8}\), and an absolute tolerance of \(10^{-10}\). In Algorithm 2, we use sparsity \(s=1000\) and extension \(N=64\) for the search space \(\Gamma\), i.e., \(\Gamma=[0,64]^{11}\) and \(|\Gamma|\approx 8.75 \cdot 10^{19}\).

The cMR1L comparison from [12] uses the same parameters. As a learning-based baseline, we train a TFNO with \(529\,233\) parameters on \(10\,000\) and \(40\,000\) training samples, respectively, resulting in \(10\,100\) and \(40\,500\) total PDE solves including validation data. We denote these approaches as TFNO\(_{10k}\) and TFNO\(_{40k}\), respectively.

Figure 4 illustrates the relative approximation error for all methods. OMP+ achieves a median error of approximately \(4\cdot 10^{-5}\), comparable to cMR1L. This improves over OMP, which attains a median error of roughly \(2\cdot 10^{-4}\). The baselines TFNO\(_{40k}\) and TFNO\(_{10k}\) yield slightly larger median errors of about \(3\cdot 10^{-4}\) and \(5\cdot 10^{-4}\), respectively.

Figure [fig:heat95runtime] shows the number of solves, the corresponding computation times, and the remaining runtime excluding PDE solves. OMP required slightly more than \(3\) hours, with the dominant portion spent on solving the PDEs for sample generation. Similarly, cMR1L required roughly \(700\) hours due to the significantly larger number of required samples. The PDE-specific OMP+ version was the fastest overall, with a total runtime of about \(1\) hour, roughly half of which was spent on solving the PDEs. The TFNO\(_{10k}\) and TFNO\(_{40k}\) required approximately \(6\) and \(23\) hours of training time, respectively.

The effect of the sparsity parameter \(s\) on the approximation quality is illustrated in Figure [fig:heat95OMP4395err] for OMP+ with extension \(N=64\). The relative error \(\text{err}_{\text{rel}}\) decreases as \(s\) increases, while the spread of errors remains within roughly one order of magnitude in each setting.

3.2 Burgers’ equation↩︎

Our second example is the non-linear 1D Burgers’ equation with homogeneous boundary conditions, namely \[\label{eq:burger951d} \begin{alignedat}{3} \partial_\tau u + u \partial_x u &= \nu \partial_{xx} u, & \qquad x \in (0,1),\, \tau &\in (0,1) \\ u(x,0) &= f(x), & \qquad x &\in (0,1)\\ u(0,\tau) = u(L,\tau) &= 0, & \qquad \tau &\in (0,1), \end{alignedat}\tag{24}\] which we solve for the viscosity parameter \(\nu = 0.05\). In contrast to the previous example, we are only interested in the solution \(u(\cdot,1)\) at the final time. Again, given \(\boldsymbol{a} \in [-1,1]^9\), we consider initial conditions \(f\) given by the truncated sine expansion \[\begin{align} \label{eq:burger95initial95approx} f(x) \approx \sum_{\ell=1}^9 a_\ell \sin (\ell \pi x), & & x \in [0,1]. \end{align}\tag{25}\] As before, we transform the spatial variable via \(\mathcal{T} x = \frac{1}{2}(x+1)\) and use solve_ivp with the Radau method, this time with a relative tolerance of \(10^{-6}\) and an absolute tolerance of \({10^{-8}}\). In Algorithm 2, we also use sparsity \(s=1000\) and extension \(N=64\). Since the target function space is one-dimensional, OMP+ would use only \(3\) space-time evaluations per PDE solve in the recommended default configuration. To improve the efficiency, we increase this to \(9\), giving \(10\) evaluations per solve, consistent with the previous example. Regarding the results in [12], already sparsity \(s=500\) and extension \(N=32\) lead to runtimes on the order of several days, making larger configurations infeasible. We highlight this mismatch in \(s\) by denoting the comparison as cMR1L\(^\ast\). The TFNO baseline with \(81\,049\) parameters was trained on \(40\,000\) samples, resulting in \(40\,500\) total PDE solves including validation data.

Figure 5 shows the relative and absolute approximation errors for all approaches. The TFNO median error is almost two orders of magnitude smaller than those of the remaining approaches, which range between \(10^{-2}\) and \(10^{-1}\). The relative error \(\text{err}_{\text{rel}}(\boldsymbol{a})\) exhibits large upper whiskers for all approaches, which are absent for the absolute error \(\text{err}_{\text{abs}}(\boldsymbol{a})\).

Figure [fig:burger95runtime] shows the computation time and number of PDE solves for each approach. cMR1L\(^\ast\) required significantly more solves than the rest, resulting in substantially longer runtime. OMP and TFNO both finished in less than \(5\) hours, while OMP+ needed only about \(1\) hour.

Lastly, Figure 6 depicts the absolute approximation error \(\text{err}_{\text{abs}}(\boldsymbol{a})\) for OMP+ with varying numbers of additional spatial evaluations per PDE solve (\(3\), \(9\), and \(15\)) and oversampling factors (\(c \in \{1, 3\}\)). The error magnitude remains similar across all six variations, with accuracy improving slightly for more spatial evaluations or higher oversampling factor \(c\).

3.3 Parametric diffusion equation↩︎

As third example, we consider the 2D diffusion equation on \(\Omega = [0,1]^2\) with homogeneous Dirichlet boundary conditions and an affine random diffusion coefficient in some domain \(\Omega_{\boldsymbol{y}}\), see [50], [51] namely \[\label{eq:rand95diff95eq} \begin{align} -\nabla \cdot (a(\boldsymbol{x},\boldsymbol{y})\nabla u(\boldsymbol{x},\boldsymbol{y})) &= f(\boldsymbol{x}), \qquad & \boldsymbol{x} \in \Omega,\;\boldsymbol{y} \in \Omega_{\boldsymbol{y}},\\ u(\boldsymbol{x},\boldsymbol{y}) &= 0, \qquad & \boldsymbol{x} \in \partial\Omega,\;\boldsymbol{y} \in \Omega_{\boldsymbol{y}}. \end{align}\tag{26}\] Here, the differential operator \(\nabla\) acts with respect to the spatial variable \(\boldsymbol{x}\). As in [52] and [53], we use the affine coefficient \[\begin{align} a(\boldsymbol{x},\boldsymbol{y}) = 1 + \sum_{j=1}^{20} y_j \psi_j(\boldsymbol{x}), \qquad \boldsymbol{x} \in \Omega,\;\boldsymbol{y} \in [-1,1]^{20}, \end{align}\] with \(\psi_j(\boldsymbol{x}) = c j^{-\mu} \cos(2\pi m_1(j)x_1)\cos(2\pi m_2(j)x_2).\) The parameter are chosen as \(c = 0.9/\zeta(2)\) and \(\mu = 2\), together with \(k(j) \mathrel{\vcenter{:}}= \lfloor -1/2 + \sqrt{1/4 + 2j} \rfloor\), \(m_1(j) \mathrel{\vcenter{:}}= j-\frac{k(j) (k(j)+1)}{2}\) and \(m_2(j) \mathrel{\vcenter{:}}= k(j)-m_1(j)\). Moreover, we set \(f \equiv 1\). As before, we transform \(\boldsymbol{x}\) componentwise by the affine map \(\mathcal{T}(z)=\frac{1}{2}(z+1)\). In total, this yields a \(22\)-dimensional approximation problem.

As in [53], we solve the PDE 26 numerically using FEniCS on a mesh with \(5000\) finite elements. In Algorithm 2, we choose \(s=1000\) and \(N=64\). Moreover, we restrict the candidate set \(\Gamma\) by imposing \(\max_{\boldsymbol{k} \in \Gamma}\|\boldsymbol{k}\|_0 = 7\). Such a restriction reflects the common observation that, in many high-dimensional parametric PDEs, higher-order interactions contribute considerably less than lower-order ones. This choice is adopted from [12] to ensure a comparable experimental setup. The cMR1L comparison uses the same parameters. The TFNO with \(529\,233\) parameters is trained on \(85\,000\) samples, resulting in a total of \(86\,000\) PDE solves including validation data.

Figure 7 compares the relative approximation error \(\mathrm{err}(\boldsymbol{y})\) for all four approaches. The TFNO achieves a median error of roughly \(2.5\cdot10^{-4}\), almost an order of magnitude below the other approaches, whose median errors are around \(10^{-3}\). On the other hand, the TFNO’s upper whisker spans almost a full order of magnitude, compared to roughly half an order of magnitude for the other approaches.

Figure [fig:diff95pde95runtime] shows the number of PDE solves and corresponding runtimes for all approaches. As in the previous example, cMR1L required significantly more solves and runtime than OMP and OMP+, which finished in roughly \(8\) and \(2\) hours, respectively. Notably, the TFNO training time is comparable to the total runtime of cMR1L.

3.4 Discussion↩︎

Now, we discuss the numerical results and the observed trade-offs between approximation quality and computational effort.

3.4.0.1 Sparse approximation

A key observation is that the approximation accuracy of OMP and OMP+ depends primarily on the sparsity of the solution \(u\) in the chosen basis. For the heat equation, \(u\) is highly sparse in the multivariate Chebyshev basis due to the linear dependence on the parameters \(a_\ell\) in 23 . Consequently, the index set \(I_{\mathrm{est}}\) detected by Algorithm 2 already performs well for small sparsity \(s = |I_{\mathrm{est}}|\), and increasing the sparsity \(s\) further improves accuracy, as seen in Figure [fig:heat95OMP4395err]. The increase saturates once the diameter of the search space \(\Gamma\) (which is restricted by \(N\)) is exhausted.

The Burgers’ equation, by contrast, showed the opposite behavior: \(u\) is not well represented by sparse Chebyshev expansions, as discussed in [12]. This results in much lower approximation accuracy at the same sparsity \(s\). Moreover, increasing the amount of training data improves accuracy only marginally, as seen in Figure 6, suggesting that the observed limitations stem from the basis rather than insufficient data for the sparse recovery step.

These results highlight that the choice of basis critically affects the approximation quality. When the detected index set \(I_{\mathrm{est}}\) exhibits little sparsity structure, experimenting with alternative bases is advisable. Indeed, OMP and OMP+ are applicable for arbitrary product bases, provided that an efficient routine for the matrix-vector multiplication with the basis matrix 17 is available.

3.4.0.2 Samples and PDE solves

Applying sparse recovery methods such as OMP to the dimension-incremental framework yields a substantial advantage in sample complexity over cubature approaches based on rank-1 lattices. While the rank-1 lattice methods in [12] enable efficient reconstruction via FFT, this benefit is outweighed by the large number of required PDE solves, which scales as \(s^2\) in the worst case. In contrast, the number of samples required by OMP scales only linearly in \(s\). Furthermore, the additional samples per PDE solve in OMP+ (which are not usable for cubature-based methods) yield a substantial additional reduction in computational cost (by up to \(80\%\)) without sacrificing approximation quality, clearly visible across all experiments. As a consequence, our methods close the gap in terms of sample complexity to learning-based approaches such as TFNOs, though the dimension-incremental framework does require samples with a particular structure.

Moreover, we want to emphasize that the TFNO models serve purely as baselines, following the generic guidelines of the neural operator library [31] without further modification. The comparison should therefore not be interpreted as a rigorous benchmark, but rather as an indication of the approximation quality and computational complexity achievable by our methods relative to off-the-shelf neural operator configurations.

3.4.0.3 Interpretability of the approximation

An important advantage of the proposed approach is its interpretability. The detected index set \(I_{\mathrm{est}}\) reveals which input variables are most influential, which parameter interactions are relevant, and to what extent higher-order interactions matter. The resulting approximation is thus not merely a black-box predictor, but also yields qualitative insight into the underlying PDE. In contrast, such structural information is far less accessible in neural operator approaches such as TFNOs, where the learned representation is distributed across a large number of trainable parameters.

References↩︎