Maximum-likelihood Estimators in Physics-Informed Neural Networks for
  High-dimensional Inverse Problems
Physics-informed neural networks (PINNs) have proven a suitable mathematical
scaffold for solving inverse ordinary (ODE) and partial differential equations
(PDE). Typical inverse PINNs are formulated as soft-constrained multi-objective
optimization problems with several hyperparameters. In this work, we
demonstrate that inverse PINNs can be framed in terms of maximum-likelihood
estimators (MLE) to allow explicit error propagation from interpolation to the
physical model space through Taylor expansion, without the need of
hyperparameter tuning. We explore its application to high-dimensional coupled
ODEs constrained by differential algebraic equations that are common in
transient chemical and biological kinetics. Furthermore, we show that
singular-value decomposition (SVD) of the ODE coupling matrices (reaction
stoichiometry matrix) provides reduced uncorrelated subspaces in which PINNs
solutions can be represented and over which residuals can be projected.
Finally, SVD bases serve as preconditioners for the inversion of covariance
matrices in this hyperparameter-free robust application of MLE to
``kinetics-informed neural networks''.
Gabriel S. Gusmão, Andrew J. Medford
Maximum-likelihood Estimators in Physics-Informed Neural Networks for High-dimensional Inverse Problems

###### Abstract

Physics-informed neural networks (PINNs) have proven a suitable mathematical scaffold for solving inverse ordinary (ODE) and partial differential equations (PDE). Typical inverse PINNs are formulated as soft-constrained multi-objective optimization problems with several hyperparameters. In this work, we demonstrate that inverse PINNs can be framed in terms of maximum-likelihood estimators (MLE) to allow explicit error propagation from interpolation to the physical model space through Taylor expansion, without the need of hyperparameter tuning. We explore its application to high-dimensional coupled ODEs constrained by differential algebraic equations that are common in transient chemical and biological kinetics. Furthermore, we show that singular-value decomposition (SVD) of the ODE coupling matrices (reaction stoichiometry matrix) provides reduced uncorrelated subspaces in which PINNs solutions can be represented and over which residuals can be projected. Finally, SVD bases serve as preconditioners for the inversion of covariance matrices in this hyperparameter-free robust application of MLE to "kinetics-informed neural networks".

physics-informed neural network, surrogate approximator, maximum-likelihood, singular-value decomposition, dimensionality reduction, catalysis, transient, chemical kinetics

## 1 Main

Physics-informed [1] (PINNs), inspired or constrained [2] neural-networks comprise an interdisciplinary emerging area connecting machine learning to several fields of science, engineering and mathematics. It has shown that neural networks (NNs) provide a flexible structure for the solution of problems arising from fully- or partially-known differential equations, such as in heat-transfer [3, 4], compressible and incompressible, invisible and non-inviscible flows [5, 6], advection-diffusion equations [7], and other problems involving conservation laws and physical constraints.

The formulation of inverse problems based on PINNs typically entails softly-constrained cost functions with several penalty weights (hyperparameters) that need be tuned for a given set of data of known underlying physics (differential equation) [4, 5, 8, 9, 10, 11, 12]. Other strategies such as weighting based on gradient statistics, inverse Dirichlet weighting, and multiple gradient descent based on the Karush-Kuhn-Tucker (KKT) formulation have been applied [13] to satisfy initial values or boundary conditions. However, it still remains to be devised a general strategy for balancing the contribution between data interpolation and physical model satisfaction. Mean-squared errors (MSE) cost functions are analogous to maximizing the error likelihood if errors are normally distributed, uncorrelated, homoscedastic and of equal variance [14]. Consequently, the weights of MSE cost-function components in inverse PINNs problems consist of variance-weighted maximum-likelihood estimators (MLE) around uncorrelated homoscedastic normally distributed error terms. The need for supervised screening of such weighting hyperparameters adds a cumbersome layer
of complexity to the inverse problem solution, and prevents direct propagation of error between the loss function and the fitted parameters.

Although the applications of PINNs have spanned problems in time and/or spatial (ODEs and PDEs) domains (input spaces) [15, 16, 10, 17, 18], few works have assessed inverse problems of high-dimensionality in the state space (targets) [19, 20], which is an intrinsic feature of large interacting systems, such as chemical reactions and biological systems. The few exceptions are the forward solution of the Allen-Cahn, Hamilton-Jacobi-Bellman and Black-Scholes equations [21] in partial stochastic differential equations. Chemical reaction networks may include thousands of chemical species and reactions [22] in often stiff systems of coupled nonlinear differential equations. Studies involving high-dimensional inverse chemical kinetics problems using standard numerical methods have also been limited [23]. They are computationally expensive, requiring successive computation of ODE integrals to estimate gradients through adjoint or sensitivity methods [24, 25, 26]. PINNs might also serve as an auxiliary integration and derivative estimation method in conjunction with mechanism-discovery strategies based on symbolic regression. These include methods that generate many integral-form model candidates, e.g. ALAMO [27], require derivative estimates, e.g. SINDy [28], or successive numerical integration for validation [29]. Among the few studies involving NN-based approaches to chemical kinetics, De Florio et. al. used the "extreme theory of functional connection" to derive upper-bounds for generalization error in stiff forward PINNs. They tested their approach against the _POLLU_ problem [30], consisting of a stiff ODE encompassing \(20\) chemical species in homogeneous phase. Wu et. al. [31] explored the use of lasso-regularized neural ordinary differential equations (NODE) to enforce sparsity in inferring the mass-action rate law coupling matrix for homogeneous chemical systems involving six chemical species. However, studies involving PINNs for inverse complex reaction systems, especially in heterogeneous catalysis, are still scarce [32]. Such reacting systems are ubiquitous in the chemical industry and the ability to investigate the detailed mechanisms and elementary processes that underlie chemical reactions would allow the rigorous modeling, control and design of chemical reactors, and the optimization of industrial processes.

Kinetics-informed neural networks, KINNs, have been proposed by the authors and co-workers as a particular application of PINNs to high-dimensional heterogeneous chemical kinetics, which structurally tackles the issue of differential-algebraic equation (DAE) constraints and estimation of derivatives across the NN structure [32]. However, one of the main weakness of the original inverse KINNs formulation is that, as inverse PINNs in general, it consists of a multi-objective optimization problem. It required a priori estimation of a hyperparameter that conveyed the ratio between variances of state values and their derivatives. The inverse solution was shown to be a function of the chosen hyperparameter, leading to a Pareto frontiers of efficient solutions. A more effective approach is necessary to rigorously model the probability density function of residuals. The MLE formulation is advantageous since it may be applied to complex error distributions, e.g., full covariance with correlated error terms for Gaussian error [33], and its structure allows the propagation of variance between residuals. However, MLEs constructed around multivariate Gaussian distributions require the inversion of error covariance matrices, and correlated errors lead to ill-conditioning with respect to inversion. Better conditioning can be attained by finding a suitable preconditioner for the coupled differential equations [34]. The ODEs that describe the rates of change of species in a chemical system are dictated by material balances (conservation laws), which are conveyed as reaction stoichiometries [35, 36, 37]. For large chemical systems, these linear relationships between changes in concentrations of chemical species lead to inherently correlated errors in states and their derivatives, which results in ill-conditioned covariance matrices.

In this work, we show that the Taylor expansion of residuals in Gaussian MLEs allows direct error propagation between interpolation and physics models for high-dimensional ODE inverse problems, thereby removing the need for hyperparameter tuning in inverse PINNs, especially for large chemical systems (KINNs). Furthermore, we resort to the singular-value decomposition (SVD) as a framework to generate orthonormal bases for the range and nullspace of the ODE coupling matrix (reaction stoichiometry matrix). The range basis serves as preconditioner since the projection of KINNs residuals on it removes correlated terms from the error covariance matrix, which mitigate the issue of ill-conditioning with respect to inversion, allowing the estimation of precision matrices. In this robust KINNs formulation (rKINNs) we incorporate the three types of learning biases, as defined by Karniadakis et. al. [38]: (i) "observational bias" in the sampling of model error at points where observed data is not available, (ii) "inductive bias" in the conservation law defined by the stoichiometry matrix, the topological incorporation of the normalization operator that enforces DAE constraints, as well as in the representation of the ODEs in the preconditioning space obtained from SVD of the stoichiometry matrix, and (iii) "learning bias" in assuming model and interpolation residuals are normally distributed and constructing the optimization problem in terms of MLEs.
## 2 Results

### Theoretical Overview

#### 2.1.1 Mean-field chemical kinetics

In its conceptualization, KINNs was designed around mean-field kinetic models with power-law kinetics with Arrhenius temperature, \(\theta\), dependence for a set of chemical species whose concentrations (states) are given by an array \(\mathbf{c}\in\mathbb{R}_{+}^{n}\). A mean-field microkinetic model (MKM) is defined as in (2.1.1):

\[\mathbf{\dot{c}}=\frac{d}{dt}\mathbf{c}=f(\mathbf{c},\mathbf{p}(\theta))= \mathbf{M}\,\mathbf{r}(\mathbf{c},\mathbf{p}(\theta))=\mathbf{M}\,(\mathbf{k} (\mathbf{p}(\theta))\circ\psi(\mathbf{c})) \tag{2.1.1}\]

Where \(\mathbf{M}\in\mathbb{Z}^{n\times m}\) is the stoichiometry matrix (ODE coupling matrix), which can include reactions between chemical species possibly in different phases (e.g. surface and gas/liquid-phases), \(\psi(\cdot):\mathbb{R}_{+}^{n}\to\mathbb{R}_{+}^{m}\) represents the power-law kinetics mapping, and \(\ \mathbf{k}:=\{\mathbf{k}=\exp(-\mathbf{p}(\theta))\ \in\ \mathbb{R}_{m}^{m},\ \theta\in \mathbb{R}_{+},\ \mathbf{p}(\theta)\in\mathbb{R}^{m}\}\) denotes Arrhenius-like temperature dependent rate constant with \(\mathbf{p}\) as a vector-valued function that includes entropic and enthalpic contributions to a transition state.

#### 2.1.2 Inverse KINNs residuals

"Surrogate approximator" (SA) were introduced in KINNs [32] as a composite set of NNs that serve as basis functions with \(\mathbf{\omega_{\mathrm{s}}}\) weights for the solution of chemical kinetics ODEs, Eq. (2.1.1), with \(\mathbf{c}(t)\approx\mathbf{x}(t,\mathbf{\omega_{\mathrm{s}}})\), with the ability to structurally represent DAE constraints through an operator \(\mathbf{C_{N}}\) (section 4.3). Such an approach builds upon the standard PINNs concept [39, 40, 41], which is a direct result of the universal approximation theorem [42, 43, 44, 45]. The inverse problem is defined in terms of the multi-objective optimization function involving interpolation and physics-model residuals, \(\varepsilon_{\mathbf{x}_{i}}\) and \(\varepsilon_{\mathbf{x}_{i}}\in\mathbb{R}^{n}\) for \(i=1,2,...,d\) datapoints, respectively. The interpolation error, Eq. (2.1.2), is defined by the difference between measured states values, \(\tilde{\mathbf{x}}_{i}\in\mathbb{R}^{n}\), and estimated state values, \(\mathbf{x}(t_{i},\mathbf{\omega_{\mathrm{s}}})\in\mathbb{R}^{n}\), for every \(t_{i}\):

\[\varepsilon_{\mathbf{x}_{i}}(\mathbf{\omega_{\mathrm{s}}})=\mathbf{x}(t_{i},\bm {\omega_{\mathrm{s}}})-\tilde{\mathbf{x}}_{i} \tag{2.1.2}\]

such that if \(\tilde{\varepsilon}_{\mathbf{x}_{i}}(\mathbf{\omega_{\mathrm{s}}})=0\) then the SA goes exactly through data point \(i\). The physics-model residual, \(\varepsilon_{\mathbf{x}_{i}}\), relates the SA derivative and the output of differential equation (e.g. the kinetic model) in Eq. (2.1.3).

\[\varepsilon_{\mathbf{x}_{i}} =\hat{\mathbf{x}}(t_{i},\mathbf{\omega_{\mathrm{s}}})-f(\mathbf{x}_{i },\mathbf{p}) \tag{2.1.3}\] \[=\hat{\mathbf{x}}(t_{i},\mathbf{\omega_{\mathrm{s}}})-\mathbf{M}\,( \mathbf{k}(\mathbf{p}(\theta))\circ\psi(\mathbf{x}(t_{i},\mathbf{\omega_{\mathrm{ s}}})))\]

such that if \(\varepsilon_{\tilde{\mathbf{x}}_{i}}=0\) then the differential equation in Eq. (2.1.1) is exactly satisfied at time \(t_{i}\).

In the original regularized KINNs loss function, a hyperparameter \(\alpha\) weights the interpolation in relation to the model-fitting MSE, Eq. (2.1.4).

\[j_{t}(\mathbf{\omega_{\mathrm{s}}},\mathbf{p})\,\mathrm{MSE}(\mathbf{\varepsilon_{ \mathrm{x}}})+\alpha\,\mathrm{MSE}(\mathbf{\varepsilon_{\mathrm{x}}})=\frac{1}{d} \sum_{i}^{d}\mathbf{\varepsilon}_{\mathbf{x}_{i}}^{T}\mathbf{\varepsilon_{\tilde{ \mathbf{x}}_{i}}}+\alpha\frac{1}{d}\sum_{i}^{d}\mathbf{\varepsilon}_{\mathbf{x}_{i }}^{T}\mathbf{\varepsilon_{\mathbf{x}_{i}}} \tag{2.1.4}\]

The loss function, \(j_{t}\) is then minimized with respect to \(\mathbf{\omega_{\mathrm{s}}}\) and \(\mathbf{p}\) for a fixed \(\alpha\). The error in kinetic model parameters, \(\mathbf{\varepsilon_{\mathrm{p}}}\), is defined as:

\[\mathbf{\varepsilon_{\mathrm{p}}}^{*}=\mathbf{p}-\mathbf{p}^{*} \tag{2.1.5}\]

where \(\mathbf{p}\) represents the parameters of the physical model at any point in the optimization and \(\mathbf{p}^{*}\) represents the true model parameters. Unlike \(\varepsilon_{\mathbf{x}_{i}}\) and \(\mathbf{\varepsilon_{\tilde{\mathbf{x}}_{i}}}\), \(\mathbf{\varepsilon_{\mathrm{p}}^{*}}\) is not dependent on time, and cannot be evaluated explicitly from kinetic data since \(\mathbf{p}^{*}\) is not known in general. Rather, it is implicitly assumed that minimizing the KINNs loss function also minimizes \(\mathbf{\varepsilon_{\mathrm{p}}^{*}}\).

### Derivation of maximum-likelihood estimator

#### 2.2.1 Error probability density function

To circumvent the need for a regularization hyperparameter (\(\alpha\)) and enable rigorous error propagation, the original KINNs problem, Eq. (2.1.4), and most or all related PINN problems, can be reformulated in terms of maximum likelihood estimators, MLEs. The key idea of MLE is to assume an underlying probability distribution function, \(\pi(\varepsilon)\), for the error terms in the regression problem and to maximize the likelihood of observing the errors given the data. Let \(\mathbf{\varepsilon_{\mathrm{x}}}\), \(\mathbf{\varepsilon_{\mathrm{k}}}\), \(\mathbf{\varepsilon_{\mathrm{p}}}\) denote the errors related to the interpolation, physical model and its parameters, respectively. Assuming
errors are independent between timepoints \(i\), and that \(\mathbf{\varepsilon}_{\mathbf{x}}\) is normally distributed with covariance \(\mathbf{\Sigma}_{\mathbf{x}}\), the associated probability distribution function conditioned to the observed data \(\tilde{\mathbf{x}}\) and its covariance \(\mathbf{\Sigma}_{\mathbf{x}}\) is given in Eq. (2.2.1).

\[\pi(\cap_{i=1}^{n}\mathbf{\varepsilon}_{\mathbf{x}_{i}},\mathbf{\varepsilon}_{\dot{ \mathbf{x}}_{i}},\mathbf{\varepsilon}_{\mathbf{p}}|\mathbf{\Sigma}_{\mathbf{x}},\tilde{ \mathbf{x}}_{i})=\prod_{i}^{n}\pi(\mathbf{\varepsilon}_{\mathbf{x}_{i}},\mathbf{ \varepsilon}_{\dot{\mathbf{x}}_{i}},\mathbf{\varepsilon}_{\mathbf{p}}|\mathbf{\Sigma}_ {\mathbf{x}},\tilde{\mathbf{x}}_{i}) \tag{2.2.1}\]

The main advantage of the MLE formulation to KINNs is that it allows direct connection between errors associated with model-fitting and those related to data interpolation. With the independent variables \(\mathbf{x}\) and \(\mathbf{p}\), and assuming \(\hat{\mathbf{x}}(\mathbf{x})\) through \(\mathbf{\omega}_{\mathbf{s}}\), a variation in \(\mathbf{\varepsilon}_{\hat{\mathbf{x}}}\) can be represented by its total derivative:

\[\delta\mathbf{\varepsilon}_{\dot{\mathbf{x}}}=\underbrace{\partial_{\mathbf{x}} \hat{\mathbf{x}}\delta\mathbf{x}}_{\text{representation}}-\underbrace{(\partial_{ \mathbf{x}}f(\mathbf{x},\mathbf{p})\delta\mathbf{x}+\partial_{\mathbf{p}}f( \mathbf{x},\mathbf{p})\delta\mathbf{p})}_{\text{structural}}+\mathcal{O}(\delta \mathbf{x}^{2},\delta\mathbf{p}^{2}) \tag{2.2.2}\]

which represents the sensitivity of the model residuals, \(\delta\mathbf{\varepsilon}_{\dot{\mathbf{x}}}\), in a Taylor expansion as a function of variations in model parameter, \(\delta\mathbf{p}\), and states, \(\delta\mathbf{x}\). In Eq. (2.2.2), \(\partial_{\mathbf{x}}f(\mathbf{x},\mathbf{p})\) and \(\partial_{\mathbf{p}}f(\mathbf{x},\mathbf{p})\) are the kinetic model gradients with respect to states and physical model parameters, respectively, about \(\mathbf{x}\) and \(\mathbf{p}\), and they represent the contribution of the physical model to the variations in the residuals of the SA. The partial derivatives between the time derivatives of the SA with respect to its states, \(\partial_{\mathbf{x}}\hat{\mathbf{x}}\), are implicit functions of the NN architecture, i.e., depth, width, connectivity, choice of activation functions and, ultimately, their weights, \(\mathbf{\omega}_{\mathbf{s}}\), as in Eq. (2.2.3).

\[\partial_{\mathbf{x}}\hat{\mathbf{x}}=\frac{\partial\hat{\mathbf{x}}}{ \partial\mathbf{\omega}_{\mathbf{s}}}\frac{\partial\mathbf{\omega}_{\mathbf{s}}}{ \partial\mathbf{x}}=\frac{\partial\hat{\mathbf{x}}}{\partial\mathbf{\omega}_{ \mathbf{s}}}\left(\frac{\partial\mathbf{x}}{\partial\mathbf{\omega}_{\mathbf{s}}} \right)^{-1} \tag{2.2.3}\]

The universal approximation theorem extends the approximation of functions not only to their values but their gradients [46, 47]. Therefore, in the limit of sufficiently complex NNs, the models states, \(\mathbf{x}\), and the model derivatives, \(\hat{\mathbf{x}}\), will be independent. This is analogous to the case of explicit Galerkin methods that employ bases which enforce orthogonality between states and associated derivatives [48, 49]. However, for smaller NNs there should exist some correlation between states and their derivatives, with the ideal scenario being \(\partial_{\mathbf{x}}\hat{\mathbf{x}}\propto\partial_{\mathbf{x}}f\), where SA behavior mimics that of the physical model making a PINN approximate a NODE. Nonetheless, in this work we neglect the representation error, \(\partial_{\mathbf{x}}\hat{\mathbf{x}}\), based on the assumption that the feed-forward NN basis is large enough to allow states and their derivatives to be represented independently. With the assumption of normality distribution of residuals, the derivative residuals covariance matrix at a timepoint \(i\), \(\mathbf{\Sigma}_{\hat{\mathbf{x}}_{i}}\) can be expressed as projections of the interpolation covariance onto the gradients of the model at the respective timepoints, as in section 4.1. The joint probability density function (PDF) for the MLE approach is then given by the product of independent distributions:

\[\pi\left(\cap_{i=1}^{n}\mathbf{\varepsilon}_{\mathbf{x}_{i}},\mathbf{ \varepsilon}_{\dot{\mathbf{x}}_{i}},\mathbf{\varepsilon}_{\mathbf{p}}|\mathbf{\Sigma} _{\mathbf{x}},\tilde{\mathbf{x}}_{i}\right) \approx\ \pi(\cap_{i=1}^{n}\mathbf{\varepsilon}_{\mathbf{x}_{i}},\mathbf{ \varepsilon}_{\dot{\mathbf{x}}_{i}}|\mathbf{\Sigma}_{\mathbf{x}},\mathbf{\Sigma}_{ \dot{\mathbf{x}}_{i}},\tilde{\mathbf{x}}_{i},\mathbf{p}) \tag{2.2.4}\] \[=\prod_{i}^{n}\pi(\mathbf{\varepsilon}_{\dot{\mathbf{x}}_{i}}|\mathbf{ \Sigma}_{\dot{\mathbf{x}}_{i}})\pi(\mathbf{\varepsilon}_{\mathbf{x}_{i}}|\mathbf{ \Sigma}_{\mathbf{x}}) \tag{2.2.5}\]

In Eq. (2.2.5), we imply that model derivative and interpolation covariances are assumed constant for a given sample of evaluated model and model derivative errors. In section 2.2.2, we show how this assumption about the statistical model can be used within a maximum-likelihood estimator framework.

#### 2.2.2 Maximum likelihood estimation

Maximizing likelihood estimators are equivalent to minimizing the negative log-likelihood of error PDFs [50, 51]. In MLE form, Eq. (2.2.5) represents a summation, as in Eq. (2.2.6), where \(\pi\) are assumed to be Gaussian PDFs and \(\mathbf{\Omega}=\mathbf{\Sigma}^{-1}\) is the precision matrix, which is held constant within parameter optimization and, therefore, does not depend on the suppressed variables \(\mathbf{\omega}_{\mathbf{s}}\) and \(\mathbf{p}\). For brevity, we lump all parameters being optimized into a single variable, \(\mathbf{\omega}_{t}=[\mathbf{\omega}_{\mathbf{s}},\mathbf{p}]\), yielding.

\[\begin{split}\min_{\mathbf{\omega}_{t}}&\ell_{\mathrm{t}}= j_{\mathrm{t}}=\frac{1}{d}\sum_{i}^{d}\left(\mathbf{\varepsilon}_{\dot{\mathbf{x}}_{i}}^{T} \mathbf{\Omega}_{\dot{\mathbf{x}}_{i}}\mathbf{\varepsilon}_{\dot{\mathbf{x}}_{i}}+ \mathbf{\varepsilon}_{\dot{\mathbf{x}}_{i}}^{T}\mathbf{\Omega}_{\mathbf{x}}\mathbf{ \varepsilon}_{\mathbf{x}_{i}}\right)\\ \text{s.t.}&\mathbf{\Omega}_{\dot{\mathbf{x}}_{i}}= \left(\mathbf{\Sigma}_{\dot{\mathbf{x}}_{i}}^{\mathbf{x}}+\mathbf{\Sigma}_{\dot{ \mathbf{x}}_{i}}^{\mathbf{p}}\right)^{-1}\\ &\mathbf{\Omega}_{\mathbf{x}}\ =\mathbf{\Sigma}_{\mathbf{x}}^{-1}\\ &\mathbf{p}\in\mathbb{R}^{\dim(\mathbf{p})},\ \mathbf{\omega}_{\mathbf{s}}\in\mathbb{R}^{\dim(\mathbf{\omega}_{\mathbf{s}})}\end{split} \tag{2.2.6}\]
### Range-Nullspace Decomposition

Coupled ODEs describing rank-deficient dynamical system can be reduced by defining proper bases for states and their derivatives. The algebraic analysis of the stoichiometry coupling matrix for information recovery and data reconstructions has been focus of investigation in chemical and biological systems [52, 53]. Here we generalize the algebraic analysis of the coupling matrix, \(\mathbf{M}\) through it singular-value decomposition (SVD), which is determined by Eq. (2.3.1).

\[\mathbf{M}\triangleq\mathbf{U}\mathbf{S}\mathbf{V}^{\dagger} \tag{2.3.1}\]

Where \(\mathbf{U}\in\mathbb{R}^{n\times n}\), \(\mathbf{V}\in\mathbb{R}^{m\times m}\) are orthonormal basis of \(\mathbb{R}^{n}\) and \(\mathbb{R}^{m}\), respectively, and \(\mathbf{S}\in\mathbb{R}^{n\times m}\) is a matrix with singular values of \(\mathbf{M}\) along its diagonal entries and zeroes for the off-diagonal ones. Therefore, \(\mathbf{U}\) is a basis for concentrations and their derivatives over time, which can be broken down into \(\mathbf{U}^{\mathrm{R}}=\{u_{ij}\mid s_{jj}>0,\ \forall i\}\), which is a basis for the range of \(\mathbf{M}\), and \(\mathbf{U}^{\mathrm{N}}=\{u_{ij}\mid s_{jj}=0,\ \forall i\}\), which is a basis for the null-space or kernel of \(\mathbf{M}^{T}\).

\[\mathbf{U}=\begin{bmatrix}\mathbf{U}^{\mathrm{R}}&\mathbf{U}^{\mathrm{N}} \end{bmatrix} \tag{2.3.2}\]

Any net non-zero rate of reaction must be in \(\mathbf{U}^{\mathrm{R}}\), while the boundary conditions and inner fluxes within the chemical system lie in the null-space \(\mathbf{U}^{\mathrm{N}}\), allowing the concentration vector to be decomposed into a range and nullspace component:

\[\mathbf{c}_{(t)} =\mathbf{c}_{(t)}^{\mathrm{R}}+\mathbf{c}^{\mathrm{N}} \tag{2.3.3}\] \[\mathbf{c}_{(t)}^{\mathrm{R}} =\mathbf{U}^{\mathrm{R}}(\mathbf{U}^{\mathrm{R}})^{T}\mathbf{c}_{ (t)} =\mathbf{U}^{\mathrm{R}}\mathbf{z}_{(t)}^{\mathrm{R}}\] (2.3.4) \[\mathbf{c}^{\mathrm{N}} =\mathbf{U}^{\mathrm{N}}(\mathbf{U}^{\mathrm{N}})^{T}\mathbf{c}_{ (t)} =\mathbf{U}^{\mathrm{N}}\mathbf{z}^{\mathrm{N}} \tag{2.3.5}\]

where \(\mathbf{z}_{(t)}^{\mathrm{N}}\) and \(\mathbf{z}_{(t)}^{\mathrm{R}}\) are the subspace projections of \(\mathbf{c}\) onto the null-space and range, respectively. As for concentration derivatives, they exclusively reside in the range of \(\mathbf{M}\), and can be decomposed using the projections as follows.

\[\hat{\mathbf{c}}_{(t)} =\hat{\mathbf{c}}_{(t)}^{\mathrm{R}}=\mathbf{U}^{\mathrm{R}}( \mathbf{U}^{\mathrm{R}})^{T}\hat{\mathbf{c}}_{(t)}=\mathbf{U}^{\mathrm{R}} \hat{\mathbf{z}}_{(t)}^{\mathrm{R}} \tag{2.3.6}\] \[\mathbf{0} =(\mathbf{U}^{\mathrm{N}})^{T}\hat{\mathbf{c}}_{(t)}=\hat{ \mathbf{z}}^{\mathrm{N}} \tag{2.3.7}\]

The non-zero singular values from SVD are associated with the subspace over which all the model derivative variance is confined. In section 4.7, we show that the SA can be built such that the time dependence is only contained in the span of \(\mathbf{U}^{\mathrm{R}}\), with a time-independent shift in \(\mathbf{U}^{\mathrm{N}}\). Hence, the projection of Eq. (2.1.2) onto \(\mathbf{U}^{\mathrm{N}}\) will lead to a time-independent error \(\hat{\mathbf{c}}_{\mathbf{z}}^{\mathrm{N}}\in\mathbf{U}^{\mathrm{N}}\) between the SA estimates and the observed data. Assuming that \(\mathbf{z}^{\mathrm{N}}\) can be estimated from the observed data, the nullspace solution can be set a priori as \(\mathbf{z}^{\mathrm{N}}=\left\langle\left(\mathbf{U}^{\mathrm{N}}\right)^{T} \tilde{\mathbf{x}}\right\rangle\). In Eq. (2.2.6), it is assumed that \(\mathbf{\Omega}_{\mathbf{x}}\) can be defined as \(\mathbf{\Sigma}_{\mathbf{x}}^{-1}\); however, when \(\mathbf{M}\) has zero singular values, \(\mathbf{\Sigma}_{\mathbf{x}}\) is rank-deficient and, hence, ill-conditioned with respect to inversion. Projecting \(\mathbf{\Sigma}_{\mathbf{x}}\) onto \(\mathbf{U}^{\mathrm{R}}\) provides a route to preserve the total variance and allows the inversion of the projected matrix, i.e, if \(\left(\mathbf{U}^{\mathrm{R}}\right)^{T}\mathbf{\Sigma}_{\mathbf{x}}\mathbf{U} ^{\mathrm{R}}=\mathbf{\Sigma}_{\mathbf{z}}\) is full-rank, \(\mathbf{\Omega}_{\mathbf{z}}=\left(\mathbf{\Sigma}_{\mathbf{z}}\right)^{-1}\) can be computed. The same rationale holds for \(\mathbf{\Sigma}_{\mathbf{\hat{x}}}\), since derivative residuals lie only in \(\mathbf{U}^{\mathrm{R}}\).

### Case Studies and Discussions

#### 2.4.1 Representative 'Deep' Microkinetic Model

Here, we utilize the most complex of the anecdotal test cases from the prior KINNs study [32], which includes latent species and reactions between them. The elementary steps are given in Eq. (2.4.1), and the corresponding matrix representation is shown in Supplementary Material, S5.
Heterogeneous Reaction Network

\[\begin{split} A&+*\frac{k_{\ast}^{d}}{k_{\ast}^{d}}A*\\ B&+*\frac{k_{\ast}^{d}}{k_{\ast}^{d}}B*\\ C&+*\frac{k_{\ast}^{d}}{k_{\ast}^{d}}C*\\ A&*+*\frac{k_{\ast}^{e}}{k_{\ast}^{e}}2D*\\ B&*+*\frac{k_{\ast}^{e}}{k_{\ast}^{e}}2E*\\ D&*\;+E*\frac{k_{\ast}^{t}}{k_{\ast}^{t}}F*\;+*\\ F&*\;+E*\frac{k_{\ast}^{e}}{k_{\ast}^{e}}C*\;+* \end{split} \tag{2.4.1}\]

An overview of the generated synthetic data is shown in Fig. 1, where each initial condition is referred to as an "experiment", and data is separated by gas-phase concentrations (**x\({}_{\text{o}}\)**) that can be directly measured and surface concentrations (**x\({}_{\text{*}}\)**) that can only be measured indirectly, which is a common challenge in the dynamics of heterogeneous catalytic systems, i.e. the existence of "latent" states arising from adsorbed surface species. The calibrated signals for surface concentrations were obtained from the raw synthetic data by using the closed-form solution for surface reconstruction defined in Eq. (S3.5), section S3, with an eigenvalue cutoff of \(5\times 10^{-3}\) to find the calibration coefficients \(\gamma\), which are shown in the parity plot in Fig. S3.1. The coverages are normalized by construction, according to Eq. (S3.3), and satisfy the nullspace time invariance, which expressed by the differential equation in \(z\)-space, Fig. 1. Therefore, for the synthetic experiments under study, the bulk- and surface-phases dynamics can be reconstructed directly rather than learned through NN training.

Figure 1: Bulk phase (a and e) and calibrated surface coverage fraction (b and f); **z**-Space: range, **z\({}^{\text{R}}\)** (c and g), and nullspace, **z\({}^{\text{N}}\)** (d and h), projections, where all variance is contained in the range projection and time invariance in the nullspace, for IC1, \(\{x_{A},x_{B},x_{C},x_{\ast}\}_{t_{0}}=\{0.6,0.4,0.0,1.0\}\), and IC2, \(\{x_{A},x_{B},x_{C},x_{\ast}\}_{t_{0}}=\{0.2,0.3,0.5,1.0\}\).
#### 2.4.2 Application and Comparison of naive KINNs and robust KINNs

The range-nullspace decomposition in Fig. 1 leads to the dimensionality reduction in the inverse problem, from a total of initially nine degrees of freedom, i.e. three bulk-phase species, and six surface species (adsorbates), while the free-sites balance "\(\ast\)" degree-of-freedom is accounted for by the normalization operator \(\mathrm{C_{\mathbf{N}}}\), only the seven range components are taken into consideration in the optimization. Therefore, the KINNs output layer included three linear components for bulk phases and seven DAE-constrained components with underlying NN with six outputs for surface species. rKINNs' SA only maps the seven components of the range space, while the three nullspace components are a priori estimated.

In the naive approach, KINNs, the hyperparameter \(\alpha\) controls the expected ratio between interpolation (data) and model (physics) variances. From the MLE perspective, this represents a (naive) assumption that the underlying interpolation and model residuals are normally distributed with covariances \(\mathbf{\Sigma_{x}}=\sigma_{\mathbf{x}}\mathbf{I}\) and \(\mathbf{\Sigma_{\hat{x}}}=\sigma_{\hat{x}}\mathbf{I}\), with \(\alpha=\sigma_{\hat{x}}/\sigma_{\mathbf{x}}\), such that smaller values of \(\alpha\) denote lower variance for the model residual as compared to the interpolation residuals, and the opposite for larger values of \(\alpha\). In Fig. 2 (a), the inverse naive-KINNs are applied to the case-study system. The data- and model-residuals evolution as a function of \(\alpha\) define a Pareto-set of efficient solutions. A parity plot for the kinetic model parameters at the inflection points depicted by open (tightening) and closed (relaxation) diamonds in Fig. 2 (a) is shown in (b). Optimal kinetic model parameters are associated with the "elbows" or inflection points in the relaxation direction. Increasing \(\alpha\) to larger values leads to SA overfitting to data, whereas reducing it indefinitely leads to SA "forgetting" the data, which leads to a scenario of producing degenerate physics solutions [32]. The process of optimizing \(\alpha\) through multi-objective optimization is computationally expensive, since each iteration requires 300 epochs of NN training, and leads to ambiguity in the solution due to hysteresis in the parameters as a function of \(\alpha\).

Fig. 2 (c) illustrates the main advantage of rKINNs over KINNs. Interpolation (\(\ell_{\mathbf{x}}\)) and model (\(\ell_{\mathbf{\hat{x}}}\)) negative log likelihoods were computed from the interpolation residuals projected onto the final (stationary) covariance matrix obtained by rKINNs, and model residuals per time point according to Eq. (4.7.3), respectively. In addition to the faster per-epoch convergence, rKINNs reaches a stable or stationary point (Fig. 2 (e), black-filled circle). Such a behavior is attributed to the stationarity in the SA interpolation and, thus, of the sample covariance, \(\hat{\Sigma}_{\mathbf{x}}\), and its propagation through the model gradients, \(\hat{\Sigma}_{\hat{\mathbf{x}}_{i}}\), over time points \(t_{i}\). Different choices of hidden-layers activation functions sets have led to indistinguishable results for stationary model parameters and confidence intervals for rKINNs, e.g., \(\{20_{\times 3}\}\) full swish \(\{swish_{\times 3}\}\), \(\{100_{\times 1}\}\) with \(\{rbf\}\), thus highlighting the robustness of rKINNs to not only weight the variances between interpolation and physics, but also to converge to similar convex basins irrespective of the choice of activation functions (basis functions) and SA architecture.

In Fig. 3, synthetic data (solid circles), SA estimates (open circles, dashed line) and physical model integral based on the initial conditions estimated from SA (solid lines) are plotted for bulk species (a and d) and adsorbates (b and e). Additionally, the physical model, \(f(\mathbf{x})\), and SA algorithmic differentiation, \(\hat{\mathbf{x}}\), parity plots for all species are shown in Fig. 3 (c) and (f). The results illustrate a strong agreement between the numerical integration of the model with recovered parameters (solid lines) and the synthetic experimental data for both experiments.

## 3 Discussion

In this work, we introduce a hyperparameter-free approach for inverse PINNs problems applied to high-dimensional coupled non-linear ODEs. We express the probability density function for the residuals of PINNs applied to heterogeneous chemical kinetics systems, KINNs [32], as normal distributions. We reformulate the KINNs regularized MSE cost functions as MLEs with direct propagation of SA interpolation error through the physical (kinetic) model gradients. The model gradients act as projectors to connect the cost functions for data interpolation and physical derivatives to eliminate the hyperparameter that is often required in the in inverse PINNs. In addition, we address the issue of ill-conditioning of covariance matrices by using SVD of the coupling matrix to find uncorrelated subspaces that function as pre-conditioners over which the residuals and covariance matrices may be projected. SVD also allows the decomposition of the SA representation into the range (time variant) and nullspace (time invariant) components of the stoichiometry matrix, thus leading to an intrinsic dimensionality reduction at the representation level, where only the range components are mapped by SAs' NNs. Such a decomposition also allows a closed-form solution for the calibration factors needed to relate semi-quantitative signals of latent adsorbate populations to explicit adsorbate concentrations.

The novel rKINNs framework with MLE-estimators and the range-nullspace decomposition is applied to a complex case study in heterogeneous catalysis. It efficiently converges to a well-defined solution that is consistent with the results of expensive multi-objective optimization with the prior naive-KINN approach. Moreover, the rKINN formalism allows for rigorous propagation of error between the model and parameter spaces, enabling calculation of error bars on parameter estimates that are robust to noise in the data and the architecture of the underlying NNs The direct connection
between interpolation and physics done through the construction of "physically-coupled" residuals probability density functions. This approach provides a general self-regularized framework for inverse PINNs problems based on error propagation through MLE forms, which are applicable to high-dimensional coupled ODEs. Even in the case of highly correlated outputs, SVD - or other orthogonal decomposition methods - can be used to define orthogonal minimal subspaces for the residuals projection and SA representation and robust evaluation of MLEs as NNs weights are trained. This framework applications span other types of reactor models, such as: homogeneous continuous flow stirred-tank reactor equations would include information about the flows at the boundary in the nullspace and, hence, may not be time invariant yet still separable from the kinetics, and uniform plug-flow reactors differential equations, which requires a change in variable over which it is integrated. Finally, additional effort would be necessary to expand the application of the concepts here devised to non-uniform transient systems, where partial differential equations would need to be solved.

One limitation of the current rKINN approach is that semi-quantitative signals for all latent species are required. However, it common in heterogeneous catalysis the lack of information about some or all latent species. Additional work is required to address this scenario, where some assumptions about the covariance structure of unobservable species will be required. This is complicated by the fact that the uniqueness of the solution to the inverse problem with incomplete information is not guaranteed, and further studies are needed to establish the limits of solutions obtained from inverse problems in the limit of missing information or extreme stiffness. Moreover, in the case study presented

Figure 2: Inverse KINNs vs rKINNs - (a) sensitivity analysis of naive KINNs regularization parameter \(\alpha\) with Pareto set estimate as a function of interpolation and model MSEs; highlighted tightening (increasing \(\alpha\), open diamond) and relaxation (decreasing \(\alpha\) closed diamond) inflection points; (b) naive KINNs model parameters parity plot at inflection points in highlighted in (a); (c) naive KINNs residuals represented in terms of negative log-likelihoods (\(\alpha\)-colored circles) and rKINNs (open circles) log-likelihoods reconstructed in terms of the final stable point (black-filled circle in the inset of (c)); (d) rKINNs model parameters at the stable point (error bars defined as \(2\) standard deviations based on the Fischer information estimated through Hessian analysis with algorithmic differentiation). One epoch is defined as \(100\) parameter-update iterations.
it is assumed the physical model structure is fully known, which is often not the case in heterogeneous catalysis. This can potentially be overcome by coupling mechanism discovery algorithms [31, 54, 55] with (r)KINNs or other frameworks for solving the inverse problem, or by the use of NODEs that can learn surrogate functions for the system dynamics [56, 57, 58, 59].

## 4 Methods

### Probability density function

We assume that variations about \(\mathbf{x}\), and incumbent parameters, \(\mathbf{p}\), in Eq. (2.2.2) are random variables that follow independent multivariate normal distributions:

\[\begin{split}\delta\mathbf{x}\in\hat{\boldsymbol{\varepsilon}}_ {\mathbf{x}}&\sim\mathcal{N}(0,\boldsymbol{\Sigma}_{\mathbf{x}}) \\ \delta\mathbf{p}\in\hat{\boldsymbol{\varepsilon}}_{\mathbf{p}}& \sim\mathcal{N}(0,\boldsymbol{\Sigma}_{\mathbf{p}})\end{split} \tag{4.1.1}\]

With these assumptions, we can show that the variations about \(\hat{\boldsymbol{\varepsilon}}_{\hat{\mathbf{x}}_{i}}\) also comprise a multivariate normal distribution, whose covariance matrix can be estimated by the outer product of Eq. (2.2.2) at timepoint \(i\). Neglecting \(\mathcal{O}(\boldsymbol{\varepsilon}^{2})\) terms and assuming orthogonality between \(\boldsymbol{\varepsilon}_{\mathbf{x}}\) and \(\boldsymbol{\tilde{\varepsilon}}_{\mathbf{p}}\) yields:

\[\begin{split}\boldsymbol{\Sigma}_{\hat{\mathbf{x}}_{i}}& =\{\boldsymbol{\Sigma}_{\hat{\mathbf{x}}_{i}}|\mathbf{x}_{i}, \hat{\mathbf{x}}_{i},\mathbf{p}\}=\langle\delta\boldsymbol{\varepsilon}_{ \hat{\mathbf{x}}_{i}}\otimes\delta\boldsymbol{\varepsilon}_{\hat{\mathbf{x}}_ {i}}\rangle\\ &=\langle(-\partial_{\mathbf{x}}f_{i}\hat{\boldsymbol{\varepsilon }}_{\mathbf{x}}-\partial_{\mathbf{p}}f_{i}\boldsymbol{\tilde{\varepsilon}}_{ \mathbf{p}})\otimes(-\partial_{\mathbf{x}}f_{i}\boldsymbol{\tilde{\varepsilon }}_{\mathbf{x}}-\partial_{\mathbf{p}}f_{i}\boldsymbol{\tilde{\varepsilon}}_{ \mathbf{p}})\rangle\\ &=\partial_{\mathbf{x}}f_{i}\boldsymbol{\hat{\varepsilon}}_{ \mathbf{x}}\otimes\boldsymbol{\varepsilon}_{\mathbf{x}}\left(\partial_{ \mathbf{x}}f_{i}\right)^{T}+\partial_{\mathbf{p}}f_{i}\boldsymbol{\tilde{ \varepsilon}}_{\mathbf{p}}\otimes\boldsymbol{\tilde{\varepsilon}}_{\mathbf{p}} \left(\partial_{\mathbf{p}}f_{i}\right)^{T}\text{ for }\hat{\boldsymbol{\varepsilon}}_{\mathbf{x}} \perp\hat{\boldsymbol{\varepsilon}}_{\mathbf{p}}\\ &=\partial_{\mathbf{x}}f_{i}\boldsymbol{\Sigma}_{\mathbf{x}} \left(\partial_{\mathbf{x}}f_{i}\right)^{T}+\partial_{\mathbf{p}}f_{i} \boldsymbol{\Sigma}_{\mathbf{p}}\left(\partial_{\mathbf{p}}f_{i}\right)^{T}\\ &=\boldsymbol{\Sigma}_{\hat{\mathbf{x}}_{i}}^{\mathbf{x}}+ \boldsymbol{\Sigma}_{\hat{\mathbf{x}}_{i}}^{\mathbf{p}}\end{split} \tag{4.1.2}\]

where the arguments of the physical model \(f\) are suppressed for brevity, and \(f_{i}=f(\mathbf{x}_{i},\mathbf{p})\). This implies that the perturbations to the model residual can also be considered as a random variable that follows a multivariate normal

Figure 3: KINNs results: synthetic data (solid markers), surrogate approximator (open markers) down-sampled to every two datapoints, and numerical integral with initial conditions estimated from SA at \(t=0\) for bulk-phase (a and b) and adsorbates (b and e). Parity plot of standardized derivatives for all chemical species (c and f) from the kinetic model \(f(\mathbf{x})\) and from algorithmic differentiation of the surrogate approximator (\(\hat{\mathbf{x}}\)). Physical model and SA derivatives were standardized based on the mean standard deviation per state to facilitate the visualization of their correlation.
distribution:

\[\delta\epsilon_{\mathbf{\hat{x}}_{i}}\in\hat{\mathbf{\varepsilon}}_{\mathbf{\hat{x}}_{ i}}\sim\mathcal{N}(0,\mathbf{\Sigma}_{\mathbf{\hat{x}}_{i}}) \tag{4.1.3}\]

Notably, the covariance matrix for the underlying model parameter error distribution is generally unavailable, since the true parameters, \(\mathbf{p}\), are unknown. However, Eq. (2.2.2) can be re-arranged to provide a least-squares approximation for \(\epsilon_{\mathbf{p}}\), where the current evaluated states and derivative errors, \(\epsilon_{\mathbf{x}}\) and \(\epsilon_{\mathbf{\hat{x}}}\) are considered as samples drawn from the underlying random variables defined in Eq. (4.1.1) and Eq. (4.1.3), and the representation contribution in \(\partial_{\mathbf{x}}\mathbf{\hat{x}}\) is neglected since model parameter errors should not depend on the SA:

\[\epsilon_{\mathbf{p}_{i}}=\left(\partial_{\mathbf{p}}f_{i}^{T}\partial_{ \mathbf{p}}f_{i}\right)^{-1}\partial_{\mathbf{p}}f_{i}^{T}\left(\epsilon_{ \mathbf{\hat{x}}_{i}}-\partial_{\mathbf{x}}f_{i}\epsilon_{\mathbf{x}_{i}}\right) \tag{4.1.4}\]

We assume that the covariance of the data interpolation can be obtained from the residuals in \(\mathbf{x}\). Letting \(\epsilon_{\mathbf{x}}\) be a sample of the random variable \(\hat{\mathbf{\varepsilon}}_{\mathbf{x}}\) (where \(\epsilon_{\mathbf{x}}\) is defined by Eq. (2.1.2)) yields:

\[\mathbf{\Sigma}_{\mathbf{x}}=\langle\epsilon_{\mathbf{x}}-\langle\epsilon_{ \mathbf{x}}\rangle\otimes\epsilon_{\mathbf{x}}-\langle\epsilon_{\mathbf{x}} \rangle\rangle \tag{4.1.5}\]

In this case, "recentering" (subtracting the sample mean when estimating sample covariances) is not strictly necessary; however, to avoid numerical instabilities from large initial residuals, intermediate estimated covariances are recentered during the KINNs training, which does not affect the final results given that the fundamental structural assumption is that \(\langle\epsilon_{\mathbf{x}}\rangle\rightarrow\mathbf{0}\) as the optimization proceeds. Notably, Eq. (4.1.4) yields a time-dependent parameter error, \(\epsilon_{\mathbf{p}_{i}}\), based on the time-dependence of \(\epsilon_{\mathbf{x}}\) and \(\epsilon_{\mathbf{\hat{x}}}\), and is hence different from the true time-independent parameter error, \(\epsilon_{\mathbf{p}}^{*}\) (Eq. (2.1.5)). However, since we are treating \(\epsilon_{\mathbf{x}}\) and \(\epsilon_{\mathbf{\hat{x}}}\) as independent samples from underlying multivariate normal distributions, the interpretation is that \(\epsilon_{\mathbf{p}_{i}}\) is also a sample from \(\hat{\mathbf{\varepsilon}}_{\mathbf{p}}\) (Eq. (4.1.1)) that obeys the structure imposed by Eq. (4.1.4). Thus, the model parameter covariance matrix can be estimated as \(\mathbf{\Sigma}_{\mathbf{p}}=\langle\epsilon_{\mathbf{p}}-\langle\epsilon_{ \mathbf{p}}\rangle\otimes\epsilon_{\mathbf{p}}-\langle\epsilon_{\mathbf{p}} \rangle\rangle\), with \(\langle\epsilon_{\mathbf{p}}\rangle\rightarrow\mathbf{0}\) as the optimization reaches a critical point.

### Covariance inversion stabilization

In the MLE formulation, a point of attention is that often the data- and model-parameter covariance structures are unknown, either requiring a priori assumptions or the utilization of methods for covariance inference from the data itself [60, 61] or from sample residuals [62, 63, 64]. In this work, we assume a lack of prior knowledge, and thus that the interpolation error is homoscedastic. The precision matrix is then updated periodically throughout the optimization at each epoch. To stabilize successive covariance inversions and account for the fact that residuals are not initially centered, the mean squared error is added as uncorrelated additive Gaussian noise to the covariances, Eq. (4.2.1), which is similar to the strategy adopted to stabilize Gaussian Process Regression [65]. As \(\langle\epsilon_{\mathbf{x}}\rangle\to 0\) at the critical points of \(\ell_{\mathrm{t}}\), this approach will not affect the final solution of the inverse problem.

\[\mathbf{\Sigma}_{\mathbf{x}/\hat{\mathbf{x}}}\leftarrow\mathbf{\Sigma}_{ \mathbf{x}/\hat{\mathbf{x}}}+\mathrm{diag}(|\langle\epsilon_{\mathbf{x}} \rangle|) \tag{4.2.1}\]

As the model converges this stabilization term reduces, and can be removed prior to error estimation to provide more accurate error estimates.

### Neural-network representation of heterogeneous systems

Heterogeneous chemical systems comprise unbound (gas-phase), adsorbed molecules and intermediates. Unbound species (observable) and adsorbate or intermediate - latent - species are represented by individual NNs \(\mathbf{\hat{x}}_{\mathbf{o}}\) and \(\mathbf{\hat{x}}_{\mathbf{*}}\), respectively, where a normalization operator \(\mathrm{C}_{\mathrm{N}}\) was introduced to enforce normalization of surface species (DAE constraint) by construction, i.e., \(0\leq\mathbf{x}_{\mathbf{*}}=\mathrm{C}_{\mathbf{N}}[\mathbf{\hat{x}}_{\mathbf{* }}]\leq 1\) and \(|\mathbf{x}_{\mathbf{*}}|_{1}=1\ \forall\ \mathbf{\hat{x}}_{\mathbf{*}}\in\mathbb{R}^{\mathrm{dim}( \mathbf{x}_{\mathbf{*}})-1}\). In this work, we have reformulated the normalization operator \(\mathrm{C}_{\mathbf{N}}\) from its previous trigonometric description to a simpler structure in terms of products of logistic functions that enforces bijection (Supplementary Material, S1).

\[\mathbf{x}=\begin{bmatrix}\mathbf{\hat{x}}_{\mathbf{o}}\\ \mathrm{C}_{\mathbf{N}}[\mathbf{\hat{x}}_{\mathbf{*}}]\end{bmatrix} \tag{4.3.1}\]

### KINNs architecture

In this work, derivatives are obtained through algorithmic differentiation (AD) with just-in-time (JIT) compilation using the JAX library [66, 67]. In addition, SA training and model fitting are trained using the _Adam_ optimization algorithm. The code necessary for reproducing the results of this work can be found on GitHub ([https://github.com/gusmaogabriels/kinn/tree/rkinns](https://github.com/gusmaogabriels/kinn/tree/rkinns)).

The MLE-formulation was utilized to build the cost function to be minimized for learning the ODE solutions using rKINNs. The JAX library provided the routines for the forward AD, allowing the obtainment of the gradients with respect to the SAs' parameters, which were optimized with the _Adam_ algorithm [68]. For both KINNs and rKINNs the underlying SA NNs consists of three-hidden layers, each of which with \(20\) neurons \(\{20_{\times 3}\}\) and corresponding activation functions \(\{tanh,\,swish,\,tanh\}\).
### Data Generation

Two sets of initial conditions (ICs) were used to generate synthetic data, which were built to represent scenarios of initially high and low concentration of reactants and products: IC1, \(\left\{x_{A},x_{B},x_{C},x_{*}\right\}_{t_{0}}=\left\{0.6,0.4,0.0,1.0\right\}\), and IC2, \(\left\{x_{A},x_{B},x_{C},x_{*}\right\}_{t_{0}}=\left\{0.2,0.3,0.5,1.0\right\}\), respectively. Additive homoscedastic noise with \(\sigma=0.025\) was added to the generated data. Signals for surface coverages were re-normalized to simulate unknown calibration factors, although it is assumed that semi-quantitative signals for all surface states are known. Numerical integrations for synthetic data generation were carried out with the stiff-nonstiff algorithm LSODA from the FORTRAN ODEPACK [69, 70]. Forward and reverse rate constants conveyed in Table 1 were chosen to depict over two orders-of-magnitude variance in rate constants that typical of reaction networks in heterogeneous catalysis. A total of \(100\) time-points were sampled in logarithmic-space, i.e., more points were sampled at the beginning of the experiments, where derivatives exhibit higher variance.

### Latent states

In many complex dynamical systems the species can be decomposed into measurable states and "latent" states that cannot be directly observed. This is the case in heterogeneous catalysis, where unbound gas phase species are typically observable, and denoted by \(\mathbf{c_{o}}\), while \(\mathbf{c_{a}}\) denotes the concentration of the latent variables, i.e., adsorbed and intermediate surface species. In this case, the range and null-space basis can be distinguished as:

\[\mathbf{U}=\left[\begin{bmatrix}\mathbf{U_{o}}\\ \mathbf{U_{\star}}\end{bmatrix}^{\mathrm{R}}\quad\begin{bmatrix}\mathbf{U_{o} }\\ \mathbf{U_{\star}}\end{bmatrix}^{\mathrm{N}}\right] \tag{4.6.1}\]

This leads to partitioned concentration vectors in Eq. (4.6.2).

\[\begin{bmatrix}\mathbf{c_{o}}\\ \mathbf{c_{\star}}\end{bmatrix}_{(t)}=\begin{bmatrix}\mathbf{c_{o}}\\ \mathbf{c_{\star}}\end{bmatrix}^{\mathrm{R}}_{(t)}+\begin{bmatrix}\mathbf{c_{o }}\\ \mathbf{c_{\star}}\end{bmatrix}^{\mathrm{N}} \tag{4.6.2}\]

Latent dynamics are characterized by incomplete information about the evolution of a subset of states, which can be partially "hidden" or indirectly inferable through the calibration of intensity measurements in conjunction with known associated conservation laws. For example, in heterogeneous catalysis, operando analytical methods, such as temperature-programmed reaction spectroscopy, may be utilized to track catalyst surface composition over time-on-stream [71]. The raw spectroscopic data (signal), \(\mathbf{y}\) is linear correlated with the latent states \(\mathbf{c_{\star}}\), as in Eq. (4.6.3), with calibration factors, \(\boldsymbol{\gamma}\), based on the conservation of total active sites constraints.

\[\mathbf{c_{\star}}=\mathbf{y_{\star}}\circ\boldsymbol{\gamma},\quad\mathrm{ s.t.}\ \mathbf{1}^{T}\mathbf{c_{\star}}=1 \tag{4.6.3}\]

In section S3, we use the null-space equations to a derive a closed-form solution for the calibration factor \(\boldsymbol{\gamma}\) in the scenario where \(\mathbf{y_{\star}}\) signal can be retrieved.

### Robust KINNs - MLE & Range-Nullspace decomposition

In Eq. (2.1.1), the underlying ODE can be decomposed into its components in the nullspace (\(\mathrm{N}(\mathbf{M})\)) and range (\(\mathrm{R}(\mathbf{M})\)), and SVD provides a basis for the representation of states and its derivatives in those spaces. Estimated variances can be propagated and separated into nullspace and range components. Such a decomposition eliminates linearly-dependent terms, creates denser covariance matrices, and mitigates the propagation of numerical errors to the precision matrices due to ill-conditioned covariance matrix-inversions. The SA representation is constructed to structurally separate the orthogonal components of the time-independent (nullspace) and -dependent (range) solutions, Fig. 4. By combining Eqs. (2.3.3) to (2.3.5) and Eq. (4.6.1), the SA is built as a function of the nullspace and range coverage solutions, i.e. \(\mathrm{C_{\mathbf{N}}[\cdot]}\) in Eq. (4.7.1) and Eq. (4.7.2), respectively. The time-invariant solution \(\mathbf{z}^{\mathrm{N}}\) consists of the sum of a solution in the

\begin{table}
\begin{tabular}{l c c c c c c c} \hline \hline  & \multicolumn{6}{c}{Rate Constants} \\ \cline{2-9} Type & \(\frac{k_{1}^{d}}{k_{1}^{d}}\) & \(\frac{k_{2}^{d}}{k_{2}^{d}}\) & \(\frac{k_{3}^{d}}{k_{3}^{d}}\) & \(\frac{k_{3}^{c}}{k_{3}^{c}}\) & \(\frac{k_{4}^{c}}{k_{4}^{c}}\) & \(\frac{k_{3}^{c}}{k_{3}^{c}}\) & \(\frac{k_{1}^{c}}{k_{1}^{c}}\) \\ \hline forward & \(20\) & \(24\) & \(16\) & \(640\) & \(160\) & \(560\) & \(640\) \\ reverse & \(8\) & \(12\) & \(40\) & \(960\) & \(80\) & \(160\) & \(240\) \\ \hline \hline \end{tabular}
\end{table}
Table 1: Forward and reverse reactions rate constants for the complex reaction mechanism (_dcs_) synthetics data generation.
convex-hull of \((\mathbf{U}_{\bullet}^{\mathrm{N}})^{T}\), since \(\mathrm{C}_{\mathrm{N}}[\mathbf{\hat{x}}_{\bullet}^{\mathrm{N}}]\) implies a convex-affine combination of its basis vectors, translated by solutions in the positive linear-span of \((\mathbf{U}_{\bullet}^{\mathrm{N}})^{T}\), which parameterized by \(\mathbf{\hat{x}}_{\mathbf{o}}^{\mathrm{N}}\) through \(\mathbf{U}_{\bullet}^{\mathrm{N}:\mathrm{N}}\). A visual representation of how the all subspaces are interconnected is shown in Fig. 4.

\[\mathbf{z}^{\mathrm{N}}=\left(\mathbf{U}_{\bullet}^{\mathrm{N}}\right)^{+} \mathrm{C}_{\mathrm{N}}[\mathbf{\hat{x}}_{\bullet}^{\mathrm{N}}]+\mathbf{U}_{ \bullet}^{\mathrm{N}:\mathrm{N}}\mathbf{\hat{x}}_{\mathbf{o}}^{\mathrm{N}} \tag{4.7.1}\]

The range-component solution depends on the nullspace estimates for the latent variables, \(\mathbf{x}_{\bullet}^{\mathrm{N}}=\mathbf{U}_{\bullet}^{\mathrm{N}}\mathbf{z} ^{\mathrm{N}}\). Similarly to the nullspace solution, to enforce normalization, \(\mathbf{x}_{\bullet}=\mathrm{C}_{\mathrm{N}}[\mathbf{\hat{x}}_{\bullet(t)}^{ \mathrm{R}}]\) is first calculated and variations in the latent species solutions \(\mathbf{z}^{\mathrm{R}}\) are parameterized by \(\mathbf{U}_{\bullet}^{\mathrm{R}:\mathrm{N}}\in\mathrm{N}[\mathbf{U}_{\bullet }^{\mathrm{R}}]\), as in Eq. (4.7.2). Importantly, components in \(\mathbf{U}_{\bullet}^{\mathrm{N}:\mathrm{N}}\) and \(\mathbf{U}_{\bullet}^{\mathrm{R}:\mathrm{N}}\) are orthogonal to the latent space solution and, therefore, represent variations in the gas phase that do not affect the dynamics of the surface phase. The latent space solution is computed as

\[\mathbf{z}^{\mathrm{R}}=\left(\mathbf{U}_{\bullet}^{\mathrm{R}}\right)^{+} \left(\mathrm{C}_{\mathrm{N}}[\mathbf{\hat{x}}_{\bullet(t)}^{\mathrm{R}}]- \mathbf{x}_{\bullet}^{\mathrm{N}}\right)+\mathbf{U}_{\bullet}^{\mathrm{R}: \mathrm{N}}\mathbf{\hat{x}}_{\mathbf{o}\left(t\right)}^{\mathrm{R}} \tag{4.7.2}\]

where \(\mathbf{\hat{x}}_{\mathbf{o}\left(\mathbf{t}\right)}^{\mathrm{R}}\in\mathbb{R}^ {d}\) for \(d=\mathrm{nullity}[\mathbf{U}_{\bullet}^{\mathrm{R}}]\), and \(\mathbf{\hat{x}}_{\mathbf{o}\left(\mathbf{t}\right)}^{\mathrm{N}}\in\mathbb{R} ^{d}\) for \(d=\mathrm{nullity}[\mathbf{U}_{\bullet}^{\mathrm{N}}]\). With such a structure, all the variance is constrained to the \(\mathrm{R}[\mathbf{M}]\), which is parameterized by the time-invariant nullspace. We refer to the reduced optimization problem formulation as rKINN (Robust KINN), which is conveyed in Eq. (4.7.3), where \(\mathbf{\varepsilon}_{\mathbf{z}}^{\mathrm{R}}=(\mathbf{U}^{\mathrm{R}})^{T} \mathbf{\tilde{x}}-\mathbf{z}^{\mathrm{R}}=(\mathbf{U}^{\mathrm{R}})^{T} \mathbf{\varepsilon}_{\mathbf{x}}\) and \(\mathbf{\varepsilon}_{\mathbf{z}}^{\mathrm{R}}=(\mathbf{U}^{\mathrm{R}})^{T} \mathbf{\tilde{}}\mathbf{\varepsilon}_{\mathbf{x}}\) are the state- and state-derivative errors projected onto the \(\mathrm{R}[\mathbf{M}]\) subspace, respectively. This leads to a revised definition of the optimization problem,

\[\min_{\mathbf{\omega}_{t}} \ell_{\mathrm{t}}=\frac{1}{n}\sum_{i}^{n}\left((\mathbf{ \varepsilon}_{\mathbf{z}_{i}}^{\mathrm{R}})^{T}\mathbf{\Omega}_{\mathbf{ \delta}_{i}}^{\mathrm{R}}\mathbf{\varepsilon}_{\mathbf{\delta}_{i}}^{\mathrm{R }}+(\mathbf{\varepsilon}_{\mathbf{z}_{i}}^{\mathrm{R}})^{T}\mathbf{\Omega}_{ \mathbf{\delta}_{i}}^{\mathrm{R}}\mathbf{\varepsilon}_{\mathbf{\delta}_{i}}^{ \mathrm{R}}\right)\] (4.7.3) s.t. \[\mathbf{\Omega}_{\mathbf{z}_{i}}^{\mathrm{R}}=\left((\mathbf{U}^{ \mathrm{R}})^{T}\left(\mathbf{\Sigma}_{\mathbf{\hat{x}}_{i}}^{\mathrm{x}}+ \mathbf{\Sigma}_{\mathbf{\delta}_{i}}^{\mathrm{p}}\right)\mathbf{U}^{\mathrm{ R}}\right)^{-1}\] \[\mathbf{\Omega}_{\mathbf{z}}^{\mathrm{R}}=\left((\mathbf{U}^{ \mathrm{R}})^{T}\mathbf{\Sigma}_{\mathbf{x}}\mathbf{U}^{\mathrm{R}}\right)^{-1}\] \[\mathbf{p}\in\mathbb{R}^{\mathrm{dim}(\mathbf{p})},\ \mathbf{\omega}_{\mathbf{s}}\in\mathbb{R}^{ \mathrm{dim}(\mathbf{\omega}_{\mathbf{s}})}\]

which is solved with the rKINNs algorithm 1.

```
Input: Observed data (\(\mathbf{t}\), \(\mathbf{\tilde{x}}\)), stoichiometry matrix \(\mathbf{M}\), max. iterations \(n_{max}\), no. of epochs \(n_{e}\), tolerance \(tol\)
2: Build NN, initialize SA parameters \(\mathbf{\omega}_{\mathbf{s}}\), kinetic model parameters \(\mathbf{p}\) Evaluate \(\mathbf{\Omega}_{\mathbf{z}}^{\mathrm{R}}\) and \(\mathbf{\Omega}_{\mathbf{z}}^{\mathrm{R}}\).
4:while epoch \(<n_{e}\) or \(\ell_{\mathrm{t}}(\mathbf{\omega}_{t})\) < \(tol\)do whileiter < \(n_{max}\)do
6: Predict states \(\mathbf{x}(t,\mathbf{\omega}_{\mathbf{s}})\) from SA over \(\mathbf{t}\) Predict states derivatives \(\mathbf{\hat{x}}(t,\mathbf{\omega}_{\mathbf{s}})\) from SA over \(\mathbf{t}\)
8: Compute projected errors \(\mathbf{\varepsilon}_{\mathbf{z}}^{\mathrm{R}}\) and \(\mathbf{\varepsilon}_{\mathbf{z}}^{\mathrm{R}}\).  Compute objective function \(\ell_{\mathrm{t}}(\mathbf{\omega}_{t})\) over \(\mathbf{t}\)
10: Compute gradients of \(\ell_{\mathrm{t}}(\mathbf{\omega}_{t})\) with AD; \(\partial_{\mathbf{\omega}_{t}}\ell_{\mathrm{t}}(\mathbf{\omega}_{\mathbf{s}},\mathbf{p})\) Update SA parameters \(\mathbf{\omega}_{\mathbf{s}}\) and kinetic model parameters \(\mathbf{p}\) using Adam
12: Update \(\partial_{\mathbf{x}}f(\mathbf{x},\mathbf{p})\) and \(\partial_{\mathbf{p}}f(\mathbf{x},\mathbf{p})\) Update \(\mathbf{\Sigma}_{\mathbf{\hat{x}}_{i}}^{\mathrm{x}}\) and \(\mathbf{\Sigma}_{\mathbf{\hat{x}}_{i}}^{\mathrm{p}}\)
14: Update \(\mathbf{\Sigma}_{\mathbf{x}}\) and \(\mathbf{\Sigma}_{\mathbf{\hat{x}}}\) with uncorrelated additive noise \(|\langle\mathbf{\varepsilon}_{\mathbf{x}}\rangle|\) Update covariance matrices \(\mathbf{\Sigma}_{\mathbf{x}}^{\mathrm{R}}\) and \(\mathbf{\Sigma}_{\mathbf{\hat{x}}_{i}}^{\mathrm{R}}\)
16: Evaluate precision matrices \(\mathbf{\Omega}_{\mathbf{z}}^{\mathrm{R}}\) and \(\mathbf{\Omega}_{\mathbf{\hat{x}}_{i}}^{\mathrm{R}}\) Output: Kinetic model parameters \(\mathbf{p}\)
```

**Algorithm 1** Robust KINNs

## Acknowledgements

G.S.G. is a PhD research fellow with International Business Machine Corporation (IBM) since mid-2021 through the IBM Global University Program Awards. Acknowledgement is made to the donors of the American Chemical Society Petroleum Research Fund for partial support of this research. This work was supported by the U.S. Department of Energy (USDOE), Office of Energy Efficiency and Renewable Energy (EERE), Industrial Efficiency and Decarbonization
Office (IEDO) Next Generation R&D Project DE-FOA-0002252-1775 under contract no. DE-AC07-05ID14517. We acknowledge Dr. John Kitchin for his formative work on the utilization of neural networks as bases for the solution of simple coupled forward kinetics ODEs (kitchingroup.cheme.cmu.edu/blog). The authors are grateful to Dr. M. Ross Kunz and Dr. Adam Yonge for the enlightening discussions involving the maximum-likelihood formulation of PINNs problems and uncertainty quantification, to Ms. Sushree Jagriti Sahoo for her comments on the algebraic decomposition of ODE coupling matrices, and to Mr. Dingqi Nai for general assessments of the application of KINNs to experimental data. G.S.G. thanks especially Dr. William Bradley and Mr. Zachary Kilwein for the extensive debate-like discussions on the application of PINNs and NODEs to general Chemical Engineering problems.

Figure 4: Surrogate approximator structure based on range-nullspace decomposition using SVD. Initial-values and inner fluxes are contained in the time-invariant nullspace component (bottom). Time variations are represented by the fully-connected feed-forward neural network (top) with \(k\) hidden layers, which includes both unbound, \(\hat{\mathbf{x}}_{\bullet}^{\mathrm{R}}\), and bound, \(\hat{\mathbf{x}}_{\bullet}^{\mathrm{R}}\) latent representations that, through \(C_{\mathrm{N}}\) and linear operations, leads to \(\mathbf{z}^{\mathrm{R}}\). The MLE framework is constructed around the \(\mathbf{z}^{\mathrm{R}}\) lower dimension representation assuming that \(\mathbf{z}^{\mathrm{N}}\) can be obtained directly from measured data.
### Conflict of Interest

The authors declare the nonexistence of conflicts of interest.

### Nomenclature

**Ordinary Differential Equation**

\(\mathbf{M}\): Stoichiometry matrix

\(x\): State variable

\(\mathbf{x}\): State variable array

\(\mathbf{\dot{x}}\): State variable time derivatives

\(\mathbf{c}\): Concentrations vector

\(\mathbf{r}\): Rate of reaction vector

\(\mathbf{\Sigma}\): Covariance matrix

\(\mathbf{\Omega}\): Precision matrix

\(\mathbf{\tilde{y}}\): Unscaled measurements

\(\mathbf{\tilde{x}}\): Scaled measurements

\(\mathbf{\gamma}\): Calibration factor

**Surrogate Models**

\(\phi\): Activation function

\(\mathrm{C_{N}}\): Normalization-constraint operator

\(\ell\): Objective function errors: \(\ell_{\mathrm{t}},\ \ell_{\mathbf{x}},\ \ell_{\mathbf{\hat{x}}}\) data and model likelihoods, respectively

**Residuals**

\(\boldsymbol{\varepsilon}\): Residuals array

\(\boldsymbol{\varepsilon}_{i}\): Residuals array at time point \(t_{i}\)

\(\boldsymbol{\hat{\varepsilon}}\): Residuals random variable

\(\delta\): Sensitivity through total derivative

**Parameters**

\(\mathbf{p}\): Kinetic model parameters

\(\boldsymbol{\omega_{\mathrm{s}}}\): Surrogate model parameters

**Dimensionality Reduction**

\(\mathbf{U^{R}}\): Orthogonal basis for the range of \(\mathbf{M}\), i.e., \(\mathrm{row}[\mathbf{M}]\)

\(\mathbf{U^{N}}\): Orthogonal basis for the nullspace of \(\mathbf{M}\), i.e., \(\mathrm{null}[\mathbf{M}]\)

\(\mathbf{V^{R}}\): Orthogonal basis for the range of \(\mathbf{M}^{T}\), i.e., \(\mathrm{row}[\mathbf{M}^{T}]\)

\(\mathbf{V^{N}}\): Orthogonal basis for the nullspace of \(\mathbf{M}^{T}\), i.e., \(\mathrm{null}[\mathbf{M}^{T}]\)

\(\mathbf{S}\): Diagonal matrix with singular values of \(\mathbf{M}\)

**Abbreviations**

SVD: Singular-value Decomposition

## References

* [1] M. Raissi, P. Perdikaris, and G. E. Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. _Journal of Computational Physics_, 378:686-707, 2 2019.
* [2] Dehao Liu and Yan Wang. Multi-Fidelity Physics-Constrained Neural Network and Its Application in Materials Modeling. _Journal of Mechanical Design_, 141(12), 12 2019.
* [3] Shengze Cai, Zhicheng Wang, Sifan Wang, Paris Perdikaris, and George Em Karniadakis. Physics-informed neural networks for heat transfer problems. _Journal of Heat Transfer_, 143(6), 6 2021.
* [4] Siddhartha Mishra and Roberto Molinaro. Physics informed neural networks for simulating radiative transfer. _Journal of Quantitative Spectroscopy and Radiative Transfer_, 270:107705, 8 2021.
* [5] Ameya D. Jagtap, Ehsan Kharazmi, and George Em Karniadakis. Conservative physics-informed neural networks on discrete domains for conservation laws: Applications to forward and inverse problems. _Computer Methods in Applied Mechanics and Engineering_, 365:113028, 6 2020.
* [6] Shengze Cai, Zhiping Mao, Zhicheng Wang, Minglang Yin, and George Em Karniadakis. Physics-informed neural networks (PINNs) for fluid mechanics: a review. _Acta Mechanica Sinica/Lixue Xuebao_, 37(12):1727-1738, 12 2021.
* [7] Guofei Pang, Lu Lu, and George Em Karniadakis. fPINNs: Fractional Physics-Informed Neural Networks. _SIAM Journal on Scientific Computing_, 41(4):A2603-A2626, 1 2019.
* [8] Lu Lu, Raphael Pestourie, Wenjie Yao, Zhicheng Wang, Francesc Verdugo, and Steven G. Johnson. Physics-Informed Neural Networks with Hard Constraints for Inverse Design. _SIAM Journal on Scientific Computing_, 43(6):B1105-B1132, 1 2021.
* [9] Siddhartha Mishra and Roberto Molinaro. Estimates on the generalization error of physics-informed neural networks for approximating a class of inverse problems for PDEs. _IMA Journal of Numerical Analysis_, 42(2):981-1022, 4 2022.
* [10] Jeremy Yu, Lu Lu, Xuhui Meng, and George Em Karniadakis. Gradient-enhanced physics-informed neural networks for forward and inverse PDE problems. _Computer Methods in Applied Mechanics and Engineering_, 393:114823, 4 2022.
* [11] Xiaoli Chen, Liu Yang, Jinqiao Duan, and George Em Karniadakis. Solving Inverse Stochastic Problems from Discrete Particle Observations Using the Fokker-Planck Equation and Physics-Informed Neural Networks. _SIAM Journal on Scientific Computing_, 43(3):B811-B830, 1 2021.
* [12] Ravi G. Patel, Indu Manickam, Nathaniel A. Trask, Mitchell A. Wood, Myoungkyu Lee, Ignacio Tomas, and Eric C. Cyr. Thermodynamically consistent physics-informed neural networks for hyperbolic systems. _Journal of Computational Physics_, 449:110754, 1 2022.
* [13] Xiaocong Du, Shreyas Kolala Venkataramanaiah, Zheng Li, Chih-Wei Chang, Yuan Gao, Tonghe Wang, Suryanarayana Maddu, Dominik Sturm, Christian L Muller, and Ivo F Sbalzarini. Inverse Dirichlet weighting enables reliable training of physics informed neural networks. _Machine Learning: Science and Technology_, 3(1):015026, 2 2022.
* [14] Fumio Hayashi. Finite-Sample Properties of OLS. In _Econometrics_, pages 47-59. 2000.
* [15] Han Gao, Matthew J. Zahr, and Jian-Xun Wang. Physics-informed graph neural Galerkin networks: A unified framework for solving PDE-governed forward and inverse problems. _Computer Methods in Applied Mechanics and Engineering_, 390:114502, 2 2022.
* [16] Xuhui Meng, Zhen Li, Dongkun Zhang, and George Em Karniadakis. PPINN: Parareal physics-informed neural network for time-dependent PDEs. _Computer Methods in Applied Mechanics and Engineering_, 370, 10 2020.
* [17] Liu Yang, Xuhui Meng, and George Em Karniadakis. B-PINNs: Bayesian physics-informed neural networks for forward and inverse PDE problems with noisy data. _Journal of Computational Physics_, 425:109913, 1 2021.
* [18] Kirill Zubov, Zoe McCarthy, Yingbo Ma, Francesco Calisto, Valerio Pagliarino, Simone Azeglio, Luca Bottero, Emmanuel Lujan, Valentin Sulzer, Ashutosh Bharambe, Nand Vinchhi, Kaushik Balakrishnan, Devesh Upadhyay, and Chris Rackauckas. NeuralPDE: Automating Physics-Informed Neural Networks (PINNs) with Error Approximations. 7 2021.
* [19] Son Ich Ngo and Young-Il Lim. Solution and Parameter Identification of a Fixed-Bed Reactor Model for Catalytic CO2 Methanation Using Physics-Informed Neural Networks. _Catalysts 2021, Vol. 11, Page 1304_, 11(11):1304, 10 2021.
* [20] Alireza Yazdanii, Lu Luid, Maziar Raissid, and George Em Karniadakisid. Systems biology informed deep learning for inferring parameters and hidden dynamics. 2020.
* [21] Jiequn Han, Arnulf Jentzen, and E. Weinan. Solving high-dimensional partial differential equations using deep learning. _Proceedings of the National Academy of Sciences of the United States of America_, 115(34):8505-8510, 8 2018.
* [22] J.W. Thybaut and G.B. Marin. Single-Event MicroKinetics: Catalyst design for complex reaction networks. _Journal of Catalysis_, 308:352-362, 12 2013.
* [23] Rob J. Berger, Freek Kapteijn, Jacob A. Moulijn, Guy B. Marin, Juray De Wilde, Maria Olea, De Chen, Anders Holmen, Luca Lietti, Enrico Tronconi, and Yves Schuurman. Dynamic methods for catalytic kinetics. _Applied Catalysis A: General_, 342(1-2):3-28, 6 2008.
* [24] Patricia Rubert-Nason, Manos Mavrikakis, Christos T. Maravelias, Lars C. Grabow, and Lorenz T. Biegler. Advanced solution methods for microkinetic models of catalytic reactions: A methanol synthesis case study. _AIChE Journal_, 60(4):1336-1346, 4 2014.
* [25] Srinivas Rangarajan, Christos T. Maravelias, and Manos Mavrikakis. Sequential-Optimization-Based Framework for Robust Modeling and Design of Heterogeneous Catalytic Systems. _Journal of Physical Chemistry C_, 121(46):25847-25863, 11 2017.
* [26] Guy B. Marin, Vladimir V. Galvita, and Gregory S. Yablonsky. Kinetics of chemical processes: From molecular to industrial scale. _Journal of Catalysis_, 404:745-759, 12 2021.
* [27] Zachary T Wilson and Nikolaos V Sahinidis. The ALAMO approach to machine learning. _Computers and Chemical Engineering_, 106:785-795, 2017.
* [28] Steven L. Brunton, Joshua L. Proctor, J. Nathan Kutz, and William Bialek. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. _Proceedings of the National Academy of Sciences of the United States of America_, 113(15):3932-3937, 4 2016.
* Methodological Frameworks. 1 2023.
* [30] J. G. Verwer. Gauss-Seidel Iteration for Stiff ODES from Chemical Kinetics. _[https://doi.org/10.1137/0915076_](https://doi.org/10.1137/0915076_), 15(5):1243-1250, 7 2006.
* [31] Qin Wu, Talin Avanesian, Xiaohui Qu, and Hubertus Van Dam. PolyODENet: Deriving mass-action rate equations from incomplete transient kinetics data. _The Journal of Chemical Physics_, 157(16):164801, 10 2022.
* [32] Gabriel S. Gusmao, Adhika P. Retnanto, Shashwati C. da Cunha, and Andrew J. Medford. Kinetics-informed neural networks. _Catalysis Today_, 4 2022.
* [33] L T Biegler, J J Damiano, and G E Blau. Nonlinear parameter estimation: A case study comparison. _AIChE Journal_, 32(1):29-45, 1986.
* [34] A. J. Wathen. Preconditioning. _Acta Numerica_, 24:329-376, 5 2015.
* [35] J. A. Dumesic. The Microkinetics of heterogeneous catalysis. _ACS professional reference book_, 1993.
* [36] D. Constales, G. S. Yablonsky, and G. B. Marin. The C-matrix: Augmentation and reduction in the analysis of chemical composition and structure. _Chemical Engineering Science_, 110:164-168, 2014.
* [37] Evgeniy A Redekop, Gregory S Yablonsky, Denis Constales, Palghat A Ramachandran, John T Gleaves, and Guy B Marin. Elucidating complex catalytic mechanisms based on transient pulse-response kinetic data. _Chemical Engineering Science_, 110:20-30, 5 2014.
* [38] George Em Karniadakis, Ioannis G. Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang. Physics-informed machine learning. _Nature Reviews Physics_, 3(6):422-440, 6 2021.
* [39] A.J. Meade and A.A. Fernandez. The numerical solution of linear ordinary differential equations by feedforward neural networks. _Mathematical and Computer Modelling_, 19(12):1-25, 6 1994.
* [40] A.J. Meade and A.A. Fernandez. Solution of nonlinear ordinary differential equations by feedforward neural networks. _Mathematical and Computer Modelling_, 20(9):19-44, 11 1994.
* [41] I.E. E. Lagaris, Aristidis Likas, and D.I. I. Fotiadis. Artificial neural networks for solving ordinary and partial differential equations. _IEEE Transactions on Neural Networks_, 9(5):987-1000, 5 1998.
* [42] G. Cybenko. Approximation by superpositions of a sigmoidal function. _Mathematics of Control, Signals, and Systems_, 2(4):303-314, 12 1989.
* [43] Kurt Hornik, Maxwell Stinchcombe, and Halbert White. Universal approximation of an unknown mapping and its derivatives using multilayer feedforward networks. _Neural Networks_, 3(5):551-560, 1 1990.
* [44] Kurt Hornik. Approximation capabilities of multilayer feedforward networks. _Neural Networks_, 4(2):251-257, 1991.
* [45] Moshe Leshno, Vladimir Ya Lin, Allan Pinkus, and Shimon Schocken. Multilayer feedforward networks with a nonpolynomial activation function can approximate any function. _Neural Networks_, 6(6):861-867, 1 1993.
* [46] Allan Pinkus. Approximation theory of the MLP model in neural networks. pages 143-195, 1999.
* [47] Nam Mai-Duy and Thanh Tran-Cong. Approximation of function and its derivatives using radial basis function networks. _Applied Mathematical Modelling_, 27(3):197-220, 3 2003.
* [48] Philip W Livermore. Galerkin orthogonal polynomials. 2009.
* [49] Philip W Livermore, Glenn R Ierley, P W Livermore, and \(\cdot\) G R Ierley. Quasi-L p norm orthogonal Galerkin expansions in sums of Jacobi polynomials Orthogonal expansions. 54:533-569, 2010.
* [50] Anders Hald. On the history of maximum likelihood in relation to inverse probability and least squares. _Statistical Science_, 14(2):214-222, 5 1999.
* [51] In Jae Myung. Tutorial on maximum likelihood estimation. _Journal of Mathematical Psychology_, 47(1):90-100, 2 2003.
* [52] F. H. MacDougall. Thermodynamic Theory of Affinity. By Th. De Donder and Pierre Van Rysselberghe. _The Journal of Physical Chemistry_, 41(5), 1937.
* [53] Evgeniy A. Redekop, Gregory S. Yablonsky, Denis Constales, Palghat A. Ramachandran, Cathryn Pherigo, and John T. Gleaves. The Y-Procedure methodology for the interpretation of transient kinetic data: Analysis of irreversible adsorption. _Chemical Engineering Science_, 66(24):6441-6452, 12 2011.
* [54] Weiqi Ji and Sili Deng. Autonomous Discovery of Unknown Reaction Pathways from Data by Chemical Reaction Neural Network. _The Journal of Physical Chemistry A_, 125(4):1082-1092, 2 2021.
* [55] Weiqi Ji, Franz Richter, Michael J. Gollner, and Sili Deng. Autonomous kinetic modeling of biomass pyrolysis using chemical reaction neural networks. _Combustion and Flame_, 240:111992, 6 2022.
* [56] Christopher Rackauckas, Yingbo Ma, Julius Martensen, Collin Warner, Kirill Zubov, Rohit Supekar, Dominic Skinner, Ali Ramadhan, and Alan Edelman. Universal Differential Equations for Scientific Machine Learning. 1 2020.
* [57] William Bradley, Gabriel S. Gusmao, Andrew J. Medford, and Fani Boukouvala. Training Stiff Dynamic Process Models via Neural Differential Equations. volume 49, pages 1741-1746. Elsevier, 1 2022.
* [58] Jun Yin, Jiali Li, Iftekhar A Karimi, and Xiaonan Wang. Generalized reactor neural ODE for dynamic reaction process modeling with physical interpretability. _Chemical Engineering Journal_, 452:139487, 2023.
* [59] Felix A. Doppel and Martin Votsmeier. Efficient machine learning based surrogate models for surface kinetics by approximating the rates of the rate-determining steps. _Chemical Engineering Science_, 262:117964, 11 2022.
* [60] Murali R. Rajamani and James B. Rawlings. Estimation of the disturbance structure from data using semidefinite programming and optimal weighting. _Automatica_, 45(1):142-148, 1 2009.
* [61] Raf Roelant, Denis Constales, Gregory S. Yablonsky, Roger Van Keer, Michael A. Rude, and Guy B. Marin. Noise in temporal analysis of products (TAP) pulse responses. _Catalysis Today_, 121(3-4):269-281, 3 2007.
* [62] Peter J. Rousseeuw and Katrien Van Driessen. A Fast Algorithm for the Minimum Covariance Determinant Estimator. _Technometrics_, 41(3):212-223, 8 1999.
* [63] Olivier Ledoit and Michael Wolf. A well-conditioned estimator for large-dimensional covariance matrices. _Journal of Multivariate Analysis_, 88(2):365-411, 2 2004.
* [64] Yilun Chen, Ami Wiesel, Yonina C. Eldar, and Alfred O. Hero. Shrinkage Algorithms for MMSE Covariance Estimation. _IEEE Transactions on Signal Processing_, 58(10):5016-5029, 10 2010.
* [65] Sivaram Ambikasaran, Daniel Foreman-Mackey, Leslie Greengard, David W. Hogg, and Michael O'Neil. Fast Direct Methods for Gaussian Processes. _IEEE Transactions on Pattern Analysis and Machine Intelligence_, 38(2):252-265, 2 2016.
* [66] Roy Frostig, Matthew Johnson, and Chris Leary. Compiling machine learning programs via high-level tracing. In _SysML_, 3 2018.
* [67] James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, and Skye Wanderman-Milne. {JAX}: composable transformations of {P}ython+{N}um{P}y} programs. _[http://github.com/googlefjax_](http://github.com/googlefjax_), 2018.
* Conference Track Proceedings_. International Conference on Learning Representations, ICLR, 12 2015.
* [69] Alan C. Hindmarsh. LSODE and LSODI, two new initial value ordinary differential equation solvers. _ACM SIGNUM Newsletter_, 15(4):10-11, 12 1980.
* [70] Linda Petzold. Automatic Selection of Methods for Solving Stiff and Nonstiff Systems of Ordinary Differential Equations. _SIAM Journal on Scientific and Statistical Computing_, 4(1):136-148, 3 1983.
* [71] Robert J. Madix. The application of flash desorption spectroscopy to chemical reactions on surfaces: Temperature programmed reaction spectroscopy. _Critical Reviews in Solid State and Materials Sciences_, 7(2):143-152, 1978.
Maximum-likelihood Estimators in Physics-Informed Neural Networks for High-dimensional Inverse Problems

Gabriel S. Gusmao

School of Chemical & Biomolecular Engineering

Georgia Institute of Technology

Atlanta, GA 30332

{gusmaogabriels, ajm}@gatech.edu

Andrew J. Medford

School of Chemical & Biomolecular Engineering

Georgia Institute of Technology

Atlanta, GA 30332

{gusmaogabriels, ajm}@gatech.edu

###### Abstract

We propose a
\[x_{i}(\mathbf{\omega_{s}},t) =(1-\sigma\left(\hat{x}_{i}(\mathbf{\omega_{s}},t)\right))\prod_{j<i} \sigma\left(\hat{x}_{j}(\mathbf{\omega_{s}},t)\right)\ \forall\ i<p;\ i,j\in\mathbb{N}\] (S1.2) \[x_{p}(\mathbf{\omega_{s}},t) =\prod_{j<p}\sigma\left(\hat{x}_{j}(\mathbf{\omega_{s}},t)\right)\]

## S2 MLE Optimality Conditions

The optimality condition for cost function of the minimization problem in eq. (2.2.6) can be represented as expected values in eq. (S2.1).

\[\nabla\left\langle\mathbf{\varepsilon}_{\mathbf{x}}^{T}\mathbf{\Omega}_{\mathbf{k}} \mathbf{\varepsilon}_{\mathbf{\hat{x}}}+\mathbf{\varepsilon}_{\mathbf{x}}^{T}\mathbf{ \Omega}_{\mathbf{x}}\mathbf{\varepsilon}_{\mathbf{x}}\right\rangle=\mathbf{0}\] (S2.1)

Which can be expanded over the independent variables \(\mathbf{x}\) and \(\mathbf{p}\) as in

\[\nabla\left\langle\mathbf{\varepsilon}_{\mathbf{x}}^{T}\mathbf{\Omega}_{ \mathbf{k}}\mathbf{\varepsilon}_{\mathbf{\hat{x}}}+\mathbf{\varepsilon}_{\mathbf{x}}^ {T}\mathbf{\Omega}_{\mathbf{x}}\mathbf{\varepsilon}_{\mathbf{x}}\right\rangle =\left[\partial_{\mathbf{x}},\partial_{\mathbf{x}},\partial_{ \mathbf{p}}\right]\left\langle\mathbf{\varepsilon}_{\mathbf{x}}^{T}\mathbf{\Omega}_{ \mathbf{k}}\mathbf{\varepsilon}_{\mathbf{\hat{x}}}+\mathbf{\varepsilon}_{\mathbf{x}}^ {T}\mathbf{\Omega}_{\mathbf{x}}\mathbf{\varepsilon}_{\mathbf{x}}\right\rangle\] (S2.2) \[\propto\left\langle\left[\partial_{\mathbf{x}}\mathbf{\hat{x}}+ \mathbf{I}-\partial_{\mathbf{x}}f(\mathbf{x},\mathbf{p}),\right.\right.- \partial_{\mathbf{p}}f(\mathbf{x},\mathbf{p})]^{T}\,\mathbf{\Omega}_{\mathbf{k}} \mathbf{\varepsilon}_{\mathbf{\hat{x}}}\right\rangle-\left\langle\left[\partial_ {\mathbf{\hat{x}}}\mathbf{x},\ \mathbf{I}\right]^{T}\mathbf{\Omega}_{\mathbf{x}}\mathbf{\varepsilon}_{ \mathbf{x}}\right\rangle\] (S2.3) \[=-\left\langle\begin{bmatrix}\partial_{\mathbf{x}}f(\mathbf{x}, \mathbf{p})^{T}-(\partial_{\mathbf{x}}\mathbf{\hat{x}})^{T}&\mathbf{I}\\ \partial_{\mathbf{p}}f(\mathbf{x},\mathbf{p})^{T}&\mathbf{0}\\ -\mathbf{I}&(\partial_{\mathbf{x}}\mathbf{x})^{T}\end{bmatrix}\begin{bmatrix} \mathbf{\Omega}_{\mathbf{\hat{x}}}&\mathbf{0}\\ \mathbf{0}&\mathbf{\Omega}_{\mathbf{x}}\end{bmatrix}\begin{bmatrix}\mathbf{\varepsilon}_{ \mathbf{\hat{x}}}\\ \mathbf{\varepsilon}_{\mathbf{x}}\end{bmatrix}\right\rangle=\mathbf{0}\] (S2.4)

With the assumption that \(\partial_{\mathbf{x}}\mathbf{\hat{x}}\) and \(\partial_{\mathbf{\hat{x}}}\mathbf{x}=(\partial_{\mathbf{x}}\mathbf{\hat{x}}) ^{-1}\) do not contribute to structural variance, the condition defined in eq. (S2.4) can be simplified to

\[\left\langle\begin{bmatrix}\partial_{\mathbf{x}}f(\mathbf{x}, \mathbf{p})^{T}&\mathbf{I}\\ \partial_{\mathbf{p}}f(\mathbf{x},\mathbf{p})^{T}&\mathbf{0}\end{bmatrix} \begin{bmatrix}\mathbf{\Omega}_{\mathbf{\hat{x}}}&\mathbf{0}\\ \mathbf{0}&\mathbf{\Omega}_{\mathbf{x}}\end{bmatrix}\begin{bmatrix}\mathbf{\varepsilon }_{\mathbf{\hat{x}}}\\ \mathbf{\varepsilon}_{\mathbf{x}}\end{bmatrix}\right\rangle=\mathbf{0}\] (S2.5)

From the final results in eq. (4.1.2), eq. (S2.5) can be expressed as a function about \(\mathbf{\Omega}_{\mathbf{x}}\) and \(\mathbf{\Omega}_{\mathbf{p}}\) only, but expressing \(\partial_{\mathbf{x}}f(\mathbf{x}_{i},\mathbf{p})\) and \(\partial_{\mathbf{p}}f(\mathbf{x}_{i},\mathbf{p})\) in terms of their singular value decompositions.

\[\partial_{\mathbf{x}}f(\mathbf{x},\mathbf{p}) =\mathbf{U}_{\mathbf{x}}\mathbf{S}_{\mathbf{x}}\mathbf{V}_{\mathbf{ x}}^{T}\] (S2.6) \[\partial_{\mathbf{p}}f(\mathbf{x},\mathbf{p}) =\mathbf{U}_{\mathbf{p}}\mathbf{S}_{\mathbf{p}}\mathbf{V}_{\mathbf{ p}}^{T}\] \[\text{with }\mathbf{U}^{T}\mathbf{U} =\mathbf{I},\ \mathbf{V}^{T}\mathbf{V}=\mathbf{I}\] \[\text{and }\mathbf{S} =\mathrm{diag}(\mathbf{s}),\ \text{singular values}\]

Which leads to eq. (S2.7) for the general case where \(\dim(\mathbf{p})>\dim(\mathbf{x})\).

\[\left\langle\begin{bmatrix}\mathbf{V}_{\mathbf{x}}\mathbf{S}_{ \mathbf{x}}\mathbf{U}_{\mathbf{x}}^{T}\mathbf{U}_{\mathbf{x}}\mathbf{S}_{ \mathbf{x}}^{-1}\mathbf{V}_{\mathbf{x}}^{T}\mathbf{\Omega}_{\mathbf{x}}\mathbf{V}_{ \mathbf{x}}\mathbf{S}_{\mathbf{x}}^{-1}\mathbf{U}_{\mathbf{x}}^{T}&\mathbf{\Omega}_ {\mathbf{x}}\\ \mathbf{V}_{\mathbf{p}}\mathbf{S}_{\mathbf{p}}\mathbf{U}_{\mathbf{p}}^{T} \mathbf{U}_{\mathbf{p}}\mathbf{S}_{\mathbf{p}}^{-1}\mathbf{V}_{\mathbf{p}}^{T} \mathbf{\Omega}_{\mathbf{p}}\mathbf{V}_{\mathbf{p}}\mathbf{S}_{\mathbf{p}}^{-1} \mathbf{U}_{\mathbf{p}}^{T}&\mathbf{0}\end{bmatrix}\begin{bmatrix}\mathbf{\varepsilon }_{\mathbf{x}}\\ \mathbf{\varepsilon}_{\mathbf{x}}\end{bmatrix}\right\rangle =\left\langle\begin{bmatrix}\mathbf{\Omega}_{\mathbf{x}}\left( \left(\partial_{\mathbf{x}}f(\mathbf{x},\mathbf{p})^{T}\right)^{+}\mathbf{ \varepsilon}_{\mathbf{\hat{x}}}+\mathbf{\varepsilon}_{\mathbf{x}}\right)\\ \left\langle\mathbf{V}_{\mathbf{p}}\mathbf{V}_{\mathbf{p}}^{T}\mathbf{\Omega}_{ \mathbf{p}}\left(\partial_{\mathbf{p}}f(\mathbf{x},\mathbf{p})^{T}\right)^{+} \mathbf{\varepsilon}_{\mathbf{\hat{x}}}\right\rangle\end{bmatrix}\right\rangle\] (S2.7) \[=\begin{bmatrix}\mathbf{\Omega}_{\mathbf{x}}\left(\left(\partial_{ \mathbf{x}}f(\mathbf{x},\mathbf{p})\right)^{T}\right)^{+}\mathbf{\varepsilon}_{ \mathbf{\hat{x}}}+\mathbf{\varepsilon}_{\mathbf{x}}\right)\\ \left\langle\mathbf{V}_{\mathbf{p}}\mathbf{V}_{\mathbf{p}}^{T}\mathbf{\Omega}_{ \mathbf{p}}\left(\partial_{\mathbf{p}}f(\mathbf{x},\mathbf{p})^{T}\right)^{+} \mathbf{\varepsilon}_{\mathbf{\hat{x}}}\right\rangle\end{bmatrix}=\mathbf{0}\]

In which \(\mathbf{A}^{+}\) is the Moore-Penrose pseudo-inverse of matrix \(\mathbf{A}\)[2, 3]. For \(\mathbf{\Omega}_{\mathbf{x}}\) full-rank, the optimality conditions simplify to:

\[\left\langle\mathbf{\varepsilon}_{\mathbf{x}}\right\rangle =\mathbf{0}\] (S2.8) \[\left\langle\left(\partial_{\mathbf{x}}f(\mathbf{x},\mathbf{p})^{T} \right)^{+}\mathbf{\varepsilon}_{\mathbf{\hat{x}}}\right\rangle =\mathbf{0}\] (S2.9) \[\left\langle\mathbf{V}_{\mathbf{p}}\mathbf{V}_{\mathbf{p}}^{T}\mathbf{ \Omega}_{\mathbf{p}}\left(\partial_{\mathbf{p}}f(\mathbf{x},\mathbf{p})^{T} \right)^{+}\mathbf{\varepsilon}_{\mathbf{\hat{x}}}\right\rangle =\mathbf{0}\] (S2.10)

Condition eq. (S2.8) follows naturally from the assumption of centered interpolation errors. Condition eq. (S2.9) conveys the expected value for orthogonality between derivative residuals and the range of \(\partial_{\mathbf{x}}f\), and condition eq. (S2.10) denotes the expected value for the relationship between derivative residuals and gradients of the physical model with respect to the incumbent model parameter. Importantly, \(\mathbf{V}_{\mathbf{p}}\) is also a function of \(\mathbf{x}\) and \(\mathbf{p}\) since it comprises an orthogonal basis for the row space of \(\partial_{\mathbf{p}}f(\mathbf{x},\mathbf{p})\).
## S3 Latent state reconstruction

With a mean-field approximation, and assuming the spectroscopic data includes the dynamics of all kinetically relevant reaction intermediates in an NPT-system, calibration coefficients can be retrieved by minimizing the raw signal projection onto the coupling matrix nullspace. With the nomenclature in eq. (2.1.2), i.e., taking \(\tilde{\mathbf{x}}\) as observed values, the underlying material balance can be expressed in terms of the nullspace equations as in eq. (S3.1).

\[(\mathbf{U}_{\mathbf{o}}^{\mathrm{N}})^{T}\left(\tilde{\mathbf{x}}_{\mathbf{o} (t_{i})}-\tilde{\mathbf{x}}_{\mathbf{o}(t_{j})}\right)+(\mathbf{U}_{\star}^{ \mathrm{N}})^{T}\left(\tilde{\mathbf{x}}_{\star(t_{i})}-\tilde{\mathbf{x}}_{ \star(t_{j})}\right)=\varepsilon_{\mathbf{x}_{ij}}^{\mathrm{N}}\] (S3.1)

In the absence of statistical noise, the nullspace error, \(\varepsilon_{\mathbf{x}_{i}}^{\mathrm{N}},\,\forall i\), would be zero. Latent states can be represented in terms of raw signal, \(\mathbf{y}\in\mathbb{R}_{+}^{n}\), and a calibration factor per chemical intermediate, \(\boldsymbol{\gamma}\in\mathbb{R}_{+}^{\mathrm{dim}(\star)}\), as in eq. (S3.2).

\[\mathbf{x}_{\star}=\mathbf{y}_{\star}\circ\boldsymbol{\gamma}\] (S3.2)

Coverage concentrations also need to be normalized, which can be expressed in terms of a parameterized set of solutions for \(\boldsymbol{\gamma}\), as in the following minimization problem, eq. (S3.3).

\[\boldsymbol{\gamma}(\boldsymbol{\beta})=\boldsymbol{\gamma}^{\mathrm{R}}+ \boldsymbol{\gamma}^{\mathrm{N}}(\boldsymbol{\beta})=\left(\sum_{i}^{n}\tilde{ \mathbf{y}}_{\star(t_{i})}\tilde{\mathbf{y}}_{\star(t_{i})}^{T}\right)^{+} \sum_{i}^{n}\tilde{\mathbf{y}}_{\star(t_{i})}+\mathbf{U}_{\boldsymbol{\gamma} }^{\mathrm{N}}\boldsymbol{\beta}\] (S3.3)

Where \(\mathbf{U}_{\boldsymbol{\gamma}}^{\mathrm{N}}\) is a basis for the nullspace of \(\sum_{i}^{n}\tilde{\mathbf{y}}_{\star(t_{i})}\tilde{\mathbf{y}}_{\star(t_{i})} ^{T}\), which can be computed from its eigendecomposition. By replacing eq. (S3.2) with the partial solution of eq. (S3.3) in eq. (S3.1), the reduced form in eq. (S3.4) is obtained. For realistic scenarios, statistical noise will lead to non-zero eigenvalues, which are associated with the nullspace of \(\sum_{i}^{n}\tilde{\mathbf{y}}_{\star(t_{i})}\tilde{\mathbf{y}}_{\star(t_{i})} ^{T}\) and, therefore, a cutoff criterion must be established to obtain the nullspace basis by selecting the associated eigenvectors.

\[\begin{split}&(\mathbf{U}_{\mathbf{o}}^{\mathrm{N}})^{T}\left( \tilde{\mathbf{x}}_{\mathbf{o}(t_{i})}-\tilde{\mathbf{x}}_{\mathbf{o}(t_{j})} \right)+(\mathbf{U}_{\star}^{\mathrm{N}})^{T}\operatorname{diag}\left(\tilde {\mathbf{y}}_{\star(t_{i})}-\tilde{\mathbf{y}}_{\star(t_{j})}\right) \boldsymbol{\gamma}(\boldsymbol{\beta})\\ &=\mathbf{v}_{\mathbf{o}_{ij}}+\mathbf{V}_{\mathbf{\ast}_{ij}} \boldsymbol{\gamma}(\boldsymbol{\beta})\end{split}\] (S3.4)

The minimization of the MSE of eq. (S3.4) in terms of \(\boldsymbol{\beta}\), considering only the asymmetrical pairs \((i,j)\,:\,j>i\), leads to the closed-form solution in eq. (S3.5):

\[\boldsymbol{\gamma}=-\mathbf{U}_{\boldsymbol{\gamma}}^{\mathrm{N}}\left(( \mathbf{U}_{\boldsymbol{\gamma}}^{\mathrm{N}})^{T}\left(\sum_{i}^{n}\sum_{j}^{ i-1}\mathbf{V}_{\mathbf{\ast}_{ij}}^{T}\mathbf{V}_{\mathbf{\ast}_{ij}}\right) \mathbf{U}_{\boldsymbol{\gamma}}^{\mathrm{N}}\right)^{+}(\mathbf{U}_{ \boldsymbol{\gamma}}^{\mathrm{N}})^{T}\sum_{i}^{n}\sum_{j}^{i-1}\mathbf{V}_{ \mathbf{\ast}_{ij}}^{T}\left(\mathbf{v}_{\mathbf{o}_{ij}}+\mathbf{V}_{\mathbf{ \ast}_{ij}}\boldsymbol{\gamma}^{\mathrm{R}}\right)+\boldsymbol{\gamma}^{ \mathrm{R}}\] (S3.5)

With \(\boldsymbol{\beta}\) from eq. (S3.5) replaced in eq. (S3.3), multiplicative or calibration factors \(\boldsymbol{\gamma}\) can be found to satisfy the chemical system nullspace condition in eq. (2.3.5) under coverage normalization.

## S4 Model Parameters Conditional Covariance - Calibration Factor

The calibration factors, \(\boldsymbol{\gamma}\), can be included in the rKINNs objective function by expressing them error terms for interpolation and model in the original space, eqs. (2.1.2) and (2.1.3). We need first note that for reconstructed/calibrated data, the SAs have predefined nullspace invariances, \(\mathbf{z}^{\mathrm{N}}\), which is estimated directly from the unscaled data, \(\tilde{\mathbf{y}}\), and optimal calibration factors, \(\boldsymbol{\gamma}^{\star}\), eq. (S3.5), such that,

\[\begin{split}\mathbf{x}(t_{i},\boldsymbol{\omega}_{\mathbf{s}}, \boldsymbol{\gamma})&=\mathbf{U}^{\mathrm{R}}\mathbf{z}^{ \mathrm{R}}(t_{i},\boldsymbol{\omega}_{\mathbf{s}})+\mathbf{U}^{\mathrm{N}} \mathbf{z}^{\mathrm{N}}(\boldsymbol{\gamma})\\ &=\mathbf{U}^{\mathrm{R}}\mathbf{z}^{\mathrm{R}}(t_{i}, \boldsymbol{\omega}_{\mathbf{s}})+\mathbf{U}^{\mathrm{N}}\mathbf{z}^{\mathrm{N}} \\ &=\mathbf{U}^{\mathrm{R}}\mathbf{z}^{\mathrm{R}}(t_{i}, \boldsymbol{\omega}_{\mathbf{s}})+\mathbf{U}^{\mathrm{N}}\frac{1}{n}\sum_{j=1} ^{n}\left((\mathbf{U}^{\mathrm{N}})^{T}\tilde{\mathbf{x}}_{j}\odot \boldsymbol{\gamma}\right)\\ &=\mathbf{U}^{\mathrm{R}}\mathbf{z}^{\mathrm{R}}(t_{i}, \boldsymbol{\omega}_{\mathbf{s}})+\frac{1}{n}\sum_{j=1}^{n}\left(\mathbf{U}^ {\mathrm{N}}(\mathbf{U}^{\mathrm{N}})^{T}\operatorname{diag}(\tilde{\mathbf{x}}_ {j})\right)\boldsymbol{\gamma}\end{split}\] (S4.1)

Such a formulation conveys the effect of \(\boldsymbol{\gamma}\) in the nullspace solution, i.e., linear shifts in the data and, consequently, on the SAs. However, to allow the propagation of errors to the range equations, the error expression need to be formulated based on the unscaled mapping \(\mathbf{y}=\mathbf{x}\odot(\boldsymbol{\gamma}^{\star})^{-1}\).

\[\begin{split}\boldsymbol{\varepsilon}_{\mathbf{x}_{i}}(\boldsymbol{ \omega}_{\mathbf{s}},\boldsymbol{\gamma})&=\mathbf{x}(t_{i}, \boldsymbol{\omega}_{\mathbf{s}},\boldsymbol{\gamma})-\tilde{\mathbf{y}}_{i} \odot\boldsymbol{\gamma}\\ \therefore\,\,\boldsymbol{\varepsilon}_{\mathbf{y}_{i}}=\boldsymbol{ \varepsilon}_{\mathbf{x}_{i}}(\boldsymbol{\omega}_{\mathbf{s}},\boldsymbol{ \gamma})\odot(\boldsymbol{\gamma}^{\star})^{-1}&=\mathbf{x}(t_{i}, \boldsymbol{\omega}_{\mathbf{s}},\boldsymbol{\gamma})\odot(\boldsymbol{ \gamma}^{\star})^{-1}-\tilde{\mathbf{y}}_{i}\end{split}\] (S4.2)
Let \(\mathbf{y}(t_{i},\mathbf{\omega}_{\mathbf{s}},\mathbf{\gamma})=\mathbf{x}(t_{i},\mathbf{\omega }_{\mathbf{s}},\mathbf{\gamma})\odot(\mathbf{\gamma}^{*})^{-1}\), \(\mathbf{\varepsilon}_{\mathbf{x}_{i}}\) can be expressed in terms of \(\mathbf{y}(t_{i},\mathbf{\omega}_{\mathbf{s}},\mathbf{\gamma})\) as follows,

\[\begin{split}\mathbf{\varepsilon}_{\mathbf{x}_{i}}&= \mathbf{y}(t_{i},\mathbf{\omega}_{\mathbf{s}},\mathbf{\gamma})\odot\mathbf{\gamma}-\tilde{ \mathbf{y}}_{i}\\ \mathbf{\varepsilon}_{\mathbf{\dot{x}_{i}}}&=\dot{ \mathbf{x}}(t_{i},\mathbf{\omega}_{\mathbf{s}},\mathbf{\gamma})-\mathbf{M}\left( \mathbf{k}(\theta,\mathbf{p})\circ\psi(\mathbf{x}(t_{i},\mathbf{\omega}_{\mathbf{s }}))\right)\\ &=\dot{\mathbf{y}}(t_{i},\mathbf{\omega}_{\mathbf{s}},\mathbf{\gamma}) \odot\mathbf{\gamma}-\mathbf{M}\left(\mathbf{k}(\theta,\mathbf{p})\circ\psi( \dot{\mathbf{y}}(t_{i},\mathbf{\omega}_{\mathbf{s}},\mathbf{\gamma})\odot\mathbf{\gamma })\right)\end{split}\] (S4.3)

With \(\mathbf{\varepsilon}_{\mathbf{x}_{i}}^{\text{R/N}}=(\mathbf{U}^{\text{R/N}})^{T} \mathbf{\varepsilon}_{\mathbf{x}_{i}}\) and \(\mathbf{\varepsilon}_{\mathbf{z}_{i}}^{\text{R/N}}=(\mathbf{U}^{\text{R/N}})^{T} \mathbf{\varepsilon}_{\mathbf{\dot{x}_{i}}}\), the negative log-likelihood after range-nullspace decomposition can be evaluated as in eq. (4.7.3), and the Fisher information matrix, \(\mathcal{I}\), which is equivalent to the Hessian matrix of the MLE, \(\ell_{\text{t}}\), can be estimated as a function of \(\mathbf{p}\) and \(\mathbf{\gamma}\).

\[\mathcal{I}(\mathbf{p},\mathbf{\gamma})=\nabla_{[\mathbf{p}\;\mathbf{\gamma}]}^{2} \ell_{\text{t}}=\begin{bmatrix}\partial_{\mathbf{p}}^{2}\ell_{\text{t}}& \partial_{\mathbf{p}}\partial_{\mathbf{\gamma}}\ell_{\text{t}}\\ \partial_{\mathbf{\gamma}}\partial_{\mathbf{p}}\ell_{\text{t}}&\partial_{ \mathbf{\gamma}}^{2}\ell_{\text{t}}\end{bmatrix}\] (S4.4)

The corresponding asymptotic covariance matrix if given by:

\[\mathbf{\Sigma}_{[\mathbf{p}\;\mathbf{\gamma}]}=\frac{1}{n}(\mathcal{I}(\mathbf{p}, \mathbf{\gamma}))^{-1}=\begin{bmatrix}\mathbf{\Sigma}_{\mathbf{p}\mathbf{p}}&\mathbf{ \Sigma}_{\mathbf{p}\mathbf{\gamma}}\\ \mathbf{\Sigma}_{\mathbf{\gamma}\mathbf{p}}&\mathbf{\Sigma}_{\mathbf{\gamma}\mathbf{\gamma}} \end{bmatrix}\] (S4.5)

Where the diagonal terms represent the covariance matrices for model parameters and calibration factors, and off-diagonal terms are the covariances between them. The conditional covariance for the kinetic model parameters, \(\mathbf{\Sigma}_{\mathbf{p}}\) is evaluated as:

\[\mathbf{\Sigma}_{\mathbf{p}}=\mathbf{\Sigma}_{\mathbf{p}\mathbf{p}}+\mathbf{\Sigma}_{ \mathbf{p}\mathbf{\gamma}}\left(\mathbf{\Sigma}_{\mathbf{\gamma}\mathbf{\gamma}}\right)^{-1} \mathbf{\Sigma}_{\mathbf{\gamma}\mathbf{p}}\] (S4.6)

Which can be from Taylor expansion of the residuals assuming \(\mathbf{\gamma}\) and \(\mathbf{p}\) uncorrelated and normally distributed, and using Cholesky decomposition of the covariance matrices to generated perturbations around the optimal values.

## S5 Reaction Network Type _dcs_

The shape of the stoichiometry matrices reflect the elementary reactions of the system such that the rows are associated with the mass balance of each species, and the columns are associated with reaction rates, therefore having the same size as kinetic parameters k. The \(s\)-type elementary steps embody elementary reactions between surface intermediates
exclusively. In the following, \(D*\), \(E*\) and \(F*\) are reaction intermediates for which there are no stable desorbed counterparts, which cannot be directly measured or inferred from ordinary analytical chemistry techniques.

\[A+*\frac{k_{1}^{k_{1}^{d}}}{k_{1}^{d}}A*\] ( **d.1**) \[B+*\frac{k_{2}^{d}}{k_{2}^{d}}B*\] ( **d.2**) \[C+*\frac{k_{3}^{d}}{k_{3}^{d}}C*\] ( **d.3**) \[A*+*\frac{k_{5}^{d}}{k_{3}^{d}}2D*\] ( **c.1**) \[B*+*\frac{k_{4}^{d}}{k_{4}^{d}}2E*\] ( **c.2**) \[D*+E*\frac{k_{1}^{d}}{k_{1}^{d}}F*+*\] ( **s.1**) \[F*+E*\frac{k_{5}^{d}}{k_{5}^{d}}C*+*\] ( **c.3**) \[\ln(\mathbf{k}_{0})^{T}=[3.00\text{ }2.08\text{ }3.18\text{ }2.48\text{ }2.77\text{ }3.69\text{ }6.46\text{ }6.87\text{ }5.08\text{ }4.38\text{ }6.46\text{ }5.48\text{ }6.33\text{ }5.08]\]

\[M=\left[\begin{array}{cccccccccccccccc}-1&1&0&0&0&0&0&0&0&0&0&0&0&0&0\\ 0&0&-1&1&0&0&0&0&0&0&0&0&0&0\\ 0&0&0&0&-1&1&0&0&0&0&0&0&0&0\\ 1&-1&0&0&0&0&-1&1&0&0&0&0&0&0\\ 0&0&1&-1&0&0&0&-1&1&0&0&0&0&0\\ 0&0&0&0&1&-1&0&0&0&0&0&0&1&-1\\ 0&0&0&0&0&0&2&-2&0&0&-1&1&0&0\\ 0&0&0&0&0&0&0&0&2&-2&-1&1&-1&1\\ 0&0&0&0&0&0&0&0&0&0&1&-1&-1&1\\ -1&1&-1&1&-1&1&-1&1&1&-1&1&-1&1&-1\end{array}\right]\] (S5.1)

\[\mathbf{x}^{T}=[x_{A}\text{ }x_{B}\text{ }x_{C}\text{ }x_{A*}\text{ }x_{B*} \text{ }x_{C*}\text{ }x_{D*}\text{ }x_{E*}\text{ }x_{F*}\text{ }x_{*}]\]
