# COMBINED GALERKIN AND REGRESSION BASED ALGORITHM FOR PARAMETER DEPENDANT PDE'S

Hassan Ballout, Yvon Maday, Christophe Prud'Homme

► To cite this version:

Hassan Ballout, Yvon Maday, Christophe Prud'Homme. COMBINED GALERKIN AND REGRESSION BASED ALGORITHM FOR PARAMETER DEPENDANT PDE'S. 2025. hal-05143191

HAL Id: hal-05143191

https://cnrs.hal.science/hal-05143191v1

Preprint submitted on 4 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# COMBINED GALERKIN AND REGRESSION BASED ALGORITHM FOR PARAMETER DEPENDANT PDE'S*

HASSAN BALLOUT†, YVON MADAY‡, AND CHRISTOPHE PRUD'HOMME†

**Abstract.** Solving parameterized partial differential equations (PDEs) efficiently remains a major challenge in computational science, especially when the parameter-to-solution map is highly nonlinear and the (linear) Kolmogorov N-width decays slowly. Although linear reduced basis methods (RBM) perform well in low-complexity regimes, their efficiency deteriorates in cases where the solution manifold cannot be approximated accurately by (low dimensional) linear subspaces. Non-linear approaches, particularly those based on autoencoders and neural networks, offer enhanced representational power but often lack robustness, interpretability, and error control. In this work, we extend and deepen the new framework for the reduction of nonlinear models known as the Non-linear Compressive Reduced Basis Method (NLCRBM). Our approach combines nonlinear manifold compression with Galerkin projection to retain physical structure while achieving efficient dimensionality reduction. The method applies to a broad range of parameterized PDEs and supports adaptive mechanisms for error monitoring and manifold refinement. We illustrate its advantages through representative numerical experiments and discuss its potential to overcome key limitations of both linear RBM and data-driven reduction techniques.

**Key words.** Model reduction, parameterized PDEs, nonlinear manifold, reduced basis, Galerkin projection, autoencoder, neural networks, Kolmogorov width, sensing numbers, nonlinear approximation.

**AMS subject classifications.** 65M60, 41A65, 68T07

**1. Introduction.** Physical and industrial phenomena are frequently modeled by parameter-dependent partial differential equations (PDEs), whose repeated solution simulations for various parameter values can become extremely costly. To address this challenge, model reduction techniques—especially those based on Reduced Basis Methods (RBM)—have proven effective in accelerating the resolution of such PDEs. RBM techniques construct a reduced-dimensional space (the reduced space) from a small number of precomputed solutions during an offline stage and then solve the problem rapidly in this reduced space during an online stage for new parameter values. This framework enables near-real-time simulations in many settings, such as optimal control or uncertainty quantification, provided that the solution manifold can be well approximated in a low-dimensional linear subspace. The goal is not to offer a universal method for all problems, but rather to design a reduction strategy adapted to the particular class of parameter-dependent PDEs of interest.

A formal understanding of RBM performance is grounded in the concept of the solution manifold $\mathcal{K}$, the set of all PDE solutions as the parameters vary. This set forms a manifold in a high-dimensional function space (such as a Hilbert space), and the complexity of the problem is related to the effective dimension of this manifold. Traditional RBM seeks to approximate $\mathcal{K}$ by a low-dimensional vectorial subspace, with the optimal error quantified by the Kolmogorov N-width, which measures the

<footnote>*Submitted to the editors DATE.</footnote>

<footnote>**Funding:** This work has been partially funded by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant No 810367), project EMC2 (YM), and the France 2030 NumPEx Exa-MA (ANR-22-EXNU-0002) project managed by the French National Research Agency (ANR).</footnote>

<footnote>† Cemosis, IRMA UMR 7501, University of Strasbourg, CNRS, (hassan.ballout@math.unistra.fr, prudhomme@unistra.fr)</footnote>

<footnote>‡Sorbonne Université, CNRS, Université Paris Cité, Laboratoire Jacques-Louis Lions (LJLL), F-75005 Paris, France (yvon.maday@sorbonne-universite.fr).</footnote>
best achievable error when projecting $\mathcal{K}$ onto an $N$-dimensional linear subspace. If the $N$-width decays rapidly with $N$, the solution can be accurately represented with a small number of well-selected basis functions, and the Galerkin-projected RBM provides accurate, certifiable, and computationally efficient approximations.

However, in many practical cases, the Kolmogorov $N$-width decays slowly with $N$, reflecting the failure of any low-dimensional linear subspace to approximate $\mathcal{K}$ effectively. This so-called *Kolmogorov barrier* arises in problems with strong nonlinear parameter dependence, such as those with convection-dominated regimes or moving features such as shocks or interfaces. In such cases, linear RBMs require an excessively large number of modes to achieve acceptable accuracy, revealing a fundamental limitation of linear approximation. Geometrically, while $\mathcal{K}$ may have a low intrinsic (nonlinear) dimension, its linear approximation requires a high-dimensional subspace due to curvature or folding in the ambient space.

Several approaches have been explored to overcome this Kolmogorov barrier, including problem-specific transformations (see e.g. [16] and the references therein) and *machine learning* techniques (see e.g. [8] and the references therein). On the one hand, deterministic approaches, such as moving-frame methods or domain mapping, exploit expert knowledge to reduce complexity, but they are difficult to generalize. On the other hand, *neural network-based* methods aim to learn nonlinear, low-dimensional representations of $\mathcal{K}$ through *autoencoders*, which map high-dimensional solutions to a low-dimensional latent space via an encoder and reconstruct approximations via a decoder. Linear encoders and decoders correspond to traditional projection-based methods, whereas nonlinear maps offer greater flexibility and compression power. If the solution manifold can be unfolded in the latent space, accurate approximations can be obtained with a latent dimension comparable to the number of parameters.

Recent works have demonstrated this potential using deep learning for model reduction in parameterized PDEs. While these *data-driven* methods have shown promising results, they suffer from limitations: insufficient accuracy, lack of error guarantees, dependence on training data, and sensitivity to neural network architecture. These issues make purely data-driven approaches difficult to trust in applications that require robustness and accuracy.

To address these drawbacks, a promising direction is to retain the Galerkin framework while introducing *nonlinear compression* of the solution space. The idea is to define a nonlinear approximation manifold (through a decoder) and project the PDE onto it, thereby benefiting from both physical fidelity and nonlinear expressiveness. The *Nonlinear Compressive Reduced Basis Method* (NLCRBM) proposed in [5, 7] exemplifies this strategy, introducing a nonlinear encoder to better capture the intrinsic degrees of freedom of the problem. This approach achieves dimensionality reduction governed by the parameter space, rather than by the Kolmogorov $N$-width, and demonstrates significant efficiency gains in a test case.

NLCRBM has currently been validated on rather simple problems, and its application to more complex PDEs remains limited. In this article, we explore new aspects of the method's implementation, highlighting both its positive and limited features, and propose a new nonlinear reduced basis method that extends and generalizes NLCRBM. Our method remains Galerkin-based but introduces improved strategies for constructing the nonlinear manifold and its projection. It applies to a wider class of parameterized PDEs, including time-dependent and nonlinear cases, and offers adaptive mechanisms for controlling the approximation error. We present several numerical examples to demonstrate the robustness and efficiency of the method and discuss its advantages and limitations in comparison to both traditional RBM and
purely data-driven techniques.

The outline of the paper is as follows: in section 2, we recall the various notions of complexity, some of which were initially introduced in functional analysis as part of approximation theory, and which took on an important role in signal theory in the 2000s with the notions of parsimonious compression and compressed sensing theory. In section 3 we introduce the two test cases that will serve as illustrations for our analysis and allow us to compare our approach to more classical Reduced Basis Methods. In particular, we use a test case introduced in [15]. After recalling in section 4 the basics of the reduced basis method based on the linear notion of Kolmogorov N-width, we introduce in the next section 5 the different nonlinear approaches designed to overcome the Kolmogorov barrier. Our approach, based on the NLCRBM originally introduced in [5, 7], is then presented in detail in Section 6. We discuss its advantages, potential limitations, and the enhancements we propose to overcome these issues, leading to the improved NLCRBM++ method. Finally, Section 7 presents numerical results that support our conclusions.

**2. Notions of widths.** Let us consider a general Hilbert space $Y$, with a norm denoted as $||\cdot||_Y$, and a compact subset $\mathcal{K} \subset Y$, typically representing the solution manifold of a parametric PDE,

$$ (2.1) \qquad \mathcal{K} = \{u(\mu) \mid \mu \in \mathcal{D}\}, $$

where $u: \mathcal{D} \to Y$ is the solution map, and $\mathcal{D} \subset \mathbb{R}^p$ is a compact parameter set. To evaluate the complexity of $\mathcal{K}$, the classical approach relies, for a given dimension $n$, on the data of two continuous operators

• Encoder: $E: Y \to \mathbb{R}^n$

• Decoder: $D: \mathbb{R}^n \to Y$

An approximation of $u \in \mathcal{K}$ is then given by the composition:

$$ u \approx D(E(u)), $$

the quality of which can be measured by its maximum distortion over $\mathcal{K}$:

$$ \mathcal{E}_n(\mathcal{K}; E, D) = \max_{v \in \mathcal{K}} \|v - D(E(v))\|_Y. $$

This leads naturally to the concept of widths, by optimizing, for a given $n$, over a class of admissible encoder-decoder pairs $(E, D)$:

$$ \inf_{E,D} \mathcal{E}_n(\mathcal{K}; E, D). $$

Depending on the assumptions on $E$ and $D$, this yields several notions of widths:

• **Kolmogorov n-width:** If the decoder $D$ is restricted to be linear and $E$ is the orthogonal projection onto an $n$-dimensional subspace $\mathcal{V} \subset X$, we recover the "linear" Kolmogorov width:

$$ (2.2) \qquad d_n(\mathcal{K}) := \inf_{\mathcal{V} \subset Y, \dim(\mathcal{V})=n} \max_{v \in \mathcal{K}} \min_{w \in \mathcal{V}} \|v - w\|_Y, $$

which measures how well $\mathcal{K}$ can be approximated by a linear subspace of dimension $n$.

• **Sensing numbers:** If the encoder $E$ is restricted to be linear, but the decoder $D$ is allowed to be nonlinear, we obtain the notion of sensing numbers:

$$ (2.3) \qquad s_n(\mathcal{K}) := \inf_{\substack{D \in \mathcal{C}(\mathbb{R}^n, Y) \\ \lambda_1, \dots, \lambda_n \in Y'}} \max_{v \in \mathcal{K}} \|v - D(\lambda_1(v), \dots, \lambda_n(v))\|_Y. $$
• **Nonlinear Kolmogorov n-width:** If E and D are only required to be continuous (or potentially Lipschitz), we obtain the nonlinear Kolmogorov width (or manifold width):

$$ (2.4) \qquad \delta_n(\mathcal{K}) := \inf_{\substack{E \in \mathcal{C}(Y, \mathbb{R}^n) \\ D \in \mathcal{C}(\mathbb{R}^n, Y)}} \max_{v \in \mathcal{K}} \|v - D(E(v))\|_Y. $$

While linear approximation methods based on the Kolmogorov n-width $d_n$ are limited by the Kolmogorov barrier, as discussed in the introduction, nonlinear model reduction techniques based on the notions of $\delta_n$ and $s_n$ are designed to overcome this limitation. In particular, approaches based on the sensing numbers $s_n$ can be viewed as intermediate strategies, aiming to match the performance of $\delta_n$ while employing simpler, linear encoders.

One is then naturally led to consider the smallest dimension that allows for a full reconstruction of the elements in $\mathcal{K}$. This is equivalent to the existence of an encoder-decoder pair $(E, D)$ such that the identity on $\mathcal{K}$ is recovered:

$$ \mathrm{Id}_{\mathcal{K}} = D \circ E. $$

This motivates the definition of minimal dimensions associated with different types of approximation widths. For the nonlinear Kolmogorov width, we define:

$$ n_{\delta}(\mathcal{K}) := \min \{n \in \mathbb{N} \mid \delta_n(\mathcal{K}) = 0\}. $$

Similarly, for the sensing number, we define:

$$ n_s(\mathcal{K}) := \min \{n \in \mathbb{N} \mid s_n(\mathcal{K}) = 0\}. $$

While the Kolmogorov width $d_n(\mathcal{K})$ may decay slowly, $n_{\delta}(\mathcal{K})$ and $n_s(\mathcal{K})$ are directly linked to the intrinsic dimension of the solution manifold $\mathcal{K}$, i.e., the number of parameters $p$. In [8], focusing on the nonlinear width $\delta_n$, rigorous bounds for $n_{\delta}$ were established:

$$ p \le n_{\delta}(\mathcal{K}) \le 2p + 1. $$

In particular, if the solution map $u: \boldsymbol{\mu} \mapsto u(\boldsymbol{\mu})$ is continuous and injective, and thus bijective onto its image $\mathcal{K}$, then $n_{\delta}(\mathcal{K}) = p$, since $u^{-1}$ can serve as a continuous encoder<fnref n="1" />.

If we restrict to linear encoders, the dimension $n_s(\mathcal{K})$ satisfies the lower bound $n_s(\mathcal{K}) \ge n_{\delta}(\mathcal{K}) \ge p$. However, no general upper bound is known. In [3], it was shown numerically for a specific test case that $n_s(\mathcal{K}) = p$ using a linear encoder based on projection onto SVD modes. In [2], a more rigorous analysis proves that, at least locally, $n_s(\mathcal{K}) = p$: let us consider a solution $u^* = u(\boldsymbol{\mu}^*) \in \mathcal{K}$, and for all $\boldsymbol{\mu} \in B(\boldsymbol{\mu}^*, r) \cap \mathcal{D}$ let us express the solution $u(\boldsymbol{\mu})$ in an SVD-like approximation,

$$ u(\boldsymbol{\mu}) \approx u^* + \sum_{i=1}^{p} \tilde{u}_i^r(\boldsymbol{\mu}) \phi_i^r, $$

with $\tilde{u}_i^r(\boldsymbol{\mu}) = \langle u(\boldsymbol{\mu}) - u^*, \phi_i^r \rangle$, for $i = 1, \dots, p$. It is proven in [2] that there exists $R > 0$ such that the mapping

$$ (2.5) \qquad \boldsymbol{\mu} \mapsto (\tilde{u}_i^r(\boldsymbol{\mu}))_{i=1,\dots,p} $$

<footnote n="1">¹In [8], it is also stated that, under some regularity assumption, an improved error bounds $n_{\delta}(\mathcal{K}) \le 2p$ or even $n_{\delta}(\mathcal{K}) \le 2p - 1$ can be obtained</footnote>
is bijective for all $r \le R$. Hence, the complete information is encoded in the first $p$ coefficients $\tilde{u}_i^r(\mu)$, $i = 1, \dots, p$, and we have $n_s(K^{*,r}) = p$, where $K^{*,r} = \{u(\mu) | \mu \in B(\mu^*, r) \cap D\}$.

This indicates that even with simple linear encoders, one can achieve performance close to, or possibly equal to, the optimal dimension $n_\delta(K)$.

In this work, we present model reduction techniques grounded in the framework of nonlinear widths, with a particular focus on the NLCRBM based on the sensing numbers width (see [7, 3]). We highlight its strengths and limitations and propose an enhanced version of the method. These approaches are applied to different test cases to assess and compare their performances with other existing reduction methods.

**3. Test cases.** We introduce two test cases. In Section 3.1, we study a multi-parameter elliptic problem to highlight the relation between $n_s(K)$ and the parameter dimension $p$, as well as the challenges arising in multi-dimensional parameter spaces i.e. $p > 1$. In Section 3.2, we present a transport problem that serves as a benchmark for comparing various model order reduction techniques.

**3.1. Test case 1: A multiparameter problem.** We consider the problem of designing a thermal fin to remove heat from a surface effectively; see e.g. [13]. The two-dimensional fin, shown in Figure 1, consists of a vertical central **post** and $N_f$ horizontal **fins**. The thermal conductivity $k^i$, for $i = 1, \dots, N_f$ of the $i$-th fin (normalized with respect to the post conductivity $k^0 = 1$) may depend on a parameter, e.g. $k^i = \mu_i \in [0.1, 10]$, for $i = 1, \dots, p-1$, while $k^i = 1$ for $i = p, \dots, N_f$. The Biot number, which quantifies convective heat transfer to ambient air is denoted as Bi (larger Bi implies better convection) and may also vary $Bi = \mu_p \in [0.01, 1]$. The design parameters $\mu = (\mu_1, \dots, \mu_p)$ may take values in the compact design domain $D = ([0.1, 10]^{p-1} \times [0.01, 1]) \subset \mathbb{R}^p$.

The fin conducts heat from a prescribed uniform flux source at the root, $\Gamma_{\text{root}}$, through the large-surface-area fins to surrounding flowing air.

The full PDE formulation and boundary conditions are given in [3]. The natural functional space is here $Y = X = H^1(\Omega)$ We recall here only the corresponding variational formulation: for a given $\mu \in D$, find $u(\mu) \in X$ such that

$$ (3.1) \qquad a(u(\mu), v; \mu) = f(v), \quad \forall v \in X, $$

where the bilinear and linear forms are given by

$$ a(u, v; \boldsymbol{\mu}) = \sum_{i=0}^{N_f} k^i \int_{\Omega_i} \nabla u \cdot \nabla v \, d\boldsymbol{x} + \text{Bi} \int_{\Gamma_{\text{ext}}} uv \, d\boldsymbol{x}, \quad f(v) = \int_{\Gamma_{\text{root}}} v \, d\boldsymbol{x}. $$

The domains $\Omega_i$ and boundaries $\Gamma_{\text{root}}$ and $\Gamma_{\text{ext}}$ are shown in Figure 1, which illustrates the geometry and labeling used in this setup.

A second-order finite element (FE) method is used as the high-fidelity solver. The solution $u(\mu) \in X$ is approximated by $u_N(\mu) \in X_N \subset X$, where $X_N$ is a conforming FE space of dimension $N$. The corresponding high-fidelity solution manifold is

$$ \mathcal{K}_N = \{u_N(\mu) | \mu \in D\} \subset X_N. $$

After this discretization, there is thus a natural shift from $Y = X$ to $Y = X_N$ still provided with the same $H^1$-norm.

In this work, we fix $N_f = 5$, and the FE discretization is based on approximately $N = 90,000$ degrees of freedom (DoFs). The accuracy of this approximation was
<figure><img src="image_2.png" /><figcaption>FIG. 1. Schematic of thermal fin geometry with labeled subdomains $\Omega_i$ and boundaries $\Gamma_{\text{root}}$, $\Gamma_{\text{ext}}$.</figcaption></figure>

evaluated: the average $L^2$-norm error between the second-order FE solution and a third-order FE solution on the same mesh is approximately $10^{-4}$. This FE error may become dominant when the reduced basis error is smaller. All FE computations are performed using Feel++ library [12].

In Figure 2, we show representative high-fidelity solutions for fixed thermal conductivities $k^i = 1$ for all $i$, and varying Biot numbers $Bi \in \{0.01, 0.1, 1\}$. As expected, the average temperature near the root decreases with increasing $Bi$, reflecting enhanced convective cooling.

<figure><img src="image_3.png" /><figcaption>FIG. 2. High-fidelity solutions for $Bi = 0.01$ (left), $Bi = 0.1$ (middle), and $Bi = 1$ (right)</figcaption></figure>

## 3.2. Test case 2: Linear transport equation.
We follow [9, 15] and consider the one-dimensional linear transport equation:

$$ \partial_t u(x,t) + c \partial_x u(x,t) = 0, \quad x \in \mathbb{R}, \quad t \in (0, \infty), $$

whose exact solution is given by $u(x,t) = u_0(x - ct)$, where $c$ is a constant advection velocity, set to $c = 10$ and $u_0(x)$ is the specified initial condition at $t = 0$. The quantity $u(x,t)$ may represent, for example, the concentration of a pollutant being advected in a one-dimensional flow with constant velocity $c$. We consider the initial condition of the form:

$$ u_0(x) := u(x, 0) = \frac{1}{\sqrt{0.0002 \pi}} \exp \left( -\frac{(x - 0.1)^2}{0.0002} \right), \quad x \in \mathbb{R}. $$
We consider a spatial domain $x \in [0, 1] \subset \mathbb{R}$ and a time interval $t \in \mathcal{D} = [0, 0.1] \subset \mathbb{R}$.
Note that in this problem, we can treat time $t$ as the only parameter, i.e., $\mu = t \in \mathcal{D}$ (hence $p = 1$). The high-fidelity representation of the exact solution $u$ is obtained by evaluating it on $\mathcal{N} = 2^{12}$ equidistant points in the interval $[0, 1]$. Figure 3 shows the initial condition $u_0$ and the direction of the transport. Note that in [15], the authors consider the set $\mathcal{K}$ composed of vectors of $Y = \mathbb{R}^N$, i.e., they work directly with the vectors of degrees of freedom $u_{\mathcal{N}} \in \mathbb{R}^N$ associated with a high-fidelity approximation $u_{\mathcal{N}}$ of $u$.

<figure><img src="image_4.png" /><figcaption>FIG. 3. The initial condition and the direction of the transport.</figcaption></figure>

**4. Linear approximation: traditional RBM.** Let $\mathcal{V} \subset Y$ be some appropriate $N$-dimensional subspace of $Y$. For simplicity, and as in Section 3.2, we now consider $Y = \mathbb{R}^N$, and identify each high-fidelity approximation $u_{\mathcal{N}} \in X_{\mathcal{N}} \subset X$ with its vector of degrees of freedom $u_{\mathcal{N}} \in \mathbb{R}^N$ (e.g., the coordinates of $u_{\mathcal{N}}$ in a finite element basis).

Let $\mathcal{V} = [\varphi_1, \dots, \varphi_N] \in \mathbb{R}^{N \times N}$ be a basis matrix of $\mathcal{V}$ with columns orthonormal with respect to the Euclidean inner product $\langle \cdot, \cdot \rangle$. The traditional RBM, grounded in the concept of the Kolmogorov $N$-width $d_N(\mathcal{K})$, can be interpreted as an orthogonal projection onto the reduced space $\mathcal{V}$. Specifically, the projection of a high-fidelity vector $u_{\mathcal{N}} \in \mathbb{R}^N$ is given by:

$$P_{\mathcal{V}}(u_{\mathcal{N}}) = VV^{\top}u_{\mathcal{N}},$$

where $P_{\mathcal{V}}: \mathbb{R}^N \to \mathbb{R}^N$ denotes the orthogonal projection operator onto $\mathcal{V}$.

This projection can equivalently be expressed as the composition of a linear encoder-decoder pair $(E_V, D_V)$, defined as:

$$\bullet E_V : u_{\mathcal{N}} \mapsto \tilde{u}_N = V^{\top}u_{\mathcal{N}} \in \mathbb{R}^N,$$

$$\bullet D_V : \tilde{u}_N \mapsto u_N = V\tilde{u}_N \in \mathbb{R}^N.$$

Both $E_V$ and $D_V$ are linear, and the full projection is given by $P_V = D_V \circ E_V$.

The determination of the optimal subspace $\mathcal{V}$ realizing the minimal approximation (projection) error is generally an open problem. In practice, this issue has been addressed using two main approaches within the framework of traditional RBM. On one hand, the Singular Value Decomposition (SVD), or Proper Orthogonal Decomposition (POD) based on SVD, provides an optimal $N$-dimensional subspace in the mean squared error sense (i.e. based on $L^2(\mathcal{D})$-norm). On the other hand, greedy
algorithms construct quasi-optimal subspaces with respect to the worst-case error (i.e. based on $L^\infty(D)$-norm). In this work, we focus on the SVD/POD framework, and for completeness, we briefly recall its basic notions and properties.

We introduce a training set $\mathcal{D}_M = \{\mu_1, \mu_2, \dots, \mu_M\} \subset \mathcal{D}$ of $M < N$ parameter values, chosen large enough to sufficiently represent the parameter domain $\mathcal{D}$. For each $\mu_i \in \mathcal{D}_M$, we compute a high-fidelity solution vector $\mathbf{u}^{(i)} = \mathbf{u}_N(\mu_i) \in \mathbb{R}^N$, and then we assemble the snapshot matrix $S = [\mathbf{u}^{(1)}, \mathbf{u}^{(2)}, \dots, \mathbf{u}^{(M)}] \in \mathbb{R}^{N \times M}$. The SVD of $\mathbf{S}$ reads:

$$ \mathbf{S} = \mathbf{U} \mathbf{\Sigma} \mathbf{\Psi}^{\top}, $$

where $\mathbf{U} = [\zeta_1, \dots, \zeta_N] \in \mathbb{R}^{N \times N}$, $\mathbf{\Sigma} = [\psi_1, \dots, \psi_M] \in \mathbb{R}^{M \times M}$ are orthogonal matrices, and $\mathbf{\Sigma} \in \mathbb{R}^{N \times M}$ is a rectangular diagonal matrix containing the singular values of $\mathbf{S}$ in decreasing order $\sigma_1 \ge \sigma_2 \ge \dots \ge \sigma_M \ge 0$ on its main diagonal. By the Schmidt-Eckart-Young-Mirsky theorem, the best rank-$N$ approximation of $\mathbf{S}$ in the Frobenius norm is given by:

$$ (4.1) \qquad \tilde{\mathbf{S}}_N = \sum_{i=1}^{N} \sigma_i \zeta_i \psi_i^{\top} = \mathbf{V} \mathbf{V}^{\top} \mathbf{S}, $$

where $\mathbf{V} = [\zeta_1, \dots, \zeta_N] \in \mathbb{R}^{N \times N}$ is the matrix formed by the first $N$ left singular vectors of $\mathbf{S}$. This corresponds to the orthogonal projection of $\mathbf{S}$ onto the subspace $\mathcal{V} = \text{span}\{\zeta_1, \dots, \zeta_N\}$, that is, $\tilde{\mathbf{S}}_N = P_{\mathcal{V}}(\mathbf{S})$. This approximation satisfies the optimality property:

$$ (4.2) \qquad \sum_{i=1}^{M} \|\mathbf{u}^{(i)} - P_{\mathcal{V}}(\mathbf{u}^{(i)})\|_2^2 = \|\mathbf{S} - \tilde{\mathbf{S}}_N\|_F^2 = \min_{\mathbf{A} \in \mathbb{R}^{N \times M} \atop \text{rank}(\mathbf{A}) \le N} \|\mathbf{S} - \mathbf{A}\|_F^2 = \sum_{i=N+1}^{M} \sigma_i^2. $$

In particular, the subspace $\mathcal{V}$ spanned by the first $N$ left singular vectors minimizes the total projection error of the snapshots over all $N$-dimensional subspaces of $\mathbb{R}^N$. In addition, from (4.2), we have an explicit expression for the error committed using this projection, which can be considered as a lower bound for all approximations of $\mathbf{S}$ on an $N$-dimensional subspace of $\mathbb{R}^N$.

We note that the POD is based on the SVD and offers a computational advantage by working with the correlation matrix $\mathbf{C} = \mathbf{S}^{\top}\mathbf{S} \in \mathbb{R}^{M \times M}$, where typically $M \ll N$, thereby avoiding the high-fidelity dimension $N$ during the eigenvalue problem. Additionally, POD provides the possibility of recovering the optimality property (4.2) in more relevant functional norms, rather than just in the Euclidean norm. This can be done by working with $\mathbf{C} = \mathbf{S}^{\top}\mathbf{M}\mathbf{S}$, where $\mathbf{M} \in \mathbb{R}^{N \times N}$ is a high-fidelity matrix associated to a given inner product. For further details on the POD algorithm, we refer the reader to [3] or [14].

After constructing the reduced basis $\mathbf{V} = [\varphi_1, \dots, \varphi_N]$, one common approach to approximate, for a new value of the parameter $\mu$ the generalized coordinates vector $\hat{\mathbf{u}}_N \approx \tilde{\mathbf{u}}_N = \mathbf{V}^{\top}\mathbf{u}_N(\mu) \in \mathbb{R}^N$ – without evaluating the high fidelity solution $\mathbf{u}_N(\mu)$ – is to apply a Galerkin projection method. To formulate this, we now switch back to the functional setting. Each reduced basis vector $\varphi_i \in \mathbb{R}^N$ corresponds to a high-fidelity function $\varphi_i \in X_N \subset X$, and we define the reduced functional space as

$$ X_N := \text{span}\{\varphi_1, \dots, \varphi_N\} \subset X_N. $$
We then seek a reduced basis approximation $u_N(\boldsymbol{\mu}) \in X_N$ expressed as:

$$ (4.3) \qquad u_N(\boldsymbol{\mu}) = \sum_{i=1}^{N} \hat{u}_i^{\boldsymbol{\mu}} \varphi_i $$

such that:

$$ (4.4) \qquad a(u_N(\boldsymbol{\mu}), v_N; \boldsymbol{\mu}) = f(v_N; \boldsymbol{\mu}), \quad \forall v_N \in X_N. $$

For linear problems, by testing with the basis functions $\{\varphi_i\}_{i=1}^N$, we obtain the following linear system:

$$ (4.5) \qquad \mathbf{A}_N^{\boldsymbol{\mu}} \hat{\mathbf{u}}_N^{\boldsymbol{\mu}} = \mathbf{f}_N^{\boldsymbol{\mu}}, $$

where $(\mathbf{A}_N^{\boldsymbol{\mu}})_{ij} = a(\varphi_i, \varphi_j; \boldsymbol{\mu})$, and $(\mathbf{f}_N^{\boldsymbol{\mu}})_i = f(\varphi_i; \boldsymbol{\mu})$, for $i, j = 1, \dots, N$.

The first test case in 3.1, being elliptic, allows the traditional linear RBM based on the Kolmogorov N-width $d_N(\mathcal{K})$ to perform very well, i.e., $d_N(\mathcal{K})$ decays rapidly. The reduced space is generated by applying a POD in the norm induced by the bilinear form $a(\cdot, \cdot; \bar{\boldsymbol{\mu}})$, where $\bar{\boldsymbol{\mu}} = \{1, 1, 1, 1, 1, 0.1\}$ is a reference parameter, i.e., on $C = S^{\top} \mathbf{A}^{\boldsymbol{\mu}} S$.

In table 1, we show the number of POD modes, $N$, needed to reach a given accuracy $\varepsilon = 10^{-5}$ as a function of the parameter space dimension $p$. The accuracy is defined from the notion of Relative Information Content (RIC) i.e. $N$ is chosen as the smallest integer satisfying:

$$ I(N) = \frac{\sum_{i=1}^{N} \lambda_i}{\sum_{i=1}^{M} \lambda_i} \geq 1 - \varepsilon^2 $$

, where $\lambda_1, \dots, \lambda_M$ are the eigenvalues of $C$.

<table><caption>TABLE 1
$N$ vs $p$ for $\varepsilon = 10^{-5}$</caption><tr><td><i>p</i></td><td>1</td><td>2</td><td>3</td><td>4</td><td>5</td><td>6</td></tr><tr><td><i>N</i></td><td>8</td><td>21</td><td>33</td><td>45</td><td>54</td><td>64</td></tr></table>

Even though $N$ increases almost linearly with $p$, the decay can still be described as fast. Solving the reduced linear system (4.5) of complexity $O(N^3)$ can be performed in real-time, even for $p = 6$, corresponding to $N = 64$ modes. However, the choice of this simple multiparametric problem is mainly to illustrate the relation between $n_s(\mathcal{K})$, and $p$ for different values of $p$ up to 6.

For test case 3.2, the decay of the Kolmogorov N-width is significantly slower than in test case 3.1, even for $p = 1$, which is typical of transport-dominated and hyperbolic problems. This behavior is illustrated in figure 4, which shows the slow decay of the singular values. Note that for this test case, the reduced basis is constructed by applying the SVD directly (and not POD) to the snapshot matrix $\mathbf{S}$, to remain consistent with [15].
<figure><img src="image_5.png" /><figcaption>FIG. 4. Singular value decay.</figcaption></figure>

**5. Nonlinear decoder-based methods.** As discussed in the introduction, many nonlinear model order reduction techniques have been developed to mitigate the Kolmogorov barrier. Most of these approaches are based on the notion of the nonlinear Kolmogorov n-width $\delta_n(K)$, meaning they employ both a nonlinear encoder and a nonlinear decoder aiming at a complexity that depends on $n_\delta \ll N$, where $N$ is chosen such that $d_N(K) \le \varepsilon$, for a given accuracy $\varepsilon$.

Other methods based on the notion of sensing numbers $s_n(K)$ rely on a simpler linear encoder $E_V$ and a nonlinear decoder $D_{V,H}: \mathbb{R}^n \to \mathbb{R}^N$ defined as:

$$D_{V,H} = V + H,$$

where $V \in \mathbb{R}^{N \times n}$, and $H: \mathbb{R}^n \to \mathbb{R}^N$ is a nonlinear correction map. This formulation aims to approximate the solution manifold $K$ by the manifold:

$$K_n(V, H) = \{D_{V,H}(\tilde{u}_n) | \tilde{u}_n \in E_V(K)\}.$$

In this work, we focus on formulations where $H$ takes a particular form based on a feature map $F: \mathbb{R}^n \to \mathbb{R}^r$, such that:

$$H(\tilde{u}_n) = WF(\tilde{u}_n),$$

with $W \in \mathbb{R}^{N \times r}$, where $r$ is the dimension of the nonlinear correction term. The decoder then becomes:

$$(5.1) \qquad D_{V,W}(\tilde{u}_n) = V\tilde{u}_n + WF(\tilde{u}_n),$$

and we denote the approximation manifold $K_n(V, H)$ as:

$$(5.2) \qquad K_n(V, W) = \{D_{V,W}(\tilde{u}_n) | \tilde{u}_n \in E_V(K)\}.$$

We note that even though $K_n(V, W) \not\subset V$, we have $K_n(V, W) \subset V+W$, where $V, W$ denote the linear subspaces spanned by the columns of $V$ and $W$ respectively.

Many methods have been developed based on different choices of encoder-decoder pairs, i.e., different choices of $V, W$, and feature maps $F$. The objective is to reduce the complexity from $N$ to $n$, $n_s \le n \ll N$. In [4, 9, 15], the feature map $F$ is chosen to be a homogeneous quadratic function

$$F_q: \mathbb{R}^n \to \mathbb{R}^{r_q}, \quad \tilde{u}_n = (\tilde{u}_1, \dots, \tilde{u}_n) \mapsto (\tilde{u}_1\tilde{u}_1, \tilde{u}_1\tilde{u}_2, \dots, \tilde{u}_i\tilde{u}_j, \dots, \tilde{u}_n\tilde{u}_n), \text{ with } i \le j$$
where the dimension of the correction term $r_q = n(n+1)/2$. In [4, 9], the basis functions in the linear part, i.e., the columns of the matrix $V$, are chosen to be the first $n$ left singular vectors of a snapshot matrix $S \in \mathbb{R}^{N \times M}$. The matrix $W$ is then obtained by minimizing the approximation error in the Frobenius norm through the following least-squares problem:

$$W = \arg \min_{W' \in \mathbb{R}^{N \times r_q}} \left\| VV^T S + W'F_q(V^T S) - S \right\|_F^2,$$

where $F_q(V^T S) \in \mathbb{R}^{r_q \times M}$ denotes the column-wise evaluation of the quadratic feature map applied to the reduced coordinates. We refer to this method as the Quadratic SVD-based Method (QSVDM) and refer the reader to [4] for the complete algorithm and additional practical details.

The SVD/POD modes, while being quasi-optimal for linear approximation, are not necessarily optimal for a specific nonlinear decoder or feature map $F$. The feature map operates on the encoded representation—in this case, the projection onto the column space of $V$. However, this projection step can lead to a loss of information that may be essential for the nonlinear correction performed by $F$, especially if the basis $V$ is not tailored to the structure of $F$. In the standard SVD/POD approach, the basis $V$ is computed independently of the choice of $F$, which can limit performance. This issue is studied in detail in [15], where a new method is introduced to construct the basis $V$ specifically for the quadratic feature map $F_q$. The proposed approach selects $n$ basis vectors from the first $m \gg n$ SVD/POD modes using a greedy algorithm, rather than simply taking the first $n$ modes. We refer to this method as the Quadratic Greedy-based Method (QGM), and we refer the reader to [15] for the complete algorithm and further implementation details.

The QGM demonstrates a significant improvement in accuracy over the QSVDM. However, it remains suboptimal in the sense that it requires a reduced dimension $n \gg n_s$. This limitation stems from the use of a fixed quadratic feature map $F_q$, which lacks the expressiveness needed to fully exploit the intrinsic structure of the data.

To overcome this issue, the NLCRBM replaces the fixed feature map with a trainable nonlinear map $F_\theta$. This map is designed to reconstruct the coordinates associated with the SVD/POD basis functions of an index greater than $n$, using only the first $n$ coefficients. By allowing $F_\theta$ to be sufficiently flexible, the method achieves a complexity that depends only on $n \approx n_s$.

**6. Nonlinear compressive RBM.** In this section, we present the original NL-CRBM and its formulation using a Galerkin approach. We then discuss its limitations and introduce an enhanced variant, denoted NLCRBM++, which aims to address these issues. Finally, we provide an error analysis to characterize the expected accuracy of this method.

**6.1. NLCRBM.** Unlike the quadratic methods introduced in Section 5, the NLCRBM seeks to reproduce the same approximation as the traditional RBM by retaining the same reduced basis $\{\phi_i\}_{i=1}^N$ corresponding to a target accuracy $\epsilon$, while significantly reducing the computational complexity. This is achieved by relying only on $n \approx n_s \ll N$ degrees of freedom.

This reduction is justified by the need for only $n_s$ linear measurements when combined with a suitable nonlinear decoder. In the definition (2.3) of the sensing numbers $s_n(K)$, the optimal linear measurements $\lambda_1, \dots, \lambda_{n_s}$ are not available in
practice. However, the moments associated with the first $n \approx n_s$ SVD/POD modes $\{\langle \cdot, \varphi_i \rangle\}_{i=1}^n$ provide good, though not necessarily optimal, candidates.

This approach is supported by a local analysis in [2], where the bijection (2.5) maps the first $p = n_s$ coefficients $\{\langle u(\boldsymbol{\mu}), \varphi_i \rangle\}_{i=1}^p$ and the parameter $\boldsymbol{\mu}$. Over larger parameter domains $\mathcal{D}$, a slightly larger value $n > p$ may be required to capture the necessary information.

Once the first $n$ coefficients are known, the remaining $r = N-n$ coefficients can be reconstructed using a nonlinear feature map $F: \mathbb{R}^n \to \mathbb{R}^r$, which is approximated by a trainable nonlinear regression model $F_\theta$. The error introduced by this approximation is referred to as the “regression error”.

The NLCRBM can be described using a linear encoder $E_V$ constructed for example from the $n$ leading SVD/POD modes, where $V = [\varphi_1, \dots, \varphi_n] \in \mathbb{R}^{N \times n}$, a non-linear decoder $D_{V,W}$ that is of the form (5.1), where $W = [\varphi_{n+1}, \dots, \varphi_N] \in \mathbb{R}^{N \times r}$ contains the remaining $r$ SVD modes and a regression model $F_\theta$.

Similarly to what is done in Section 4, once the reduced basis $\{\varphi_i\}_{i=1}^N$ is computed, we adopt a Galerkin approach to approximate the vector of generalized coordinates $\hat{u}_N^\mu = [\hat{u}_n^\mu, \hat{u}_r^\mu] \approx \tilde{u}_N^\mu \in \mathbb{R}^N$, where $\hat{u}_n^\mu \in \mathbb{R}^n$ stands now for an approximation of the first $n$ coefficients: $V^\top u_N(\boldsymbol{\mu})$ and $\hat{u}_r^\mu \in \mathbb{R}^r$ stands for an approximation of the remaining ones: $W^\top u_N(\boldsymbol{\mu})$. Additionally $\hat{u}_r^\mu$ can be expressed in terms of $\hat{u}_n^\mu$, through a nonlinear map $F$, $\hat{u}_r^\mu = F(\hat{u}_n^\mu)$. Therefore, it suffices to solve for $\hat{u}_n^\mu$ and from Problem (4.4), it follows that:

$$ (6.1) \qquad a(u_n(\boldsymbol{\mu}), \varphi_i; \boldsymbol{\mu}) + a(u_r(\boldsymbol{\mu}), \varphi_i; \boldsymbol{\mu}) = f(\varphi_i; \boldsymbol{\mu}), \quad \forall i \in \{1, \dots, n\}, $$

where $u_n(\boldsymbol{\mu}) = \sum_{i=1}^n \hat{u}_i^\mu \varphi_i$ and $u_r(\boldsymbol{\mu}) = \sum_{k=n+1}^N \hat{u}_k^\mu \varphi_k$. This leads to the following nonlinear system in matrix form:

$$ (6.2) \qquad A_n^\mu \hat{u}_n^\mu = f_n^\mu - B^\mu F(\hat{u}_n^\mu), $$

where $(A_n^\mu)_{ij} = a(\varphi_i, \varphi_j; \boldsymbol{\mu})$, $(f_n^\mu)_i = f(\varphi_i; \boldsymbol{\mu})$, and $(B^\mu)_{ik} = a(\varphi_i, \varphi_k; \boldsymbol{\mu})$ for $i, j = 1, \dots, n$ and $k = n+1, \dots, N$.

In practice, $F$ is approximated by $F_\theta$ and we look for an approximation $\tilde{u}_N^\mu = [\tilde{u}_n^\mu, \tilde{u}_r^\mu] \approx \hat{u}_N^\mu$ by imposing:

$$ \tilde{u}_r^\mu \equiv F_\theta(\tilde{u}_n^\mu) \approx F(\tilde{u}_n^\mu), $$

the nonlinear system (6.2) is thus replaced by:

$$ (6.3) \qquad A_n^\mu \tilde{u}_n^\mu = f_n^\mu - B^\mu F_\theta(\tilde{u}_n^\mu), $$

Note that assembling the matrix $B^\mu$ similarly to $A_n^\mu$ and $f_n^\mu$ can benefit from affine decomposition or hyperreduction techniques to avoid operations involving the high-fidelity dimension $\mathcal{N}$. Note also that instead of solving the set of equations (6.3) one may rather minimize a residual over the set of all vectors $\tilde{v}_N = [\tilde{v}_n, F_\theta(\tilde{v}_n)] \in \mathbb{R}^N$ with unknown $\tilde{v}_n \in \mathbb{R}^n$.

**6.2. NLCRBM limitations.** To address the Kolmogorov barrier, NLCRBM reduces the model to $n \approx n_s$ degrees of freedom with a nonlinear decoder. This reduction, however, comes with limitations, discussed below before introducing an improved method.
**6.2.1. Loss of Linearity in the Reduced System.** A key limitation of the NLCRBM approach, also encountered in quadratic methods, lies in the loss of linearity induced by the use of a nonlinear decoder. Even if the original problem is linear, the introduction of a nonlinear reconstruction map renders the reduced problem nonlinear. As shown when adopting a Galerkin approach, the linear system (4.5) of dimension $N$ in the traditional RBM is transformed into a nonlinear system (6.2) of dimension $n$ in the NLCRBM, where $n \ll N$.

While this dimension reduction is expected to offer computational benefits, solving a nonlinear system introduces new challenges, particularly related to well-posedness.

The existence of a solution $\hat{u}_n^\mu \in \mathbb{R}^n$ to the reduced nonlinear problem (6.2) is guaranteed by the existence of the traditional RBM solution $\hat{u}_N^\mu \in \mathbb{R}^N$ of (4.5), since one can simply extract the first $n$ components of $\hat{u}_N^\mu$.

However, uniqueness is not guaranteed. In Section 6.2.2, we construct a toy example that explicitly demonstrates the non-uniqueness of the reduced nonlinear system. This problem is experienced in practice when solving (6.3), even when the reduced nonlinear problem is well-posed within the training domain of $F_\theta$. However, it is generally difficult to control the extrapolation behavior of $F_\theta$ in the multidimensional case $n > 1$. Certain regression models, such as artificial neural networks (ANN), may introduce irrelevant solutions outside the training region, and the solver may converge to these depending on the initial guess.

A second limitation concerns the computational efficiency. While the NLCRBM formulation avoids the $O(N^3)$ complexity of solving system (4.5), the complexity of solving system (6.3) depends on several factors, such as the number of iterations $K$ required by the solver and the inference cost $R(n, N)$ of the regression model $F_\theta$. Since $F_\theta$ must be evaluated at each iteration, the total online complexity includes a term of the form $K \cdot R(n, N)$, which can dominate the overall cost—particularly when $R(n, N)$ is non-negligible, as is often the case when using complex neural network architectures.

**6.2.2. Non-Uniqueness.** In this part, we illustrate an important limitation of the NLCRBM (present even for quadratic type methods: QSVDM or QGM): the lack of uniqueness in the solution of the resulting nonlinear reduced system. Unlike the traditional RBM, where the Galerkin projection leads to a well-posed linear system for linear problems, the nonlinear structure introduced by the regression model may give rise to multiple solutions.

Let us construct a toy example in $Y = \mathbb{R}^2$ to highlight this behavior and introduce $e_1 = (1, 0)$, $e_2 = (0, 1)$ as being the canonical basis set. The full-order solution $u(\mu) \in \mathbb{R}^2$ is defined as:

$$u(\mu) = c_1(\mu)e_1 + \gamma c_2(\mu)e_2,$$

where $\mu \in \mathcal{D} = [-\mu_{\max}, \mu_{\max}]$, $c_1(\mu) = \mu/\alpha$, and $c_2(\mu) = c_1(\mu)^2/\beta$. By construction, the two components $c_1$ and $c_2$ are orthogonal over $\mathcal{D}$. The constants $\alpha$ and $\beta$ are normalization factors chosen to normalize $c_1$ and $c_2$ respectively. Furthermore, $\gamma < 1$ is selected to preserve the singular value ordering after applying the SVD to the set $S = \{u(\mu), \mu \in \mathcal{D}\}$, i.e. $\varphi_i = e_i$, for $i = 1, 2$.

We then define a bilinear form $a(u, v; \mu) = u^\top A(\mu)v$, and a linear form $f(v; \mu) = v^\top b(\mu)$ where:

$$A(\mu) = \begin{bmatrix} 1 & \mu \\ \mu & 1 \end{bmatrix}, \text{ and } b(\mu) = \begin{bmatrix} c_1(\mu) + \gamma c_2(\mu) \mu \\ c_1(\mu) \mu + \gamma c_2(\mu) \end{bmatrix}.$$
By construction, $u(\mu)$ satisfies the following variational problem:

$$a(u(\mu), v; \mu) = f(v; \mu) \quad \text{for all } v \in \mathbb{R}^2.$$

Since the SVD provides the canonical basis vectors $\varphi_1 = e_1$ and $\varphi_2 = e_2$ as mentioned before, the traditional RBM approximation $u_N(\mu)$ with $N = 2$ yields:

$$u_2(\mu) = c_1(\mu)\varphi_1 + \gamma c_2(\mu)\varphi_2 = u(\mu),$$

i.e., it exactly recovers the full-order solution.

In contrast, the NLCRBM approximation $u_n(\mu)$ with $n = 1$ and $N = 2$, takes the form:

$$u_1(\mu) = a_1(\mu)\varphi_1 + \gamma/\beta a_1(\mu)^2\varphi_2,$$

which assumes that the second coefficient is expressed via a nonlinear map from the first (that has to be learned beforehand and is obviously quadratic here). This approximation is then bound to satisfy:

$$a(u_1(\mu), \varphi_1; \mu) = f(\varphi_1; \mu).$$

Substituting $u_1(\mu)$ into the bilinear form leads to a quadratic equation in $a_1$:

$$\gamma/\beta a_1(\mu)^2 a(\varphi_2, \varphi_1; \mu) + a_1(\mu) a(\varphi_1, \varphi_1; \mu) = f(\varphi_1; \mu),$$

that admits two distinct solutions $\forall \mu \neq 0$:

$$a_1^*(\mu) = c_1(\mu), \quad \text{and} \quad a_1^{**}(\mu) = -\frac{\beta}{\gamma\mu} - c_1(\mu) \neq c_1(\mu).$$

This example illustrates that the nonlinear reduced system can have multiple solutions, and depending on the initialization or numerical strategy, the solver may converge to an irrelevant solution $a_1^{**}$, which does not correspond to the true reduced coefficient.

In Figure 5, we illustrate this phenomenon: for $\mu_{max} = 0.9$, the left panel shows the exact coefficients $c_1(\mu)$ as a function of $\mu$, along with their corresponding traditional RBM approximations. The right panel displays the spurious coefficients $a_1^{**}(\mu)$.

<figure><img src="image_6.png" /><figcaption>FIG. 5. Left: exact coefficients $c_1(\mu)$ and traditional RBM approximations. Right: irrelevant NLCRBM approximations $a_1^{**}(\mu)$.</figcaption></figure>
### 6.3. NLCRBM++
We now present an enhanced version of the NLCRBM method, referred to as NLCRBM++, which aims to address the limitations discussed in Section 6.2.

The fundamental assumption behind the original NLCRBM is that the reduced $N-n$ coefficients $\hat{u}_r^\mu$ can be expressed as a function of the first $n$ coefficients $\hat{u}_n^\mu$, these ones being the true and only degrees of freedom of the approximation. This assumption implicitly relies on the ability to recover the parameter $\mu$ from $\hat{u}_n^\mu$, and thus recover $\hat{u}_r^\mu$.

However, in a broad class of parametric PDE problems, including the type considered in this work, the parameter $\mu$ is known and is precisely the one that is used during the online stage, to approximate the solution $u(\mu)$ directly.

Based on this observation, NLCRBM++ proposes to bypass the intermediate step of using $\hat{u}_n^\mu$, and instead predicts directly the remaining coefficients from the parameter $\mu$. More precisely, we replace the original nonlinear map $F: \mathbb{R}^n \to \mathbb{R}^r$ by a new map:

$$F^+: \mathcal{D} \in \mathbb{R}^p \to \mathbb{R}^r, \quad \mu \mapsto \hat{u}_r^\mu.$$

By doing so, the input dimension is independent of $n_s$, associated with the sensing numbers of $S$. This tuning in the input dimension can be beneficial for the accuracy of the regression model $F_\theta^+$ used to approximate $F^+$, particularly if the model is sensitive to the input dimension. In addition, the number of degrees of freedom $n$ becomes a user-defined parameter of the discretization method. While the choice of $n$ does not affect the input $\mu \in \mathcal{D}$ of $F_\theta^+$, it directly impacts the output dimension $N-n$ and the regression error. This aspect is discussed in more detail in Section 6.4.2.

The first $n$ coefficients, that are not predicted by $F_\theta^+$, can then be computed via a Galerkin approach, and the nonlinear system (6.3) is replaced by the following linear system:

$$ (6.4) \qquad A_n^\mu \hat{u}_n^\mu = f_n^\mu - B^\mu F_\theta^+(\mu). $$

where the approximation of $F^+$ by $F_\theta^+$, introduces a regression error as discussed in Section 6.1.

In Section 6.2, we showed that representing the higher coefficients $\hat{u}_r^\mu$ as a non-linear function of the lower ones $\hat{u}_n^\mu$ significantly reduces the number of degrees of freedom. However, this representation introduces challenges due to the nonlinearity embedded in the reduced system. The NLCRBM++ approach achieves a similar reduction in dimension while avoiding these difficulties.

The linear system (6.4) is well-posed and thus resolves the non-uniqueness issue present in the NLCRBM formulation. Furthermore, the inference cost $R(n, N)$ associated with evaluating $F_\theta^+$ is incurred only once, leading to a total online complexity of $O(n^3) + R(n, N)$.

From an accuracy standpoint, when using general-purpose regression models such as artificial neural networks (ANNs), decision trees, or random forests, we expect both NLCRBM and NLCRBM++ to achieve comparable accuracy—particularly when the first $p$ SVD coefficients are sufficient to reconstruct the parameter $\mu$, making the input space for both $F$ and $F^+$ effectively identical up to the bijection (2.5). In cases where this is not true, NLCRBM++ is expected to provide better accuracy due to its direct use of the parameter $\mu$.

### 6.4. Error analysis
We conclude this section with an error analysis of the proposed methods, specifically NLCRBM and NLCRBM++, to better understand their performance.
**6.4.1. Lower bound.** We first establish a lower bound for approximations on an $n$-dimensional nonlinear manifold $\mathcal{K}_n(V, W)$ of the form (5.2), based on the decoder (5.1).

As already mentioned in Section 5, we have $\mathcal{K}_n(V, W) \subset V+W$, where $\dim(V+W) = n+r = N$. Moreover, as seen in Section 4, the best rank-$N$ approximation $\tilde{S}_N$ of a snapshot matrix $S$ in the mean squared error sense is given by the projection onto the linear subspace spanned by the first $N$ left singular vectors of $S$, associated with the singular values $\sigma_1, \dots, \sigma_N$.

Consequently, any approximation on the nonlinear manifold $\mathcal{K}_n(V, W)$ cannot outperform this best linear approximation, regardless of the choice of $V$ and $W$. From the optimality property of the singular value decomposition (see (4.2)), the following lower bound holds:

$$ (6.5) \qquad \sum_{i=1}^{M} \min_{\hat{u}^{(i)} \in \mathcal{K}_n(V,W)} \left\| \hat{u}^{(i)} - u^{(i)} \right\|_2^2 \ge \sum_{i=1}^{M} \min_{\hat{u}^{(i)} \in V+W} \left\| \hat{u}^{(i)} - u^{(i)} \right\|_2^2 \ge \sum_{j=N+1}^{M} \sigma_j^2, $$

and the equality stands in the last part of the above inequality with the choice made in the NLCRBM.

**6.4.2. Regression Error and Choice of n.** We now focus on the NLCRBM (respectively, NLCRBM++) method, and recall that both methods aim to reproduce the same approximation as the traditional RBM. This is achieved by constructing the matrices $V$ and $W$ from the left singular vectors of the snapshot matrix. As indicated above, the error of these methods could, in case of perfect training (no regression error), coincide with the lower bound stated in (6.5), which we will refer to as the “reduced basis error” in the context of these two methods.

In actual practice, approximating the unknown feature map $F$ (respectively, $F^+$) with a learned model $F_\theta$ (respectively, $F_\theta^+$) introduces a “regression error” defined as:

$$ (6.6) \qquad \epsilon_n(u) = \sum_{i=n+1}^{N} |\hat{u}_i - \tilde{u}_i| $$

where $\hat{u}_N = (\hat{u}_1, \dots, \hat{u}_N)$ is the full set of generalized coordinates approximation, solution of (4.5), and $(\tilde{u}_{n+1}, \dots, \tilde{u}_N)$ are the predicted coefficients by $F_\theta$ (respectively, $F_\theta^+$). This error adds to the reduced basis error. Therefore, the total error of the NLCRBM and NLCRBM++ methods consists of two components: the reduced basis error and the regression error.

For the reduced basis error component, one can reduce it by increasing the total dimension $N = n+r$. The advantage of NLCRBM and NLCRBM++ is that the complexity does not scale with $O(N^3)$ as in the traditional RBM. Notably, this can be done by increasing $r$ while keeping $n$ fixed, offering greater flexibility. This is in contrast to quadratic methods, where increasing $N$ inherently requires increasing $n$, since the number of quadratic correction terms is given by $r = \frac{n(n+1)}{2}$, and thus also depends on $n$. Additionally, the quadratic methods QGM and QSVDM do not achieve the same level of accuracy as the traditional RBM with a reduced dimension $N = n+r$ — i.e., they do not reach the lower bound (6.5) — since they do not rely on the optimal SVD basis.

On the other hand, the regression error (6.6) depends on the approximation quality of the learned model $F_\theta$ (or $F_\theta^+$). In practice, for high-dimensional parameter
spaces (i.e., $p > 1$), achieving very low regression error often requires a large and diverse training dataset as well as a complex model architecture, which may significantly increase the inference cost.

In the NLCRBM++ method, the input dimension of the regression model $F_{\theta}^{+}$ is fixed and depends solely on the parameter dimension $p$; it is independent of the number $n$ of coefficients computed by projection. In contrast, in the standard NLCRBM approach, $n$ also determines the input dimension of the regression model $F_{\theta}$, since the first $n$ coefficients are used as input to predict the remaining $r = N - n$ coefficients.

Let us denote $e_i^n = |\hat{u}_i - \tilde{u}_i|$, for $i = n+1, \dots, N$, the regression error on the $i$-th coefficient when using a model with $n$ input components. As discussed earlier, for the NLCRBM, the essential information about the parameter $\mu$ is typically encoded in the first few coefficients. Let $n_{\min}$ be the smallest number of such coefficients sufficient to recover $\mu$ and enable meaningful regression.

We then distinguish two regimes:

* **When $n < n_{\min}$:** the regression problem is ill-posed—the model lacks sufficient input information to approximate the target coefficients accurately. In this case, it is generally impossible to make $e_i^n$ small for $i > n$.

* **When $n \ge n_{\min}$:** the behavior depends on the type of regression model:

- For structured models such as polynomials, increasing $n$ effectively enriches the model (e.g., enabling access to higher-degree terms), enhancing its approximation power. In this case, we often observe that for $n_{\min} \le n_1 < n_2$, the regression error satisfies $e_i^{n_2} \le e_i^{n_1}$ for $i > n_2$. This phenomenon is supported by results in [3], and extended in [6] where an adaptive strategy incrementally expands the model's input by including predicted coefficients.

- For general-purpose models such as artificial neural networks (ANNs) or random forests, increasing the number of input coefficients beyond $n_{\min}$ does not improve the regression accuracy. Once the necessary information is available, adding more inputs has negligible effect, i.e., $e_i^n \approx e_i^{n_{\min}}$ for $i > n$.

This is the regime considered in this work. In all numerical experiments, the input dimension of the regression model $F_{\theta}$ is fixed to $n_{\min}$, and $n$ refers solely to the number of coefficients computed via projection. The remaining $r = N - n$ coefficients are predicted using the learned model.

Interestingly, even in this case, it is possible to mitigate the regression error without necessarily improving the regression model itself. This can be achieved just by increasing the number $n$ of coefficients that are computed by projection rather than regression. Indeed, as a consequence of the decay of the singular values, the averaged behavior of the sequence $\hat{u}_N = (\hat{u}_1, \dots, \hat{u}_N)$, is a monotonically decreasing in modulus. Note now that the regression models generally exhibit a fixed relative error; this implies again, in average, a decreasing absolute error, i.e.,

$$|\hat{u}_i - \tilde{u}_i| \ge |\hat{u}_{i+1} - \tilde{u}_{i+1}| \quad \text{for } i \ge n+1,$$

The dominant contribution to the regression error then comes from the first few predicted coefficients. This explains why, increasing $n$ allows to decrease the overall regression error.

This strategy introduces a trade-off: increasing $n$ improves accuracy but also increases computational cost. For NLCRBM, the choice of $n$ is constrained by the requirement $n \ge n_{\min}$ to make the regression functional. In contrast, NLCRBM++ allows more flexibility—even the extreme case $n = 0$ is admissible, leading to a fully
regression-based method that predicts all *N* coefficients directly from the parameter (see e.g.[10] for a similar approach using Gaussian process regression).

As a final remark for this section, it is worth noting that when the decay of the singular values is slow, increasing *n* may not significantly reduce the error but still adds computational overhead. In such scenarios, minimizing *n* while improving the regression model may be more beneficial.

**7. Numerical Results.** To validate the performance of the proposed NLCRBM and NLCRBM++ methods, as well as the error analysis discussed in Section 6.4, we assess their effectiveness on the test cases in 3.

**7.1. Test case 1.:** As mentioned earlier, this test case aims, first, to demonstrate the practical link between $n_s$ (and $n_{\min}$) and the parameter dimension $p$, and, second, to highlight the challenges encountered when applying the NLCRBM method to multi-parametric problems, underline the need for the improved NLCRBM++ formulation.

In Figure 6, we illustrate the mean absolute prediction error $|\hat{u}_i - \tilde{u}_i|$ for $i = n+1, \dots, N$ over a test set. We observe that $n_s \le n_{\min} \le p+1$ for $1 \le p \le 6$. When using an ANN for $F_\theta$, good accuracy is obtained with $n_{\min} = p+1$ (for $p \ge 2$), and increasing $n$ (not illustrated) beyond this point does not improve the accuracy of the regression model. This indicates that the essential information is primarily encoded in the first $n_{\min}$ components.

The prediction accuracy can, however, be further enhanced either by increasing the training dataset size $n_{\text{train}}$, or by using a deeper or wider ANN architecture. In this experiment, we use $n_{\text{train}} = 9 \times 10^5$, and for the most complex setting ($p=6$), we employ a 5-layer ANN with architecture (32–64–128–256–512) in TensorFlow [1]. For other regression models, we use the Scikit-learn library [11].

<figure><img src="image_7.png" /><figcaption>FIG. 6. Regression error for different regression models and parameter dimension $p$.</figcaption></figure>
As a comparison between different approaches, we implement the online phase consisting of solving problems (4.5), (6.3) or (6.4) with different methods.

In Tables 2–4, we illustrate the full reconstruction errors with respect to the high-fidelity solution $u_N$ in both the $H^1$ and $L^2$ norms for $p = 1, 3, 6$.

In Table 2, we observe that when the regression model (spline, in this case) is sufficiently accurate—i.e., the regression error is smaller than the reduced basis error—both NLCRBM and NLCRBM++ (even for $n = 0$) attain the lower bound corresponding to the reduced basis error achieved by the traditional RBM.

<table><caption>TABLE 2
Mean $L^2$ and $H^1$ errors $P = 1$.</caption><thead><tr><th>Model</th><th>Mean $L^2$ Error</th><th>Mean $H^1$ Error</th></tr></thead><tbody><tr><td>Traditional RBM ($N = 1$)</td><td>6.50 × 10<sup>-1</sup></td><td>8.24 × 10<sup>-1</sup></td></tr><tr><td>Spline NLCRBM ($n = 1$)</td><td>5.93 × 10<sup>-6</sup></td><td>1.08 × 10<sup>-5</sup></td></tr><tr><td>Spline NLCRBM++ ($n = 0$)</td><td>5.93 × 10<sup>-6</sup></td><td>1.08 × 10<sup>-5</sup></td></tr><tr><td>Traditional RBM ($N = 8$)</td><td>5.93 × 10<sup>-6</sup></td><td>1.08 × 10<sup>-5</sup></td></tr></tbody></table>

In Table 3, the ANN achieves the best performance. However, the regression error remains dominant and prevents reaching the reduced basis error limit. In addition, we can already observe the gain obtained by choosing $n = 4$ compared to $n = 0$ following the error analysis in 6.4.2.

<table><caption>TABLE 3
Mean $L^2$ and $H^1$ errors $P = 3$.</caption><thead><tr><th>Model</th><th>Mean $L^2$ Error</th><th>Mean $H^1$ Error</th></tr></thead><tbody><tr><td>Traditional RBM ($N = 4$)</td><td>1.21 × 10<sup>-1</sup></td><td>1.81 × 10<sup>-1</sup></td></tr><tr><td>ANN NLCRBM ($n = 4$)</td><td>2.37 × 10<sup>-5</sup></td><td>4.50 × 10<sup>-5</sup></td></tr><tr><td>ANN NLCRBM++ ($n = 4$)</td><td>2.76 × 10<sup>-5</sup></td><td>4.68 × 10<sup>-5</sup></td></tr><tr><td>ANN NLCRBM++ ($n = 0$)</td><td>1.45 × 10<sup>-4</sup></td><td>2.10 × 10<sup>-4</sup></td></tr><tr><td>Traditional RBM ($N = 33$)</td><td>4.28 × 10<sup>-6</sup></td><td>2.12 × 10<sup>-5</sup></td></tr></tbody></table>

The NLCRBM method suffers from the issue of non-uniqueness, as discussed in Section 6.2.2. This limitation is particularly evident in practice for $p = 6$ when using ANN regressors, where it becomes difficult to ensure convergence across the test set. Therefore, in Table 4, we focus on the NLCRBM++ method, which exhibits approximately the same accuracy for a given $n$, while providing significantly better stability and reduced online computation time.

<table><caption>TABLE 4
Mean $L^2$ and $H^1$ errors $P = 6$.</caption><thead><tr><th>Model</th><th>Mean $L^2$ Error</th><th>Mean $H^1$ Error</th></tr></thead><tbody><tr><td>Traditional RBM ($N = 7$)</td><td>9.27 × 10<sup>-2</sup></td><td>1.55 × 10<sup>-1</sup></td></tr><tr><td>ANN NLCRBM++ ($n = 0$)</td><td>1.38 × 10<sup>-3</sup></td><td>1.50 × 10<sup>-3</sup></td></tr><tr><td>ANN NLCRBM++ ($n = 7$)</td><td>1.59 × 10<sup>-4</sup></td><td>2.52 × 10<sup>-4</sup></td></tr><tr><td>Traditional RBM ($N = 64$)</td><td>3.00 × 10<sup>-6</sup></td><td>3.33 × 10<sup>-5</sup></td></tr></tbody></table>

From the results above, we observed that for $p > 1$, achieving a regression model with perfect accuracy becomes challenging. As a consequence, the total reconstruction
error is often dominated by the regression component. As discussed in Section 6.4.2, this issue can be mitigated by increasing $n$, not as an input dimension to the model, but to restrict the regression task on predicting only the last $N-n$ coefficients.

This effect is illustrated in Figure 7 for $p = 3$ and $N = 33$, where we show the decay of the reconstruction error in the $L^2$-norm as $n$ increases. The curve corresponds to the NLCRBM++ method using an ANN for the regression model. For comparison, we also include the lower bound (corresponding to the reduced basis error) and the fully regression-based case with $n = 0$, which highlights the gain in accuracy achieved by increasing $n$.

<figure><img src="image_8.png" /><figcaption>FIG. 7. $L^2$ error as a function of $n$ for NLCRBM++ based on ANN model, compared with the RB lower bound and full-regression case.</figcaption></figure>

**7.2. Test case 2.** In this section, we evaluate the performance of the proposed NLCRBM and NLCRBM++ methods on the test case 3.2. For comparison, we include the results of standard quadratic approaches based on SVD and greedy selection, reproducing the benchmarks reported in [15]. This allows us to assess the improvements brought by these methods over existing reduced-order modeling techniques. Note that the reduced approximations are obtained by directly projecting the high-fidelity solution vectors $u_N(\mu) \in \mathbb{R}^N$ onto the approximation spaces defined by each method. No Galerkin projection is performed, in consistency with the analysis presented [15].

Figure 8 illustrates the points discussed in Section 6.4. The curve labeled *Linear* corresponds to the traditional RBM with $n$ basis functions (i.e., the first $n$ left singular vectors), for $n = 1, \dots, 20$. The *Lower Bound* curve corresponds to the error bound given in Equation (6.5), which is the error of a traditional RBM with reduced dimension $N = n + \frac{n(n+1)}{2}$. As previously discussed, this serves as a lower bound for any approximation within a subspace of dimension $N$.

The curves corresponding to the quadratic methods QSVDM and QGM are reproduced from [15]. We observe that, due to their use of different bases than the SVD/POD modes, these methods do not reach the lower bound. The QGM method
performs better than QSVDM but is known to be sensitive to preprocessing and may
face challenges in wider parameter spaces.

For the NLCRBM method, the number of independent coefficients is $n_{\min} = p + 1 = 2$. The regression model $F_{\theta}$ is chosen to be a random forest with $n_{\min}$ input features and $\frac{n(n+1)}{2}$ output targets, to allow a fair comparison with other approaches. As seen in the figure, the NLCRBM matches the lower bound until the regression error from $F_{\theta}$ becomes dominant.

In the NLCRBM++ setting, the regression model $F_{\theta}^{+}$ takes only $p = 1$ input feature, allowing us to use a spline model that achieves high accuracy and enables the method to reach the lower bound. Additionally, as previously discussed, the total number of basis functions $N$ in both NLCRBM and NLCRBM++ is independent of $n$. Therefore, one can choose a large value $N = N_{\max} = 230$, and achieve high accuracy even with $n = n_{\min}$ in NLCRBM or $n = 0$ in NLCRBM++, subject to the regression model's accuracy.

In this particular test case, increasing $n$ does not significantly reduce the regression error, as the decay of the Kolmogorov width $d_n(K)$, and hence the singular values, is extremely slow.

<figure><img src="image_9.png" /><figcaption>FIG. 8. Reconstruction error curves for various reduced-order models.</figcaption></figure>

**Acknowledgments.** As part of the “France 2030” initiative, this work has benefited from a national grant managed by the French National Research Agency (Agence Nationale de la Recherche), attributed to the ExaMA project of the NumPEx PEPR program under reference ANR-22-EXNU-0002, and has also been funded in part by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant No. 810367), project EMC2 (Y.M.). We acknowledge the constructive discussions with Dr Joubine Aghili on the NLCRBM method and for his valuable comments throughout this work.
[1] M. ABADI, et al., TensorFlow: Large-scale machine learning on heterogeneous systems, 2015, http://tensorflow.org/. Software available from tensorflow.org.

[2] J. AGHILI, H. BALLOUT, C. PRUD'HOMME, AND Y. MADAY, Nonlinear compressive reduced basis approximation: when Taylor meets Kolmogorov, (2025). In preparation.

[3] H. BALLOUT, Y. MADAY, AND C. PRUD'HOMME, Nonlinear compressive reduced basis approximation for multi-parameter elliptic problem, Springer Nature Switzerland, Cham, 2024, pp. 55-73, https://doi.org/10.1007/978-3-031-75802-7-4, https://doi.org/10.1007/978-3-031-75802-7-4.

[4] J. BARNETT AND C. FARHAT, Quadratic approximation manifold for mitigating the Kolmogorov barrier in nonlinear projection-based model order reduction, Journal of Computational Physics, 464 (2022), p. 111348, https://doi.org/https://doi.org/10.1016/j.jcp.2022.111348, https://www.sciencedirect.com/science/article/pii/S0021999122004107.

[5] J. BARNETT, C. FARHAT, AND Y. MADAY, Neural-network-augmented projection-based model order reduction for mitigating the Kolmogorov barrier to reducibility, Journal of Computational Physics, 492 (2023), p. 112420.

[6] A. BENSALAH, A. NOUY, AND J. SOFFO, Nonlinear manifold approximation using compositional polynomial networks, 2025, https://arxiv.org/abs/2502.05088, https://arxiv.org/abs/2502.05088.

[7] A. COHEN, C. FARHAT, Y. MADAY, AND A. SOMACAL, Nonlinear compressive reduced basis approximation for PDE's, Comptes Rendus. Mécanique, 351 (2023), pp. 357-374, https://doi.org/10.5802/crmeca.191.

[8] N. FRANCO, A. MANZONI, AND P. ZUNINO, A deep learning approach to reduced order modelling of parameter dependent partial differential equations, Mathematics of Computation, 92 (2022), p. 483-524, https://doi.org/10.1090/mcom/3781, http://dx.doi.org/10.1090/mcom/3781.

[9] R. GEELEN, S. WRIGHT, AND K. WILLCOX, Operator inference for non-intrusive model reduction with quadratic manifolds, Computer Methods in Applied Mechanics and Engineering, 403 (2023), p. 115717, https://doi.org/https://doi.org/10.1016/j.cma.2022.115717, https://www.sciencedirect.com/science/article/pii/S0045782522006727.

[10] M. GUO AND J. S. HESTHAVEN, Reduced order modeling for nonlinear structural analysis using Gaussian process regression, Computer Methods in Applied Mechanics and Engineering, 341 (2018), pp. 807-826, https://doi.org/https://doi.org/10.1016/j.cma.2018.07.017, https://www.sciencedirect.com/science/article/pii/S0045782518303487.

[11] PEDREGOSA ET AL., Scikit-learn: Machine learning in Python, Journal of Machine Learning Research, 12 (2011), pp. 2825-2830.

[12] C. PRUD'HOMME ET AL., feelpp/feelpp: Feel++ v111 preview 9, 2024, https://doi.org/10.5281/zenodo.10837178, https://zenodo.org/records/10837178. Accessed: 2024-07-22.

[13] C. PRUD'HOMME, D. V. ROVAS, K. VEROY, L. MACHIELS, Y. MADAY, A. T. PATERA, AND G. TURINICI, Reliable Real-Time Solution of Parametrized Partial Differential Equations: Reduced-Basis Output Bound Methods, Journal of Fluids Engineering, 124 (2001), pp. 70-80, https://doi.org/10.1115/1.1448332, https://hal.science/hal-00798326.

[14] A. QUARTERONI, A. MANZONI, AND F. NEGRI, Reduced Basis Methods for Partial Differential Equations: An Introduction, vol. 92 of Unitext, Springer, 2015, https://doi.org/10.1007/978-3-319-15431-2.

[15] P. SCHWERDTNER AND B. PEHERSTORFER, Greedy construction of quadratic manifolds for non-linear dimensionality reduction and nonlinear model reduction, 2024, https://arxiv.org/abs/2403.06732, https://arxiv.org/abs/2403.06732.

[16] T. TADDEI, Compositional maps for registration in complex geometries, SIAM Journal on Scientific Computing, 47 (2025), pp. B220-B251.