Disentangled Representation Learning for Parametric Partial Differential Equations

## Abstract

Abstract Neural operators (NOs) have demonstrated remarkable success in learning mappings between function spaces, serving as efficient approximators for the forward solutions of complex physical systems governed by partial differential equations (PDEs). However, while effective as black-box solvers, they offer limited insight into the underlying physical mechanism, due to the lack of interpretable representations of the physical parameters that drive the system. To tackle this challenge, we propose a new paradigm for learning disentangled representations from neural operator parameters, thereby effectively solving an inverse problem. Specifically, we introduce DisentangO, a novel hyper-neural operator architecture designed to unveil and disentangle the latent physical factors of variation embedded within the black-box neural operator parameters. At the core of DisentangO is a multi-task neural operator architecture that distills the varying parameters of the governing PDE through a task-wise adaptive layer, coupled with a hierarchical variational autoencoder that disentangles these variations into identifiable latent factors. By learning these disentangled representations, our model not only enhances physical interpretability but also enables more robust generalization across diverse physical systems. Empirical evaluations across supervised, semi-supervised, and unsupervised learning contexts show that DisentangO effectively extracts meaningful and interpretable latent features, bridging the divide between predictive performance and physical understanding in neural operator frameworks.

## 1 Introduction

Interpretability in machine learning (ML) refers to the ability to understand and explain how models make decisions (Rudin et al., ; Molnar, ). As ML systems grow more complex, especially with the use of deep learning and ensemble methods (Sagi and Rokach, ), the reasoning behind their predictions can become opaque. Interpretability addresses this challenge by making the models more transparent, enabling users to trust the outcomes, detect biases, and identify potential flaws. It is a critical factor in applying AI responsibly, especially in fields where accountability and fairness are essential (Cooper et al., ). In physics, where models endeavor to capture the governing laws and physical principles, understanding how a model arrives at its predictions is vital for verifying that it aligns with known scientific theories. Interpretable models also allow researchers to gain insights into the learned relationships, ensuring that the predictions are physically meaningful and consistent with real-world behaviors. This transparency is key for advancing scientific discovery, validating results, and enhancing trust in models used for complex physical systems.

Discovering interpretable representations of physical parameters in learning physical systems is challenging due to the intricate nature of real-world phenomena and the often implicit relationships between variables. In physics, parameters like force, energy, and velocity are governed by well-established laws, and extracting these parameters in a way that aligns with physical intuition requires models that go beyond mere pattern recognition. Traditional ML models may fit the data but fail to represent these parameters in a physically interpretable way. To address this, recent developments include the integration of physical constraints into the learning process (Raissi et al., ), such as embedding conservation laws (Liu et al., , ), symmetries (Mattheakis et al., ), or invariances (Cohen and Welling, ) directly into model architectures. Additionally, methods like symbolic regression (Biggio et al., ) and sparse modeling (Carroll et al., ) aim to discover simple, interpretable expressions that capture the underlying dynamics. However, balancing model complexity with interpretability remains a significant challenge, as overly complex models may obscure the true physical relationships, while oversimplified models risk losing critical details of the system’s behavior.

In this work, we introduce DisentangO, a novel variational hierarchical neural operator architecture to disentangle physical factors of variation from black-box neural operator parameters for solving parametric PDEs. Neural operators (NOs) (Li et al., , ) learn mappings between infinite-dimensional function spaces in the form of integral operators, making them powerful tools for discovering continuum physical laws by manifesting the mappings between spatial and/or spatiotemporal data; see You et al. (); Liu et al. (, , ); Ong et al. (); Cao (); Lu et al. (, ); Goswami et al. (); Gupta et al. () and references therein. However, most NOs play the role of efficient forward surrogates for the underlying physical system under a supervised learning setting. As a result, they act as black-box universal approximators for a single physical system, and are lack of interpretability with respect to the underlying physical laws. In contrast, the key innovation of DisentangO lies in its use of a hierarchical architecture that distills the varying physical parameters of the governing PDE from multiple physical systems through an unsupervised learning setting. The distilled variations are further disentangled into distinct physical factors to enhance physical understanding and promote robust generalization. Consequently, DisentangO effectively extracts meaningful and interpretable physical features, thereby simultaneously solving both the forward and inverse problems. Our key contributions are:

- •
We bridge the divide between predictive accuracy and physical interpretability, and introduce a new paradigm that simultaneously performs physics modeling (i.e., as a forward PDE solver) and governing physical mechanism discovery (i.e., as an inverse PDE solver).
- •
We propose a novel variational hierarchical neural operator architecture, which we coin DisentangO. DisentangO extracts the key physical factors of variation from black-box neural operator parameters of multiple physical systems. These factors are then disentangled into distinct latent factors that enhance physical interpretation and robust generalization.
- •
We provide theoretical analysis on the component-wise identifiability of the true generative factors in physics modeling: by learning from multiple physical systems, the variability of hidden physical states in these systems promotes identifiability.
- •
We explore the practical utility of disentanglement rather than pursuing it as an end goal, and perform experiments across a broad range of settings including supervised, semi-supervised, and unsupervised learning. Our results show that DisentangO effectively extracts meaningful and interpretable physical features for enhanced interpretability.

## 2 Background and related work

Neural operators.
Learning complex physical systems directly from data is crucial in many scientific and engineering applications (Ghaboussi et al., ; Liu et al., ; Ghaboussi et al., ; Carleo et al., ; Karniadakis et al., ; Zhang et al., ; Cai et al., ; Pfau et al., ; He et al., ; Jafarzadeh et al., ). Often, the governing laws of these systems are unknown and must be inferred from data through learning models that are both *resolution-invariant* to ensure consistent performance across various levels of discretizations and *interpretable* to domain experts. Neural operators (NOs) address the first aim by learning mappings between infinite-dimensional function spaces (Li et al., , , ; You et al., ; Ong et al., ; Cao, ; Lu et al., , ; Goswami et al., ; Gupta et al., ). By capturing the relationship between spatial and spatio-temporal data, NOs offer a promising approach to discovering continuum physical surrogates, enabling more accurate and consistent predictions of the underlying responses of complex systems. Nevertheless, NOs lack an interpretable representation of the physical states driving the system, and as a result offer limited insight into the underlying physical mechanisms.

Forward and inverse PDE learning.
Existing NOs focus on providing efficient surrogates for the underlying physical systems as forward PDE solvers. They are often employed as black-box universal approximators but lack interpretability of the underlying physical principles. On the other hand, several deep learning methods have been developed as inverse PDE solvers (Fan and Ying, ; Molinaro et al., ; Jiang et al., ; Chen et al., ), aiming to reconstruct the PDE parameters from solution data. However, the inverse problem is typically more challenging due to its ill-posedness. To tackle this pathology, many NOs incorporate prior information via governing PDEs (Yang et al., ; Li et al., ), regularizers (Dittmer et al., ; Obmann et al., ; Ding et al., ; Chen et al., ), or additional operator structures (Uhlmann, ; Lai et al., ; Yilmaz, ). However, all these work rely on prior information on the operator form, which is not realistic in many real-world problems where neither the parameter field nor the model form are available. Worth noting is a recent work that leverages the attention mechanism to build a nonlinear kernel map in NO (Yu et al., ), with which a data-dependent kernel can be constructed to characterize the inverse mapping from data to the hidden governing mechanism. While effective, this approach does not disentangle the learned kernel or extract physically interpretable parameters. To the best of our knowledge, our proposed DisentangO is the first that tackles both forward (i.e., physics prediction) and inverse (i.e., physics discovery) PDE learning problems, while simultaneously inferring the governing distinct physical parameters from the learned NO.

Disentangled representation learning.
Disentangled representation learning seeks to separate data into distinct, interpretable factors, where each factor captures an independent aspect of variation in the input. This area has critical implications for transfer learning, generative modeling, and fairness in AI, as it allows models to generalize across tasks by leveraging isolated features, making them more interpretable and modular. Disentangled representation learning has seen significant advances, primarily driven by models like $\beta$-VAE (Higgins et al., ), FactorVAE (Kim and Mnih, ), and InfoGAN (Chen et al., ), who encourage disentanglement by controlling the trade-off between reconstruction accuracy and latent regularization, or introduce mutual information constraints to encourage the discovery of interpretable latent factors. Early work in disentanglement largely focus on supervised or semi-supervised approaches (Ridgeway and Mozer, ; Shu et al., ; Mathieu et al., ), requiring labeled data for factors of variation. However, more recent advancements have concentrated on unsupervised learning methods (Duan et al., ), where disentanglement is achieved without explicit labels. Despite these efforts, fully unsupervised disentanglement remains challenging (Locatello et al., ), with ongoing research focusing on improving metrics, robustness, and applicability to more complex and high-dimensional data. In this context, although a large body of research has been done in disentanglement, its applications remain largely confined to fields like computer vision and robotics, where the latent dimensions typically correspond to human-interpretable (often visually recognizable) features, with very limited exploration in areas such as physical system learning. Moreover, existing work focuses on disentangling from data, whereas our work is the first attempt that explores disentanglement from black-box neural network parameters.

## 3 DisentangO: disentangled representation learning from neural operator parameters

We consider a series of complex systems with different hidden physical parameters:

$$ $\mathcal{K}_{\bm{b}}[\bm{u}](\bm{x})=\bm{f}(\bm{x})\text{ ,}\quad\bm{x}\in{ \Omega}\text{ .}$ (3.1) $$

Here, ${\Omega}\subset\mathbb{R}^{s}$ is the domain of interest, $\bm{f}(\bm{x})$ is a function representing the loading on ${\Omega}$, and $\bm{u}(\bm{x})$ is the corresponding solution of this system. $\mathcal{K}_{\bm{b}}$ represents the unknown governing law, e.g., balance laws, which is determined by the (possibly unknown and high-dimensional) parameter field $\bm{b}$. For instance, in a material modeling problem, $\mathcal{K}_{\bm{b}}$ often stands for the constitutive law and $\bm{b}$ can be a vector ($\bm{b}\in\mathbb{R}^{d_{b}}$) representing the homogenized material parameter field or a vector-valued function ($\bm{b}\in L^{\infty}({\Omega};\mathbb{R}^{d_{b}})$) representing the heterogeneous material properties. Both scenarios are considered in our empirical experiments in Section 4.

Many physical modeling tasks can be formulated as either a forward or an inverse PDE-solving problem. In a forward problem setting, one aims to find the PDE solution when given PDE information, including coefficient functions, boundary conditions, initial conditions, and loading sources. That means, given the governing operators $\mathcal{K}$, the parameter (field) $\bm{b}$, and loading field $\bm{f}(\bm{x})$ in equation 3.1, the goal is to solve for the corresponding solution field $\bm{u}(\bm{x})$, through classical PDE solvers (Brenner and Scott, ) or data-driven approaches (Lu et al., ; Li et al., ; You et al., ). As a result, a forward map is constructed:

$$ $\mathcal{G}:(\bm{b},\bm{f})\rightarrow\bm{u}\text{ .}$ (3.2) $$

Here, $\bm{b}$ and $\bm{f}$ are input vectors/functions, and $\bm{u}$ is the output function.

Conversely, solving an inverse PDE problem involves reconstructing the underlying full or partial PDE information from PDE solutions, where one seeks to construct an inverse map:

$$ $\mathcal{H}:(\bm{u},\bm{f})\rightarrow\bm{b}\text{ .}$ (3.3) $$

Unfortunately, solving an inverse problem is typically more challenging due to the ill-posed nature of the PDE model. In general, a small number of function pairs $(\bm{u},\bm{f})$ from a single system does not suffice in inferring the underlying parameter field $\bm{b}$, making the inverse problem generally non-identifiable (Molinaro et al., ). Herein, we propose a novel neural architecture to alleviate the curse of ill-posedness without requiring prior knowledge. The key ingredients are: 1) the construction of a multi-task NO architecture, which solves both the forward (physics prediction) and inverse (physics discovery) problems simultaneously from multiple systems, and 2) a generative model to disentangle the features from NO parameters that contain critical information of $\bm{b}$.

### 3.1 Variational hierarchical neural operator as a multi-task solver

Figure: Figure 1: Illustration of DisentangO’s architecture.
Refer to caption: x1.png

Formally, we consider $S$ training datasets from different systems, each of which contains $n_{train}$ function pairs $(\bm{u},\bm{f})\in\mathcal{U}\times\mathcal{F}$:

$$ $\mathcal{D}_{\rm tr}=\{\{(\bm{u}^{\eta}_{i}(\bm{x}),\bm{f}^{\eta}_{i}(\bm{x})) \}_{i=1}^{n_{train}}\}_{\eta=1}^{S}\text{ .}$ (3.4) $$

The index $\eta$ corresponds to different hidden parameter field/set $\bm{b}^{\eta}\in\mathcal{B}$. Here, $\mathcal{U}\subset L^{2}({\Omega};\mathbb{R}^{d_{u}})$, $\mathcal{F}\subset L^{2}({\Omega};\mathbb{R}^{d_{f}})$ represent the (infinite-dimensional) Banach spaces of the solution function $\bm{u}$ and loading function $\bm{f}$, respectively, and $\mathcal{B}$ is a finite-dimensional manifold where the hidden parameter $\bm{b}$ takes values on. We further assume that the data satisfies the following generating process:

$$ $\bm{b}|\bm{z}=m(\bm{z})\text{ ,}\;\bm{z}\sim\mathbb{P}_{z}\text{ ,}$ (3.5) $$

where $\bm{z}\in\mathcal{Z}\subset\mathbb{R}^{d_{z}}$ is the embedding of $\bm{b}$ in the latent space and $m:\mathcal{Z}\rightarrow\mathcal{B}$ is a one-to-one function.
Then, we assume for each $\bm{b}^{\eta}$ that there exists an underlying parametric PDE system of equation 3.1. Hence, each function pair $(\bm{u}_{i}^{\eta},\bm{f}^{\eta}_{i})$ satisfies:

$$ $\bm{u}|\bm{f}=\mathcal{K}_{\bm{b}}^{-1}[\bm{f}]\text{ ,}\;\bm{f}\sim\mathbb{P} _{f}\text{ ,}$ (3.6) $$

where $\mathcal{K}_{\bm{b}}^{-1}$ is the ground-truth solution operator of the parametric PDE system of equation 3.1 with a given parameter $\bm{b}$.

Our first objective is to construct a surrogate operator of $\mathcal{K}_{\bm{b}}^{-1}:\mathcal{F}\rightarrow\mathcal{U}$, as the forward PDE solver for each system:

$$ $\mathcal{G}[\bm{f}^{\eta}_{i};\theta^{\eta}](\bm{x})\approx\mathcal{K}_{\bm{b} ^{\eta}}^{-1}[\bm{f}_{i}^{\eta}](\bm{x})=\bm{u}_{i}^{\eta}(\bm{x})\text{ ,}$ (3.7) $$

where $\theta^{\eta}$ denotes the corresponding trainable parameter set. Then, as the second objective, we construct an inverse map by discovering the hidden parameter field/set from $\theta^{\eta}$:

$$ $\mathcal{H}(\theta^{\eta};\Theta_{H})\approx\tilde{\bm{z}}^{\eta}\text{ .}$ (3.8) $$

Here, $\Theta_{H}$ is the trainable parameters of $\mathcal{H}$, and $\tilde{\bm{z}}^{\eta}\in\mathbb{R}^{d_{z}}$ denotes the learned latent variables, which can be transformed to the ground-truth latent variable $\bm{z}^{\eta}$ via an invertible function $h$. The goal of disentanglement is not to discover $\bm{b}^{\eta}$, since some features may not be identifiable. For instance, in a Dirichlet boundary condition problem, $\bm{u}^{\eta}_{i}(\bm{x})=0$ for all $\bm{x}\in\partial\Omega$, making $\bm{b}^{\eta}$ not learnable on $\partial\Omega$. Hence, a more realistic goal is to discover the underlying mechanism of $\bm{b}$ in the space of identifiability, i.e., $\bm{z}$, or its equivalent transformation $\tilde{\bm{z}}$.

Note that in the first objective of equation 3.7, one is interested in finding the solution operator $\mathcal{G}^{\eta}$ under a supervised setting, in the form of function-to-function mappings, for all (hidden) parameters $\bm{b}^{\eta}$ in the range of interest; while in the second objective of equation 3.8, the goal is to construct a vector-to-vector mapping from $\theta^{\eta}$ to $\tilde{\bm{z}}^{\eta}$, which is unsupervised due to the hidden physics nature. We propose to employ variational autoencoder (VAE) (Kingma and Welling, ) as the representation learning approach for the second objective, and a meta-learned NO architecture as a universal solution operator for the first objective. The key is to pair these architectures as a hierarchical neural operator, which simultaneously solves both forward and inverse modeling problems.

We now define the learning architecture. The overall objective is to maximize the log data likelihood:

$$ $\max E(\log(p(\theta,\mathcal{D})))=\max\left[E(\log(p(\bm{u}|\bm{f},\theta)+ \log(p(\theta|\bm{f}))+\log(p(\bm{f}))\right].$ $$

Note that $p(\bm{f})$ remains constant over different $\theta$, and the assumption in equation 3.6 yields $\log(p(\theta|\bm{f}))=\log(p(\theta))$. Additionally, the assumption in equation 3.5 guarantees that $\theta$ is generated from the latent space over $\bm{z}$, hence $\log(p(\theta))=\log\int_{\bm{z}}p(\theta|\bm{z})p(\bm{z})d\bm{z}$. The overall objective then becomes:

$$ $\displaystyle\max E(\log(p(\theta,\mathcal{D})))=$ $\displaystyle\max\left[E(\log(p(\bm{u}|\bm{f},\theta)))+E\left(\log\int_{\bm{z }}p(\theta|\bm{z})p(\bm{z})d\bm{z}\right)\right].$ (3.9) $$

However, the second term in this formulation is generally intractable. We use a variational posterior $q(\theta|\bm{z})$ to approximate the actual posterior $p(\theta|\bm{z})$ and maximize the evidence lower bound (ELBO):

$$ $\displaystyle L_{ELBO}=\dfrac{1}{S}\sum_{\eta=1}^{S}\left[E_{q(\bm{z}^{\eta}| \theta^{\eta})}\log p(\theta^{\eta}|\bm{z}^{\eta})-D_{KL}\left(q(\bm{z}^{\eta} |\theta^{\eta})||p(\bm{z}^{\eta})\right)\right].$ $$

Putting everything together, we obtain the loss functional:

$$ $L_{loss}=\dfrac{1}{S}\sum_{\eta=1}^{S}\left[-E(\log(p(\bm{u}|\bm{f},\theta^{ \eta})))-E_{q(\bm{z}^{\eta}|\theta^{\eta})}\log p(\theta^{\eta}|\bm{z}^{\eta}) +D_{KL}\left(q(\bm{z}^{\eta}|\theta^{\eta})||p(\bm{z}^{\eta})\right)\right].$ (3.10) $$

The above formulation naturally lends itself to a hierarchical variational autoencoder (HVAE) architecture (Vahdat and Kautz, ). The encoder aims to obtain the posterior $q_{\mu_{z},\Sigma_{z}}(\bm{z}^{\eta}|\theta^{\eta})$, which estimates the inverse map $\mathcal{H}$:

$$ $\hat{\bm{z}}\sim q_{\mu_{z},\Sigma_{z}}(\hat{\bm{z}}^{\eta}|\theta^{\eta}) \text{ .}$ (3.11) $$

Then, the first decoder, which processes the estimated latent variable $\hat{\bm{z}}$ and reconstructs the corresponding NO parameter $\hat{\theta}$:

$$ $\hat{\theta}=\hat{g}(\hat{\bm{z}})\text{ .}$ (3.12) $$

Lastly, the second decoder reconstructs the forward map $\mathcal{G}$, by further taking the estimated $\hat{\theta}$ and the loading function $\bm{f}$ and estimating the output function $\bm{u}$:

$$ $\hat{\bm{u}}=\hat{\mathcal{G}}[\bm{f};\hat{\theta}]\text{ .}$ (3.13) $$

For this purpose, we employ the implicit Fourier neural operator (IFNO)(^*^**We also point out that the proposed strategy is generic and hence also applicable to other neural operators.), which is a provably universal approximator for PDE solution operators with a relatively small number of trainable parameters (You et al., ). To provide a
universal $\hat{\mathcal{G}}$ for different domains, a multi-task learning strategy is adopted following Zhang et al. (), which applies task-wise adaptation only to the first layer of NOs. Specifically, an $L$-layer multi-task IFNO has the following form:

$$ ${\mathcal{G}}[\bm{f};{\theta}^{\eta}](\bm{x})={\mathcal{G}}[\bm{f};{\theta}_{P }^{\eta},\theta_{J},\theta_{Q}](\bm{x}):=\mathcal{Q}_{\theta_{Q}}\circ( \mathcal{J}_{\theta_{J}})^{L}\circ\mathcal{P}_{{\theta}^{\eta}_{P}}[\bm{f}]( \bm{x})\text{ ,}$ (3.14) $$

where $\mathcal{P}$, $\mathcal{Q}$ are shallow-layer MLPs that map a low-dimensional vector into a high-dimensional vector and vice versa, parameterized by $\theta_{Q}$ and ${\theta}^{\eta}_{P}$, respectively. Each intermediate layer, $\mathcal{J}$, is constructed as a mimetic of a fixed-point iteration step and parameterized by $\theta_{J}$. Supported by the universal approximator analysis in Zhang et al. (), different PDEs share the common iterative ($\theta_{J}$) and projection ($\theta_{Q}$) parameters, with all the information about parameter $\bm{b}$ embedded in the task-wise lifting parameters $\theta^{\eta}_{P}$. Hence, the encoder and first decoder structure only need to encode and reconstruct $\theta^{\eta}_{P}$, since all other NO parameters are invariant to the change of $\bm{z}$. This construction substantially reduces the degree of freedoms in $\theta^{\eta}$, making the invertibility assumption in the next section feasible. With a slight abuse of notation, in the following contents we use $\theta^{\eta}$ to denote the task-wise parameter, i.e., the lifting layer parameter $\theta^{\eta}_{P}$, and estimate $\bm{z}^{\eta}$ from it.

### 3.2 Disentangling the underlying mechanism

The proposed DisentangO aims to identify and disentangle the components of the latent representation $\bm{z}$, which serves as an inverse PDE solver. However, to capture the true physical mechanism, a natural question arises: is it really possible to identify the latent variables of interest (i.e., $\bm{z}$) with only observational data $\mathcal{D}_{\rm tr}$? We now show that, by learning the model $(p_{\hat{\bm{z}}},\hat{g},\hat{\mathcal{H}},\hat{\mathcal{G}})$ that matches the true marginal data distribution in all domains, we can indeed achieve this identifiability for the generating process proposed in equation 3.5 and equation 3.6.

First, we show that the latent variable $\bm{z}$ can be identified up to an invertible component-wise transformation: for the true latent variable $\bm{z}\in\mathbb{R}^{d_{z}}$, there exists an invertible function $h:\mathbb{R}^{d_{z}}\rightarrow\mathbb{R}^{d_{z}}$,
such that $\hat{\bm{z}}=h(\bm{z})$.

###### Theorem 1 .

We follow the data-generating process in equation 3.5 and equation 3.6, and make the following assumptions:

- •
(Density Smoothness and Positivity): $p_{\theta}$ and $p_{\bm{z}}$ are both smooth and positive.
- •
(Invertibility): The task-wise parameter $\theta$ can be generated by $\bm{z}$ through an invertible and smooth function $\mathcal{H}^{-1}$. Moreover, for each given $\bm{f}$, we denote $\mathcal{G}_{f}(\theta):=\mathcal{G}[\bm{f};\theta](\bm{x})$ as the operator mapping from $\mathbb{R}^{d_{\theta}}$ to $\mathcal{U}$. $\mathcal{G}_{f}$ is also one-to-one with respect to $\theta$.

Then, by learning $(p_{\hat{\bm{z}}},\hat{g},\hat{\mathcal{H}},\hat{\mathcal{G}})$ to achieve:

$$ $p_{\hat{\bm{u}}|\bm{f}}=p_{\bm{u}|\bm{f}}\text{ ,}$ (3.15) $$

where $\bm{u}$ and $\hat{\bm{u}}$ are generated from the true process and the estimated model, respectively, $\bm{z}$ is identifiable up to an invertible function $h$.

Proof: Please see Appendix A.

Additionally, with additional assumptions on conditional independence and datum variability, we can further obtain the following theoretical results on component-wise identifiability: for each true component $z_{i}$, there exists a corresponding estimated component $\hat{z}_{j}$ and an invertible function $h_{i}:\mathbb{R}\rightarrow\mathbb{R}$, such that $z_{i}=h_{i}(\hat{z}_{j})$.

###### Theorem 2 .

We follow the data-generating process in equation 3.5 and equation 3.6, and assume that the Density Smoothness and Positivity condition as well as the Invertibility condition in Theorem 1 are both satisfied, together with the following additional assumptions:

- •
(Conditional Independence): Conditioned on $\bm{b}$, each component of $\bm{z}$ is independent of each other: $\log p_{\bm{z}|\bm{b}}(\bm{z}|\bm{b})=\sum_{i=1}^{d_{z}}\log p_{z_{i}|\bm{b}}(
z_{i}|\bm{b})$.
- •
(Linear independence) There exists $2d_{z}+1$ values of $\bm{b}$, such that the $2d_{z}$ vectors $w(\bm{z},\bm{b}^{j})-w(\bm{z},\bm{b}^{0})$ with $j=1,\cdots,2d_{z}$ are linearly independent.
Here,
$w(\bm{z},\bm{b}^{j}):=\left(\dfrac{\partial q_{1}(z_{1},\bm{b}^{j})}{\partial z
_{1}},\cdots,\dfrac{\partial q_{d_{z}}(z_{d_{z}},\bm{b}^{j})}{\partial z_{d_{z
}}},\dfrac{\partial^{2}q_{1}(z_{1},\bm{b}^{j})}{\partial z_{1}^{2}},\cdots,
\dfrac{\partial^{2}q_{d_{z}}(z_{d_{z}},\bm{b}^{j})}{\partial z_{d_{z}}^{2}}
\right),$
(3.16)
with $q_{i}(z_{i},\bm{b}):=\log p_{z_{i}|\bm{b}}$.

Then, by learning $(p_{\hat{\bm{z}}},\hat{g},\hat{\mathcal{H}},\hat{\mathcal{G}})$ to achieve equation 3.15, $\bm{z}$ is component-wise identifiable.

Proof: Please see Appendix A.

Intuitively, the linear independence condition in Theorem 2 indicates sufficient variability between physical systems (tasks). This is a common assumption in the nonlinear ICA literature for domain adaptation (Hyvarinen et al., ; Khemakhem et al., ; Kong et al., ). To our best knowledge, this is the first time the component-wise identifiability results are discussed in the context of multi-task neural operator learning.

###### Theorem 1 .

###### Theorem 2 .

### 3.3 A generic algorithm for different supervision scenarios

In this section, we discuss the practical utility of DisentangO under three scenarios:

- •
(SC1: Unsupervised) No value or label is given for $\bm{b}^{\eta}$ for each task.
- •
(SC2: Semi-supervised) The value of $\bm{b}^{\eta}$ is not given, but a label $c(\bm{b}^{\eta})$ (e.g., a classification of $\bm{b}^{\eta}$) is given.
- •
(SC3: Supervised) The value of $\bm{b}^{\eta}$ is given.

The architecture is shown in Figure 1, with a pseudo algorithm provided in Algorithm 1 of Appendix.

To obtain the posterior $q$ in equation 3.11, we assume each latent variable satisfies a Gaussian distribution of distinct means and diagonal covariance, i.e., $q_{\mu_{z},\Sigma_{z}}(\hat{\bm{z}}|\theta):=\mathcal{N}(\mu_{z}(\theta_{P}),
\Sigma^{2}_{z}(\theta_{P}))$, then estimate its mean and covariance using an MLP. As a result, the KL-divergence term in the ELBO admits a closed form:

$$ $\text{SC1/SC2:}\;D_{KL}\left(q(\hat{\bm{z}}|\theta)||p(\bm{z})\right)=\dfrac{1 }{2}\sum_{i=1}^{d_{z}}\left((\Sigma_{z})_{i}^{2}+(\mu_{z})^{2}_{i}-2\log(( \Sigma_{z})_{i})-1\right)\text{ .}$ $$

In the supervised setting, we take $\bm{z}=\bm{b}$, and set $\mu_{z}(\theta_{P}^{\eta})=\bm{b}^{\eta}$, then:

$$ $\text{SC3:}\;D_{KL}\left(q(\hat{\bm{z}}|\theta)||p(\bm{z})\right)=\dfrac{1}{2} \sum_{i=1}^{d_{z}}\left((\Sigma_{z})_{i}^{2}+(\mu_{z}-\bm{b})^{2}_{i}-2\log(( \Sigma_{z})_{i})-1\right)\text{ .}$ $$

For the first decoder, the likelihood $p(\theta^{\eta}_{P}|\bm{z}^{\eta})$ is a factorized Gaussian with mean $\mu_{\theta}$ and covariance $\Sigma_{\theta}$, computed from another MLP.
By taking as input Monte Carlo samples once for each $\bm{z}$, the reconstruction accuracy term can be approximated as:

$$ $E_{q(\bm{z}|\theta_{P})}\log p(\theta_{P}|\bm{z})\approx-\sum_{i=1}^{d_{\theta }}\left(\dfrac{((\theta_{P})_{i}-(\mu_{\theta})_{i})^{2}}{2(\Sigma_{\theta})^{ 2}_{i}}+\log(\Sigma_{\theta})_{i}+c\right)\text{ ,}$ $$

where ${d_{\theta}}$ is the dimension of $\theta_{P}$, and $c$ is a constant.

For the second decoder, we parameterize it as a multi-task NO in equation 3.14, and model the discrepancy between $\hat{\mathcal{G}}_{[\hat{\theta}^{\eta}_{P},\theta_{I},\theta_{Q}]}[\bm{g}^{
\eta}_{j}](\bm{x}_{k})$ and the ground truth $\bm{u}(\bm{x}_{k})$ as an additive independent unbiased Gaussian random noise $\epsilon$, with

$$ $\hat{\mathcal{G}}_{[\hat{\theta}^{\eta}_{P},\theta_{I},\theta_{Q}]}[\bm{g}^{ \eta}_{j}](\bm{x}_{k}))=\bm{u}_{j}^{\eta}(\bm{x}_{k})+\epsilon_{\eta,j,k}\text { ,}\quad\epsilon_{\eta,j,k}\sim\mathcal{N}(0,\varpi^{2})\text{ ,}$ $$

where $\varpi$ is a tunable hyperparameter. Then, with a uniform spatial discretization of size $\Delta x$ in a domain ${\Omega}\subset\mathbb{R}^{d_{u}}$,
the log likelihood after eliminating the constant terms can be written as:

$$ $\dfrac{1}{2\varpi^{2}\Delta x^{2d_{u}}}\sum_{\eta=1}^{S}\sum_{j=1}^{n_{train}} {\left|\left|\hat{\mathcal{G}}_{[\hat{\theta}^{\eta}_{P},\theta_{I},\theta_{Q} ]}[\bm{g}^{\eta}_{j}](\bm{x}_{k})-\bm{u}_{j}^{\eta}(\bm{x}_{k})\right|\right|} _{L^{2}({\Omega})}^{2}.$ $$

In our empirical tests in Section 4, for the unsupervised and semi-supervised scenarios we choose $q$ as a standard Gaussian distribution following the independence assumption of Theorem 2. Additionally, to avoid overparameterization, we take $\Sigma_{\theta}=\sigma^{2}_{\theta}\mathbf{I}$, with $\sigma_{\theta}$ being a tunable hyperparameter. In this case, the total objective can be further simplified as:

$$ $\displaystyle L_{loss}=$ $\displaystyle\dfrac{1}{S}\sum_{\eta=1}^{S}\left(\beta_{d}\sum_{j=1}^{n_{train} }{\left|\left|\hat{\mathcal{G}}_{[\hat{\theta}^{\eta}_{P},\theta_{I},\theta_{Q }]}[\bm{g}^{\eta}_{j}]-\bm{u}^{\eta}_{j}\right|\right|}^{2}_{L^{2}({\Omega})}+ {\left|\left|\hat{\theta}^{\eta}_{P}-\mu^{\eta}_{\theta}\right|\right|}^{2}+ \beta_{KL}{\left|\left|\mu^{\eta}_{\bm{z}}\right|\right|}^{2}\right)\text{ ,}$ $$

for the unsupervised scenario, and

$$ $\displaystyle L_{loss}=$ $\displaystyle\dfrac{1}{S}\sum_{\eta=1}^{S}\left(\beta_{d}\sum_{j=1}^{n_{train} }{\left|\left|\hat{\mathcal{G}}_{[\hat{\theta}^{\eta}_{P},\theta_{I},\theta_{Q }]}[\bm{g}^{\eta}_{j}]-\bm{u}^{\eta}_{j}\right|\right|}^{2}_{L^{2}({\Omega})}+ {\left|\left|\hat{\theta}^{\eta}_{P}-\mu^{\eta}_{\theta}\right|\right|}^{2}+ \beta_{KL}{\left|\left|\mu^{\eta}_{\bm{z}}-\bm{b}^{\eta}\right|\right|}^{2} \right)\text{ ,}$ $$

for the fully supervised scenario. Here, $\beta_{d}:=\frac{\sigma^{2}_{\theta}}{\varpi^{2}\Delta x^{2d_{u}}}$ and $\beta_{KL}:=\sigma^{2}_{\theta}$ are treated as tunable hyperparameters. In the semi-supervised scenario, we incorporate partial supervision by equipping the above loss with a constraint following Locatello et al. ():

$$ $\displaystyle L_{loss}=\dfrac{1}{S}\sum_{\eta=1}^{S}$ $\displaystyle\left(\beta_{d}\sum_{j=1}^{n_{train}}{\left|\left|\hat{\mathcal{G }}_{[\hat{\theta}^{\eta}_{P},\theta_{I},\theta_{Q}]}[\bm{g}^{\eta}_{j}]-\bm{u} ^{\eta}_{j}\right|\right|}^{2}_{L^{2}({\Omega})}+{\left|\left|\hat{\theta}^{ \eta}_{P}-\mu^{\eta}_{\theta}\right|\right|}^{2}\right.+\left.\beta_{KL}{\left |\left|\mu^{\eta}_{\bm{z}}\right|\right|}^{2}+\beta_{cls}L_{c}(\bm{z}^{\eta},c (\beta^{\eta}))\right)$ $$

where $\beta_{cls}$ is a tunable hyperparameter, and $L_{c}(\bm{z}^{\eta},c(\beta^{\eta}))$ is the corresponding loss term depending on the provided information. For instance, when taking $c$ as the classification label of $\beta^{\eta}$, $L_{c}$ can be taken as the cross-entropy loss.

We note that the coefficient $\beta_{KL}$ is closely connected to the adjustable hyperparameter in $\beta-$VAE (Higgins et al., ). Taking a larger $\beta_{KL}$ encourages disentanglement (Locatello et al., ), while it also adds additional constraint on the implicit capacity of the latent bottleneck that results in information loss (Burgess et al., ). On the other hand,
the data reconstruction term forces the latent factors to contribute to the “reconstruction” of the more complicated solution field in a global perspective, and hence increasing $\beta_{d}$ is anticipated to alleviate information loss.
In our empirical studies, we demonstrate the interplay between these two tunable hyperparameters.

## 4 Experiments

We assess the performance of DisentangO across a variety of physics modeling and discovery datasets. In the following, we denote a DisentangO model with a latent dimension of $n$ as DisentangO-$n$. Our evaluation focuses on the following key aspects. Firstly, we demonstrate the capability of DisentangO in forward PDE learning, and compare its performance with a number of relevant baselines. In particular, we select three NO-related baselines, i.e., FNO (Li et al., ), NIO (Molinaro et al., ) and MetaNO (Zhang et al., ), as well as two non-NO baselines, i.e., VAE (Kingma and Welling, ) and $\beta$-VAE (Higgins et al., ).
Secondly, we showcase the merits of DisentangO in inverse modeling for interpretable physics discovery.
Thirdly, we perform parametric studies on the associated disentanglement hyperparameters. Lastly, we provide interpretations of the disentangled parameters and point out future directions.
Further details are provided in Appendix C.

### 4.1 Supervised forward and inverse PDE learning

**Table 1: Test errors and the number of trainable parameters. Bold number highlights the best method.**
| Models | #param | Test errors |  |
| --- | --- | --- | --- |
|  |  | data | latent supervision |
| DisentangO | 697,383 | 1.65% | 4.63% |
| MetaNO | 295,842 | 1.59% | - |
| NIO | 709,260 | - | 15.16% |
| FNO | 698,284 | 2.45% | 14.55% |
| VAE | 697,672 | 49.97% | 16.34% |
| $\beta$-VAE | 697,672 | 51.16% | 16.47% |

We start by investigating DisentangO’s capability in solving both forward and inverse PDE problems in a fully supervised setting. Specifically, we consider the constitutive modeling of anisotropic fiber-reinforced hyperelastic materials governed by the Holzapfel–Gasser–Ogden (HGO) model, where the data-generating process is controlled by sampling the governing material parameter set $\{E,\nu,k_{1},k_{2},\alpha\}$ and the latent factors can be learned consequently in a supervised fashion.
In this setting, the model takes as input the padded traction field and learns to predict the displacement.

As the number of latent factors is fixed for supervised learning in this example, we defer our ablation study to the second experiment of semi-supervised learning. We report in Table 1 our experimental results from DisentangO as well as its comparison with the selected baselines. Note that only FNO, VAE and $\beta$-VAE are able to handle both forward and inverse problems, whereas MetaNO and NIO can only solve the forward problem and the inverse problem, respectively. With this caveat in mind, MetaNO achieves the best performance in forward PDE learning, with DisentangO performing comparably well and beating the third best model by 32.7% in accuracy. In terms of inverse modeling with full latent supervision, DisentangO is the only method that can hold the error below 5%, outperforming the second best method by 68.2%.

### 4.2 Semi-supervised mechanical MNIST benchmark

We next explore in the context of semi-supervised learning and apply DisentangO and relevant baselines to the benchmark Mechanical MNIST (MMNIST) dataset (Lejeune, ). MMNIST comprises 70,000 heterogeneous material specimens undergoing large deformation, each governed by a material model of the Neo-Hookean type with a varying modulus converted from the MNIST bitmap images. In our experiment, we take 500 images (with a 420/40/40 split for training/validation/test) and generate 200 loading/response data pairs per sample on a $29\times 29$ structured grid, simulating uniaxial extension, shear, equibiaxial extension, and confined compression load scenarios. Since only partial knowledge is available for each image (i.e., the corresponding digit), we apply semi-supervised learning to the latent factors to classify the digits.

Ablation study. Firstly, we investigate DisentangO’s predictability in forward PDE learning by comparing its performance to MetaNO (i.e., the base meta-learnt NO model without the hierarchical structure). As seen in Table 2, MetaNO achieves a forward prediction error of 2.68%, which serves as the optimal bound for DisentangO. As we increase the latent dimension in DisentangO from 2 to 15, the prediction error drops from 11.49% to 5.48% and converges to the optimal bound.
Next, we study the role of the data loss term in disentanglement by gradually varying $\beta_{d}$ from $\beta_{d}=1$ to $\beta_{d}=1000$. In Table 2 we observe a consistent improvement in accuracy with the increase in $\beta_{d}$, where the boost in accuracy becomes marginal or slightly deteriorates beyond $\beta_{d}=100$. We thus choose $\beta_{d}=100$ as the best DisentangO model in this case. Besides offering an increased accuracy in forward modeling, the data loss term also enhances disentanglement as discussed in Section 3.2. This is evidenced by the illustration in Figure 2, where the unsupervised mutual information (MI) score that measures the amount of MI across latent factors consistently decreases in both DisentangO-2 and DisentangO-15 as we increase $\beta_{d}$ from $\beta_{d}=1$ (solid lines) to $\beta_{d}=100$ (dashed lines). On the contrary, the classification term poses a negative effect on disentanglement, as indicated by the increase in $\beta_{\text{cls}}$ leading to an increase in MI scores. This is reasonable because the classifier linearly combines all the latent factors to make a classification. The more accurate the classification, the stronger the correlation across latent factors. We then move on to study the effect of semi-supervised learning by comparing the model’s performance with and without latent partial supervision. While the models without semi-supervision is slightly more accurate in that the forward prediction accuracy with $\beta_{d}=1$ reaches 12.60% and 6.16% in DisentangO-2 and DisentangO-15, respectively, the latent scatterplot in Figure 3 reveals the inability of the DisentangO model without latent semi-supervision to acquire the partial knowledge of the embedded digits from data. In contrast, although the accuracy of DisentangO with latent semi-supervision slightly deteriorates due to the additional regularization effect from the classification loss term, it is able to correctly recognize the embedded digits and leverage this partial knowledge in disentangling meaningful latent factors.

**Table 2: Test errors and the number of trainable parameters in MMNIST benchmark. Models are abbreviated (0=DisentangO-2, 1=DisentangO-5, 2=DisentangO-10, 3=DisentangO-15, 4=VAE, 5=$\beta$-VAE, 6=MetaNO). Bold number highlights the best method that can handle both forward and inverse settings.**
| Models |  | 0 | 1 | 2 | 3 | 4 | 5 | 6 |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
| #param (M) |  | 0.66 | 0.97 | 1.49 | 2.02 | 2.02 | 2.02 | 0.35 |
| Error | $\beta_{d}=1$ | 12.82% | 9.56% | 7.36% | 6.29% | 16.34% | 17.13% | 2.68% |
| $\beta_{d}=10$ | 11.51% | 9.16% | 6.62% | 5.95% | - | - | - |  |
| $\beta_{d}=100$ | 11.49% | 8.43% | 6.65% | 5.48% | - | - | - |  |
| $\beta_{d}=1000$ | 11.62% | 8.22% | 6.50% | 5.80% | - | - | - |  |

Figure: Figure 2: MMNIST unsupervised scores against $\beta_{d}$ with DisentangO-2 (left) and DisentangO-15 (right). By comparing $\beta_{d}=1$ (solid lines) with $\beta_{d}=100$ (dashed lines), increasing $\beta_{d}$ forces the latent factors to maximize the contained information and in turn decreases MI, thus encouraging disentanglement.
Refer to caption: x2.png

Figure: Figure 3: MMNIST scatterplot with DisentangO-2 and $\beta_{d}=1$: left: $(\beta_{kl}=1,\beta_{cls}=100)$ (data test error 18.81%), middle: $(\beta_{kl}=10,\beta_{cls}=10)$ (data test error 16.94%), right: fully unsupervised DisentangO-2 without classification loss $(\beta_{kl}=100)$ (data test error 12.65%).
Refer to caption: x4.png

Comparison with additional baselines. We compare the performance of DisentangO with relevant additional baselines, i.e., VAE and $\beta$-VAE. We abandon the three NO-related baselines in this and the following experiments because they are black-box approximators and do not possess any mechanism to extract meaningful information without supervision. In this context, even the least accurate DisentangO-2 model (i.e., with $\beta_{d}=1$) with significantly fewer parameters outperforms the two selected baselines in Table 2 by 21.5% and 25.2%, respectively, with the best performing DisentangO-15 model beating the two baselines by 66.5% and 68.0% in accuracy, respectively. We do not proceed to compare the physics discovery capability between them as the baseline models are considerably inaccurate in forward prediction.

Interpretable physics discovery. To showcase the interpretability of DisentangO, we visualize its latent variables through latent traversal based on a randomly picked loading field as input in DisentangO-2, as displayed in Figure 6 of Appendix B.3. One can clearly see that the digit changes from “6” to “0” and then “2” from the top left moving down, and from “6” to “1” and then “7” moving to the right. Other digits are visible as well such as “7”, “9”, “4” and “8” in the right-most column. This corresponds well to the distribution of the latent clustering in Figure 3. More discussions on the interpretability of DisentangO-15 can be found in Figures 7 and 8 in Appendix B.3.

### 4.3 Unsupervised heterogeneous material learning

We now demonstrate the applicability of DisentangO in unsupervised learning in the context of learning synthetic human tissues that exhibit highly organized structures, with collagen fiber arrangements varying spatially. In this subject, it is critical to understand the underlying low-dimensional disentangled properties in the hidden latent space of this complex, high-dimensional microstructure, as inferred from experimental mechanical measurements. We generate two datasets, each containing 500 specimens and 100 loading/displacement pairs. The first dataset features variations in fiber orientation distributions using a Gaussian Random Field (GRF) and the second differs fiber angles in two segmented regions, separated by a centerline with a randomly rotated orientation.

Figure: Figure 4: Latent traversal of DisentangO-3 in unsupervised heterogeneous material learning, where the three latent dimensions control the border rotation between the two segments (top), the relative fiber orientation between the two segments (middle), and the fiber orientation of the top segment (bottom), respectively. Legend indicates fiber orientation ranging from 0 to $\pi$.
Refer to caption: x7.png

Table: Table 3: Test errors and the number of trainable parameters in experiment 3. Models are abbreviated (0=DisentangO-2, 1=DisentangO-5, 2=DisentangO-10, 3=DisentangO-15, 4=DisentangO-30, 5=VAE, 6=$\beta$-VAE, 7=MetaNO). Bold number highlights the best method that can handle both forward and inverse settings.

| Models | #param (M) | Error |  |
| --- | --- | --- | --- |
|  |  | $\beta_{d}=1$ | $\beta_{d}=100$ |
| 0 | 0.63 | 26.33% | 25.18% |
| 1 | 0.94 | 17.51% | 15.75% |
| 2 | 1.46 | 10.76% | 10.01% |
| 3 | 1.98 | 7.11% | 7.02% |
| 4 | 3.55 | 5.33% | 5.28% |
| 5 | 3.55 | 61.10% | - |
| 6 | 3.55 | 57.04% | - |
| 7 | 0.32 | 2.67% | - |

Figure: Figure 5: Unsupervised MI score against $\beta_{kl}$ with DisentangO-2 in heterogeneous material learning.

We report in Table 3 our experimental results of DisentangO with different latent dimensions and data loss strength $\beta_{d}$, along with comparisons with the baseline model results.
Consistent with the findings in the second experiment, increasing $\beta_{d}$ results in a boost in accuracy by comparing the rows in the table. The model’s predictive performance also converges to the optimal bound of MetaNO as we increase the latent dimension from 2 to 30, where the model’s prediction error improves from 25.18% to 5.28%. This significantly outperforms the best baseline model by 90.7%. On the other hand, the effect of the data loss term on disentanglement is further proved in Figure 5, where increasing $\beta_{d}$ leads to a decrease in MI between the latent factors, thus encouraging disentanglement. Lastly, we interpret the mechanism of the learned latent factors in DisentangO-3 via learning a mapping between the learned latent factors and the underlying material microstructure and subsequently performing latent traversal in each dimension. The results are shown in Figure 4, where the three latent factors manifest control on the border rotation between the two segments, the relative fiber orientation between the two segments, and the fiber orientation of the top segment, respectively.

## 5 Conclusion

We present DisentangO, a novel hyper-neural operator architecture designed to disentangle latent physical factors embedded within black-box neural operator parameters. DisentangO leverages a hierarchical neural operator architecture that extracts varying parameters of the governing PDE through a task-wise adaptive layer, then further disentangles these variations into distinct latent factors. By learning these disentangled representations, DisentangO not only enhances physical interpretability but also enables more robust generalization across diverse physical systems. We perform empirical evaluations across supervised, semi-supervised, and unsupervised learning contexts.

Limitations: Due to limited computational resource, our experiments are focused on learning a small-to-medium number ($<500$). Expanding the coverage to encompass a broader range of physical systems and enabling learning across more diverse parameters would be beneficial.

## Acknowledgments and Disclosure of Funding

L. Zhang would like to acknowledge support by the AFOSR grant FA9550-22-1-0197. Y. Yu would like to acknowledge the support by the National Science Foundation under award DMS-2436624 and the support from the National Institute of Health under award 1R01GM157589-01. Portions of this research were conducted on Lehigh University’s Research Computing infrastructure partially supported by NSF Award 2019035.

## Appendix A Proofs of identifiability analysis

We first provide the proof of Theorem 1:

Proof:
With equation 3.15, we have:

$$ $p_{\mathcal{G}[\bm{f};\theta]|\bm{f}}=p_{\hat{\mathcal{G}}[\bm{f};\hat{\theta} ]|\bm{f}}\Leftrightarrow p_{\mathcal{G}_{f}(\theta)|\bm{f}}=p_{\hat{\mathcal{G }}_{f}(\hat{\theta})|\bm{f}}\Leftrightarrow p_{\theta|\bm{f}}=p_{\mathcal{G}_{ f}^{-1}\circ\hat{\mathcal{G}}_{f}(\hat{\theta})|\bm{f}}.$ $$

Note that the parameter $\theta$ varies with the change of $\bm{b}$. Per the data generating process in equation 3.6, the distribution of $\bm{b}$ is invariant to $\bm{f}$. Therefore, the distribution of $\theta$ is also invariant to $\bm{f}$:

$$ $p_{\theta}=p_{\theta|\bm{f}}=p_{\mathcal{G}_{f}^{-1}\circ\hat{\mathcal{G}}_{f} (\hat{\theta})|\bm{f}}=p_{\mathcal{G}_{f}^{-1}\circ\hat{\mathcal{G}}_{f}(\hat{ \theta})},\;\forall\bm{f}\in\mathcal{F}.$ $$

Denoting $r:={\mathcal{G}}_{f}\circ\hat{\mathcal{G}}^{-1}_{f}$, it is the transformation between the true $\theta$ and the estimated one, and it is invertible and invariant with respect to $\bm{f}$.

We proceed to derive the relation between $\bm{z}$ and $\hat{\bm{z}}$: since $\theta=r(\hat{\theta})$, with the invertibility assumption $\theta=\mathcal{H}^{-1}(\bm{z})$ and $\hat{\theta}=\hat{\mathcal{H}}^{-1}(\hat{\bm{z}})$, we obtain:

$$ $\bm{z}=\mathcal{H}(\theta)=\mathcal{H}\circ r(\hat{\theta})=\mathcal{H}\circ r \circ\hat{\mathcal{H}}^{-1}(\hat{\bm{z}}).$ $$

Denoting $h:=\mathcal{H}\circ r\circ\hat{\mathcal{H}}^{-1}$, it is the transformation between the true latent variable and the estimated one, and it is invertible because $r$, $\mathcal{H}$ and $\hat{\mathcal{H}}$ are all invertible.∎

We now show the proof of Theorem 2.

Proof: With the independence relation assumption, we have

$$ $p_{\bm{z}|\bm{b}}(\bm{z}|\bm{b})=\prod_{i}p_{z_{i}|\bm{b}}(z_{i}),\;p_{\hat{ \bm{z}}|\bm{b}}(\hat{\bm{z}}|\bm{b})=\prod_{i}p_{\hat{z}_{i}|\bm{b}}(\hat{z}_{ i})\text{ .}$ $$

Denoting $\hat{q}_{i}:=\log p_{\hat{z}_{i}|\bm{b}}$, it yields:

$$ $\log p_{\bm{z}|\bm{b}}(\bm{z}|\bm{b})=\sum_{i}q_{i}(z_{i},\bm{b}),\;\log p_{ \hat{\bm{z}}|\bm{b}}(\hat{\bm{z}}|\bm{b})=\sum_{i}\hat{q}_{i}(\hat{z}_{i},\bm{ b})\text{ .}$ $$

With the change of variables we have

$$ $p_{\bm{z}|\bm{b}}=p_{h({\hat{\bm{z}}})|\bm{b}}=p_{{\hat{\bm{z}}}|\bm{b}}\cdot| J_{h^{-1}}|\Leftrightarrow\sum_{i}q_{i}(z_{i},\bm{b})+\log|J_{h}|=\sum_{i}\hat {q}_{i}(\hat{z}_{i},\bm{b})\text{ ,}$ $$

where $|J_{h^{-1}}|$ stands for the absolute value of the Jacobian matrix determinant of $h^{-1}$. Differentiating the above equation twice with respect to $\hat{z}_{k}$ and $\hat{z}_{q}$, $k\neq q$, yields

$$ $\sum_{i}\left(\dfrac{\partial^{2}q_{i}(z_{i},\bm{b})}{\partial z_{i}^{2}} \dfrac{\partial z_{i}}{\partial\hat{z}_{k}}\dfrac{\partial z_{i}}{\partial\hat {z}_{q}}+\dfrac{\partial q_{i}(z_{i},\bm{b})}{\partial z_{i}}\dfrac{\partial^{ 2}z_{i}}{\partial\hat{z}_{k}\partial\hat{z}_{q}}\right)+\dfrac{\partial^{2} \log|J_{h}|}{\partial\hat{z}_{k}\partial\hat{z}_{q}}=0\text{ .}$ (A.1) $$

To show the identifiability, one can rewrite the Jacobian $J_{h}$ as:

$$ $J_{h}=\left[\dfrac{\partial\bm{z}}{\partial\hat{\bm{z}}}\right].$ $$

The invertibility results shown in Theorem 1 indicates that it is full rank. Next, we will use the linear independence assumption to show that there exists one and only one non-zero component in each row of $\dfrac{\partial\bm{z}}{\partial\hat{\bm{z}}}$.

Taking $\bm{b}=\bm{b}^{0},\cdots,\bm{b}^{2d_{z}}$ in equation A.1 and subtracting them from each other, we have

$$ $\sum_{i=1}^{d_{z}}\left(\left(\dfrac{\partial^{2}q_{i}(z_{i},\bm{b}^{j})}{ \partial z_{i}^{2}}-\dfrac{\partial^{2}q_{i}(z_{i},\bm{b}^{0})}{\partial z_{i} ^{2}}\right)\dfrac{\partial z_{i}}{\partial\hat{z}_{k}}\dfrac{\partial z_{i}}{ \partial\hat{z}_{q}}+\left(\dfrac{\partial q_{i}(z_{i},\bm{b}^{j})}{\partial z _{i}}-\dfrac{\partial q_{i}(z_{i},\bm{b}^{0})}{\partial z_{i}}\right)\dfrac{ \partial^{2}z_{i}}{\partial\hat{z}_{k}\partial\hat{z}_{q}}\right)=0\text{ ,}$ $$

where $j=1,\cdots,2d_{z}$. With the linear independence condition for $w$, this is a $2d_{z}\times 2d_{z}$ linear system, and therefore the only solution is

$$ $\dfrac{\partial z_{i}}{\partial\hat{z}_{k}}\dfrac{\partial z_{i}}{\partial\hat {z}_{q}}=0\text{ ,}\quad\;\dfrac{\partial^{2}z_{i}}{\partial\hat{z}_{k} \partial\hat{z}_{q}}=0\text{ ,}$ $$

for $i=1,\cdots,d_{z}$. The first part implies that, for the $i-$th row of the Jacobian matrix $J_{h}$, we have $\dfrac{\partial z_{i}}{\partial\hat{z}_{k}}\neq 0$ for at most one element $k\in\{1,\cdots,d_{z}\}$, hence $\bm{z}$ is identifiable up to permutation and component-wise invertible transformation. ∎

## Appendix B Additional experimental details, results and discussion

Herein, we provide additional details in training and baseline models, as well as more results as a supplement of Section 4.

### B.1 Training details

In all experiments, we adopt the Adam optimizer for optimization and use a 16-layer DisentangO model. For fair comparison across different models, we tune the hyperparameters, including the learning rates, the decay rates, and the regularization parameters, to minimize the validation loss. Experiments are conducted on a single NVIDIA Tesla A100 GPU with 40 GB memory. A pseudo algorithm of all three scenarios is summarized in Algorithm 1.

Figure: Algorithm 1 A pseudo algorithm of DisentangO.

### B.2 Supervised forward and inverse PDE learning

The parameter of each baseline is given in the following, where the parameter choice of each model is selected by tuning the number of layers and the width (channel dimension), keeping the total number of parameters on the same magnitude.

- •
MetaNO: We use a 16-layer IFNO model with the lifting layer as the adaptive layer. We keep the total number of parameters in MetaNO the same as the number of parameters used in forward PDE learning in DisentangO.
- •
NIO: We closely follow the setup in Molinaro et al. () and use four convolution blocks as the encoder for the branch net and a fully connected neural network with two hidden layers of 256 neurons as the trunk net, with the number of basis functions set to 50. For the FNO part, we use one Fourier layer with width 32 and modes 8, as suggested in Molinaro et al. ().
- •
FNO: We adopt a 4-layer FNO with width 26 and modes 8. For the forward model, in addition to the loading field and coordinates as input, we also concatenate the ground-truth material properties to form the final input. For the inverse model, we employ an iterative gradient-based optimization to solve for the optimal material parameters for physics discovery.
- •
VAE: We use a 2-layer MLP of size [1681, 136, 30] as the encoder, and another 2-layer MLP of size [30, 136, 1681] as the decoder, with the size of the bottleneck layer being 30.
- •
$\beta$-VAE: The parameter choice of the $\beta$-VAE baseline is the same as the VAE baseline, except that we tune the $\beta$ hyperparameter.

### B.3 Mechanical MNIST benchmark

Figure: Figure 6: MMNIST latent traversal based on a randomly picked loading field as DisentangO-2 input.
Refer to caption: x9.png

Herein, we demonstrate the latent traversal based on a randomly picked loading field as DisentangO-2 input. One can clearly see that the digit changes from “6” to “0” and then “2” from the top left moving down, and from “6” to “1” and then “7” moving to the right. Other digits are visible as well such as “7”, “9”, “4” and “8” in the right-most column. This corresponds well to the distribution of the latent clustering in Figure 3. We also provide two exemplary MMNIST latent interpretations of the learned DisentangO-15 in Figures 7 and 8.

Figure: Figure 7: Exemplary MMNIST latent interpretation of DisentangO-15: dimension 10 encodes the information for digit ‘8’, as is evidenced by the fact that all $y<-0.4$ regions on the scatterplots contain only digit ‘8’.
Refer to caption: x10.png

Figure: Figure 8: Exemplary MMNIST latent interpretation of DisentangO-15: dimension 12 encodes the information for digit ‘3’, as is evidenced by the fact that all $y<-0.45$ regions on the scatterplots contain only digit ‘3’.
Refer to caption: x24.png

## Appendix C Data Generation

### C.1 Example 1 - Supervised forward and inverse PDE learning

we consider the constitutive modeling of anisotropic fiber-reinforced hyperelastic materials governed by the Holzapfel–Gasser–Ogden (HGO) model, whose strain energy density function can be written as:

$$ $\displaystyle\eta$ $\displaystyle=\frac{E}{4(1+\nu)}(\bar{I}_{1}-2)-\frac{E}{2(1+\nu)}\ln J$ $\displaystyle+\frac{k_{1}}{2k_{2}}(e^{k_{2}\langle S(\alpha)\rangle^{2}}+e^{k_ {2}\langle S(-\alpha)\rangle^{2}}-2)+\frac{E}{6(1-2\nu)}(\frac{J^{2}-1}{2}-\ln J )\text{ ,}$ (C.1) $$

where $\langle\cdot\rangle$ indicates the Macaulay bracket, $\alpha$, $k_{1}$ and $k_{2}$ are the fiber angle, modulus and exponential coefficient, respectively, $E$ denotes the Young’s modulus of the matrix, $\nu$ is Poisson’s ratio, and $S(\alpha)$ describes the fiber strain of the two fiber groups, $S(\alpha)=\frac{\bar{I}_{4}(\alpha)-1+|\bar{I}_{4}(\alpha)-1|}{2}$. $\bar{I}_{i}$ is the $i^{\text{th}}$ invariant of the right Cauchy-Green tensor $\bm{C}$, $\bar{I}_{1}=tr(\bm{C})$ and $\bar{I}_{4}=\bm{n}^{T}(\alpha)\bm{C}\bm{n}(\alpha)$, with $\bm{n}(\alpha)=[\cos(\alpha),\sin(\alpha)]^{T}$. In this context, the data generation process is controlled by sampling the material set $\{E,\nu,k_{1},k_{2},\alpha\}$, and the latent factors can be learned consequently in a supervised fashion.
The physical parameters are sampled from $\frac{k_{1}}{k_{2}}\sim\mathcal{U}[90,100]$ ,$k_{2}\sim\mathcal{U}[0.001,0.1]$, $E\sim\mathcal{U}[0.5001,0.6001]$, $\nu\sim\mathcal{U}[0.2,0.3]$, and $\alpha\sim\mathcal{U}[\pi/10,\pi/2]$. To generate the high-fidelity (ground-truth) dataset, we sample 220 material sets, which are split into 200/10/10 for training/validation/test, respectively. For each material set, we sample $50$ different vertical traction conditions $T_{y}(\bm{x})$ on the top edge from a random field, following the algorithm in Lang and Potthoff (); Yin et al. (). $T_{y}(\bm{x})$ is taken as the restriction of a 2D random field, $\phi(\bm{x})=\mathcal{F}^{-1}(\gamma^{1/2}\mathcal{F}(\Gamma))(\bm{x})$, on the top edge. Here, $\Gamma(\bm{x})$ is a Gaussian white noise random field on $\mathbb{R}^{2}$, $\gamma=(w_{1}^{2}+w^{2}_{2})^{-\frac{2}{4}}$ represents a correlation function, and $w_{1}$, $w_{2}$ are the wave numbers on $x$ and $y$ directions, respectively. Then, for each sampled traction loading, we solve the displacement field on the entire domain by minimizing the potential energy using the finite element method implemented in FEniCS (Alnæs et al., ). Sample data of the obtained dataset is illustrated in Figure 9. In this setting, the model takes as input the padded traction field and learns to predict the resulting displacement field.

Figure: Figure 9: Illustration of the HGO data, loading and displacement pairs of three materials $\bm{b}^{\eta}$, $\eta=1,2,3$. Top: three instances of different loadings $T_{y}(x)$, which corresponds to different $\bm{f}_{i}$. Bottom: corresponding displacement solutions $\bm{u}^{\eta}_{i}$, illustrating the impacts of system ($\bm{b}^{\eta}$) variability in solution operators.
Refer to caption: extracted/5897506/example1_hvaeno.png

### C.2 Example 2 - Mechanical MNIST benchmark

Mechanical MNIST is a benchmark dataset of heterogeneous material undergoing large deformation, modeled by the Neo-Hookean material with a varying modulus converted from the MNIST bitmap images (Lejeune, ). It contains 70,000 heterogeneous material specimens, and each specimen is governed by the Neo-Hookean material with a varying modulus converted from the MNIST bitmap images. We illustrate samples from the MMNIST dataset in Figure 10, including the underlying microstructure, two randomly picked loading fields, and the corresponding displacement fields.

Figure: Figure 10: Illustration of exemplar MMNIST samples in the semi-supervised scenario. (a): material parameter field corresponding to different $\bm{b}^{\eta}$. (b): displacement fields (second row) $\bm{u}^{\eta}_{38}$ corresponding to the same loading field (first row) $\bm{f}_{38}$. (c): displacement fields (second row) $\bm{u}^{\eta}_{138}$ corresponding to the same loading field (first row) $\bm{f}_{138}$.
Refer to caption: x38.png

### C.3 Example 3 - Unsupervised heterogeneous material learning

We generate two sets of datasets in this case, varying the material microstructure in the following two ways.
Varying fiber orientation distribution.
We generate the samples by controlling the predefined parameters of Gaussian Random Field (GRF). With the GRF sharpened by the thresholding, the values to binary field represent two distinct fiber orientations. The binary field is smoothed using a windowed convolution. To address boundary conditions, the matrix is padded with replicated edge values, ensuring that the convolution works uniformly across the entire grid. After the smoothing process, the padded sections are removed, and the remaining field is used to construct the fiber field. We use two fixed fiber angles $\frac{\pi}{3},\frac{2\pi}{3}$ for the corresponding binary field. We generate 300 material sets, each with 500 loading/displacement pairs, and divide these into training, validation, and test sets in a 200/50/50 split. Exemplar Samples from this dataset are illustrated in Figure 11. These samples demonstrate the variability of $\bm{b}$, which is of critical for the latent variable identifiability in unsupervised learning settings, as proved in Theorem 2.

Varying fiber orientation magnitude and segmentation line rotation.
Instead of controlling the fiber orientation angles on the binary field as two constant values, we sample the orientation distribution consisting of two segments with orientations $\alpha_{1}$ and $\alpha_{2}$ on each side, respectively, separated by a line passing through the center. The values of $\alpha_{1}$ and $\alpha_{2}$ are independently sampled from a uniform distribution over $[0,2\pi]$, and the centerline’s rotation is sampled from $[0,\pi]$. We generate 300 material sets, each with 500 loading/displacement pairs, and divide these into training, validation, and test sets in a 200/50/50 split.

Figure: Figure 11: Illustration of fiber orientation magnitude and segmentation line rotation. Upper: three different loading instances of $\bm{f}_{i}$. Bottom left: three different microstructure instances of $\bm{b}^{\eta}$. Bottom right: corresponding solution fields $\bm{u}^{\eta}_{i}$.
Refer to caption: extracted/5897506/example3_micro.png

Loading and displacement pairs for one microstructure.
After generating the specimens with varying fiber orientations, we feed the $\bm{b}^{\eta}(x)$ as the $\alpha(x)$ to the HGO model, and keep the material property set $\{E,\nu,k_{1},k_{2}\}$ as constant. For each microstructure, we randomly sample loading and displacement pairs for each microstructure from the previous step.

The loading in this example is taken as the body load, $\bm{f}(\bm{x})$. Each instance is generated as the restriction of a 2D random field, $\phi(\bm{x})=\mathcal{F}^{-1}(\gamma^{1/2}\mathcal{F}(\Gamma))(\bm{x})$. Here, $\Gamma(\bm{x})$ is a Gaussian white noise random field on $\mathbb{R}^{2}$, $\gamma=(w_{1}^{2}+w^{2}_{2})^{-\frac{5}{4}}$ represents a correlation function, $w_{1}$ and $w_{2}$ are the wave numbers on $x$ and $y$ directions, respectively, and $\mathcal{F}$ and $\mathcal{F}^{-1}$ denote the Fourier transform and its inverse, respectively. This random field is anticipated to have a zero mean and covariance operator $C=(-\Delta)^{-2.5}$, with $\Delta$ being the Laplacian with periodic boundary condition on $[0,2]^{2}$, and we then further restrict it to $\Omega$. For the detailed implementation of Gaussian random field sample generation, we refer interested readers to Lang and Potthoff (). Then, for each sampled loading $\bm{f}_{i}(\bm{x})$ and microstructure field $\bm{b}^{\eta}(\bm{x})$, we solve for the displacement field $\bm{u}^{\eta}_{i}(\bm{x})$ on the entire domain.

## Peer Reviews

confidence: 3
contribution: 2
initial_rating: 5
presentation: 2
questions: 1. Theorem 1 addresses the identifiability issue up to an invertible transformation, which may not depend on the differential equation parameters. Could the authors clarify the source of this identifiability issue, and specify what guarantees the theorem provides regarding distribution uniqueness despite this transformation? (While Theorem 2 confirms identifiability even without such a transformation, further elaboration could be beneficial.)

2. In the proposed neural operator model, shared parameters are introduced to reduce degrees of freedom. While the architecture appears sufficient to support the presented theorems, could the authors clarify whether this model structure is indeed minimal? For instance, in line 314, the authors parameterize the posterior distribution $q$ with projection weights $\mathcal{P}$, possibly because the weights in $\mathcal{P}$ are fewer than those in $\mathcal{Q}$ or the intermediate layer. Is this reduction the primary rationale behind the model’s design? We would greatly appreciate further details on the trade-offs in using additional weights in $\mathcal{Q}$ or the intermediate layer for parameterization and how this choice impacts the model's overall performance and interpretability.

3. In line 210, the statement that parameters in the equation cannot be learned under zero Dirichlet boundary conditions is intriguing. Could the authors briefly elaborate on this point?

4. Is it possible to empirically validate one of the assumptions in Theorems 1 and 2, such as linear independence, within the simulation? Although these assumptions may be reasonable based on literature in line 302, empirical support could further strengthen the argument.

5. Tables 2 and 3 indicate that increasing the latent dimension improves prediction results. Is there a specific reason, such as computational constraints, for selecting a latent dimension smaller than 30? The authors might consider including experiments with higher latent dimensions to explore this effect further.

6. In Table 1, MetaNO, which is applied solely to forward problems, utilizes fewer parameters than other models. Could the authors clarify the rationale behind this choice? Additionally, could the author provide the inference time or total training time for each model with the specified number of parameters to enable a more thorough comparison?

7. Could the authors briefly explain how parameter classification information was utilized in NIO and FNO within the latent supervision setting in Table 1?
soundness: 2
strengths: The paper effectively presents a clear motivation and a well-organized introduction. 

To the best of my knowledge, the theoretical contribution addressing the identifiability issue is both novel and significant.
summary: This paper presents a model based on a neural operator architecture, integrating an encoder-decoder structure to derive interpretable representations and address multiphysics problems by disentangling data into distinct factors.
weaknesses: The methodology section would benefit from additional clarification in the following areas.

1. In lines 189 and 204, the term “parameters” referred to both the parameters within the differential equation and those within the neural operator, which may introduce ambiguity. While the authors have distinguished these parameters through notation, we kindly request that the authors revise the manuscript to adopt distinct terminology for each set. This revision would further enhance clarity in distinguishing between the two. 

2. In Theorems 1 and 2, separating the assumptions from the main theorem statements could enhance readability. The current formulation may be somewhat challenging to follow, and presenting the assumptions outside each theorem could help clarify their interpretation. For instance, the authors might consider defining all relevant assumptions -- such as density smoothness, invertibility, conditional independence, and linear independence -- at the beginning of Section 3.2. These assumptions could then be referenced as needed throughout the related theorems. Additionally, the authors may wish to discuss the importance and validity of each assumption to support the theoretical framework. 

3. The term "latent supervision" appears to denote a semi-supervised learning context (SC2) in the paper. Using consistent terminology throughout the paper would enhance readability and facilitate understanding. If our interpretation is accurate, we kindly request that the authors consistently use "semi-supervised learning when referring to SC2 throughout the paper.

---

confidence: 3
contribution: 2
initial_rating: 5
presentation: 2
questions: - Interpretability: is there any way to systematically validate that the latent factors learned by DisentangO correspond to actual physical parameters, especially in the unsupervised setting with large latent dimension?
- In Table 1, it is shown that a MetaNO with half as many parameters as DisentangO matches on forward problem error. How does the comparison look with a parameter-matched MetaNO? Similarly with Table 3.
- How does the choice of neural operator architecture (IFNO for the DisentangO vs. other architectures as baselines) affect the results? How would IFNO compare on the supervised setting, and could the model be improved by using MetaNO's architecture for DisentangO?
soundness: 2
strengths: - The paper investigates an important question, interpretability in PDE learning, which has been underexplored in the literature.
- The authors provide interesting theoretical results regarding when physical parameters can be uniquely identified.
- The empirical results are promising:
  - Especially in the supervised and semi-supervised settings, the proposed method is shown to match or outperform existing SOTA neural operator models.
  - Furthermore, the distinction of supervised vs. semi-supervised vs. unsupervised learning in the setting of neural operators is interesting, and the proposed architecture is one of the few which can engage with all three regimes.
summary: The authors introduce DisentangO, a novel neural operator architecture based on a variational framing that aims to perform both physics modeling (e.g. forward PDE solving) and governing physics discovery (e.g. inverse PDE solving). The authors provide theoretical analysis showing when the true physical factors can be identified from the learned representations. The effectiveness of the method is demonstrated in supervised, semi-supervised, and unsupervised learning settings. Additionally, the authors demonstrate that the learned VAE embedding space corresponds to physically meaningful properties on the Mechanical MNIST benchmark.
weaknesses: - The main weakness of the paper is that the empirical results regarding the claims of interpretability are somewhat limited. The authors do not propose clear metrics for evaluating interpretability:
  - MMNIST: the authors show that the latent space clusters data based on digit (Fig 3, 7, 8) and that the latent space can be interpolated (Fig 6). However, these seem to be standard properties of VAEs, and it's unclear whether meaningful physical interpretability is being achieved.
  - Heterogeneous material learning: the authors show that a DisentangO with a latent dimension of 3 is able to associate a different physically-meaningful characteristic with each of its dimensions. However, this problem is quite toy, and it's unclear whether these insights transfer to larger problems/networks.
- The paper's writing clarity could be improved in places:
  - The proposed architecture is described in the intro/abstract/conclusion as a "hyper-neural operator" and "hierarchical variational autoencoder", but these terms are not clearly defined in the pape. In particular, the term "hierarchical VAE" seems to be used differently from the typical definition in the literature [1], which involves partitioning latent variables into disjoint groups and modeling conditionally ($\prod_l p(z_l \mid z_{< l})$). The model used in the paper seems to be a standard VAE, if I understand correctly.
  - Similarly, the intro mentions that DisentangO learns from "multiple physical systems", when it seems the paper focuses on different instances of the same PDEs/physical systems.
  - In general, the model architecture could be more clearly described in Section 3.1. A clearer caption in Figure 1 could also be helpful.
  - Tables 2 and 3 would be clearer if model names were included in the table rather than the abbreviation number.
- The existing empirical results could be clarified: see questions.

1. NVAE: A Deep Hierarchical Variational Autoencoder. Arash Vahdat, Jan Kautz.

---

confidence: 3
contribution: 2
initial_rating: 3
presentation: 1
questions: 1. Presentation
- The overall architecture presented in Figure 1 is not sufficiently informative. The meaning of each arrow in the figure is unclear, as some appear to represent feed-forward computations while others do not.
- Additionally, every notation should be clearly defined upon its first introduction; this paper lacks clarity in this regard.
- What are the first decoder and second decoder in Figure 1?
- How did you computed $\theta_P^{\eta}$ from the input? This is not explained in the manuscript. 
- How did you reconstructed $b^\{\eta}$ from $z$ for the semi-supervised and unsupervised settings? 

2. Experiments
- The authors conducted experiments on several benchmark datasets. What baseline algorithms, aside from neural operators and variational autoencoders (VAEs), were used for these datasets? Additionally, could you provide the test errors for the traditional methods?
soundness: 2
strengths: - The paper presents an innovative architecture that combines a hierarchical variational autoencoder (VAE) with a Fourier neural operator (FNO) to address inverse problems.
- The authors also provide a theoretical analysis on the component-wise identifiability of the model.
- The proposed method outperforms the baselines in the supervised settings.
summary: This paper introduces DisentangO, a novel variational hierarchical neural operator designed to solve forward and inverse problems simultaneously. The authors present theoretical results to substantiate the method, which is further evaluated through numerical experiments on various benchmark datasets.
weaknesses: - The presentation lacks clarity.
- The experiments conducted on semi-supervised and unsupervised datasets are insufficient, with a notable absence of comparisons to baseline algorithms.
- The learning process does not adequately incorporate information from the partial differential equations (PDEs).

---

confidence: 3
contribution: 2
initial_rating: 5
presentation: 3
questions: 1) Why do you use simple MLPs for VAE architectures as baselines and not more advanced architectures?

2) Why are relatively small models, like the 0.7M FNO, being used?
soundness: 2
strengths: 1) The authors present a novel approach based on VAEs that addresses both forward and inverse problems simultaneously.

2) They offer theoretical guarantees to support the effectiveness of their method.

3) The approach is evaluated against MetaNO and basic VAEs.

4) The paper is well-written and mathematically rigorous.
summary: The authors introduce DisentangO, a novel architecture that disentangles latent physical factors from neural operator parameters. The model is capable of solving forward and inverse problems at the same time. By combining a multi-task model with a VAE-based approach, DisentangO improves physical interpretability and generalization across PDE-driven systems. The paper offers theoretical guarantees to support the effectiveness of their method.
weaknesses: (1) The work has the potential to be highly impactful. However, several papers that take a similar approach are not cited in this work [1][2][3]. I believe these papers should be acknowledged and their contributions recognized.

(2) The networks used in the benchmarks are relatively small. For example, the FNO architecture that you used in supervised settings has only 8 modes and a width of 26, which is quite modest.

(3) The architecture used for baseline VAEs' encoders and decoders are simple MLPs. I believe that more appropriate and strnoger baselines would be convolutional VAEs, for example.

(4) The model is not compared to the state-of-the-art approaches such as [1][2][3]. The work done in [1][3] could be easily adapted to experiments in the paper. 

----

[1] Lingsch, L. E., Grund, D., Mishra, S., & Kissas, G. (2024). FUSE: Fast Unified Simulation and Estimation for PDEs. arXiv preprint arXiv:2405.14558.

[2] Fotiadis, S., Valencia, M. L., Hu, S., Garasto, S., Cantwell, C. D., & Bharath, A. A. (2023, July). Disentangled generative models for robust prediction of system dynamics. In International Conference on Machine Learning (pp. 10222-10248). PMLR.

[3] Tong, G. G., Long, C. A. S., & Schiavazzi, D. E. (2024). InVAErt networks: A data-driven framework for model synthesis and identifiability analysis. Computer Methods in Applied Mechanics and Engineering, 423, 116846.

Meta-review: This paper proposes a new variational hyper-operator that aims to disentangle physical factors of variation from the neural operator parameters used for solving parametric PDEs. The approach can handle both forward and inverse problems simultaneously. The experimental results demonstrate its effectiveness across various settings, including supervised, semi-supervised, and unsupervised learning.

The reviewer feedback is mixed, with one reviewer giving a score of 8 and another giving a score of 3. During the final discussion, reviewers recognized the relevance of interpretable learning in scientific contexts. While reviewer gJXj strongly advocated for the paper and highlighted its strengths, others indicated that further revisions would be beneficial. In particular, it would be valuable to (i) provide additional examples where the latent embedding offers a clear advantage, and (ii) discuss more thoroughly the practical benefits of interpretability. This includes studying scenarios where it may be preferable to use the interpretable model rather than a black-box alternative. Moreover, there are concerns about the effectiveness of the proposed approach. For instance, how does training and inference times compare for DisentangO and FUSE.

Overall, I found the paper interesting and well-written. Although the authors attempted to address many of the concerns during the rebuttal phase, not all reviewers were convinced. The paper seems borderline in its current state, and I believe another round of revisions to fully address the remaining issues would substantially improve it. Therefore, I recommend rejecting the paper at this time.

Decision: Reject