# Residual correction for POD-Galerkin reduced order models via conditional VAEs: Application to the bias-aware solution of forward and inverse problems

Dawen Wu, Ludovic Chamoin

## ► To cite this version:

Dawen Wu, Ludovic Chamoin. Residual correction for POD-Galerkin reduced order models via conditional VAEs: Application to the bias-aware solution of forward and inverse problems. Computer Methods in Applied Mechanics and Engineering, 2026, 449, pp.118528. 10.1016/j.cma.2025.118528. hal-05359296

HAL Id: hal-05359296

https://hal.science/hal-05359296v1

Submitted on 11 Nov 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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
<figure><img src="image_2.png" /></figure>

# Residual correction for POD-Galerkin reduced order models via conditional VAEs: Application to the bias-aware solution of forward and inverse problems

Dawen Wu ID a, b, Ludovic Chamoin ID c, *

## ARTICLE INFO

### Keywords:

- Reduced order models
- Proper orthogonal decomposition
- Conditional variational autoencoder
- Residual correction
- Parameterized PDEs

## ABSTRACT

POD-Galerkin reduced order models (PGROMs) significantly accelerate the solution of parameterized partial differential equations (PDEs) but can suffer from accuracy limitations across the parameter space. We propose a hybrid framework that augments a PGROM with a data-driven residual corrector based on a conditional variational autoencoder (cVAE). In the offline stage, we (i) construct a standard PGROM from high-fidelity snapshots and (ii) train a cVAE on PGROM residuals defined as the differences between PGROM approximations and the corresponding high-fidelity solutions. The cVAE learns a mapping from time-parameter pairs to residual fields. In the online stage, for a new parameter instance, the PGROM provides a baseline approximation while the cVAE decoder predicts the residual; adding the two yields a corrected, high-accuracy solution. The effectiveness of the proposed hybrid approach is demonstrated on three parameterized PDE problems: a heat dissipation problem with varying thermal conductivity, a pollutant transport problem with varying convection velocity, and a plane-strain elastodynamics problem with a circular inclusion whose stiffness and density are parameterized.

## 1. Introduction

Partial differential equations (PDEs) are ubiquitous in engineering and physics. They describe a wide range of physical phenomena - such as heat conduction, fluid flow, electromagnetic wave propagation and structural deformation - by describing relationships involving unknown functions and their partial derivatives. In particular, parameterized PDEs treat certain key factors in the governing equations or boundary conditions (such as material constants, geometric configurations, or source intensities) as parameters. This allows for a more flexible and realistic representation of physical systems and enables researchers to analyze system behavior across the entire parameter space.

Traditional numerical methods, including finite element, finite difference, and finite volume methods, treat the PDE corresponding to each parameter as an independent problem. For every new parameter value, these methods must start the solution process from scratch. This lack of information sharing across different parameter instances leads to significant computational expense. The challenge becomes particularly acute in multi-query scenarios - such as engineering design optimization or uncertainty quantification - where hundreds or even thousands of PDE solves may be required. Furthermore, if the PDE is three-dimensional, nonlinear, or requires fine spatial and temporal resolution, solving even a single parameter instance can be computationally intensive, making a comprehensive

<footnote>a CNRS@CREATE, 1 Create Way, #08-01 Create Tower, 138602, Singapore</footnote>

<footnote>b School of Computing, National University of Singapore, 13 Computing Drive, 117417, Singapore</footnote>

<footnote>c Université Paris-Saclay, CentraleSupélec, ENS Paris-Saclay, CNRS, Laboratoire de Mécanique Paris-Saclay (LMPS), 91190, Gif-sur-Yvette, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: dawen.wu@cnrsatcreate.sg (D. Wu), ludovic.chamoin@ens-paris-saclay.fr (L. Chamoin).</footnote>
parameter study impractical in many cases. Therefore, developing methods that are both highly accurate and computationally efficient for solving parameterized PDEs is crucial.

An effective approach to dealing with parameterized PDEs involves model order reduction techniques, such as Proper Orthogonal Decomposition (POD). This approach starts by collecting a series of representative snapshots, each corresponding to the spatially discretized solution of the PDE at different parameters and times. Applying POD, often implemented via singular value decomposition, to these snapshots extracts a low-dimensional reduced basis that captures the dominant features of the solution. The original high-dimensional model is then projected onto this reduced basis using Galerkin projection, yielding a reduced order model (ROM). The low-dimensional ROM solution is projected back into the original high-dimensional space to obtain the full approximate solution field. This approach, called POD-Galerkin (PGROM), significantly reduces the computational cost, as it has much fewer dimensions than the original full-order model. For more details on model order reduction, see [1–3].

Despite its computational advantages, maintaining high accuracy across the entire parameter space with a PGROM can be challenging, especially for problems with complex parameter dependencies or sharp solution features. The goal of this paper is to enhance the accuracy of the PGROM by employing a complementary approach-specifically, the variational autoencoder.

An autoencoder (AE) is a classical neural network (NN) model that learns a low-dimensional representation of data by training an encoder to compress the input into a latent space and a decoder to reconstruct the original data from this representation [4]. The variational autoencoder (VAE), introduced by [5], significantly extends this concept by incorporating probabilistic modeling. VAEs treat the latent space variables as following a probability distribution, typically Gaussian, optimizing a lower bound on the data log-likelihood. This probabilistic framework enhances generative capabilities and robustness compared to deterministic AEs, although challenges like posterior collapse and blurry reconstructions have spurred further research [6].

Many variants have been developed to refine the performance of the VAE and tailor it to specific tasks. For instance, the β-VAE [7] adjusts the weight of the Kullback-Leibler (KL) divergence term in the loss function to encourage better disentanglement of latent factors, leading to more interpretable representations [8]. The conditional VAE (cVAE) [9] incorporates conditioning variables (e.g., labels or parameters) into both the encoder and decoder, enabling targeted generation and prediction conditioned on specific inputs. Other notable advancements include improvements in hierarchical structures [10], integration with normalizing flows for more flexible posterior distributions [11], and exploration of alternative divergence measures [12]. These developments continue to expand the capabilities and applications of VAEs.

Consequently, VAEs and their variants have found widespread applications across diverse domains. Examples include high-fidelity image generation [10,13], anomaly detection in industrial systems and time series data [14], drug discovery through the generation of novel molecular structures [15], and even in areas like reinforcement learning and robotics [16]. Furthermore, recent studies [17,18] have demonstrated the potential of AEs and VAEs as standalone solvers for parameterized PDEs. These studies show that models trained solely on snapshot data-without explicit knowledge of the underlying governing equations-can effectively predict solutions for new parameter configurations.

Our primary contribution is the development of a hybrid method to improve the PGROM. This improvement is achieved by using a cVAE to correct the residuals. Our approach consists of three main stages:

*   **Offline – PGROM construction:** We gather a set of snapshots from simulations conducted on different PDE parameters and times. Using these snapshots, we construct a PGROM.

*   **Offline – cVAE training:** We utilize the PGROM to recompute approximate solutions for the same set of training parameters. By comparing these approximate solutions against the original high-fidelity snapshots, we calculate residuals (errors) that serve as training data for the cVAE.

*   **Online solution for a new parameter:** For a new, unseen PDE parameter, we first generate an initial approximate solution using the PGROM. Concurrently, we input the new parameter into the trained cVAE decoder to estimate the residual error. Combining these two, we derive a hybrid solution as:

$$ \text{Hybrid Solution} = \text{PGROM approximation} + \text{cVAE residual correction} $$

There has been significant research on learning nonlinear reduced models from data, which often involves using NNs to define a nonlinear mapping from high-dimensional state fields to a compact latent representation [19–21]. In contrast to these non-intrusive approaches, other methods like [22] employ NNs in an intrusive manner: to learn a nonlinear manifold onto which the governing physical equations are projected.

Data-driven ROMs typically use an AE or VAE to compress and reconstruct snapshots and advance temporal dynamics in the latent space with a separate NN. For example, [23] employed a β-VAE to learn the latent representation and a transformer to forecast the next-step latent code from a preceding latent history; future flow snapshots are then reconstructed by the decoder. This approach performs well for a single PDE configuration, but it cannot be directly applied to parameterized PDEs. For parameterized PDEs, [24] used an AE to build a latent space and use a fully-connected NN to regress from the time-parameter pair (t, λ) to the latent code. Building on this idea, [25] extended the approach by applying POD to effect a substantial dimensionality reduction, thereby shrinking the network and accelerating training. [26] explicitly incorporated physics in the loss-introducing a backbone (bias) term and residual components, together with pretraining and finetuning-to improve robustness in regimes with limited data.

Unlike the above methods, this research does not propose a new purely data-driven ROM. Instead, we adopt a hybrid residual-learning-correction strategy based on the classical equation-driven ROM. First, a PGROM delivers a baseline solution. Then, a cVAE receives the input of the time-parameter pair (t, λ) and predicts the residual between the PGROM approximation and the true solution. At inference, the predicted residual is added to the PGROM output to produce a hybrid solution with greater accuracy. In contrast to
<figure><img src="image_3.png" /><figcaption>Fig. 1. Overview of the hybrid ROM-cVAE framework. Offline preparation involves constructing the PGROM and training the cVAE with the residuals. The online solution combines the PGROM approximation with the cVAE residual prediction for enhanced accuracy.</figcaption></figure>

the data-driven ROMs, our hybrid approach significantly reduces the NN's workload. Instead of predicting the entire field, the network is tasked with a much simpler objective: learning only the residual correction. This focused approach renders the learning task more tractable and enhances robustness. Fig. 1 provides an overview of our hybrid method, contrasting it with the equation-driven and data-driven ROMs.

Our hybrid solution strategy is an instance of what is increasingly referred to as Hybrid AI in recent literature. The following discussion reviews some key developments in this area. For instance, recent research has explored various hybrid approaches that combine physics-based models with data-driven techniques to improve prediction accuracy and handle model discrepancies [27]. Some methods focus on correcting physical model parameters or discovering missing governing operators [28,29], while others directly infer and correct for model errors or biases [30–32]. Further approaches learn the dynamics of the deviation between model predictions and observed data, sometimes having to deal with unmeasured variables or ensure temporal stability [33,34]. Distinctly, our work retains the traditional PGROM and introduces a cVAE specifically trained to predict and correct the residual error inherent in the PGROM approximation across the parameter space. Thus, rather than modifying the governing equations of the model or its construction, our contribution lies in a novel hybrid framework that exploits the generative power of cVAE to improve the accuracy of established PGROMs by directly addressing their prediction residuals.

The remainder of the paper is organized as follows: Section 2 introduces the problem formulation and details the PGROM. Section 3 describes the cVAE, including its architecture, training procedure, and its role in correcting the PGROM. Section 4 presents the solution to parameterized PDEs and compares the accuracy of the standalone PGROM against the proposed hybrid model in forward and inverse problems. Finally, Section 5 concludes the paper and outlines potential directions for future research.

## 2. Equation-driven method: PGROM

This section presents the model order reduction method based on POD for parameterized PDEs. For the sake of exposition, this study considers a parabolic PDE; however, it should be noted that the proposed method can be readily extended to other types of PDEs.

### 2.1. Parameterized PDE

Consider the following parameterized parabolic PDE:

$$ 
\begin{aligned}
\rho(x; \lambda) \frac{\partial u}{\partial t} - \nabla \cdot \sigma(u, \nabla u; \lambda) &= f(x, t; \lambda) & x \in \Omega, t \in (t_0, T], \\
u(x, t = t_0; \lambda) &= u_0(x) & x \in \Omega, \\
u(x, t; \lambda) &= g(x, t; \lambda) & x \in \partial\Omega, t \in (t_0, T].
\end{aligned} 
\quad (1) 
$$
Here, $u(x, t; \lambda)$ denotes the solution defined on the spatial domain $\Omega \subset \mathbb{R}^d$ and over the time interval $t \in [t_0, T]$. The function $\rho(x; \lambda)$ represents the density or other related parameters, $\sigma$ corresponds to the stress or diffusion term, and $f(x, t; \lambda)$ is an external forcing term. The Dirichlet boundary conditions prescribe the solution on $\partial\Omega$, while the initial condition defines the state of the system at $t_0$.

## 2.2. Full order model (FOM)

By applying spatial discretization techniques such as the finite element method, Eq. (1) can be transformed into a system of first-order nonlinear ordinary differential equations. Let $y(t) \in \mathbb{R}^N$ denote the full order discrete solution vector, where $N$ is the number of degrees of freedom. The semi-discrete FOM is expressed as

$$M(\lambda) \dot{y}(t) + g(y(t); \lambda) = f(t; \lambda), \quad t \in (t_0, T], \quad (2)$$

where

* $M(\lambda) \in \mathbb{R}^{N \times N}$ is the mass matrix arising from the discretization of $\rho(x; \lambda)$.

* $g(y(t); \lambda) \in \mathbb{R}^N$ represents the discretized internal force, obtained from discretizing $-\nabla \cdot \sigma(u, \nabla u; \lambda)$, and may include nonlinear diffusion or convection effects.

* $f(t; \lambda) \in \mathbb{R}^N$ is the discretized external forcing term.

The initial condition is obtained by discretizing the PDE's initial condition $u_0(x)$, i.e., $y(t_0) = y_0$.

## 2.3. POD-Galerkin ROM (PGROM)

The FOM, as presented in Eq. (2), is usually highly dimensional, making it computationally expensive to solve, particularly when solutions involving many parameters are required. To address this issue, we use the PGROM approach, which significantly reduces dimensionality. The construction process is given below.

First, snapshots of the full order solution are collected from numerical simulations. Specifically, assume that the full order solutions are obtained at various time instants $t_i$ and for different parameters $\lambda_j$:

$$(y(t_i; \lambda_j))_{i=1,...,N_i, j=1,...,N_j}$$

These snapshots are arranged as columns to form the snapshot matrix

$$S = [y(t_1; \lambda_1), y(t_2; \lambda_1), \dots, y(t_{N_1}; \lambda_1), y(t_1; \lambda_2), \dots, y(t_{N_2}; \lambda_2), \dots, y(t_{N_i}; \lambda_i)] \in \mathbb{R}^{N \times m}, \quad (3)$$

where $m = N_i \times N_j$ is the total number of snapshots.

By performing a singular value decomposition on the matrix $S$, one can extract a reduced basis matrix $V \in \mathbb{R}^{N \times n}$ (with $n \ll N$), i.e.,

$$V = [\phi_1, \phi_2, \dots, \phi_n], \quad (4)$$

where $\phi_1, \phi_2, \dots, \phi_n$ are the left singular vectors corresponding to the $n$ largest singular values of the matrix $S$.

The full-order solution can be approximated within the subspace spanned by the columns of $V$, as follows:

$$y(t; \lambda) \approx \hat{y}(t; \lambda) = V q(t; \lambda), \quad (5)$$

where $q(t) \in \mathbb{R}^n$ represents the reduced coordinates. By projecting the FOM of Eq. (2) onto the subspace span{V} using a Galerkin projection, the ROM is obtained:

$$M^{(n)}(\lambda) \dot{q}(t; \lambda) + V^T g(V q(t; \lambda); \lambda) = V^T f(t; \lambda) \quad (6)$$

where the reduced mass matrix defined as $M^{(n)}(\lambda) = V^T M(\lambda) V$, and the initial condition is

$$q(t = t_0; \lambda) = V^T y_0. \quad (7)$$

## 3. Hybrid method: PGROM + conditional VAE

**Model Architecture.** Fig. 2 illustrates the architecture of the proposed cVAE. In this framework, the encoder receives the PGROM residual as input, while the decoder is designed to reconstruct the residual. Unlike a conventional VAE, the cVAE decoder accepts not only a sample $z$ from the latent distribution but also the conditioning variables: the PDE parameter $\lambda$ and the time $t$. Specifically, $z, \lambda$, and $t$ are typically concatenated and fed into the decoder's initial layers. Both the encoder and decoder consist of a series of convolutional and fully connected layers arranged in a mirror symmetric fashion. Once the model is trained, the cVAE decoder can predict the PGROM residual for a new parameter and time instance.

**Dataset.** The cVAE is trained on a dataset composed of PGROM residuals and associated parameters. The dataset is constructed through the following steps: First, we use the PGROM to approximate the snapshots $S$ defined in Eq. (3). The PGROM approximations are expressed as:

$$S_{\text{PGROM}} = [\hat{y}(t_1; \lambda_1), \hat{y}(t_2; \lambda_1), \dots, \hat{y}(t_{N_i}; \lambda_i)] \in \mathbb{R}^{N \times m}, \quad (8)$$
<figure><img src="image_4.png" /><figcaption>Fig. 2. Architecture of the conditional Variational Autoencoder (cVAE) used for residual prediction. The encoder maps the PGROM residual (re-shaped as an image) to a latent distribution (μ, σ). The decoder reconstructs the residual conditioned on the latent sample z, the PDE parameter κ (representing λ in general), and time t.</figcaption></figure>

where ŷ(t; λ) denotes the PGROM solution at time t and parameter λ. The residual matrix, representing the discrepancy between the true solution and the PGROM approximations, is calculated as:

$$R = S - S_{\text{PGROM}}. \qquad (9)$$

Each residual vector in R is reshaped into a two-dimensional array before being provided as input to the cVAE encoder. Additionally, label information comprising time t and parameter λ is required. The labels are formulated as:

$$\{(t_1, \lambda_1), (t_2, \lambda_1), \dots, (t_{N_t}, \lambda_{N_t})\}. \qquad (10)$$

Each pair $(t_i, \lambda_i)$ serves as input to the cVAE decoder.

**Loss Function.** For a PGROM residual $r$ corresponding to the time $t$ and parameter $\lambda$, the loss function is defined as:

$$\ell(w; r, t, \lambda) = \underbrace{\|r - p_{\omega}(z, t, \lambda)\|_{2}^{2}}_{\text{Reconstruction Loss (MSE)}} + D_{\text{KL}}(\mathcal{N}(\mu, \sigma) | \mathcal{N}(0, I)), \qquad (11)$$

where $z = \mu + \sigma \cdot \epsilon$ represents the latent variable, with $\mu$ and $\sigma$ being the outputs of the encoder, and $\epsilon$ a sample drawn from the standard normal distribution. $p_{\omega}$ denotes the decoder parameterized by $\omega$. The term $D_{\text{KL}}$ represents the KL divergence between two probability distributions. $\mathcal{N}(\mu, \sigma)$ denotes a multivariate normal distribution with mean vector $\mu$ and a diagonal covariance matrix with diagonal entries $\sigma_i^2$, where $\mu$ and $\sigma$ (vector of standard deviations) are outputs of the encoder parameterized by $\phi$. $\mathcal{N}(0, I)$ is the standard Gaussian prior.

**Residual Correction by cVAE.** After training, the cVAE is employed to enhance the accuracy of the PGROM. We refer to this approach as PGROM+cVAE, and its detailed procedure is outlined as follows:

* For a new, unknown parameter $\lambda$, the PGROM is first used to obtain the approximate solution at the corresponding times, i.e., $\{\hat{y}(t_1; \lambda), \hat{y}(t_2; \lambda), \dots, \hat{y}(t_{N_t}; \lambda)\}$.

* The parameter $\lambda$ is combined with the corresponding time instances to form the label set $\{(t_1, \lambda), (t_2, \lambda), \dots, (t_{N_t}, \lambda)\}$. This set is then fed into the cVAE decoder to generate predictions of the PGROM residuals, denoted as $\{r(t_1; \lambda), r(t_2; \lambda), \dots, r(t_{N_t}; \lambda)\}$.

* The final solution is the PGROM solution corrected by the cVAE residual prediction:

$$\{\hat{y}(t_1; \lambda) + r(t_1; \lambda), \hat{y}(t_2; \lambda) + r(t_2; \lambda), \dots, \hat{y}(t_{N_t}; \lambda) + r(t_{N_t}; \lambda)\}. \qquad (12)$$

Fig. 3 illustrates how the cVAE corrects the PGROM.

We note that the proposed cVAE framework relies on a structured mesh for the problem's spatial discretizations. This is a consequence of reshaping the residual field into a 2D image-like array, a format necessary for the convolutional layers which are highly effective at learning spatial patterns from grid-based data. Future extension of this methodology to unstructured meshes could be pursued either by interpolating the residuals onto a uniform background grid or, more fundamentally, by replacing the convolutional neural networks with graph neural networks that can operate directly on irregular mesh topologies.

## 4. Numerical tests

**Experimental Setup.** The cVAE used in this study is structured as follows: The encoder consists of three convolutional layers with output channels 16, 32, and 64, respectively, ending with a fully connected layer that outputs the latent space parameters $(\mu, \sigma)$. Each convolutional layer has a kernel size of 3, a stride of 2, and a padding of 1. The ReLU activation function is used after each
<figure><img src="image_5.png" /><figcaption>Fig. 3. Online prediction phase using the hybrid PGROM + cVAE model. For a new PDE parameter $\lambda$ (represented as $\kappa$ in the diagram for generality), the pre-computed PGROM generates a fast approximation. Concurrently, the trained cVAE decoder predicts the corresponding residual, conditioned on $\lambda$, time $t$ (implicitly handled for multiple steps), and a randomly sampled latent variable $z \sim \mathcal{N}(0, 1)$. The final, more accurate solution is the sum of the PGROM approximation and the predicted residual.</figcaption></figure>

convolutional layer. The latent space dimensionality is set to 5. The decoder mirrors the architecture of the encoder, employing ReLU similarly after its convolutional layers. The mini-batch size is set to 512, and we use the Adam optimizer for training with a learning rate of 0.001.

**Evaluation Metrics.** The performance of the model is evaluated using the following metrics:

* cVAE Loss: Defined by Eq. (11), which includes a reconstruction loss term and a KL divergence term.

* Mean Absolute Error: Defined as

$$E_{\text{mae}} = \frac{1}{m} \sum_{i=1}^{m} \| y_{\text{true}}^{(i)} - y_{\text{predict}}^{(i)} \|_{1}, \qquad (13)$$

where $y_{\text{true}}^{(i)}$ is the true snapshot, $y_{\text{predict}}^{(i)}$ the approximate solution from the model being evaluated (PGROM or PGROM + cVAE), and $m$ is the total number of snapshots in the test set. $\| \cdot \|_1$ represents the $L_1$ norm.

* Mean Maximum Error: This metric characterizes the maximum error in a snapshot and is defined as

$$E_{\text{max}} = \frac{1}{m} \sum_{i=1}^{m} \| y_{\text{true}}^{(i)} - y_{\text{predict}}^{(i)} \|_{\infty}, \qquad (14)$$

where $\| \cdot \|_{\infty}$ denotes the infinity norm.

* mean Relative L2 Error (Gap Error): Defined as

$$E_{\text{gap}} = \frac{1}{m} \sum_{i=1}^{m} \frac{\| y_{\text{true}}^{(i)} - y_{\text{predict}}^{(i)} \|_{2}}{\| y_{\text{true}}^{(i)} \|_{2}}, \qquad (15)$$

where $\| \cdot \|_2$ is the $L_2$ norm.

We test the methodology on three parameterized PDE problems. Section 4.1 considers a heat conduction problem with varying thermal conductivity, Section 4.2 investigates a convection-diffusion problem with varying convection velocity, and Section 4.3 deals with a parameterized elastodynamics problem with a circular inclusion. Finally, Section 4.4 applies the proposed framework to an inverse problem for the heat equation.

## 4.1. Parameterized heat dissipation problem

In this experiment, we consider a heat dissipation problem on a two-dimensional plane. The problem is described by the following heat equation:

$$\frac{\partial T}{\partial t} - \nabla \cdot (\kappa(x) \nabla T) = f(x), \quad x \in \Omega, t \in [0, T], \qquad (16)$$

where $T(x, t)$ is the temperature field. The spatial domain is the unit square $\Omega = [0, 1]^2$, and the maximum simulation time is $T = 0.3$.

The thermal conductivity of the material $\kappa(x)$ is a spatially dependent piecewise function defined as

$$\kappa(x) = \begin{cases} \kappa_{\text{in}}, & (x_1 - 0.5)^2 + (x_2 - 0.5)^2 \le r_{\text{in}}^2, \\ 1, & \text{otherwise,} \end{cases} \qquad (17)$$
<figure><img src="image_6.png" /><figcaption>Fig. 4. Evolution of cVAE reconstruction loss (left, log scale) and KL divergence loss (right) during training for the heat dissipation problem. Both training (orange) and testing (blue) losses are shown. The reconstruction loss converges, and the KL loss quickly drops near zero, indicating successful training without significant overfitting. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

where $\kappa_{in}$ is the thermal conductivity inside a circular region centered at (0.5, 0.5) with radius $r_{in}$, and outside the circle the conductivity is 1. Both $\kappa_{in}$ and $r_{in}$ serve as parameters in this experiment.

The heat source term $f(x)$ is defined as:

$$f(x) = \begin{cases} 2e^{-10((x_1-0.8)^2+(x_2-0.8)^2)}, & (x_1-0.8)^2+(x_2-0.8)^2 \le 0.1^2, \\ 0, & \text{otherwise,} \end{cases} \qquad (18)$$

centered on (0.8, 0.8).

The initial condition at $t = 0$ is set as:

$$T(x, t = 0) = \begin{cases} 5, & 0.3 < x_1 < 0.7 \text{ and } 0.3 < x_2 < 0.7, \\ 0, & \text{otherwise.} \end{cases} \qquad (19)$$

The boundary condition is homogeneous Dirichlet:

$$T(x, t = 0) = 0, \quad x \in \partial\Omega, t \in [0, T]. \qquad (20)$$

The snapshot dataset is constructed by iterating over different values of parameters $\kappa_{in}$ and $r_{in}$. For the training dataset, $\kappa_{in}$ takes values from 3.0 to 10.0 with a step of 0.5, and $r_{in}$ ranges from 0.1 to 0.3 with a step of 0.02. For the testing dataset, the parameters are shifted so that $\kappa_{in}$ takes values from 3.25 to 9.75 (step 0.5) and $r_{in}$ from 0.11 to 0.29 (step 0.02), ensuring that the test parameters do not coincide with the training parameters. Solution snapshots are collected at discrete time instances $t$ ranging from 0.01 to 0.30 with a step of 0.01. For each parameter pair ($\kappa_{in}$, $r_{in}$), we solve the governing equation using the finite element method to obtain the temperature field $T(x, t)$ at these time instances, which constitute the snapshots.

Fig. 4 shows the evolution of the cVAE loss values during training. It can be seen that both the reconstruction loss and the KL divergence loss decrease with increasing epochs. After 50,000 epochs, the training and testing reconstruction losses converge from an initial level on the order of $10^{-1}$ to $1.2 \times 10^{-5}$ and $7.8 \times 10^{-6}$, respectively, indicating that the model effectively reconstructs the PGROM residual. Furthermore, the close agreement between the training and testing losses demonstrates that the model does not suffer from overfitting. The KL loss converges to nearly zero after only 1000 epochs, indicating that the latent variables are effectively constrained to follow a standard normal distribution.

Fig. 5 shows the simulation results of this heat dissipation problem for a specific parameter pair, ($\kappa_{in} = 4.0$, $r_{in} = 0.12$) at time $t = 0.11$. From top-left to bottom-right, the panels show: the FOM solution, the approximate solution from the PGROM, the true residual (FOM-PGROM), the predicted residual from cVAE, and the solution obtained by combining PGROM with cVAE. Clearly, the PGROM solution deviates significantly from the true solution, and the residual plot highlights regions with high error. The cVAE predicted residual closely matches the true residual. The hybrid model (PGROM + cVAE) significantly improves the standalone PGROM result, demonstrating that cVAE effectively compensates for the shortcomings of PGROM.

Fig. 6 shows the performance of the error metrics on the test set, including $E_{mae}$, $E_{max}$, and $E_{gap}$. It is evident that the hybrid model, PGROM + cVAE, has lower errors compared to the standalone PGROM. Specifically, the PGROM + cVAE model achieves $E_{mae} = 1.57 \times 10^{-3}$, $E_{max} = 1.27 \times 10^{-2}$, and $E_{gap} = 4.00\%$, all of which are lower than the PGROM errors of $E_{mae} = 2.15 \times 10^{-2}$, $E_{max} = 6.60 \times 10^{-2}$, and $E_{gap} = 23.85\%$. These statistical results confirm the effectiveness of cVAE in error correction.

## 4.2. Parameterized pollutant transport problem

This subsection considers a pollutant transport problem, governed by the following convection-diffusion equation:

$$\frac{\partial c}{\partial t} + v \cdot \nabla c - \nabla \cdot (D(x)\nabla c) = f(x), \quad x \in \Omega, t \in [0, T], \qquad (21)$$
<figure><img src="image_7.png" /><figcaption>Fig. 5. Solution for $\kappa_{in} = 4.0$, $r_{in} = 0.12$ at time $t = 0.11$. From top-left to bottom-right: the FOM solution, the PGROM solution, the true residual (FOM-PGROM), the cVAE predicted residual, and the combined PGROM+cVAE solution.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 6. Quantitative comparison of error metrics (mean ± std dev over the test set) for the heat dissipation problem. The hybrid PGROM+cVAE model consistently achieves lower mean absolute error ($E_{mae}$), mean maximum error ($E_{max}$), and mean relative L2 error ($E_{gap}$) compared to the standalone PGROM.</figcaption></figure>

where $c(x, t)$ is the pollutant concentration, $\Omega = [0, 1]^2$ is the spatial domain, and the maximum simulation time is $T = 0.5$. The velocity field is defined as

$$v = (v_{mag} \cos \theta, v_{mag} \sin \theta), \qquad (22)$$

where the parameter $v_{mag}$ represents the magnitude of the velocity and the angle $\theta$ is fixed at $\pi/4$. The diffusion coefficient $D(x)$ is a spatially varying function defined as

$$D(x) = \begin{cases} D_{base} \times D_{ratio}, & \sqrt{(x_1 - 0.5)^2 + (x_2 - 0.5)^2} \le 0.3, \\ D_{base} \times (1 + w(x) \times (D_{ratio} - 1)), & \text{otherwise,} \end{cases} \qquad (23)$$

where the base diffusion coefficient is fixed at $D_{base} = 0.01$, $D_{ratio}$ defines the ratio of the diffusion coefficient inside the central circle (radius 0.3) to the base value $D_{base}$. The weighting function $w(x)$ is defined as

$$w(x) = \max \left\{ 0, \frac{0.4 - \sqrt{(x_1 - 0.5)^2 + (x_2 - 0.5)^2}}{0.1} \right\}. \qquad (24)$$

The pollutant source $f(x)$ is defined as a localized Gaussian function centered at (0.2, 0.2) within a radius of 0.1, and zero elsewhere:

$$f(x) = \begin{cases} s \times e^{-50((x_1 - 0.2)^2 + (x_2 - 0.2)^2)}, & (x_1 - 0.2)^2 + (x_2 - 0.2)^2 \le 0.01, \\ 0, & \text{otherwise,} \end{cases} \qquad (25)$$

where $s = 2.0$ is the strength of the source centered at (0.2, 0.2).
<figure><img src="image_9.png" /><figcaption>Fig. 7. Evolution of cVAE reconstruction loss (left, log scale) and KL divergence loss (right) during training for the pollutant transport problem. Both training (orange) and testing (blue) losses converge well, indicating effective learning. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 8. Visual comparison of solutions for the pollutant transport problem with $v_{mag} = 0.65$, $D_{ratio} = 1.95$ at $t = 0.5$. Panels arranged as in Fig. 5. The PGROM exhibits errors, particularly near sharp gradients, which are well-approximated by the cVAE residual prediction, resulting in an improved PGROM + cVAE solution.</figcaption></figure>

The initial condition is specified as

$$c(x, 0) = 0, \quad x \in \Omega, \tag{26}$$

and the boundary condition is homogeneous Dirichlet:

$$c(x, t) = 0, \quad x \in \partial\Omega, \quad t \in [0, T]. \tag{27}$$

The velocity $v_{mag}$ and the diffusion ratio $D_{ratio}$ are considered as parameters in this experiment. The snapshot dataset is constructed by varying these parameters. Specifically, for the training dataset, both $v_{mag}$ and $D_{ratio}$ take values from 0.5 to 2.0 with a step of 0.1, while for the testing dataset, both vary from 0.55 to 1.95 with a step of 0.1, ensuring that the test parameters are offset from those in the training set. Solution snapshots are collected at discrete time instances $t$ ranging from 0.01 to 0.50 with a step of 0.01. For each parameter pair ($v_{mag}$, $D_{ratio}$), we solve the governing equation using the finite element method to obtain the pollutant concentration field $c(x, t)$ at these time instances, which constitute the snapshots.

Fig. 7 shows the training of the cVAE model. It can be observed that both the reconstruction loss and the KL loss decrease steadily as the training progresses. Specifically, the reconstruction loss decreases from an initial level of approximately $5 \times 10^{-3}$ to nearly $3 \times 10^{-7}$ after 50,000 epochs, indicating that the cVAE model accurately reconstructs the PGROM residuals. In addition, the close match between the training and testing losses suggests that the model does not suffer from overfitting. The KL loss decreases rapidly and approaches zero within 1000 epochs, indicating that the latent variables have been successfully constrained to follow a standard normal distribution.
<figure><img src="image_11.png" /><figcaption>Fig. 9. Quantitative comparison of error metrics (mean ± std dev over the test set) for the pollutant transport problem. The PGROM + cVAE model demonstrates significantly lower errors ($E_{mae}$, $E_{max}$, $E_{gap}$) compared to the standalone PGROM.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 10. Convergence with respect to the number of modes for the pollutant transport instance with $v_{mag} = 0.65$ and $D_{ratio} = 1.95$ at $t = 0.5$. The top rows display the PGROM residuals, the cVAE-predicted residuals, and the residuals of the hybrid PGROM + cVAE solution. The bottom panel reports the corresponding quantitative error metrics. Smaller residual magnitudes (closer to 0) indicate higher predictive accuracy.</figcaption></figure>

Fig. 8 shows a comparison of the solutions of the convection-diffusion equation for a specific parameter combination ($v_{mag} = 0.65$, $D_{ratio} = 1.95$) at time $t = 0.5$. From top-left to bottom-right, the panels display: the FOM solution, the approximate solution from the PGROM, the true residual (FOM-PGROM), the predicted residual from cVAE, and the hybrid model (PGROM + cVAE). It is evident that the PGROM error is primarily concentrated in regions where the pollutant concentration changes rapidly. The cVAE predicted residual aligns closely with the true residual in both location and magnitude, indicating that the model effectively captures the deficiencies of the PGROM.

Fig. 9 compares the errors of the PGROM and PGROM + cVAE over the test dataset. The results show that the hybrid model, PGROM + cVAE, outperforms the standalone PGROM on all metrics. Specifically, the PGROM + cVAE model achieves $E_{mae} = 2.75 \times 10^{-4}$, $E_{max} = 6.19 \times 10^{-3}$, and $E_{gap} = 2.2\%$, compared to the standalone PGROM results of $E_{mae} = 6.93 \times 10^{-4}$, $E_{max} = 8.88 \times 10^{-3}$, and $E_{gap} = 5.2\%$. These results show that cVAE successfully corrects the PGROM error, thereby improving the overall accuracy and reliability of the approximate solutions.

Fig. 10 demonstrates the convergence of the proposed PGROM + cVAE as the number of POD modes increases. The top rows contrast the PGROM residuals with the cVAE-predicted residuals and the residuals after PGROM + cVAE correction, while the bottom
<figure><img src="image_13.png" /><figcaption>**Fig. 11.** Space-time solution of the plane-strain elastodynamics problem at $Y_{\text{ratio}} = 5.0$ and $\rho_{\text{ratio}} = 1.0$. Columns show snapshots at $t = 0.01, 0.02, 0.03, 0.04, 0.05$ (left to right). Rows (top to bottom): $u_x$ (horizontal displacement), $u_y$ (vertical displacement), $\epsilon_{xx}$ (normal strain in x), $\epsilon_{yy}$ (normal strain in y), and $\epsilon_{xy}$ (shear strain). The red circle marks the circular inclusion. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

panel reports the corresponding error metrics. As the mode count grows, the standalone PGROM error decreases and PGROM + cVAE follows the same monotone trend. With 30 modes, the MSE of PGROM + cVAE reaches $2.28 \times 10^{-7}$.

### 4.3. Parameterized elastodynamics with a circular inclusion (plane strain)

We consider two-dimensional linear elastodynamics under plane strain on the unit square $\Omega = [0, 1]^2$ over a final time $T = 0.5$. The unknown is the displacement field $u : \Omega \times [0, T] \to \mathbb{R}^2$, $u(x, t) = (u_x, u_y)^\top$. The governing equation reads

$$ \rho(x; \lambda) \partial_t u - \nabla \cdot \sigma(u; x, \lambda) = 0 \quad x \in \Omega, t \in (0, T], \qquad (28) $$

with the isotropic (plane-strain) constitutive law

$$ \sigma(u; x, \lambda) = \lambda_L(x; \lambda) \text{tr}(\varepsilon(u)) I + 2 \mu(x; \lambda) \varepsilon(u), \quad \varepsilon(u) = \frac{1}{2} (\nabla u + (\nabla u)^\top), \qquad (29) $$

where $I$ is the $2 \times 2$ identity. In terms of Young's modulus $Y(x; \lambda)$ and Poisson's ratio $\nu = 0.30$, the Lamé parameters are

$$ \mu(x; \lambda) = \frac{Y(x; \lambda)}{2(1 + \nu)}, \quad \lambda_L(x; \lambda) = \frac{\nu Y(x; \lambda)}{(1 + \nu)(1 - 2\nu)}. \qquad (30) $$
<figure><img src="image_14.png" /><figcaption>Fig. 12. Interpolation test for the parameter $Y_{ratio} = 0.975$, $\rho_{ratio} = 1.625$ at time $t = 0.1$. Left to right: Ground truth, PGROM+cVAE, direct cVAE, POD-DL-ROM.</figcaption></figure>

Material stiffness and density are parameterized by a circular inclusion of center $c = (0.5, 0.5)$ and radius $r_{in} = 0.20$. With parameter vector $\lambda = (Y_{ratio}, \rho_{ratio})$ and baselines $Y_{base} = 10^3$, $\rho_{base} = 1$,

$$ Y(x; \lambda) = \begin{cases} Y_{\text{base}} Y_{\text{ratio}}, & \|x - c\| \le r_{\text{in}}, \\ Y_{\text{base}}, & \text{elsewhere,} \end{cases} \qquad \rho(x; \lambda) = \begin{cases} \rho_{\text{base}} \rho_{\text{ratio}}, & \|x - c\| \le r_{\text{in}}, \\ \rho_{\text{base}}, & \text{elsewhere.} \end{cases} \qquad (31) $$

Boundary and initial conditions are

$$ u = 0 \qquad \text{on } \Gamma_D := \{x = (x_1, x_2) \in \partial\Omega : x_1 = 0\}, \qquad (32) $$

$$ \sigma(u) \mathbf{n} = \tau(x, t) \qquad \text{on } \Gamma_N^\tau := \{(x_1, x_2) \in \partial\Omega : 0.4 \le x_1 \le 0.6, x_2 = 1\}, \qquad (33) $$

$$ \sigma(u) \mathbf{n} = 0 \qquad \text{on } \Gamma_N^0 := \partial\Omega \setminus (\Gamma_D \cup \Gamma_N^\tau), \qquad (34) $$

$$ u(x, 0) = 0, \quad \partial_t u(x, 0) = 0 \qquad \text{in } \Omega, \qquad (35) $$

where $\mathbf{n}$ denotes the outward unit normal on $\partial\Omega$, and the prescribed vertical traction is $\tau(x, t) = (0, s(t))^{\top}$ with

$$ s(t) = \begin{cases} \sin(\pi t/T_p), & 0 \le t \le T_p, \\ 0, & t > T_p, \end{cases} \qquad T_p = 0.10. $$

As an example, Fig. 11 displays $u_x$, $u_y$, and the strain components $\epsilon_{xx}$, $\epsilon_{yy}$, and $\epsilon_{xy}$ for $Y_{ratio} = 5.0$ and $\rho_{ratio} = 1.0$.

We denote the parameter by $\lambda = (Y_{ratio}, \rho_{ratio})$ and construct three Cartesian-product grids. The training set uses a uniform step of 0.05 on [0.60, 2.00],

$$ \Lambda_{\text{train}} = \{(Y_{\text{ratio}}, \rho_{\text{ratio}}) : Y_{\text{ratio}} \in \{0.60, 0.65, \dots, 2.00\}, \rho_{\text{ratio}} \in \{0.60, 0.65, \dots, 2.00\}\}. $$

The interpolation test adopts an offset grid strictly within the training range to avoid overlap,

$$ \Lambda_{\text{test}}^{\text{int}} = \{(Y_{\text{ratio}}, \rho_{\text{ratio}}) : Y_{\text{ratio}} \in \{0.625, 0.675, \dots, 1.975\}, \rho_{\text{ratio}} \in \{0.625, 0.675, \dots, 1.975\}\}, $$

and the extrapolation test probes generalization beyond the training bounds,

$$ \Lambda_{\text{test}}^{\text{ext}} = \{(Y_{\text{ratio}}, \rho_{\text{ratio}}) : Y_{\text{ratio}} \in \{0.40, 2.20, 2.30, 2.40, 2.50\}, \rho_{\text{ratio}} \in \{0.40, 2.20, 2.30, 2.40, 2.50\}\}. $$

This subsection assesses the proposed PGROM + cVAE against two data-driven ROMs.

* Direct cVAE: A cVAE that directly maps $(t, \lambda)$ to the full snapshot without using POD-Galerkin; its network and training follow the cVAE setup in Section 3, except that the model's input and output are the entire snapshot rather than the PGROM residual.

* POD-DL-ROM [25]: This approach first utilizes POD as a pre-processing step to compress the high-fidelity solution snapshots into a set of intrinsic coordinates of a linear trial manifold. An AE then learns a further dimensionality reduction of these POD generalized coordinates to a low-dimensional latent space. A fully-connected NN is then trained to predict the latent code from the time-parameter pair.

All methods are trained on $\Lambda_{\text{train}}$ and evaluated on both $\Lambda_{\text{test}}^{\text{int}}$ and $\Lambda_{\text{test}}^{\text{ext}}$. Errors are reported over all $(t, \lambda)$ for $u_x$ and $u_y$.
<figure><img src="image_15.png" /><figcaption>**Fig. 13.** Extrapolation test for the parameter $Y_{\text{ratio}} = 2.4$, $\rho_{\text{ratio}} = 2.2$ at time $t = 0.05$. Left to right: Ground truth, PGROM + cVAE, direct cVAE, POD-DL-ROM.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 14.** Comparison of the $E_{\text{gap}}$ error for $u_x$ and $u_y$ on the interpolation test set ($\Lambda_{\text{test}}^{\text{int}}$) and the extrapolation test set ($\Lambda_{\text{test}}^{\text{ext}}$).</figcaption></figure>

Figs. 12–13 show examples of interpolation and extrapolation, respectively. In interpolation (Fig. 12), PGROM + cVAE is closest to the ground truth. POD-DL-ROM displays moderate phase lag and amplitude bias, while direct cVAE degrades notably via oversmoothing. In the extrapolation case (Fig. 13), PGROM + cVAE still maintains high accuracy and acceptable error, whereas both direct cVAE and POD-DL-ROM exhibit markedly larger discrepancies outside the training range.

Fig. 14 gives the $E_{\text{gap}}$ error on the full interpolation and extrapolation test sets. In interpolation, PGROM + cVAE attains the lowest average errors: 0.07 ($u_x$) and 0.02 ($u_y$), outperforming POD-DL-ROM (0.52 for $u_x$, 0.15 for $u_y$) and direct cVAE (3.65 for $u_x$, 2.56 for $u_y$). In extrapolation, PGROM + cVAE increases modestly to 0.11 ($u_x$) and 0.03 ($u_y$), whereas POD-DL-ROM rises to 16.7 and 5.85, and direct cVAE remains high at 3.65 and 2.55.

We conclude with a brief comparison of the proposed hybrid PGROM + cVAE against data-driven ROMs:

* **Advantages: improved accuracy and robustness.** As evidenced by the quantitative results above, our hybrid approach attains consistently lower errors across the tests. The gap is most pronounced in the extrapolation setting, where data-driven ROMs tend to deteriorate sharply outside the training range, whereas PGROM + cVAE remains reliable. This stems from the division of labor: the PGROM supplies a physics-consistent baseline governed by the equations, while the cVAE only learns the residual correction rather than the full state, simplifying the learning task and improving generalization.

* **Limitation: higher computational cost.** For fairness we matched the number of training iterations, making the offline training effort comparable to that of data-driven ROMs. In the online phase, however, the data-driven ROMs produce a new solution via a single forward pass with negligible cost, whereas PGROM + cVAE must still integrate the reduced ODE system to obtain the PGROM baseline before adding the cVAE correction. Consequently, its online cost is higher than an end-to-end data-driven ROM.

### 4.4. Application to an inverse problem in the heat equation

This subsection demonstrates the application of the proposed PGROM + cVAE framework to an inverse problem in PDEs. Specifically, we consider the heat equation as defined in Section 4.1. To simplify the problem setup, we focus on exploring the parameter $\kappa_{\text{in}}$ while keeping the other parameter, $r_{\text{in}}$, fixed at a constant value of 0.17. The objective of the inverse problem is to recover the unknown thermal conductivity parameter, $\kappa_{\text{in}}$, from a sparse set of temperature observations at specific spatial locations and time instances. Specifically, we consider nine distinct spatial observation points: [0.25, 0.25], [0.25, 0.5], [0.25, 0.75], [0.5, 0.25], [0.5, 0.5],
<figure><img src="image_17.png" /><figcaption>Fig. 15. Comparison of temperature field snapshots for different $κ_{in}$ parameters in the heat equation inverse problem. Top: True temperature field for the ground truth parameter $κ_{in}$ = 5.80. Middle: Temperature field predicted by the standard PGROM based on its identified optimal parameter $κ_{in}$ = 3.25. Bottom: Temperature field predicted by the PGROM+cVAE method based on its identified optimal parameter $κ_{in}$ = 5.75. Green dots indicate the spatial observation points used for parameter identification. Columns correspond to times $t$ = 0.1, 0.2, and 0.3, respectively.</figcaption></figure>

[0.5, 0.75], [0.75, 0.25], [0.75, 0.5], [0.75, 0.75], and three time instances: $t$ = 0.1, 0.2, and 0.3. These observation points (marked points) are visualized in Fig. 15. The ground truth value for the unknown parameter we aim to identify is $κ_{in}$ = 5.80.

We use a simple parameter sweep strategy to solve this inverse problem. A series of candidate values for $κ_{in}$ are selected over the interval [3.0, 10.0]. For each candidate value $κ_{in}$, a forward simulation of the heat equation is performed using the respective model (either the standard PGROM or the PGROM+cVAE) to generate a set of predicted snapshots. These predicted snapshots are then compared against the “true” snapshots (generated using the known ground truth parameter $κ_{in}$ = 5.80), and the Mean Squared Error (MSE) is computed. The candidate $κ_{in}$ value that yields the minimum MSE is considered the optimal estimate identified by that model. This parameter identification strategy is applied independently using both the standard PGROM and the PGROM+cVAE framework to compare their efficacy.

Fig. 16 compares the MSE profiles obtained for different candidate $κ_{in}$ values using both the PGROM and the PGROM+cVAE approaches. The standard PGROM identifies an optimal value of $κ_{in}$ = 3.25, where the MSE is approximately 0.122. In contrast, the PGROM+cVAE framework identifies an optimal value of $κ_{in}$ = 5.75. This estimate is significantly closer to the true parameter value ($κ_{in}$ = 5.80) and achieves a substantially lower MSE of approximately 0.011.

Fig. 15 visually compares the true temperature field with the predicted temperature fields obtained by the PGROM and the PGROM+cVAE. A comparison reveals that the temperature field predicted by the standard PGROM exhibits significant discrepancies from the true temperature field. In contrast, the temperature field predicted by the PGROM+cVAE method shows a much closer agreement with the true temperature field. This result demonstrates that the PGROM+cVAE framework offers higher accuracy in solving such inverse problems.
<figure><img src="image_18.png" /><figcaption>**Fig. 16.** Comparison of Mean Squared Error (MSE) for varying candidate $\kappa_{in}$ values in the inverse problem of identifying the unknown parameter in the heat equation. Left: MSE obtained using the standard PGROM. Right: MSE obtained using the proposed PGROM + cVAE framework. The ground truth value $\kappa_{in} = 5.80$ is indicated by the dashed vertical line. The optimal $\kappa_{in}$ identified by each method is marked with a star symbol.</figcaption></figure>

## 5. Conclusions

This paper presented a hybrid computational framework designed to improve the accuracy of PGROM for parameterized PDEs. We proposed correcting the PGROM approximation using a data-driven approach. Specifically, a cVAE was trained during an offline phase to learn the mapping from PDE parameters and time instances to the corresponding PGROM residual errors, using snapshots generated from high-fidelity simulations. In the online phase, the trained cVAE decoder efficiently predicts the residual for new parameters, which is then added to the fast PGROM solution to yield a significantly more accurate hybrid prediction. The effectiveness of this PGROM+cVAE methodology has been demonstrated through numerical experiments on three different parameterized PDE problems: a heat dissipation problem with varying thermal conductivity, a pollutant transport problem governed by convection-diffusion with varying velocity, and a plane-strain elastodynamics problem with a circular inclusion whose stiffness and density are parameterized.

We identify three primary directions for future research. Firstly, the scope of problems can be broadened; while this study focused on parameterized PDEs, the core concept of correcting a baseline approximate solution with a data-driven model holds potential for a wider range of scientific problems beyond PDEs where approximations introduce systematic errors. Secondly, the choice of the underlying approximation method could be diversified. Although we utilized PGROMs, the correction strategy could complement other model reduction techniques, such as those based on the Reduced Basis method or Proper Generalized Decomposition, by adapting to their distinct error structures and potentially enhancing predictive capabilities. Thirdly, the machine learning component itself offers avenues for advancement. Exploring enhancements to the VAE framework-such as employing different variants (like $\beta$-VAE for potentially more disentangled latent spaces), incorporating attention mechanisms to better capture complex residual structures, or integrating physics-informed regularization into the loss function-could lead to further improvements in accuracy.

### CRediT authorship contribution statement

Dawen Wu: Writing – original draft, Software, Methodology, Formal analysis, Conceptualization; Ludovic Chamoin: Writing – review & editing, Supervision, Project administration, Methodology, Investigation, Funding acquisition, Conceptualization.

### Data availability

Data will be made available on request.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgements

This research is part of the programme DesCartes and is supported by the National Research Foundation, Prime Minister's Office, Singapore under its Campus for Research Excellence and Technological Enterprise (CREATE) programme.

### References

[1] F. Chinesta, A. Huerta, G. Rozza, K. Willcox, Model Reduction Methods, John Wiley & Sons, Ltd, 2017, pp. 1–36. https://doi.org/10.1002/9781119176817.ecm2110

[2] P. Benner, S. Gugercin, K. Willcox, A survey of projection-based model reduction methods for parametric dynamical systems, SIAM Rev. 57 (4) (2015) 483–531.
[3] A. Quarteroni, A. Manzoni, F. Negri, Reduced Basis Methods for Partial Differential Equations: an Introduction, 92, Springer, 2015.

[4] D.E. Rumelhart, G.E. Hinton, R.J. Williams, et al., Learning internal representations by error propagation, 1985.

[5] D.P. Kingma, Auto-encoding variational bayes, (2013). arXiv:1312.6114

[6] J. Lucas, G. Tucker, R. Grosse, M. Norouzi, Understanding posterior collapse in generative latent variable models, ICLR Workshop on Deep Generative Models for Highly Structured Data, (2019).

[7] I. Higgins, L. Matthey, A. Pal, C. Burgess, X. Glorot, M. Botvinick, S. Mohamed, A. Lerchner, Beta-vae: learning basic visual concepts with a constrained variational framework, in: International Conference on Learning Representations, 2017.

[8] C.P. Burgess, I. Higgins, A. Pal, L. Matthey, N. Watters, G. Desjardins, A. Lerchner, Understanding disentangling in β-VAE, 2018, arXiv:1804.03599

[9] K. Sohn, H. Lee, X. Yan, Learning structured output representation using deep conditional generative models, in: C. Cortes, N. Lawrence, D. Lee, M. Sugiyama, R. Garnett (Eds.), Advances in Neural Information Processing Systems, 28, Curran Associates, Inc., 2015. https://proceedings.neurips.cc/paper_files/paper/2015/file/8d55a249e6baa5c06772297520da2051-Paper.pdf.

[10] A. Vahdat, J. Kautz, NVAE: A deep hierarchical variational autoencoder, Adv. Neural Inf. Process. Syst. 33 (2020) 19667–19679.

[11] D.P. Kingma, T. Salimans, R. Jozefowicz, X. Chen, I. Sutskever, M. Welling, Improved variational inference with inverse autoregressive flow, in: Proceedings of the 30th International Conference on Neural Information Processing Systems, NIPS’16, Curran Associates Inc., Red Hook, NY, USA, 2016, p. 4743–4751.

[12] S. Zhao, J. Song, S. Ermon, InfoVAE: balancing learning and inference in variational autoencoders, Proc. AAAI Conf. Artif. Intelligence 33 (01) (2019) 5885–5892. https://doi.org/10.1609/aaai.v33i01.33015885

[13] A. Razavi, A. van den Oord, O. Vinyals, Generating diverse high-Fidelity images with VQ-VAE-2, in: Advances in Neural Information Processing Systems, 2019.

[14] Q. Li, X. Wang, Y. Zhang, Anomaly detection for industrial time series data using variational autoencoders, IEEE Trans. Ind. Inf. 15 (5) (2019) 2852–2861.

[15] W. Jin, R. Barzilay, T. Jaakkola, Junction tree variational autoencoder for molecular graph generation, in: International Conference on Machine Learning, PMLR, 2018, pp. 2323–2332.

[16] D. Hafner, T. Lillicrap, I. Fischer, R. Villegas, D. Ha, H. Lee, J. Davidson, Learning latent dynamics for planning from pixels, in: International Conference on Machine Learning, PMLR, 2019, pp. 2555–2565.

[17] E. Saetta, R. Tognaccini, G. Iaccarino, Uncertainty quantification in autoencoders predictions: applications in aerodynamics, J. Comput. Phys. 506 (2024) 112951.

[18] E. Saetta, R. Tognaccini, G. Iaccarino, Machine learning to predict aerodynamic stall, Int. J. Comput. Fluid Dyn. 36 (7) (2022) 641–654.

[19] B. Kramer, B. Peherstorfer, K.E. Willcox, Learning nonlinear reduced models from data with operator inference, Annu. Rev. Fluid Mech. 56 (1) (2024) 521–548.

[20] S. Lee, S. Lee, K. Jang, H. Cho, S. Shin, Data-driven nonlinear parametric model order reduction framework using deep hierarchical variational autoencoder, Eng. Comput. (2024) 1–16.

[21] X. Zhang, L. Jiang, Conditional variational autoencoder with gaussian process regression recognition for parametric models, J. Comput. Appl. Math. 438 (2024) 115532.

[22] K. Lee, K.T. Carlberg, Model reduction of dynamical systems on nonlinear manifolds using deep convolutional autoencoders, J. Comput. Phys. 404 (2020) 108973.

[23] A. Solera-Rico, C. Sanmiguel Vila, M. Gómez-López, Y. Wang, A. Almashjary, S.T.M. Dawson, R. Vinuesa, β-Variational Autoencoders and transformers for reduced-order modelling of fluid flows, Nat. Commun. 15 (1) (2024) 1361.

[24] S. Fresca, L. Dede', A. Manzoni, A comprehensive deep learning-based approach to reduced order modeling of nonlinear time-dependent parametrized PDEs, J Sci Comput 87 (2) (2021) 61.

[25] S. Fresca, A. Manzoni, POD-DL-ROM: Enhancing deep learning-based reduced order models for nonlinear parametrized PDEs by proper orthogonal decomposition, Comput. Methods Appl. Mech. Eng. 388 (2022) 114181. https://doi.org/10.1016/j.cma.2021.114181

[26] S. Brivio, S. Fresca, A. Manzoni, PTPI-DL-ROMs: Pre-trained physics-informed deep learning-based reduced order models for nonlinear parametrized PDEs, Comput. Methods Appl. Mech. Eng. 432 (2024) 117404.

[27] F. Chinesta, E. Cueto, E. Abisset-Chavanne, J.L. Duval, F.E. Khaldi, Virtual, digital and hybrid twins: a new paradigm in data-based engineering and engineered data, Arch. Comput. Methods Eng. 27 (2020) 105–134.

[28] D. Di Lorenzo, V. Champaney, C. Germoso, E. Cueto, F. Chinesta, Data completion, model correction and enrichment based on sparse identification and data assimilation, Appl. Sci. 12 (15) (2022) 7458.

[29] C. Ghnatios, F. Chinesta, Discovering PDEs corrections from data within a hybrid modeling framework, Mathematics 13 (1) (2024) 5.

[30] M. Bonavita, P. Laloyaux, Machine learning for model error inference and correction, J. Adv. Model. Earth Syst. 12 (12) (2020) e2020MS002232.

[31] F. Galarce, J. Mura, A. Caiazzo, Bias and multiscale correction methods for variational state estimation, Appl. Math. Model. 138 (2025) 115761.

[32] A. Nóvoa, A. Racca, L. Magri, Inferring unknown unknowns: regularized bias-aware ensemble kalman filter, Comput. Methods Appl. Mech. Eng. 418 (2024) 116502.

[33] C. Ghnatios, X. Kestelyn, G. Denis, V. Champaney, F. Chinesta, Learning data-driven stable corrections of dynamical systems-application to the simulation of the top-oil temperature evolution of a power transformer, Energies 16 (15) (2023) 5790.

[34] C. Ghnatios, D. Jiang, Y. Tourbier, A. Cimetière, F. Chinesta, Data-Completion and model correction by means of evanescent regularization, Appl. Sci. 13 (17) (2023) 9616.