1 Introduction↩︎

Parametrized partial differential equations (PDEs) arise in optimization, control, inverse problems, and uncertainty quantification, where the same governing model must be solved for many parameter values [1]. High-fidelity finite element, finite difference, and spectral discretizations can make such many-query studies computationally expensive [2]–[4]. Similar bottlenecks arise in computational energy science, including subsurface flow, transport, and data-driven energy-materials modeling, where simulations and high-fidelity data generation often span multiple scales [5]. Reduced-order models (ROMs) reduce this cost by separating an offline basis-construction stage from a low-dimensional online solve. This strategy is effective when the solution manifold admits a compact approximation and the offline data sufficiently identify the relevant solution components [1], [6], [7]. When only a few high-fidelity snapshots are available, however, the reduced space may fail to capture components that are essential for out-of-sample prediction.

Reduced basis methods (RBMs) and proper orthogonal decomposition (POD) are two widely used approaches for constructing reduced spaces. RBMs select parameter samples adaptively, often through greedy algorithms guided by a posteriori error estimators [8]–[13]. In contrast, POD starts from a snapshot ensemble and computes an orthonormal basis that is optimal, in the sense of minimizing the mean-square projection error over the given snapshots. Introduced in the analysis of turbulent flows [14], [15], POD is closely related to principal component analysis, empirical orthogonal functions, and the Karhunen–Loeve expansion [16]–[18]. Its strong compression capability has made POD a standard basis-construction tool in computational fluid dynamics, porous-media flow, and related applications [19]–[26].

POD depends on the empirical distribution of the available snapshots. When the snapshots are sparse or unevenly distributed, the resulting basis may accurately reproduce the sampled states but fail to represent unsampled regions of the solution manifold. Adding more POD modes within the same snapshot span cannot remedy this deficiency if that span does not contain the solution components needed for prediction. This limitation contrasts with greedy RBMs, which use an error indicator to target parameter values that are poorly represented by the current reduced space [27]. Consequently, when snapshot data are insufficient, POD may incur substantial out-of-sample errors, especially for problems involving multiple scales, localized structures, high-frequency components, or strong parameter sensitivity.

Several POD variants incorporate additional structure into the modal decomposition. Spectral POD (SPOD) extracts frequency-resolved coherent structures in statistically stationary flows [28], [29], with applications to jet dynamics [30]–[32]. Multiscale POD (mPOD) combines POD with multiresolution analysis by decomposing snapshots into prescribed frequency bands before computing local POD modes [33]–[36]. These methods exploit spectral or multiscale information primarily to improve modal analysis and flow decomposition. In contrast, we use spectral subspaces to enrich POD-Galerkin reduced spaces for predictive reduced-order modeling of parametrized PDEs in limited-data regimes.

Machine-learning approaches provide another route for parametrized PDEs. Physics-informed neural networks (PINN) and Deep Ritz methods approximate PDE solutions from residual or variational formulations [37], [38], while Deep BSDE methods target high-dimensional parabolic problems [39]–[41]. Operator-learning approaches, such as DeepONet and Fourier neural operators, learn solution maps directly from data [42], [43]. Hybrid methods further integrate traditional numerical methods with data-driven modeling techniques, including random-feature PDE solvers [44], non-intrusive ROMs for flow problems [45], nonlinear reduced-basis methods with online neural adaptation [46], POD-enhanced deep ROMs [47], neural ROM identification [48], and PINN/RBM integrations [49], [50]. Although these approaches provide powerful alternatives, they often involve substantial training data, nonlinear optimization, or online adaptation. We instead focus on a complementary ROM setting: improving POD-Galerkin reduced bases in limited-snapshot regimes while retaining the transparency and projection structure of classical ROMs.

We propose a spectral-subspace-augmented POD-Galerkin method (SS-POD) for parametrized PDEs with limited snapshot data. SS-POD starts from a problem-adapted spectral approximation space and decomposes it into mutually orthogonal spectral subspaces. The snapshots are then projected onto these subspaces, POD bases are computed locally, and the resulting local bases are assembled into an augmented reduced basis. This construction mitigates the tendency of dominant high-energy snapshot components to obscure lower-energy components that may be relevant for out-of-sample prediction, while preserving the projection structure and interpretability of classical POD-Galerkin ROMs. An energy-balancing rule determines the spectral partition by assigning comparable amounts of snapshot energy to the local POD problems. We also extend the construction to nonlinear parametrized PDEs by coupling SS-POD with the discrete empirical interpolation method (DEIM) [51]–[54]. The main contributions of this work are: (i) a spectral-subspace-augmented reduced basis for data-scarce POD-Galerkin modeling; (ii) an energy-balancing rule for selecting the spectral partition; and (iii) a DEIM extension for nonlinear terms. We validate the proposed method through numerical comparisons with POD-Galerkin and spectral-Galerkin baselines on five PDE benchmarks.

The rest of the paper is organized as follows. Section 2 reviews POD-Galerkin ROMs and the offline-online decomposition. Section 3 presents the SS-POD construction, the energy-balancing rule, the nonlinear extension, and the reconstruction estimate. Section 4 presents numerical benchmark results. Section 5 concludes the paper and discusses limitations and possible extensions.

2 POD-Galerkin Method Revisited↩︎

In this section, we briefly review the POD-Galerkin method. For simplicity, consider a linear elliptic problem with homogeneous Dirichlet boundary conditions parameterized by \(\mu\), \[\label{eqn:general95gov95eqn} \mathcal{L}(\mu) u(\boldsymbol{x}; \mu) = f(\boldsymbol{x}; \mu), \quad \boldsymbol{x} \in \Omega,\tag{1}\] where \(\mu \in \mathbb{P}\) and \(\mathbb{P}\) denotes the parameter space. After spatial discretization, we write \(\mathbb{V}\) for the finite-dimensional trial space. The space \(\mathbb{V}\) is equipped with the inner product \(\left\langle \cdot,\cdot\right\rangle _{\mathbb{V}}\) and the induced norm \(\left\Vert \boldsymbol{v}\right\Vert _{\mathbb{V}}=\sqrt{\left\langle \boldsymbol{v},\boldsymbol{v}\right\rangle _{\mathbb{V}}}\). For the homogeneous Dirichlet case, \(\mathbb{V}\) is the discrete analogue of \(H_0^1(\Omega)\), so the boundary condition is encoded in the approximation space.

The full-order variational problem associated with Eq. 1 is written as: find \(\boldsymbol{u}(\mu)\in\mathbb{V}\) such that \[\label{eqn:variational95problem} A(\boldsymbol{u}(\mu), \boldsymbol{v}; \mu) = F(\boldsymbol{v}; \mu), \quad \forall \boldsymbol{v} \in \mathbb{V}_t,\tag{2}\] where \(A(\cdot,\cdot;\mu)\) and \(F(\cdot;\mu)\) denote the discrete bilinear and linear forms, \(\boldsymbol{u}(\mu)\) denotes the discrete solution in \(\mathbb{V}\), and \(\mathbb{V}_t\) the test space. In standard Galerkin projection, \(\mathbb{V}_t = \mathbb{V}\).

Many-query applications require this full-order problem to be solved for many values of \(\mu\). Solving it directly each time is usually too expensive. POD-Galerkin reduction therefore separates the computation into an offline stage, where a reduced space is learned, and an online stage, where only the reduced coefficients are solved.

In the offline stage, snapshots \(\boldsymbol{U}=[\boldsymbol{u}_1,\dots,\boldsymbol{u}_{N_{\mathrm{s}}}]\) are collected from full-order solutions, and POD is used to extract a reduced basis \(\boldsymbol{\Phi}=\boldsymbol{\Phi}_{\mathrm{POD}}\) from the snapshot matrix. In the online formulation below, we write \(\boldsymbol{\Phi}=[\boldsymbol{\phi}_1,\dots,\boldsymbol{\phi}_{N_{\mathrm{p}}}]\) (\(N_{\mathrm{p}} < N_{\mathrm{g}}\)) for a generic reduced basis, because the same formulation will later be used for SS-POD. The reduced approximation then seeks \(\boldsymbol{u}(\mu)\) in the span of these \(N_{\mathrm{p}}\) modes, \[\label{eqn:galerkin95approx} \boldsymbol{u}(\mu) \approx \hat{\boldsymbol{u}}(\mu) = \sum_{i=1}^{N_{\mathrm{p}}} a_i(\mu)\boldsymbol{\phi}_i = \boldsymbol{\Phi}\boldsymbol{a}(\mu),\tag{3}\] where \(\boldsymbol{a}(\mu)=[a_1(\mu),\dots,a_{N_{\mathrm{p}}}(\mu)]^\mathrm{T}\) contains the reduced coefficients.

Let \(\mathbb{V}_t\) be the reduced test space spanned by \(\tilde{\boldsymbol{\Phi}}=[\tilde{\boldsymbol{\phi}}_1,\dots,\tilde{\boldsymbol{\phi}}_{N_{\mathrm{p}}}]\). Substituting Eq. 3 into Eq. 2 and testing against \(\tilde{\boldsymbol{\phi}}_j\in\mathbb{V}_t\) gives the reduced system \[\sum_{i=1}^{N_{\mathrm{p}}} A(\boldsymbol{\phi}_i,\tilde{\boldsymbol{\phi}}_j;\mu)a_i(\mu) = F(\tilde{\boldsymbol{\phi}}_j;\mu), \quad j=1,\dots,N_{\mathrm{p}}\;.\] Note that for SS-POD, the test basis not necessary coincides with the trial basis \(\boldsymbol{\Phi}\); problem-specific choices of the test space will be discussed later and presented in the numerical examples. The projection step reduces the online solve to an \(N_{\mathrm{p}}\times N_{\mathrm{p}}\) linear system.

This model reduction strategy is efficient if two requirements are met. First, the solution manifold \(\mathbb{M}=\{\boldsymbol{u}(\mu)\mid \mu\in\mathbb{P}\}\) must be well approximated by a low-dimensional space. This property is often associated with rapid decay of the Kolmogorov \(N\)-width [6], [7]. Second, the reduced operators must be evaluated without incurring full-order cost. A standard assumption used to satisfy the second requirement is an affine parameter decomposition, \[\label{eqn:affine95decomposition95operator} \mathcal{L}(\mu) = \sum_{q=1}^{Q} B_q(\mu) \mathcal{L}_q,\tag{4}\] where \(\mathcal{L}_q\) are parameter-independent operators and \(B_q(\mu)\) are scalar functions of \(\mu\). Under this affine decomposition, the reduced matrices associated with \(\mathcal{L}_q\) can be precomputed offline and assembled online through the coefficients \(B_q(\mu)\). For nonlinear or non-affine operators, this separation generally requires an additional hyper-reduction step. In this work, such terms are handled by combining SS-POD with the discrete empirical interpolation method (DEIM), which will be discussed in Sec. 3.1.3.

3 The spectral-subspace-augmented POD-Galerkin (SS-POD) method↩︎

3.1 Construction of the SS-POD Methodology↩︎

3.1.1 Hierarchy of Projection Operators and Subspaces↩︎

POD extracts coherent structures from solution snapshots and provides the optimal low-rank approximation of the snapshot matrix in the chosen norm. Its success in ROMs is closely related to the Kolmogorov \(N\)-width of the solution manifold [1], [6], [55]. This connection, however, should be interpreted with care: POD is optimal for the available snapshot ensemble, not automatically for the full parametric manifold. When the training set is sparse or poorly distributed, the empirical snapshot span may miss directions that are important for out-of-sample parameters.

Spectral methods provide an a priori approximation space independent of the sampled parameters. Such spaces tolerate sparse training data better than empirical bases, but they may require more modes than a POD basis fitted to the target solution family. SS-POD combines these two ingredients by projecting snapshots onto orthogonal spectral subspaces and then performing POD locally. The spectral prior prevents dominant empirical directions from determining the whole reduced space, while the local POD steps keep data adaptation within each spectral window. Fig. 1 summarizes the construction.

Let \(\boldsymbol{\Phi}_{\mathrm{spec}}=[\boldsymbol{e}_1,\dots,\boldsymbol{e}_{N_{\mathrm{max}}}]\) be a collection of \(N_{\mathrm{max}}\) normalized spectral basis functions. We assume that \(N_{\mathrm{max}}\) is large enough for the truncated spectral space \(\mathbb{E}=\mathrm{span}\{\boldsymbol{e}_1,\dots,\boldsymbol{e}_{N_{\mathrm{max}}}\}\) to approximate the relevant solution states to the desired accuracy. Depending on the geometry and boundary conditions, the basis \(\{\boldsymbol{e}_j\}\) may be chosen from Fourier modes, Chebyshev polynomials, spherical harmonics, or other problem-adapted spectral functions. In the discrete setting, the modes are normalized so that \(\left\langle \boldsymbol{e}_i,\boldsymbol{e}_j\right\rangle _{\mathbb{V}}=\delta_{ij}\). The space \(\mathbb{E}\) is decomposed into orthogonal subspaces, \[\mathbb{E} = \mathrm{span}\{ \boldsymbol{e}_1, \dots, \boldsymbol{e}_{N_{\mathrm{max}}} \} = \mathbb{E}_1 \oplus \mathbb{E}_2 \oplus \dots \oplus \mathbb{E}_{N_{\mathrm{sub}}},\] where \(N_{\mathrm{sub}}\) is the number of spectral windows. Each subspace is spanned by a consecutive subset of spectral modes, \[\mathbb{E}_n = \mathrm{span} \{ \boldsymbol{e}_{k_n+1}, \dots, \boldsymbol{e}_{k_{n+1}} \},\] where the decomposition indices satisfy \(k_1=0\), \(k_{N_{\mathrm{sub}}+1}=N_{\mathrm{max}}\), and \(k_i<k_j\) for \(i<j\). SS-POD selects an energy-balanced decomposition index set from the spectral coefficients of the available snapshots: \[\mathbb{I}^* = \{ k_1^*, \dots, k_{N_{\mathrm{sub}}+1}^* \},\] The procedure for selecting \(\mathbb{I}^*\) is detailed in Sect. 3.1.2.

Given snapshots \(\boldsymbol{u}_i\) (\(i=1,\dots,N_{\mathrm s}\)) and the partition \(\mathbb{I}^*\), the data are projected onto each subspace \(\mathbb{E}_n\) by the \(\mathbb{V}\)-orthogonal projection \(\boldsymbol{\Pi}_n:\mathbb{R}^{N_{\mathrm g}}\to\mathbb{R}^{N_{\mathrm g}}\): \[\boldsymbol{u}_i^{(n)} = \boldsymbol{\Pi}_n \boldsymbol{u}_i = \sum_{j= k_n +1}^{k_{n+1}} \left\langle \boldsymbol{u}_i, \boldsymbol{e}_{j}\right\rangle _{\mathbb{V}} \boldsymbol{e}_{j}, \quad i=1, \dots, N_{\mathrm{s}}.\] POD is then applied separately to the projected snapshot matrices \(\boldsymbol{U}^{(n)}=[\boldsymbol{u}_1^{(n)},\dots,\boldsymbol{u}_{N_{\mathrm s}}^{(n)}]\), yielding local bases \(\boldsymbol{\Phi}^{(n)}=[\boldsymbol{\phi}^{(n)}_1,\dots,\boldsymbol{\phi}^{(n)}_{N_{\mathrm p}^{(n)}}]\). The SS-POD basis is obtained by concatenating the local bases: \[\boldsymbol{\Phi}_{\mathrm{Aug}} = [\boldsymbol{\Phi}^{(1)}, \dots, \boldsymbol{\Phi}^{(N_{\mathrm{sub}})}].\] The SS-POD basis set construction interpolates between two limiting cases: \[\label{ROSG95property} \boldsymbol{\Phi}_{\mathrm{Aug}} = \begin{cases} \boldsymbol{\Phi}_{\mathrm{POD}}, & \text{if } N_{\mathrm{sub}}=1 \text{ and } N_{\mathrm{p}}^{(1)} = N_{\mathrm{p}}, \\ \boldsymbol{\Phi}_{\mathrm{spec}}, & \text{if } N_{\mathrm{sub}}=N_{\mathrm{max}} \text{ and } N^{(n)}_{\mathrm{p}} = 1, \, n=1, \dots, N_{\mathrm{sub}}. \end{cases}\tag{5}\] The first case recovers standard POD. The second formally recovers the truncated spectral basis when every spectral window contains one active mode and the corresponding projected snapshot matrix has nonzero energy in that mode. Thus SS-POD can be viewed as a controlled interpolation between empirical POD and spectral Galerkin spaces. The parameters \(N_{\mathrm{sub}}\) and \(N_{\mathrm p}^{(n)}\) determine the tradeoff among accuracy, online dimension, and snapshot requirements. Algorithm 2 summarizes the construction of \(\boldsymbol{\Phi}_{\mathrm{Aug}}\).

3.1.2 Energy-balanced Determination of the Decomposition Index Set \(\mathbb{I}^*\)↩︎

SS-POD depends on the decomposition index set \(\mathbb{I}^*\). At the two extremes it recovers standard POD and spectral-Galerkin methods; between those extremes, the energy-balanced partition controls how much snapshot information enters each spectral window.

The choice of \(\mathbb{I}^*\) is based on the energy distribution of the snapshots in the spectral basis. Standard POD tends to prioritize directions with the largest empirical energy. For many PDE data sets, these directions are dominated by low-frequency or large-scale components, so lower-energy spectral ranges may receive too few basis vectors even when they are important for accuracy. SS-POD mitigates this imbalance by distributing the snapshot energy captured by the finite spectral prior approximately equally among the subspaces \(\mathbb{E}_n\). Applying POD to each projected matrix \(\boldsymbol{U}^{(n)}\) then gives each spectral window its own local approximation budget.

Specifically, the energy captured by the finite spectral prior, denoted by \(\mathscr{E}\), is evaluated by summing the squared spectral coefficients of the projected snapshots: \[\label{POD32energy} \begin{align} \mathscr{E} & =\sum_{i=1}^{N_{\mathrm{s}}} \left\Vert \sum_{n=1}^{N_{\mathrm{sub}}} \boldsymbol{\Pi}_n \boldsymbol{u}_i \right\Vert _{\mathbb{V}}^2 = \sum_{i=1}^{N_{\mathrm{s}}} \left\Vert \sum_{j=1}^{N_{\mathrm{max}}} \left\langle \boldsymbol{u}_i, \boldsymbol{e}_j\right\rangle _{\mathbb{V}} \boldsymbol{e}_j \right\Vert _{\mathbb{V}}^2 = \sum_{i=1}^{N_{\mathrm{s}}} \sum_{j=1}^{N_{\mathrm{max}}} \left\langle \boldsymbol{u}_i, \boldsymbol{e}_j\right\rangle _{\mathbb{V}}^2. \end{align}\tag{6}\] The target energy for each subspace is defined as \(\mathscr E_n=\mathscr E/N_{\mathrm{sub}}\). In the idealized case, the partition would satisfy \[\label{Eqn:energy95each95subspace} \sum_{i=1}^{N_{\mathrm{s}}} \sum_{j= k_n +1}^{k_{n+1}} \left\langle \boldsymbol{u}_i, \boldsymbol{e}_j\right\rangle _{\mathbb{V}}^2 = \mathscr{E}_n.\tag{7}\] Discrete partition indices usually prevent exact equality. Algorithm 4 constructs \(\mathbb{I}^*\) by balancing energy over the ordered spectral coefficients.

3.1.3 Extension to Nonlinear PDEs via SS-POD/DEIM↩︎

Although SS-POD has been introduced for linear projection spaces, the same subspace construction can be coupled with the Discrete Empirical Interpolation Method (DEIM) to handle non-affine and nonlinear terms [51]–[54], [56], [57]. DEIM is a special case of the Empirical Interpolation Method (EIM) that selects representative interpolation indices from nonlinear basis functions, thereby reducing the cost of evaluating nonlinear reduced-order models. Related hyper-reduction strategies, including reduced over-collocation, address the same full-grid nonlinear-evaluation bottleneck from a complementary sampling perspective [58].

In a projection-based ROM, the online stage should ideally be independent of the full-order dimension \(N_{\mathrm g}\). For affine linear operators this is achieved by offline precomputation. For nonlinear terms, direct evaluation usually requires pointwise operations on the full spatial grid, causing the online cost to scale with \(N_{\mathrm g}\). DEIM addresses this issue by approximating the nonlinear term as \[\label{DEIM-approximation} \boldsymbol{f}(\mu) \approx \boldsymbol{\Psi }\boldsymbol{c}(\mu),\tag{8}\] where \(\boldsymbol{\Psi}\) is the DEIM basis. Provided that \(\boldsymbol{P}^\mathrm T\boldsymbol{\Psi}\) is nonsingular, the coefficient vector \(\boldsymbol{c}(\mu)\) is determined by solving the square interpolation system \[\boldsymbol{P}^\mathrm{T} \boldsymbol{f}(\mu) = \boldsymbol{P}^\mathrm{T} \boldsymbol{\Psi }\boldsymbol{c}(\mu).\] The matrix \(\boldsymbol{P}\) selects the interpolation indices in the discretized nonlinear term. It is defined as \[\boldsymbol{P} = [\boldsymbol{i}_{I_1}, \boldsymbol{i}_{I_2}, \dots, \boldsymbol{i}_{I_{N_{\mathrm{d}}}}],\] where \(\boldsymbol{i}_{I_k}\in\mathbb{R}^{N_{\mathrm g}}\) is the \(I_k\)-th canonical basis vector.

The DEIM basis is constructed from nonlinear snapshots \(\boldsymbol{F}=[\boldsymbol{f}_1,\dots,\boldsymbol{f}_{N_{\mathrm s}}]\). POD is applied to \(\boldsymbol{F}\), and the first \(N_{\mathrm d}\) retained modes \([\boldsymbol{\psi}_1,\dots,\boldsymbol{\psi}_{N_{\mathrm d}}]\) form \(\boldsymbol{\Psi}\). The interpolation indices \(\{I_k\}_{k=1}^{N_{\mathrm d}}\) are selected iteratively by the standard DEIM residual criterion: \[I_k = \mathop{\mathrm{argmax}}_{i=1, \dots, N_{\mathrm{g}}} |\boldsymbol{i}_i^\mathrm{T}(\boldsymbol{\psi}_k - \boldsymbol{\Psi}^{(k-1)} \boldsymbol{c}^{(k)})|, \quad k=1, \dots, N_{\mathrm{d}},\] where \(\boldsymbol{\Psi}^{(k)} = [\boldsymbol{\psi}_1, \dots, \boldsymbol{\psi}_k]\), and \(\boldsymbol{c}^{(k)}\) is the solution to: \[(\boldsymbol{P}^{(k-1)})^\mathrm{T} \boldsymbol{\Psi}^{(k-1)} \boldsymbol{c}^{(k)} = (\boldsymbol{P}^{(k-1)})^\mathrm{T} \boldsymbol{\psi}_k, \quad k=1, \dots, N_{\mathrm{d}}.\] Here, \(\boldsymbol{P}^{(k)}\) contains the first \(k\) selected columns of \(\boldsymbol{P}\): \[\boldsymbol{P}^{(k)} = [\boldsymbol{i}_{I_1}, \boldsymbol{i}_{I_2}, \dots, \boldsymbol{i}_{I_k}], \quad k=1, \dots, N_{\mathrm{d}}.\] The resulting DEIM approximation is \[\label{eqn:DEIM95approximation} \boldsymbol{f}(\mu) \approx \boldsymbol{\Psi }(\boldsymbol{P}^\mathrm{T} \boldsymbol{\Psi})^{-1} \boldsymbol{P}^\mathrm{T} \boldsymbol{f}(\mu).\tag{9}\] The term \(\boldsymbol{\Psi }(\boldsymbol{P}^\mathrm{T} \boldsymbol{\Psi})^{-1}\) can be precomputed during the offline stage. Consequently, evaluating the nonlinear term in the reduced equations only requires the entries selected by \(\boldsymbol{P}^\mathrm{T} \boldsymbol{f}(\mu)\); this removes the full-grid nonlinear evaluation from the online reduced solve, although full-state reconstruction or postprocessing still scales with \(N_{\mathrm{g}}\). The detailed DEIM procedure is summarized in Algorithm 5.

The DEIM basis is also snapshot-dependent, so standard POD/DEIM can inherit the same data-scarcity limitations as the solution basis. We apply SS-POD to the nonlinear snapshots as well. The augmented index set \(\mathbb{I}^*_{\mathrm{D}}\) follows the same energy-based construction as \(\mathbb{I}^*\), with solution snapshots \(\boldsymbol{U}\) replaced by nonlinear snapshots \(\boldsymbol{F}\). Algorithm 6 gives the resulting SS-POD/DEIM procedure.

3.2 Error and Complexity Analysis↩︎

3.2.1 Snapshot Reconstruction Error↩︎

3.2.2 Discussions about the Online Approximation Error↩︎

The snapshot reconstruction estimate above does not by itself guarantee accuracy for a new parameter value \(\mu\) outside the training set. For an out-of-sample solution \(\boldsymbol{u}=\boldsymbol{u}(\mu)\), the SS-POD approximation error can be decomposed across the spectral subspaces as \[\begin{align} \left\Vert \boldsymbol{u} - \boldsymbol{\Phi }\boldsymbol{a} \right\Vert &= \left\Vert \sum_{n=1}^{N_{\mathrm{sub}}} \boldsymbol{u}^{(n)} - \sum_{n=1}^{N_{\mathrm{sub}}} \boldsymbol{\Phi}^{(n)} \boldsymbol{a}^{(n)} \right\Vert \\ &\le \sum_{n=1}^{N_{\mathrm{sub}}} \left\Vert \boldsymbol{u}^{(n)} - \boldsymbol{\Phi}^{(n)} \boldsymbol{a}^{(n)} \right\Vert . \end{align}\] This relation is an error decomposition, not a closed a priori bound. The terms on the right depend on how well each local POD space, learned from projected training snapshots, represents the corresponding component of the new solution. The energy-balanced partition reduces the risk that a single global POD step discards low-energy spectral components. It does not remove the usual accuracy requirements: the spectral prior must approximate the relevant solution class, the snapshot set must contain representative parameter information, and the Galerkin or DEIM online solve must be accurate enough for the target tolerance.

3.2.3 Complexity Analysis and Parallel Computing for the Offline Procedure↩︎

Let \(N_{\mathrm{basis}}\) denote the total number of reduced basis functions used by a given method. For SS-POD, \(N_{\mathrm{basis}}=\sum_{n=1}^{N_{\mathrm{sub}}}N_{\mathrm p}^{(n)}\); for standard POD and spectral-Galerkin methods, it is \(N_{\mathrm p}\) and \(N_{\mathrm{max}}\), respectively.

The main computational costs are summarized below. The estimates are intended as order-of-magnitude costs for the implementations used in this work; they do not account for all possible sparsity, tensor-product, or fast-transform optimizations.

• For some properly chosen Fourier or Chebyshev bases, the spectral coefficients can be computed with fast transforms, with cost \(\mathcal{O}(N_{\mathrm{s}}N_{\mathrm{g}}\log N_{\mathrm{g}})\) when the discretization supports such transforms.

• The identification of the energy-balanced partition \(\mathbb{I}^*\) has a complexity of \(\mathcal{O}(N_{\mathrm{s}} N_{\mathrm{g}})\).

• Computing the retained POD modes by truncated SVD costs approximately \(\mathcal{O}(N_{\mathrm{s}}N_{\mathrm{g}}N_{\mathrm{basis}})\), up to implementation-dependent constants and ignoring discarded near-null modes.

• The precomputation of reduced mass and stiffness matrices scales with two transform-based terms: \(\mathcal{O}(N_{\mathrm{basis}}^2N_{\mathrm{g}})\) for reduced matrix contractions and \(\mathcal{O}(N_{\mathrm{basis}}N_{\mathrm{g}}\log N_{\mathrm{g}})\) for transform operations. Some spectral-Galerkin matrix entries can instead be evaluated analytically.

• For \(N_\mu\) new parameter values, online affine assembly scales as \(\mathcal{O}(N_\mu N_{\mathrm{basis}}^2)\), with up to \(\mathcal{O}(N_\mu N_{\mathrm{basis}}^3)\) additional cost for dense direct reduced solves.

• Reconstructing full-order solution fields from reduced coefficients requires \(\mathcal{O}(N_{\mathrm{g}}N_{\mu}N_{\mathrm{basis}})\) operations.

The local POD computations for the projected matrices \(\boldsymbol{U}^{(n)}\) are independent across subspaces and can be parallelized. We implemented this offline step in MATLAB with the Parallel Computing Toolbox.

3.3 Implementation Details and Practical Considerations↩︎

3.3.1 Boundary Conditions: Dirichlet BC and Choice of Test Functions↩︎

Boundary conditions require special care in projection-based ROMs because a basis constructed from snapshots or spectral modes does not automatically satisfy the desired test-space constraints. For inhomogeneous Dirichlet conditions, standard remedies include the control function method [59], [60] and modified basis methods [60]. In Chebyshev-Galerkin discretizations, boundary-compatible trial and test functions can be constructed following Shen [61] or Heinrichs [62]. In the Dirichlet examples below, we project the trial basis functions onto the homogeneous boundary-compatible Chebyshev subspace. Let \(\boldsymbol{\chi}_j=\boldsymbol{T}_j-\boldsymbol{T}_{j-2}\) for \(j=2,\ldots,N_{\mathrm{max}}\). The test function associated with \(\boldsymbol{\phi}_{\mathrm{Aug},i}\) is written as \[\label{eqn:test95function} \tilde{\boldsymbol{\phi}}_{\mathrm{Aug},i} = \sum_{j=2}^{N_{\mathrm{max}}} c_{ij}\boldsymbol{\chi}_j,\tag{10}\] where the coefficients solve the Gram system \[\label{eqn:test95function95gram} \sum_{\ell=2}^{N_{\mathrm{max}}} \left\langle \boldsymbol{\chi}_\ell,\boldsymbol{\chi}_j\right\rangle _{\mathbb{V}} c_{i\ell} = \left\langle \boldsymbol{\phi}_{\mathrm{Aug},i},\boldsymbol{\chi}_j\right\rangle _{\mathbb{V}}, \qquad j=2,\ldots,N_{\mathrm{max}}.\tag{11}\] If \(\{\boldsymbol{\chi}_j\}\) is first orthonormalized in the \(\mathbb{V}\)-inner product, this reduces to the simple coefficient formula given by direct inner products. The resulting \(\tilde{\boldsymbol{\phi}}_{\mathrm{Aug},i}\) satisfies the homogeneous boundary condition by construction. This test-space construction can be interpreted as a projection of the SS-POD trial basis onto \[\mathrm{span} \{ \boldsymbol{T}_j - \boldsymbol{T}_{j-2} \}_{j=2}^{N_{\mathrm{max}}}.\] For a 1D problem with Dirichlet boundary conditions, this construction leaves two boundary degrees of freedom in the reduced coefficient system. We close the system by imposing the two boundary constraints \[\sum_i a_i(\mu)\boldsymbol{\phi}_{\mathrm{Aug},i}(x)\big|_{x=\pm1}=0.\]

3.3.2 Threshold for Spurious Basis Functions↩︎

Truncation thresholds are used to remove numerically insignificant modes and to keep \(\sum_{n=1}^{N_{\mathrm{sub}}}N_{\mathrm p}^{(n)}\le N_{\mathrm{max}}\). Modes associated with very small singular values contribute little to the snapshot reconstruction and can worsen conditioning in the reduced systems. We use \(\mathrm{tol}_{\mathrm p}\) for standard POD and \(\mathrm{tol}_{\mathrm p}^{(n)}\) for the local SS-POD truncations; only modes whose singular values exceed the prescribed tolerance are retained. Analogously, thresholds \(\mathrm{tol}_{\mathrm d}\) and \(\mathrm{tol}_{\mathrm d}^{(n)}\) are used for DEIM and SS-POD/DEIM. All the parameter values used in the reported numerical experiments are listed in Table 1.

4 Numerical Results↩︎

4.1 Technical Details↩︎

4.1.1 Notations and Parameter Setup↩︎

Let \(\boldsymbol{U}_{\mathrm{test}}=[\boldsymbol{u}(\mu_1),\dots,\boldsymbol{u}(\mu_{N_\mu})]\) collect the full-order test solutions, and let \(\hat{\boldsymbol{U}}_{\mathrm{test}}\) be the corresponding ROM approximation. The error matrix and relative error, measured in the vectorized \(\ell^2\) norm equivalently the Frobenius norm, are \[\begin{align} \boldsymbol{\epsilon }&= \boldsymbol{U}_{\mathrm{test}} - \hat{\boldsymbol{U}}_{\mathrm{test}}, \\ \epsilon^{\mathrm{r}} &= \frac{\left\Vert \boldsymbol{\epsilon}\right\Vert _{\mathrm{F}}}{\left\Vert \boldsymbol{U}_{\mathrm{test}}\right\Vert _{\mathrm{F}}}. \end{align}\] When POD and SS-POD errors are reported separately, we use \(\boldsymbol{\epsilon}_{\mathrm P}\) and \(\boldsymbol{\epsilon}_{\mathrm S}\), respectively.

For nonlinear examples, we also report the DEIM approximation error for the nonlinear snapshot matrix \(\boldsymbol{F}_{\mathrm{test}}=[\boldsymbol{f}(\mu_1),\dots,\boldsymbol{f}(\mu_{N_\mu})]\): \[\epsilon_{\mathrm{D}}^{\mathrm{r}} = \frac{\left\Vert \boldsymbol{\Psi }(\boldsymbol{P}^\mathrm{T} \boldsymbol{\Psi})^{-1} \boldsymbol{P}^\mathrm{T} \boldsymbol{F}_{\mathrm{test}} - \boldsymbol{F}_{\mathrm{test}}\right\Vert _{\mathrm{F}}}{\left\Vert \boldsymbol{F}_{\mathrm{test}}\right\Vert _{\mathrm{F}}},\] The number of reduced basis functions used for the solution approximation is denoted by \(N_{\mathrm{basis}}\), as in Sect. 3.2.3. For nonlinear terms approximated by POD/DEIM or SS-POD/DEIM, the corresponding DEIM basis size is denoted by \(N_{\mathrm{basis}}^{\mathrm D}\).

For nonlinear time-dependent Galerkin systems, the reduced coefficients are obtained by solving nonlinear algebraic systems at each time step. We use the Levenberg-Marquardt algorithm, initializing \(\boldsymbol{a}^{t+1}\) with \(\boldsymbol{a}^t\), and set the stopping tolerances to \(\mathrm{tol}_{\mathrm{fun}}=10^{-11}\) and \(\mathrm{tol}_{\mathrm{grad}}=10^{-11}\).

The MATLAB function svds is used for truncated SVD computations. Since the snapshot counts are small in the reported data-scarce tests, the requested SVD rank is initially set to \(N_{\mathrm s}\), and modes with negligible singular values are subsequently removed by the tolerances in Table 1. For the nonlinear examples, we set \(N_{\mathrm{sub}}^{\mathrm D}=N_{\mathrm{sub}}\) as an empirical choice so that the nonlinear basis receives the same spectral-subspace resolution as the solution basis.

In randomized tests, snapshots are selected using MATLAB’s default random number generator. In the sphere example, the spherical harmonic transform library of [63] is used for spherical transforms and numerical integration.

4.1.2 Temporal Discretization Scheme↩︎

For time-dependent problems, the temporal interval is discretized uniformly with step size \(\Delta t=t_{i+1}-t_i\). Because the purpose of the experiments is to compare spatial reduced bases, we use sufficiently stable temporal discretizations to limit the influence of time-stepping error. Unless otherwise stated, the Galerkin system is discretized with the second-order Adams-Moulton method, i.e., the trapezoidal rule, which is A-stable for linear autonomous problems. For the 2D Allen-Cahn example, we additionally use an operator splitting scheme, since this is a standard and efficient treatment of the reaction and diffusion components.

4.2 Linear, Time-independent Cases: 2D Poisson and Helmholtz Equations↩︎

4.2.1 2D Poisson Equation↩︎

We first consider a 2D Poisson equation with a parametrized multiscale source: \[\label{eqn:poisson} \left \{ \begin{array}{ll} -\Delta u(x,y;a_k) = f(x,y;a_k) & (x,y) \in \Omega \setminus \partial\Omega \\ u(x,y;a_k) = 0 & (x,y) \in \partial\Omega \end{array} \right.\tag{12}\] where \(\Omega=[-1,1]^2\) and \[f(x,y;a_k) = \prod_{k=1}^{N} (1+\frac{1}{2}\cos(a_k \pi (x+y)))(1+\frac{1}{2}\sin(a_k \pi (x-3y))).\] The parameter \(a_k\) varies in \([2^{k-2},1.5\times2^{k-2}]\), with \(N=2\). Fig. 7 shows the source, the reference solution, and the error distributions for a limited-snapshot test. With \(N_{\mathrm s}=15\) snapshots, standard POD gives a larger relative error than SS-POD for this test.

For SS-POD, we use tensor-product Chebyshev functions \[\label{eqn:2d95cheb} T_{i,j}(x,y) = T_{i}(x) T_{j}(y).\tag{13}\] The two-dimensional indices \((i,j)\) are ordered by \(k=\sqrt{i^2+j^2}\) so that the energy-balanced partition can be applied to a one-dimensional spectral ordering. The Galerkin projection yields \[\label{eqn:gov95poisson} -\left(\left\langle \tilde{\boldsymbol{\Phi}}, \boldsymbol{\Phi}_{xx}\right\rangle _{\mathbb{V}} + \left\langle \tilde{\boldsymbol{\Phi}}, \boldsymbol{\Phi}_{yy}\right\rangle _{\mathbb{V}}\right) \boldsymbol{a}(\mu) = \left\langle \tilde{\boldsymbol{\Phi}}, \boldsymbol{f}(\mu)\right\rangle _{\mathbb{V}},\tag{14}\] subject to the boundary condition \(\boldsymbol{\Phi}(I^{\mathrm{b}},:) \boldsymbol{a}(\mu) = \boldsymbol{0}\).

Fig. 8 (a) compares the three reduced spaces. The Chebyshev-Galerkin curve gives the data-independent spectral baseline. Standard POD decreases the error at small basis sizes and then saturates in this limited-snapshot setting. SS-POD uses the spectral prior together with the snapshots and reaches \(\epsilon_\mathrm{S}^\mathrm{r}=5.92\times 10^{-7}\) for the reported configuration.

4.2.2 Helmholtz Equation↩︎

The second elliptic benchmark is the 2D Helmholtz equation \[\label{eq:Helmholtz} \left \{ \begin{array}{ll} \Delta u + k^2 u = f(x,y) & (x,y) \in \Omega \setminus \partial\Omega \\ u(x,y) = 0 & (x,y) \in \partial\Omega \end{array} \right.\tag{15}\] with wavenumber \(k\in[8,10]\) and source \(f(x,y)=\exp(-10[(y-1)^2+(x-0.5)^2])\). The interval includes values close to eigenfrequencies of the homogeneous Dirichlet problem, for example \(k=\frac{\pi}{2}\sqrt{3^2+5^2}\approx9.1592\) [64]. The reference solution at \(k=9.1515\), shown in Fig. 9, is close to the \((3,5)\) mode and therefore provides a more oscillatory test than the Poisson example.

The corresponding Galerkin ROM is \[(\left\langle \tilde{\boldsymbol{\Phi}}, \boldsymbol{\Phi}_{xx}\right\rangle _{\mathbb{V}} + \left\langle \tilde{\boldsymbol{\Phi}}, \boldsymbol{\Phi}_{yy}\right\rangle _{\mathbb{V}} + k^2\left\langle \tilde{\boldsymbol{\Phi}}, \boldsymbol{\Phi}\right\rangle _{\mathbb{V}}) \boldsymbol{a}(k) = \left\langle \tilde{\boldsymbol{\Phi}}, \boldsymbol{f}\right\rangle _{\mathbb{V}}.\]

Fig. 8 (b) shows larger error variability for standard POD in the near-resonant regime with limited snapshots (\(N_{\mathrm{s}}=15,17\)). SS-POD reduces this variability in the tested configurations by separating the snapshot information across spectral subspaces.

4.3 Linear, Time-dependent Case: 1D Heat Equation↩︎

We next solve a 1D heat equation with a multiscale source: \[\label{eq:heat} \left \{ \begin{array}{ll} \partial_t u(x,t) = \alpha \nabla^2 u(x,t) + f(x; a_k) & x \in (-1, 1), t \in [0,1] \\ u(\pm 1, t) = 0, \quad u(x,0) = 0 & \end{array} \right.\tag{16}\] where \(f(x;a_k)=\prod_{k=1}^{N}(1+\frac{1}{2}\cos(a_k\pi x))(1+\frac{1}{2}\sin(a_k\pi x))\times10^{-3}\) with \(N=5\). The parameters are sampled from \([2^{k-1},1.5\times2^{k-1}]\). In this example the reduced basis is learned from temporal snapshots of a single reference trajectory, so data scarcity refers to sparse sampling in time rather than sparse sampling of a parameter grid. Fig. 10 illustrates the transient solution and the corresponding error maps.

The temporal discretization of the reduced system leads to \[\begin{align} &\left(\left\langle \tilde{\boldsymbol{\Phi}}, \boldsymbol{\Phi}\right\rangle _{\mathbb{V}} - \alpha \frac{\Delta t}{2}\left\langle \tilde{\boldsymbol{\Phi}}, \boldsymbol{\Phi}_{xx}\right\rangle _{\mathbb{V}}\right)\boldsymbol{a}^{t+1} \\ &\qquad = \left(\left\langle \tilde{\boldsymbol{\Phi}}, \boldsymbol{\Phi}\right\rangle _{\mathbb{V}} + \alpha \frac{\Delta t}{2}\left\langle \tilde{\boldsymbol{\Phi}}, \boldsymbol{\Phi}_{xx}\right\rangle _{\mathbb{V}}\right)\boldsymbol{a}^t + \Delta t\,\left\langle \tilde{\boldsymbol{\Phi}}, \boldsymbol{f}\right\rangle _{\mathbb{V}} . \end{align}\]

Fig. 8 (c) reports the case with \(N_{\mathrm{s}}=4\) temporal snapshots. Increasing \(N_{\mathrm{sub}}\) reduces the large errors caused by temporal under-sampling for this multiscale transient.

4.4 Nonlinear Time-dependent Case: 1D Allen Cahn Equation and SS-POD/DEIM↩︎

We now consider the 1D Allen-Cahn equation \[\label{eqn:reaction32diffusion} \left \{ \begin{array}{ll} \partial_t u = \frac{1}{L^2} \nabla^2 u + u - u^3, & x\in (-1,+\infty), \;t \in [0,2], \\ u(-1,t) = 0, \quad \lim_{x \to \infty} u(x,t) = 1,& t \in [0,2],\\ u(x,0) = u_0(x), & x\in[-1,+\infty). \end{array} \right.\tag{17}\] with \(L=30\) and \(u_0(x)=\frac{a(x)-b(x)}{a(x)+b(x)+C}\), where \(a(x)=\exp(\frac{\sqrt{2}L}{2}x+\frac{\sqrt{2}L}{2})\), \(b(x)=\exp(-\frac{\sqrt{2}L}{2}x-\frac{\sqrt{2}L}{2})\), and \(C=10^4\). The reference solution is \(u(x,t)=\frac{a(x)-b(x)}{a(x)+b(x)+C\exp(-1.5t)}\). In the numerical implementation, the semi-infinite spatial domain is truncated to \([-1,1]\).

The Galerkin reduced system is \[\left\langle \tilde{\boldsymbol{\Phi}}, \boldsymbol{\Phi}\right\rangle _{\mathbb{V}} \frac{d \boldsymbol{a}(t)}{d t} = \left\langle \tilde{\boldsymbol{\Phi}}, \frac{\boldsymbol{\Phi}_{xx}}{L^2} + \boldsymbol{\Phi}\right\rangle _{\mathbb{V}} \boldsymbol{a}(t) - \left\langle \tilde{\boldsymbol{\Phi}}, \boldsymbol{f}(\boldsymbol{a})\right\rangle _{\mathbb{V}}.\] The cubic nonlinear term is approximated by SS-POD/DEIM: \[\label{eqn:DEIM95appro95nonlinear} \boldsymbol{f}(\boldsymbol{a}) \approx \boldsymbol{\Psi}(\boldsymbol{P}^\mathrm{T} \boldsymbol{\Psi})^{-1} (\boldsymbol{P}^\mathrm{T} \boldsymbol{\Phi }\boldsymbol{a} \odot \boldsymbol{P}^\mathrm{T} \boldsymbol{\Phi }\boldsymbol{a} \odot \boldsymbol{P}^\mathrm{T} \boldsymbol{\Phi }\boldsymbol{a}),\tag{18}\] where \(\odot\) denotes the Hadamard product. Fig. 8 (d) and (e) report the data-scarce cases (\(N_{\mathrm{s}}=3,5\)). SS-POD/DEIM preserves the DEIM-style reduced nonlinear evaluation and gives a lower relative error than the standard POD/DEIM baseline in these tests.

4.5 Nonlinear, Time-dependent Case: 2D Allen-Cahn Equation↩︎

The 2D Allen-Cahn equation reads: \[\label{eq:AC} \partial_t u = \epsilon^2 \Delta u + u - u^3, \quad (x,y) \in \Omega = [-1,1]^2.\tag{19}\] Periodic boundary conditions are imposed. The initial condition is generated from a Fourier series with frequencies below 4 and Gaussian random coefficients. The simulation runs to \(t=30\), from a random initial state toward a coarsened phase pattern. We use the following first-order operator splitting scheme.

Step 1: In \([t, t + \Delta t]\), solve the nonlinear reaction part, \[\begin{align} \boldsymbol{u}_1 &\longleftarrow \boldsymbol{u}^t, \\ \boldsymbol{u}_2 &= \frac{\boldsymbol{u}_1}{\sqrt{e^{-2 \Delta t} + (1 - e^{-2 \Delta t})(\boldsymbol{u}_1)^2}}. \end{align}\] Step 2: In \([t, t + \Delta t]\), solve the following linear equation with initial state \(\boldsymbol{u}_2\), \[\begin{align} \frac{d \boldsymbol{u}_3}{dt} &= \epsilon^2 \Delta \boldsymbol{u}_3, \\ \boldsymbol{u}^{t+1} &\longleftarrow \boldsymbol{u}_3. \end{align}\]

The reduced-order approximation is applied to Step 2: \[\left\langle \tilde{\boldsymbol{\Phi}}, \frac{\boldsymbol{\Phi }\boldsymbol{a}_3 - \boldsymbol{\Phi }\boldsymbol{a}_2}{\Delta t}\right\rangle _{\mathbb{V}} = \left\langle \tilde{\boldsymbol{\Phi}}, \epsilon^2 \Delta \boldsymbol{\Phi }\boldsymbol{a}_3\right\rangle _{\mathbb{V}}.\]

The error of this operator-splitting scheme is \(\mathcal{O}(\Delta t)\). If needed, other schemes such as the second-order Strang splitting scheme [65] can be used as well and combined with SS-POD. We use the first-order splitting throughout the comparison here.

Fig. 8 (f) shows larger error during rapid coarsening (\(t=9\)s) than near the later state (\(t=29.9\)s). Fig. 12 shows the interface evolution and the corresponding error distribution under limited snapshots.

4.6 Laplace-Beltrami Problem on a Sphere and Point-wise Convergence↩︎

Finally, we test SS-POD on a non-Euclidean geometry by considering a parametrized surface elliptic problem on the unit sphere: \[\label{eqn:laplace95beltrami} \Delta_{3S} u(x,y,z;\alpha) = f(x,y,z;\alpha), \quad (x,y,z) \in \mathbb{S}^2,\tag{20}\] where \(\mathbb{S}^2=\{(x,y,z)\in\mathbb{R}^3:x^2+y^2+z^2=1\}\). The source term is constructed from a localized oscillatory profile adapted from spherical benchmark functions [66], with an additional parameter-dependent modulation. It is given by \[f(x,y,z;\alpha) = F''(g) (|\vec{k}|^2 - g^2) - 2 F'(g) g.\] The surface Laplacian is written in spherical coordinates \((\theta,\phi)\) as \[\Delta_{3S} = \frac{1}{\sin \theta}\frac{\partial}{\partial \theta}\left( \sin \theta \frac{\partial}{\partial \theta}\right) + \frac{1}{\sin^2 \theta} \frac{\partial^2}{\partial \phi^2}.\] The auxiliary function \(F(g)\) and its derivatives are \[\begin{align} F(g) &= \cos(g) e^{-\sigma g^2}, \\ F'(g) &= e^{-\sigma g^2}(-\sin g - 2 \sigma g \cos g), \\ F''(g) &= e^{-\sigma g^2} [4\sigma^2 g^2 \cos g - \cos g - 2\sigma \cos g + 4\sigma g \sin g]. \end{align}\] Here \(g=\vec{k}\cdot\vec{r}\), with \(\vec{r}=(x,y,z)\) and \(\vec{k}=(3,3,\alpha)\). The parameter \(\alpha\in[1,4]\) changes the orientation and local frequency of the wave packet, and the smoothing parameter is fixed at \(\sigma=0.005\). The corresponding analytical reference solution is \[u(x,y,z;\alpha) = \cos(g) \exp(-\sigma g^2).\] The resulting solution family contains localized gradients and parameter-sensitive oscillations. Representative solutions are shown in the first row of Fig. 13.

We use spherical harmonics as the spectral prior: \[Y_{l,m}(\theta, \phi) = (-1)^{m} \sqrt{\frac{2l+1}{4\pi}\frac{(l-m)!}{(l+m)!}} P_l^m(\cos \theta) e^{im\phi},\] where \(P_l^m\) are the associated Legendre polynomials. Because the zero mode \(Y_{0,0}\) lies in the nullspace of the Laplace-Beltrami operator, the solution is determined only up to an additive constant unless a mean constraint is imposed. The reported errors are therefore computed after removing the \(Y_{0,0}\) component. For SS-POD, the two-dimensional harmonic index \((l,m)\) is reordered into a one-dimensional sequence by increasing degree \(l\).

Fig. 13 shows the error distribution for \(N_{\mathrm{s}}=5\). Standard POD yields a relative error of \(1.85 \times 10^{-4}\) with 5 basis functions, whereas SS-POD gives \(2.14 \times 10^{-6}\) with 83 basis functions. The POD error concentrates near sharper solution variation; the SS-POD error is lower and more evenly distributed in this test.

Fig. 14 gives the parameter-space error distribution. The left panel reports global convergence, and the right panel shows point-wise errors over \(\alpha \in [1,4]\) for \(N_{\mathrm{s}}=5\). For \(N_{\mathrm{sub}}=1\), the error is smallest near training snapshot locations and larger at unobserved parameter values. Increasing \(N_{\mathrm{sub}}\) produces a more uniform error profile in this example.

5 Conclusions and Outlook↩︎

This paper introduced SS-POD, a reduced-basis construction that combines POD with a problem-adapted spectral prior. The method decomposes a finite spectral approximation space into orthogonal subspaces selected by an energy-balanced partition, projects the snapshots onto these subspaces, and performs POD locally before assembling the augmented reduced basis. The SS-POD/DEIM extension applies the same construction to nonlinear snapshots.

Across the reported benchmarks, SS-POD gives lower limited-snapshot errors than standard POD in the tested configurations and uses fewer modes than a purely spectral approximation in several cases. The examples include elliptic, parabolic, nonlinear time-dependent, and spherical problems. The common pattern is that a structured spectral prior reduces the sensitivity of POD-Galerkin ROMs to sparse or uneven snapshot information.

SS-POD remains a linear-subspace ROM, so its performance is still constrained by the intrinsic approximability of the solution manifold. For transport-dominated problems, moving discontinuities, sharp traveling layers, or other regimes with slowly decaying Kolmogorov \(N\)-width, a compact linear trial space may be insufficient. Such problems likely require localization, nonlinear approximation, or adaptive basis strategies in addition to SS-POD.

SS-POD also requires a useful spectral prior. Fourier, Chebyshev, and spherical harmonic bases fit tensor-product domains and the sphere, but irregular geometries may not admit a convenient global spectral basis. Domain decomposition is a plausible extension: one could partition the physical domain into subregions with local spectral representations, following spectral element methods [67], [68] and localized reduced basis approaches [69]. This extension is not developed here; it would require separate analysis of interface treatment, local basis coupling, and offline cost. Another route is to learn POD-compatible subspaces on complex geometries: recent neural subspace POD work uses DeepONet-learned POD subspaces to support Krylov solves on unstructured meshes and CAD-derived domains [70]. This direction addresses a limitation that SS-POD does not resolve, namely the need for a convenient spectral prior on the physical domain.

Future work should test snapshot selection more carefully. The present experiments use random snapshot selection in several tests; systematic sampling, adaptive parameter exploration, or RBF-assisted interpolation may improve parameter-space coverage. For transport-dominated regimes, SS-POD may be combined with Principal Interval Decomposition [71] or with adaptive and hybrid ROM strategies developed to mitigate the Kolmogorov barrier in multiscale kinetic transport problems [72].

References↩︎