# Physics-Informed Polynomial Chaos Expansions

Lukáš Novák<EMAIL_ADDRESS>Brno University of Technology, Brno, Czech
Republic  Himanshu Sharma, Michael D. Shields  {hsharm10,
<EMAIL_ADDRESS>Johns Hopkins University, Baltimore, USA

###### Abstract

Surrogate modeling of costly mathematical models representing physical systems
is challenging since it is typically not possible to create a large
experimental design. Thus, it is beneficial to constrain the approximation to
adhere to the known physics of the model. This paper presents a novel
methodology for the construction of physics-informed polynomial chaos
expansions (PCE) that combines the conventional experimental design with
additional constraints from the physics of the model. Physical constraints
investigated in this paper are represented by a set of differential equations
and specified boundary conditions. A computationally efficient means for
construction of physically constrained PCE is proposed and compared to
standard sparse PCE. It is shown that the proposed algorithms lead to superior
accuracy of the approximation and does not add significant computational
burden. Although the main purpose of the proposed method lies in combining
data and physical constraints, we show that physically constrained PCEs can be
constructed from differential equations and boundary conditions alone without
requiring evaluations of the original model. We further show that the
constrained PCEs can be easily applied for uncertainty quantification through
analytical post-processing of a reduced PCE filtering out the influence of all
deterministic space-time variables. Several deterministic examples of
increasing complexity are provided and the proposed method is applied for
uncertainty quantification.

###### keywords:

Polynomial Chaos Expansion , Physical Constraints , Surrogate modelling ,
Uncertainty Quantification , Physics-informed Machine Learning

##  1 Introduction

Mathematical models of real-life physical systems are typically highly
computationally demanding and contain various uncertain variables, thus it is
necessary to use surrogate models as computationally cheap approximations to
perform _uncertainty quantification_ (UQ), optimization, parametric studies,
and other tasks. These surrogate models treat the original model as a black-
box and fit the model from several deterministic simulations at given data
points in the design domain. Their practical use requires a sufficient number
of data points covering the space of input random variables, which can be
prohibitively expensive. However, it is beneficial to incorporate additional
constraints, e.g., known physical principles, to ensure realistic and
physically meaningful surrogate model behavior and reduce training data
demands.

Recently, there has been a considerable interest in developing machine learned
surrogate models capable of satisfying physical constraints – spawning an
entirely new field of physics-informed machine learning [1]. Much of this work
has focused on neural network models, including the hugely popular physics-
informed neural networks (PINNs) [2] and the related physics-informed neural
operators [3]. Numerous variants of deep neural networks with physical
constraints are now available, including physics-informed autoencoders [4],
the Deep Ritz Method [5], and Deep Galerkin Method [6] to name just a few. In
fact, the recent developments and applications of PINNs and related neural
networks are so numerous we cannot possibly cover them here. But neural
networks are not the only machine learning models that benefit from physical
constraints. Recent studies have also incorporated physical constraints into
Gaussian process regression models [7, 8], which have the benefit of providing
a natural measure of uncertainty. These physically constrained Gaussian
processes have been applied for UQ of highly complex physical systems [9].

This paper focuses on imposing constraints on Polynomial Chaos Expansions
(PCE) [10], which we treat here as a machine learning regression method [11,
12]. Our interest in PCE, stems from its usefulness for UQ tasks that include
moment estimation and sensitivity index computation [13], as well as its
convenient properties for surrogate modeling in UQ that derive from their
orthogonality properties with respect to certain probability measures of input
variables. Although PCE surrogate models have distinct benefits for UQ [14,
15], they suffer from the _curse of dimensionality_ since the number of PCE
terms grows rapidly with both dimension and maximum polynomial order.
Naturally, this significantly affects the necessary number of data points
needed for regression-based PCE. As a result, recent research has focused on
exploiting the analyzed mathematical model to minimize the size of the basis
set and maximize the information obtained from each data point. This involves
the construction of optimal sparse basis by adaptive algorithms [16],
determining the optimal positions of data points for regression [17, 18, 19]
and constructing various types of localized surrogates consisting of many low-
order PCEs [20, 21, 22]. In addition to the information obtained directly from
the mathematical model, there are often specific physical characteristics and
constraints derived from the nature of the quantity of interest (QoI) that
should be incorporated into the PCE. Therefore we propose a novel and
computationally efficient method to incorporate known physics in the form of
ordinary differential equations (ODEs) and partial differential equations
(PDEs) and their boundary conditions (Dirichlet, Neumann, or mixed) in the
PCEframework, referred to as physically constrained polynomial chaos (PC2)
expansion. We enforce the constraints at discrete points in the input domain,
referred to as virtual points, to solve a constrained least square
optimization problem for the PC2 coefficients. This approach can be
generalized to incorporate any form of physical constraints.

Exploiting known physics and boundary conditions in the construction of a PCE
approximation assures higher prediction accuracy, especially in regions of
input space containing an insufficient number of training points. In the
proposed PC2 framework, we can generally adopt any numerical optimization
technique to solve the constrained optimization problem. However, as the
number of PC2 coefficients and virtual points increases, some of the numerical
optimization techniques could lead to convergence issues or their
computational cost could be significant, which restricts its use to high
dimensional problems. In this paper, we focus specifically on equality-type
constraints with the objective of minimizing the sum of squared residuals
(SSR), and thus, it is possible to use the well-known method of Lagrange
multipliers. The Karush-Kuhn-Tucker (KKT) stationarity condition yields a
system of equations that can be solved for the PC2 coefficient. This approach
is a natural extension of the ordinary least squares (OLS) approach, as used
in PCE [23], to incorporate known equality-type constraints in a
computationally efficient manner. We further show that the approach can be
combined with sparse regression methods, specifically the Least Angle
Regression (LAR) [24] to efficiently reduce the basis set. We then apply the
method to several benchmark ODEs and PDEs to illustrate its performance.
Finally, we demonstrate how it can be used for general UQ purposes.

##  2 Non-intrusive Polynomial Chaos Expansion

Assume a probability space ($\Omega,\pazocal{F},\pazocal{P}$), where $\Omega$
is an event space, $\pazocal{F}$ is a $\sigma$-algebra on $\Omega$ and
$\pazocal{P}$ is a probability measure on $\pazocal{F}$. If the input variable
of a mathematical model, $Y=u(\pazocal{X})$, is a random variable
$\pazocal{X}(\omega),\omega\in\Omega$, the model response $Y$($\omega$) is
also a random variable. Assuming that $Y$ has finite variance, PCE represents
the output variable $Y$ as a function of an another random variable $\xi$
called the _germ_ with a known distribution

$Y=u(X)=u^{\textsf{PCE}}(\xi),$ (1)

and represents the function $u(\pazocal{X})$ by an infinite expansion of
polynomials that are orthogonal with respect to the probability density
$p_{\xi}(\xi)$ of the germ. These polynomials form a basis of the Hilbert
space $L^{2}$ ($\Omega,\pazocal{F},\pazocal{P}$) of all real-valued random
variables of finite variance, where $\pazocal{P}$ takes over the meaning of
the probability distribution. The orthogonality condition is given by the
inner product on $L^{2}$ ($\Omega,\pazocal{F},\pazocal{P}$) defined for any
two functions $\psi_{j}$ and $\psi_{k}$ with respect to the weight function
$p_{\xi}(\xi)$ as:

$\langle\psi_{j},\psi_{k}\rangle=\int\psi_{j}(\xi)\psi_{k}(\xi)p_{\xi}(\xi)\;\mathrm{d}\xi=0\quad\forall
j\neq k.$ (2)

This means that there are specific orthogonal polynomials associated with the
corresponding distribution of the germ. Orthogonal polynomials corresponding
to common distributions can be chosen according to the Wiener-Askey scheme
[25]. For further processing, it is beneficial to use orthonormal polynomials,
where the inner product of polynomials is equal to the Kronecker delta
$\delta_{jk}$, i.e. $\delta_{jk}=1$ if $j=k$, and $\delta_{jk}=0$ otherwise.

In the case of $\boldsymbol{\pazocal{X}}$ and ${\bm{\xi}}$ being vectors
containing $M$ independent random variables, the polynomial basis
$\Psi({\bm{\xi}})$ is multivariate and it is built up as a tensor product of
univariate orthonormal polynomials, i.e.

$\Psi_{{\boldsymbol{\alpha}}}({\bm{\xi}})=\prod_{i=1}^{M}\psi_{\alpha_{i}}(\xi_{i}),$
(3)

where ${{\boldsymbol{\alpha}}}\in\mathbb{N}^{M}$ is a set of integers called
the _multi-index_ reflecting polynomial degrees associated to each $\xi_{i}$.
The quantity of interest (QoI), i.e. the response of the mathematical model
$Y=u(\boldsymbol{\pazocal{X}})$, can then be represented as [26]

$Y=u(\boldsymbol{\pazocal{X}})=\sum_{{\boldsymbol{\alpha}}\in\mathbb{N}^{M}}\beta_{{\boldsymbol{\alpha}}}\Psi_{{\boldsymbol{\alpha}}}({\bm{\xi}}),$
(4)

where $\beta_{{\boldsymbol{\alpha}}}$ are deterministic coefficients and
$\Psi_{{\boldsymbol{\alpha}}}$ are multivariate orthonormal polynomials.

For practical computation, the PCE expressed in Eq. (4) must be truncated to a
finite number of terms, $P$. One can generally choose any truncation rule
(e.g. tensor product of polynomials up to the selected order $p$), but the
most common truncation is achieved by retaining only terms whose total degree
$|{\boldsymbol{\alpha}}|$ is less than or equal to a given $p$, in which case
the truncated set of PCE terms is then defined as

$\pazocal
A^{M,p}=\left\\{{{\boldsymbol{\alpha}}}\in\mathbb{N}^{M}:\left|{{\boldsymbol{\alpha}}}\right|=\sum_{i=1}^{M}\alpha_{i}\leq
p\right\\}.$ (5)

The cardinality of the truncated _index set_ $\pazocal A^{M,p}$ is given by

$\mathrm{card}\>\pazocal A^{M,p}=\frac{\left(M+p\right)!}{M!\>p!}\equiv P\,.$
(6)

When the PCE is truncated to a finite number of terms, there is an error
$\varepsilon$ in the approximation such that

$Y=\displaystyle
u{(\boldsymbol{\pazocal{X}})}=\sum_{{\boldsymbol{\alpha}}\in\pazocal
A}\beta_{{\boldsymbol{\alpha}}}\Psi_{{\boldsymbol{\alpha}}}({\bm{\xi}})+\varepsilon\,=Y^{PCE}+\varepsilon.$

From a statistical point of view, PCE is a simple linear regression model with
intercept. Therefore, it is possible to use _ordinary least squares_ (OLS)
regression to solve for the coefficients ${\boldsymbol{\beta}}$ that minimize
the error $\varepsilon$. To solve this regression problem for
${\boldsymbol{\beta}}$, we first generate ${n_{\mathrm{sim}}}$ realizations of
the input random vector $\boldsymbol{\pazocal{X}}$ and compute the
corresponding results of the model $\pazocal Y$, together called the
experimental design (ED) or training data set. Then, the vector of $P$
deterministic coefficients ${\boldsymbol{\beta}}$ can be determined by OLS as

${\boldsymbol{\beta}}=({\bm{\Psi}}^{T}{\bm{\Psi}})^{-1}\
{\bm{\Psi}}^{T}\pazocal Y,$ (7)

where ${\bm{\Psi}}$ is the data matrix

${\bm{\Psi}}=\left\\{\Psi_{ij}=\Psi_{j}({\bm{\xi}}^{(i)}),\;i=1,\ldots,{n_{\mathrm{sim}}},\;j=0,\ldots,P-1\right\\}.$
(8)

It is clear from Eq. (6) that $P$ is strongly dependent on the number of input
random variables $M$ and the maximum total degree of polynomials $p$.
Considering that estimating ${\boldsymbol{\beta}}$ by OLS requires at least
$\mathcal{O}(P\,\ln(P))$ samples for a stable solution [27, 28], the problem
can become computationally demanding in cases of large or strongly non-linear
models. To reduce computational expense, one can use advanced model selection
algorithms such as Least Angle Regression (LAR) [29, 24]to find an optimal set
of PCE terms, and thus reduce the number of samples needed to compute the
unknown coefficients. Further reduction can be obtained by incorporating
additional physical constraints.

##  3 Physically Constrained Polynomial Chaos Expansion (PC2) – Deterministic
Formulation

In this section, we propose a novel approach to perform PCE regression with
known physical constraints. The approach, referred to as Physically
Constrained Polynomial Chaos Expansion (PC2), expands the classes of
regression models that can obey physical constraints to include the PCE, which
is widely used for UQ. We begin by formulating the constrained regression
problem for the solution of deterministic PDEs.

Consider the general partial differential equation given by

$\displaystyle\mathcal{L}(\bm{x},t;u(\bm{x},t))=f(\bm{x},t),$
$\displaystyle\forall\bm{x}\in\mathcal{D},t\in\mathcal{T}$ (9)
$\displaystyle\mathcal{B}(\bm{x},t;u(\bm{x},t))=g(\bm{x},t,),$
$\displaystyle\forall\bm{x}\in\partial\mathcal{D},t\in\mathcal{T}$

where $\mathcal{T}\subset\mathbb{R}$, $\mathcal{D}\subset\mathbb{R}^{3}$ with
boundary $\partial\mathcal{D}$, $\mathcal{L}$ is a differential operator with
boundary operator $\mathcal{B}$, $u(\cdot)$ is the response of the system, and
$f,g$ are external forces/source terms. We aim to solve the PCE regression
problem described above as constrained by the general Eq. (9). That is, we
define the objective function by

$\displaystyle\mathcal{M}({\boldsymbol{\beta}})=\min_{{\boldsymbol{\beta}}}\sum_{j=1}^{n_{\mathrm{sim}}}\left[y^{j}-u^{\mathrm{PCE}}(\bm{x}^{j},t^{j})\right]^{2}=\min_{{\boldsymbol{\beta}}}\lVert\pazocal{Y}-{\bm{\Psi}}{\boldsymbol{\beta}}\rVert^{2}$
(10) s.t.
$\displaystyle\mathcal{L}(\bm{x}_{\mathrm{V}},t_{\mathrm{V}};u(\bm{x}_{\mathrm{V}},t_{\mathrm{V}}))=f(\bm{x}_{\mathrm{V}},t_{\mathrm{V}}),$
$\displaystyle\mathcal{B}(\bm{x}_{\mathrm{BC}},t_{\mathrm{BC}};u(\bm{x}_{\mathrm{BC}},t_{\mathrm{BC}}))=g(\bm{x}_{\mathrm{BC}},t_{\mathrm{BC}},)$

where we define a discrete sets of $n_{\mathrm{BC}}$ samples
$(\bm{x}_{\mathrm{BC}},t_{\mathrm{BC}})$ and $n_{\mathrm{V}}$ samples
$(\bm{x}_{\mathrm{V}},t_{\mathrm{V}})$ – referred to as boundary points and
virtual points, respectively – at which the boundary conditions and PDE will
be explicitly satisfied. Although $(\bm{x},t)$ is a deterministic vector, for
the PCE regression here it is assumed to be a random vector with uniform
distribution $\pazocal{U}$ for construction of PC2. This is important
primarily because it dictates that we use Legendre polynomials as basis
functions for both space and time according to the Wiener-Askey scheme.
Consequently, all deterministic input variables $(\bm{x},t)$ must be
transformed to standardized space, i.e.
${\bm{\xi}}\sim\pazocal{U}\left[-1,1\right]$, via an operator $\pazocal{T}$.
That is, we define the boundary points and virtual points by
${\bm{\xi}}_{\mathrm{BC}}=\pazocal{T}(\bm{x}_{\mathrm{BC}},t_{\mathrm{BC}})$
and ${\bm{\xi}}_{\mathrm{V}}=\pazocal{T}(\bm{x}_{\mathrm{V}},t_{\mathrm{V}})$,
respectively.

###  3.1 Physical Constraints & Lagrange Multipliers

The constrained optimization problem defined by Eq. 10 can be solved
efficiently using the method of Lagrange multipliers. The Lagrangian function
takes the following form:

$L({\boldsymbol{\beta}},\bm{\lambda})=\frac{1}{2}\mathcal{M}({\boldsymbol{\beta}})+\sum_{b=1}^{n_{\mathrm{BC}}}\lambda_{b}(a_{b}^{\mathrm{T}}{\boldsymbol{\beta}}-c_{b})+\sum_{v=1}^{n_{\mathrm{V}}}\lambda_{v}(a_{v}^{\mathrm{T}}{\boldsymbol{\beta}}-c_{v})$
(11)

Boundary conditions are prescribed by their type defined by boundary operator
$\mathcal{B}$, coordinates ${\bm{\xi}}_{\mathrm{BC}}$ and corresponding vector
$\mathbf{c}_{\mathrm{BC}}$ consisting of $n_{\mathrm{BC}}$ rows
$c_{b}=g({\bm{\xi}}_{\mathrm{BC}}^{(b)})$. Similarly, ODE/PDE constraints are
given by the differential operator $\mathcal{L}$, virtual points
${\bm{\xi}}_{\mathrm{V}}$, and corresponding vector $\mathbf{c}_{\mathrm{V}}$
consisting of $n_{\mathrm{V}}$ rows $c_{v}=u({\bm{\xi}}_{\mathrm{V}}^{(v)})$.
We then assemble these constraints into a matrix $\mathbf{A}$ where the
$a_{b}=\\{a_{b}^{j}=\mathcal{B}(\Psi_{j}({\bm{\xi}}_{\mathrm{BC}}^{(b)})),\quad
j=0,...,P-1\\}$ form the first $n_{\mathrm{BC}}$ rows and vectors
$a_{v}=\\{a_{v}^{j}=\mathcal{L}(\Psi_{j}({\bm{\xi}}_{\mathrm{V}}^{(v)})),\quad
j=0,...,P-1\\}$ then form the remaining $n_{\mathrm{V}}$ rows of the matrix.
Therefore, instead of solving the coefficients by ${\boldsymbol{\beta}}$ by
Eq. 7, we construct the following system of linear equations reflecting the
OLS solution with physical constraints obtained from the Karush–Kuhn–Tucker
(KKT) conditions:

$\underbracket{\begin{bmatrix}\bm{\Psi}^{T}\bm{\Psi}&\mathbf{A}^{T}\\\
\mathbf{A}&\bm{0}\end{bmatrix}}_{\mathrm{KKT\;matrix}}\begin{bmatrix}{\boldsymbol{\beta}}\\\
\bm{\lambda}\end{bmatrix}=\begin{bmatrix}\bm{\Psi}^{T}\pazocal Y\\\
\mathbf{c}\end{bmatrix}$ (12)

Construction of the $\mathbf{A}$ matrix requires derivatives of the PCE model.
We notice, however, that the PCE and its derivatives have the same
coefficients ${\boldsymbol{\beta}}$ but different basis functions. The
constraints can therefore be imposed by computing the derivatives of the PCE,
which can be done efficiently through term-wise derivatives of the basis
functions as follows:

$\frac{\partial^{n}f}{\partial
x_{i}^{n}}=\frac{\partial^{n}\left[\sum_{{\boldsymbol{\alpha}}\in\pazocal
A}\beta_{{\boldsymbol{\alpha}}}\Psi_{{\boldsymbol{\alpha}}}({\bm{\xi}})\right]}{\partial\xi_{i}^{n}}\Delta_{\Gamma}^{n}=\sum_{{\boldsymbol{\alpha}}\in\pazocal
A}\beta_{{\boldsymbol{\alpha}}}\frac{\partial^{n}\Psi_{{\boldsymbol{\alpha}}}({\bm{\xi}})}{\partial\xi_{i}^{n}}\Delta_{\Gamma}^{n}$
(13)

where $\Delta_{\Gamma}$ reflects the different size of the time-space variable
$x_{i}$ and standardized $\xi_{i}$, i.e. $\Delta_{\Gamma}=2/(x_{max}-x_{min})$
for Legendre polynomials defined on $\xi_{i}\in\left[-1,1\right]$ orthonormal
to $X_{i}\sim\pazocal{U}\left[x_{min},x_{max}\right]$. Note that the number of
terms in the PCE constraints (columns of $\mathbf{A}$) is identical to the
number in the original PCE and it is also possible to do basic arithmetic
operations on the basis to satisfy the prescribed $\mathcal{L}$ using the PCE
alone (i.e. we do not need to solve the PDEs). Although it is possible to
generate an arbitrary number of virtual points (similar to PINNs), it leads to
an over-determined system. Therefore, we use the optimal number obtained as
$n_{\mathrm{V}}=P-n_{\mathrm{BC}}$ which assures a well-determined KKT system
of equations. The total computational cost of the KKT system solution is
$\pazocal{O}\left((n_{\mathrm{V}}+n_{\mathrm{BC}}+P)^{3}\right)$. However,
additional information in the form of physical constraints incorporated in PC2
significantly reduces $n_{\mathrm{sim}}$, which could be crucial for costly
mathematical models.

###  3.2 Approximation Error Estimation

Once the PCE is constructed, it is crucial to estimate its accuracy. Further,
the PCE accuracy can be used to directly compare several PCEs in a sparse
adaptive framework to choose the best surrogate model. Ideally, the ED should
be divided into validation and training sets, but this might be extremely
computationally demanding in engineering applications with complex numerical
models. Therefore, in the field of UQ it is often preferred to estimate the
approximation error directly from the training set, without additional
sampling. A common choice of error measure is the mean squared error $R^{2}$,
which is well-known from machine learning and statistics. However, driving
$R^{2}$ down may lead to over-fitting. One of the most widely-used methods in
UQ is the leave-one-out cross-validation (LOO-CV) error $Q^{2}$, which can be
obtained analytically from a single PCE [30]. However, since the PC2 is
designed for small ED, $Q^{2}$ may not be recommended and $R^{2}$ is preferred
in most cases, possibly normalized by the variance of the QoI as commonly used
in surrogate modeling. The total mean squared approximation error consists of
three terms – error in the approximated function $R^{2}_{u}$, error associated
with failure to obey the PDE constraints $R^{2}_{\mathcal{L}}$, and error in
the boundary conditions $R^{2}_{\mathcal{B}}$, formally written as:

$R^{2}=R^{2}_{u}+R^{2}_{\mathcal{L}}+R^{2}_{\mathcal{B}}$ (14)

The tendency of $R^{2}$ to over-fitting is mitigated by the second and third
components of the error measure – the mean squared error in the given PDE and
boundary conditions.

##  4 PC2 for Uncertainty Quantification

In this section, we extend the deterministic formulation of PC2 to perform UQ
– which is the essential motivation for using PCE as opposed to other
physically constrained ML methods such as PINNs. Consider the general
stochastic partial differential equation given by

$\displaystyle\mathcal{L}(\bm{x},t,\boldsymbol{\pazocal{X}}(\omega);u(\bm{x},t,\boldsymbol{\pazocal{X}}(\omega)))=f(\bm{x},t,\boldsymbol{\pazocal{X}}(\omega)),$
$\displaystyle\forall\bm{x}\in\mathcal{D},t\in\mathcal{T},\omega\in\Omega$
(15)
$\displaystyle\mathcal{B}(\bm{x},t,\boldsymbol{\pazocal{X}}(\omega);u(\bm{x},t,\boldsymbol{\pazocal{X}}(\omega)))=g(\bm{x},t,\boldsymbol{\pazocal{X}}(\omega)),$
$\displaystyle\forall\bm{x}\in\partial\mathcal{D},t\in\mathcal{T},\omega\in\Omega$

where the meaning of the symbols is identical to Eq. (9) and
$\boldsymbol{\pazocal{X}}(\omega)\in\mathbb{R}^{d}$ is a $d$-dimensional
random vector having sample space $\Omega$. Now, we aim to solve the PCE
regression problem described above as constrained by the general Eq. (15).
That is, we define the objective function as

$\displaystyle\mathcal{M}({\boldsymbol{\beta}})=\min_{{\boldsymbol{\beta}}}\sum_{j=1}^{n_{\mathrm{sim}}}\left[y^{j}-u^{\mathrm{PCE}}(\bm{x}^{j},t^{j},\boldsymbol{\pazocal{X}}^{j})\right]^{2}=\min_{{\boldsymbol{\beta}}}\lVert\pazocal{Y}-{\bm{\Psi}}{\boldsymbol{\beta}}\rVert^{2}$
(16) s.t.
$\displaystyle\mathcal{L}(\bm{x}_{\mathrm{V}},t_{\mathrm{V}},\boldsymbol{\pazocal{X}}(\omega);u(\bm{x}_{\mathrm{V}},t_{\mathrm{V}},\boldsymbol{\pazocal{X}}(\omega)))=f(\bm{x}_{\mathrm{V}},t_{\mathrm{V}},\boldsymbol{\pazocal{X}}(\omega)),$
$\displaystyle\mathcal{B}(\bm{x}_{\mathrm{BC}},t_{\mathrm{BC}},\boldsymbol{\pazocal{X}}(\omega);u(\bm{x}_{\mathrm{BC}},t_{\mathrm{BC}},\boldsymbol{\pazocal{X}}(\omega)))=g(\bm{x}_{\mathrm{BC}},t_{\mathrm{BC}},\boldsymbol{\pazocal{X}}(\omega))$

where we define a discrete sets of $n_{\mathrm{BC}}$ boundary points
$(\bm{x}_{\mathrm{BC}},t_{\mathrm{BC}},\boldsymbol{\pazocal{X}})$ and
$n_{\mathrm{V}}$ virtual points
$(\bm{x}_{\mathrm{V}},t_{\mathrm{V}},\boldsymbol{\pazocal{X}})$. As in the
deterministic formulation, all input variables are transformed to standardized
space according to the Wiener-Askey scheme, i.e.
${\bm{\xi}}_{\mathrm{BC}}=\pazocal{T}(\bm{x}_{\mathrm{BC}},t_{\mathrm{BC}},\boldsymbol{\pazocal{X}})$
and
${\bm{\xi}}_{\mathrm{v}}=\pazocal{T}(\bm{x}_{\mathrm{V}},t_{\mathrm{V}},\boldsymbol{\pazocal{X}})$.
Note that the boundary and virtual points retain the random vector
$\boldsymbol{\pazocal{X}}$, and do not (necessarily) contain specific boundary
or virtual points for $\boldsymbol{\pazocal{X}}$. This is because physical
constraints are typically expressed in terms of the physical variables
$(\bm{x},t)$ and are not (in general) expressed in terms of the random
variables contained in $\boldsymbol{\pazocal{X}}$. That is,
$\boldsymbol{\pazocal{X}}$ does not necessarily affect the constraints
associated with either $\mathcal{L}$ and $\mathcal{B}$ (although it can). As a
result, the random vector $\boldsymbol{\pazocal{X}}$ usually serves only to
increase the dimension of the input random vector and thus the solution of the
Eq. (16) follows the same form as the deterministic formulation given in Eq.
(12).

If, on the other hand, the physical constraints associated with the PDE or BCs
are expressed in terms of the random variables in $\boldsymbol{\pazocal{X}}$,
virtual and/or boundary points can be used to enforce these constraints in the
same manner as performed above. If these are expressed as equality
constraints, the KKT equations to solve the constrained optimization do not
change. On the other hand, if these constraints are expressed as inequalities
(e.g. non-negative coefficients) then more advanced optimizers are required.
This will be the topic of a future work.

The resulting constrained PCE is expressed as a function of the complete set
of physical and random variables
${\bm{\xi}}=\pazocal{T}(\bm{x},t,\boldsymbol{\pazocal{X}})$. The form of the
PCE as a linear summation over orthonormal polynomials allows for powerful and
efficient post-processing. In particular, once a PCE approximation is created,
it is possible to directly estimate statistical moments of the output from the
expansion. The first statistical moment (the mean value) is simply the first
deterministic coefficient of the expansion
$\mu_{Y}=\big{<}Y^{1}\big{>}=\beta_{{\boldsymbol{0}}}$. The second raw
statistical moment, $\big{<}Y^{2}\big{>}$, can be estimated by

$\displaystyle\left\langle{{Y^{2}}}\right\rangle$
$\displaystyle=\int{\left[{\sum\limits_{{\boldsymbol{\alpha}}\in\pazocal{A}}{{\beta_{{\boldsymbol{\alpha}}}}{\Psi_{{\boldsymbol{\alpha}}}}\left({\bm{\xi}}\right)}}\right]}^{2}p_{{\bm{\xi}}}\left({\bm{\xi}}\right)\;\mathrm{d}{\bm{\xi}}=\sum\limits_{{\boldsymbol{\alpha}}_{1}\in\pazocal{A}}\sum\limits_{{\boldsymbol{\alpha}}_{2}\in\pazocal{A}}\beta_{{{\boldsymbol{\alpha}}}_{1}}\beta_{{{\boldsymbol{\alpha}}}_{2}}\int\Psi_{{{{\boldsymbol{\alpha}}}_{1}}}\left({\bm{\xi}}\right)\Psi_{{{{\boldsymbol{\alpha}}}_{2}}}\left({\bm{\xi}}\right)p_{\bm{\xi}}\left({\bm{\xi}}\right)\;\mathrm{d}{\bm{\xi}}$
(17)
$\displaystyle=\sum\limits_{{{\boldsymbol{\alpha}}}\in{\pazocal{A}}}{\beta_{{\boldsymbol{\alpha}}}^{2}}{\int{{\Psi_{{\boldsymbol{\alpha}}}}\left({\bm{\xi}}\right)}^{2}}{p_{\bm{\xi}}}\left({\bm{\xi}}\right)\;\mathrm{d}{\bm{\xi}}=\sum\limits_{{{\boldsymbol{\alpha}}}\in{\pazocal{A}}}{\beta_{{\boldsymbol{\alpha}}}^{2}}\left\langle{{\Psi_{{\boldsymbol{\alpha}}}},{\Psi_{{\boldsymbol{\alpha}}}}}\right\rangle=\sum\limits_{{{\boldsymbol{\alpha}}}\in{\pazocal{A}}}{\beta_{{\boldsymbol{\alpha}}}^{2}}.$

Considering the orthonormality of the polynomials, it is possible to obtain
the variance $\sigma_{Y}^{2}=\big{<}Y^{2}\big{>}-\mu_{Y}^{2}$ as the sum of
all squared deterministic coefficients except the intercept (which represents
the mean value). Note that the estimation of higher statistical central
moments, specifically skewness and kurtosis, are more complicated since they
require triple and quad products. These can be obtained analytically only for
certain polynomial families, e.g. formulas for Hermite and Legendre
polynomials (and their combination) can be found in [31]. Moreover, it can be
shown that PCE is in the form of the Hoeffding-Sobol decomposition and thus
Sobol indices can be obtained analytically by applying Eq. (17) to selected
subsets of the PCE terms [13].

As can be seen, the main advantage of PCE over other surrogate models lies in
the convenience of performing uncertainty quantification, i.e. the above
analytical moment and sensitivity estimates. The PC2 framework inherits this
convenience when appropriate conditions are established. As noted, the PCE is
expressed in terms of both the random variables $\boldsymbol{\pazocal{X}}$ and
the deterministic space-time coordinates $(\bm{x},t)$. Therefore, to properly
estimate statistical moments or sensitivity indices, we must condition the PCE
on $(\bm{x},t)$, yielding local space-time statistics. The numerical
estimation of local mean value $\mathbb{E}\left[f\rvert\bm{x},t\right]$ and
variance $\sigma^{2}_{\left[f\rvert\bm{x},t\right]}$ are based on a simple
rationale: deterministic values can be specified directly as constants in the
PCE terms. Following the concept of reduced PCE [31], the local variance
$\sigma^{2}_{\left[f\rvert x,t\right]}$ is affected only by PCE terms
containing the random variables $\boldsymbol{\pazocal{X}}$, while all other
terms are constants whose sum is equal to a local expected value
$\mathbb{E}\left[f\rvert x,t\right]$.

Formally stated, the random variables form a subset of all $M$ input variables
$\\{\bm{x},t,\boldsymbol{\pazocal{X}}\\}$, i.e.,
$\boldsymbol{\pazocal{X}}\subset\\{\bm{x},t,\boldsymbol{\pazocal{X}}\\}$. We
thus define the subset $\mathbf{u}\subseteq\mathrm{I}=\\{1,...,M\\}$ and its
complement ${\bf{u}^{c}}$ as

$\mathbf{u}=\left\\{\mathrm{i}\in\mathrm{I}:\xi_{i}\in\pazocal{T}(\boldsymbol{\pazocal{X}})\right\\},\qquad\mathbf{u}^{c}=\left\\{\mathrm{i}\in\mathrm{I}:\xi_{i}\in\pazocal{T}(\\{\bm{x},t\\})\right\\}.$
(18)

For the sake of clarity, the set of basis multivariate polynomials dependent
on $\boldsymbol{\pazocal{X}}$ and its complement are

${\pazocal{A}_{\boldsymbol{\pazocal{X}}}}=\left\\{{{\bf{\alpha}}\in{\pazocal{A}}:{\alpha_{k}}\neq
0\leftrightarrow
k\in{\bf{u}}\;}\right\\},\qquad{\pazocal{A}_{\sim\boldsymbol{\pazocal{X}}}}=\left\\{{{\bf{\alpha}}\in{\pazocal{A}}:{\alpha_{k}}\neq
0\leftrightarrow k\in{\bf{u}^{c}}\;}\right\\}.$ (19)

The local variance is further obtained by a simple post-processing of the
terms containing random variables. First, it is necessary to find a unique
basis $\pazocal{A}_{\boldsymbol{\pazocal{X}}}$ – defined by the unique set of
polynomials from $\pazocal{A}_{\boldsymbol{\pazocal{X}}}$. Although the
original $\pazocal{A}$ contains only unique terms, after the reduction there
could be identical terms, e.g. two terms differing only in polynomial degrees
associated to variables in ${\bf{u}^{c}}$. Further, coefficients of terms with
duplicate multi-indices are summed and assigned to the unique term in the
reduced PCE. The constants arising in mixed PCE terms containing deterministic
and random variables must be added to coefficients associated to the unique
basis. Finally, local mean and variance are obtained by commonly known
formulas directly from coefficients of the reduced PCE, i.e. expected value
from the coefficients associated with terms containing only deterministic
variables in complement set $\pazocal{A}_{\sim\boldsymbol{\pazocal{X}}}$ and
variance from coefficients associated to terms in
$\pazocal{A}_{\boldsymbol{\pazocal{X}}}$:

$\displaystyle\mathbb{E}\left[u\rvert\bm{x},t\right]$
$\displaystyle=\beta_{0}+\sum\limits_{{\boldsymbol{\alpha}}\in{\pazocal{A}_{\sim\boldsymbol{\pazocal{X}}}}}\beta_{{\boldsymbol{\alpha}}}$
(20) $\displaystyle\sigma^{2}_{\left[u\rvert\bm{x},t\right]}$
$\displaystyle=\sum\limits_{{\boldsymbol{\alpha}}\in{\pazocal{A}_{\boldsymbol{\pazocal{X}}}}}\beta^{2}_{{\boldsymbol{\alpha}}}$
(21)

##  5 Numerical Algorithms

In this section we present computationally efficient means to solve the
proposed approach for the PC2 coefficients. We then extend the approach to use
a reduced basis set defined through LAR. Finally, an iterative scheme is
proposed to fit the PC2 constrained by nonlinear PDEs. Implementation of the
proposed approach is simple and straightforward: instead of a common solution
by OLS, one can use the KKT system in existing non-intrusive regression-based
algorithms. Therefore, it is possible to use various adaptive sparse non-
intrusive algorithms for the construction of PCE. The main three algorithms
are presented in the following paragraphs.

The basic algorithm for the construction of the KKT normal equations builds
sub-matrix $\mathbf{A}$ containing the given physical constraints in the form
of boundary/initial conditions evaluated at ${\bm{\xi}}_{\mathrm{BC}}$ and PDE
evaluated at ${\bm{\xi}}_{\mathrm{V}}$. The sub-matrix $\mathbf{A}$ prescribes
additional constraints and initial conditions to ${\boldsymbol{\beta}}$
leading to efficient and accurate estimations for very low-size (or even
missing) ED. Finally, the unknown vector ${\boldsymbol{\beta}}$ is obtained by
OLS as a solution of the KKT matrix composed of $\mathbf{A}$ and the
information sub-matrix $\bm{\Psi}^{T}\bm{\Psi}$ prepared according to Eq. 12.
The whole algorithm is detailed in Algorithm 1.

Algorithm 1 KKT solver with virtual samples

1:Input: ED ($\\{\bm{x},t,\boldsymbol{\pazocal{X}}\\}$ and $\pazocal Y$),
$\bm{\Psi}$, boundary conditions ($\bm{\xi}_{\text{\tiny BC}}$,
$\mathbf{c}_{\mathrm{BC}}$), $\mathcal{L}$, $\mathcal{B}$, $f$

2:for $b\leftarrow 1$ to $n_{\mathrm{BC}}$ do

3: for $j\leftarrow 1$ to $P$ do

4: evaluate $a_{b}^{j}=\mathcal{B}(\Psi_{j}({\bm{\xi}}_{\mathrm{BC}}^{(b)}))$

5: add $a_{b}^{j}$ to $\mathbf{A}$

6: end for

7:end for

8:sample $n_{\mathrm{V}}=P-n_{\text{\tiny BC}}$ virtual points
${\bm{\xi}}_{\mathrm{V}}$ (MC, LHS etc.)

9:for $v\leftarrow 1$ to $n_{\mathrm{V}}$ do

10: for $j\leftarrow 1$ to $P$ do

11: evaluate $a_{v}^{j}=\mathcal{L}(\Psi_{j}({\bm{\xi}}_{\mathrm{V}}^{(v)}))$

12: add $a_{v}^{j}$ to $\mathbf{A}$

13: end for

14: evaluate $c_{v}=f({\bm{\xi}}_{\mathrm{V}}^{(v)})$

15: add $c_{v}$ to $\mathbf{c}_{\mathrm{V}}$

16:end for

17:assemble vector
$\mathbf{c}=\left[\mathbf{c}_{\mathrm{BC}},\mathbf{c}_{\mathrm{V}}\right]$

18:construct KKT normal equations according to Eq. 12

19:solve the system by OLS

20:Output: PCE coefficients ${\boldsymbol{\beta}}$

Since the KKT-based estimation of the PCE coefficients ${\boldsymbol{\beta}}$
is computationally efficient, Algorithm 1 can be easily combined with adaptive
sparse solvers such as LAR. For the sake of clarity, the whole algorithm
combining LAR and KKT solver with virtual samples is provided in Algorithm 2.

Algorithm 2 LAR-KKT algorithm

1:Input: ED ($\\{\bm{x},t,\boldsymbol{\pazocal{X}}\\}$ and $\pazocal Y$),
$\bm{\Psi}$, boundary conditions ($\bm{\xi}_{\text{\tiny BC}}$,
$\mathbf{c}_{\mathrm{BC}}$), $\mathcal{L}$, $\mathcal{B}$, $f$,
$p_{\mathrm{max}}$

2:for $p\leftarrow 1$ to $p_{\mathrm{max}}$ do

3: generate set of basis functions $\pazocal{A}$

4: identify the sequence of the most important basis functions by LAR

5: for $i\leftarrow 1$ to $P$ do

6: construct $\bm{\Psi}$

7: estimate ${\boldsymbol{\beta}}$ by KKT (Algorithm 1)

8: get approximation error $R^{2}_{\mathrm{LAR}}$

9: Optional: Check over-fitting

10: end for

11:end for

12:Output: $\\{\bm{\Psi},{\boldsymbol{\beta}}\\}$ associated to the lowest
$R^{2}_{\mathrm{LAR}}$

Figure 1: Graphical representation of LAR-KKT algorithm combining
$p$-adaptivity, LAR algorithm and KKT solver with virtual samples.

The LAR-KKT algorithm combines information from the ED and physical
constraints to identify the optimal form of the PCE approximation as
schematically depicted in Fig. 1. Information obtained from the ED is
specifically utilized to identify the most suitable basis set by standard LAR.
Physical constraints are further used to estimate ${\boldsymbol{\beta}}$
satisfying the given constraints and a suitable error measure $R^{2}$ is then
used for the selection of the most accurate PC2. However, the KKT solver can
be also used without Algorithm 2, which dramatically decreases the
computational cost of PC2 for low to moderate dimensions of the input space.
Note that the KKT solution can be generally combined with any existing
adaptive sparse regression algorithm by simply replacing standard OLS with the
solution of the KKT system.

The solution is more complicated for non-linear PDEs since it is necessary to
know the complete form of the PCE and thus ${\boldsymbol{\beta}}$. Naturally,
the PCE and its derivatives evaluated at the virtual samples are not known at
the start of the algorithm. Therefore, it is necessary to construct the PCE
(estimate coefficients ${\boldsymbol{\beta}}$) without virtual samples in the
first step and iteratively improve the estimate by evaluation of the PDE for
the virtual samples. The KKT algorithm for simple non-linear PDEs is presented
in Algorithm 3.

Algorithm 3 Iterative KKT solver for non-linear PDEs

1:Input: ED ($\\{\bm{x},t,\boldsymbol{\pazocal{X}}\\}$ and $\pazocal Y$),
$\bm{\Psi}$, boundary conditions ($\bm{\xi}_{\text{\tiny BC}}$,
$\mathbf{c}_{\mathrm{BC}}$), $\mathcal{L}$, $\mathcal{B}$

2:get $\mathbf{c}_{\mathrm{BC}}$ and ${{\boldsymbol{\beta}}}_{0}$ from
Algorithm 1 without $\mathcal{L}$

3:sample $n_{\mathrm{V}}=P-n_{\text{\tiny BC}}$ virtual points
${\bm{\xi}}_{\mathrm{V}}$ (MC, LHS etc.)

4:for $v\leftarrow 1$ to $n_{\mathrm{V}}$ do

5: for $j\leftarrow 1$ to $P$ do

6: evaluate
$a_{v}^{j}=\mathcal{L}_{{{\boldsymbol{\beta}}}_{0}}(\Psi_{j}({\bm{\xi}}_{\mathrm{V}}^{(v)}))$

7: add $a_{v}^{j}$ to $\mathbf{A}$

8: end for

9: evaluate $c_{v}=f({\bm{\xi}}_{\mathrm{V}}^{(v)})$

10: add $c_{v}$ to $\mathbf{c}_{\mathrm{V}}$

11:end for

12:assemble vector
$\mathbf{c}=\left[\mathbf{c}_{\mathrm{BC}},\mathbf{c}_{\mathrm{V}}\right]$

13:for $i\leftarrow 1$ to $n_{\mathrm{iter}}$ do

14: construct KKT normal equations

15: get ${\boldsymbol{\beta}}_{i}$ by least squares

16: for $v\leftarrow 1$ to $n_{\mathrm{V}}$ do

17: for $j\leftarrow 1$ to $P$ do

18: evaluate
$a_{v}^{j}=\mathcal{L}_{{{\boldsymbol{\beta}}_{i}}}(\Psi_{j}({\bm{\xi}}_{\mathrm{V}}^{(v)}))$

19: add $a_{v}^{j}$ to $\mathbf{A}$

20: end for

21: evaluate $c_{v}=u({\bm{\xi}}_{\mathrm{V}}^{(v)})$

22: add ${c}_{v}$ to $\mathbf{c}_{\mathrm{V}}$

23: end for

24:end for

25:Output: PCE coefficients
${\boldsymbol{\beta}}={\boldsymbol{\beta}}_{n_{\mathrm{iter}}}$

##  6 Numerical Experiments

To present the capabilities of the proposed PC2 approach, several examples are
presented for various types of PDEs. First, PC2 is applied to three 1D
examples of different nature, including an inhomogeneous ODE, an ODE with BCs
of arbitrary order, and a non-linear ODE. PC2 is then used to approximate the
solution to the 2D wave equation PDE and the heat equation PDE with an
uncertain input variable. PC2 is constructed using Algorithm 1 (KKT) and the
sparse adaptive Algorithm 2 using KKT (LAR-KKT) and the results are compared
to the standard unconstrained LAR method implemented in UQPy [32, 33]. All
compared methods are extended by $p$ adaptivity ($p\in\left[5,25\right]$)
governed by $R^{2}$ measure to obtain the highest possible accuracy. The ED is
generated by Latin Hypercube Sampling (LHS) and virtual samples are generated
by Crude Monte Carlo (MC) unless stated otherwise. The number of virtual
samples is $n_{\mathrm{V}}=P-n_{\mathrm{BC}}$. Note that Dirichlet BCs are
incorporated in the ED for LAR method for a fair comparison, though additional
information from the PDE and higher-order BCs can only be incorporated by PC2.
All numerical results were replicated for 100 independent trials and
statistical results are presented in convergence plots showing the mean
$\pm\sigma$ interval. The convergence is measured by two global quantities
characterizing the whole input domain: mean squared error (Eq. 22) and maximum
squared error (Eq. 23). Moreover, detailed insight is given also by a local
measure represented by the absolute difference between the original model $u$
and its approximation $\tilde{u}$ (Eq. 24). For numerical evaluation of
$R^{2}_{u}$ and $R^{2}_{\mathcal{L}}$ we use a validation set containing
$1000^{M}$ samples generated by MC together with $100$ samples for
$R^{2}_{\mathcal{B}}$.

$\displaystyle\epsilon_{\mathrm{mean}}\coloneqq$
$\displaystyle\mathbb{E}\Big{[}\big{(}u(\bm{x},t,\boldsymbol{\pazocal{X}})-\tilde{u}(\bm{x},t,\boldsymbol{\pazocal{X}})\big{)}^{2}\Big{]}$
(22) $\displaystyle\epsilon_{\mathrm{max}}\coloneqq$
$\displaystyle\max\Big{[}\big{(}u(\bm{x},t,\boldsymbol{\pazocal{X}})-\tilde{u}(\bm{x},t,\boldsymbol{\pazocal{X}})\big{)}^{2}\Big{]}$
(23) $\displaystyle\lvert\epsilon\rvert\coloneqq$ $\displaystyle\lvert
u(\bm{x},t,\boldsymbol{\pazocal{X}})-\tilde{u}(\bm{x},t,\boldsymbol{\pazocal{X}})\rvert$
(24)

###  6.1 1D Poisson Equation: Inhomogeneous ODE

The very first example shows the application of the PC2 to simple 1D
inhomogeneous ODE. Homogeneous and inhomogeneous ODEs/PDEs solved by PC2
differ only in vectors $\mathbf{c}_{\mathrm{V}}$ – containing specific values
for inhomogenous PDEs and zeros for homogenous PDEs. Naturally,
$\mathbf{c}_{\mathrm{V}}$ could contain constants or results of functions
$f(\bm{x}_{\mathrm{V}},t_{\mathrm{V}},\boldsymbol{\pazocal{X}}(\omega))$. The
Poisson equation and assumed BC are as follows:

$\displaystyle\Delta u(x)$ $\displaystyle=2,\qquad x\in[-1,1],$ (25)
$\displaystyle\left.\dfrac{\partial u(x)}{\partial x}\right|_{x=1}$
$\displaystyle=4,\qquad u(-1)=0.$

where $\Delta$ is the Laplace operator.

The obtained numerical results are summarized in Fig. 2. Note that although
the figure contains also $\pm\sigma$ intervals, these are not visible for the
KKT method due to extremely low variance in the obtained accuracy. The
accuracy is measured by three quantities
$R^{2}_{u},R^{2}_{\mathcal{L}},R^{2}_{\mathcal{B}}$ and their sum representing
total error ($R^{2}$). Note that while $R^{2}_{u}$ is typically a sufficient
measure in standard surrogate modeling, in physically informed surrogate
models it is also important to measure $R^{2}_{\mathcal{L}},\
R^{2}_{\mathcal{B}}$ to ensure the physical constraints are being satisfied.
Standard LAR shows convergence to the exact solution measured by $R^{2}_{u}$
from 8 samples in the domain. The KKT method is exact even for two
realizations in the domain. The KKT-LAR method combining the benefits of
sparse solution and KKT approach converges in 4 samples while having large
variance for 3 samples. This behavior of the method can be seen also in the
next linear ODE example, where the effect of sparsity is not so significant
and it is thus generally more efficient to use the standard KKT method in
simple linear ODEs. Although all compared methods converge to the exact
solution for maximum number of $n_{\mathrm{sim}}=10$, standard LAR method
leads to significant errors for lower $n_{\mathrm{sim}}$. Note that although
LAR has almost perfect accuracy measured by $R^{2}_{u}$ for
$n_{\mathrm{sim}}=4$, it still has significant errors in the remaining
criteria while KKT and LAR-KKT converge consistently in total $R^{2}$. For the
sake of completeness, the maximum errors (bottom row) complement the mean
errors (top row).

Figure 2: Obtained numerical results for 1D Poisson equation for increasing
number of samples in the domain. Top row shows mean squared errors and bottom
row shows maximum squared errors obtained by simple KKT (Algorithm 1), LAR-KKT
(Algorithm 2) and standard LAR method. Each column corresponds to a component
of total error measure $R^{2}$. Figure 3: Obtained numerical results for 1D
Euler equation for increasing number of samples in the domain. Top row shows
mean squared errors and bottom row shows maximum squared errors obtained by
simple KKT (Algorithm 1), LAR-KKT (Algorithm 2) and standard LAR method. Each
column corresponds to a component of total error measure $R^{2}$.

###  6.2 1D Euler Equation: Arbitrary BC

PC2 offers efficient $n$th derivatives of the basis functions and thus is
possible to solve PDES containing generally derivatives of arbitrary order in
$\mathcal{L}$ and $\mathcal{B}$, if the $n$th derivative of the PCE basis
exists. This example presents the 1D Euler beam equation with given BCs as
follows:

$\displaystyle\frac{\partial^{4}u(x)}{\partial x^{4}}+1$
$\displaystyle=0,\qquad x\in[0,1],$ (26) $\displaystyle u(0)=0,\qquad
u^{\prime}(0)$ $\displaystyle=0\qquad u^{\prime\prime}(1)=0,\qquad
u^{\prime\prime\prime}(1)=0$

This ODE with the prescribed BCs represents the deflection of a cantilever
beam with a uniformly distributed load having intensity $q=1$. The obtained
results are shown in Fig. 3 in the identical form as in the previous example.
The convergence trends of the compared methods are similar to the previous
example, though the differences between standard LAR and PC2 are more
significant. Consistency of PC2 convergence with respect to all criteria can
be clearly seen in convergence plots of the LAR-KKT method, while the standard
LAR approach leads to similar overfitting for $n_{\mathrm{sim}}=6$ similar to
the previous example.

###  6.3 1D Logistic Equation: Non-linear ODE

The last ODE example presenting capabilities of the proposed PC2 is a non-
linear ODE. Non-linear ODEs/PDEs solved by PC2 require iterative construction
of KKT system using Algorithm 3 with LHS sampling for virtual points, though
it is not applicable to strongly non-linear PDEs. For more complicated PDEs,
it will be necessary to use more advanced and more computationally expensive
optimizers. The example is a 1D Logistic equation commonly representing
population growth in the following form:

$\displaystyle\dfrac{\partial u(x)}{\partial x}$
$\displaystyle=u(x)\left(1-u(x)\right),\qquad x\in[-5,5],$ (27)
$\displaystyle\quad u(0)$ $\displaystyle=0.5$

The obtained numerical results are summarized in Fig. 4. The convergence trend
is similar for all methods, but PC2 achieves significantly higher accuracy.
Note that the simple KKT solver is the most accurate for a very low number of
samples, though for increasing number of samples in the domain the benefit of
the sparse solution obtained by the LAR-KKT solver leads to superior accuracy.
Interestingly, the convergence trend of standard LAR is very similar to LAR-
KKT though there is a significant difference in the absolute accuracy.
Although the total accuracy of LAR is low, the accuracy in the given BC is
perfect, possibly indicating that LAR is overfitting. In contrast to LAR, both
PC2 methods consistently converge in all presented criteria. Relatively slow
convergence of $R_{\mathcal{L}}$ and $R_{\mathcal{B}}$ obtained by the KKT
solution clearly show the necessity of sparse solvers in this example, since
some basis functions associated with high $p$ are the most suitable basis of
an accurate PCE approximation. However, the cardinality of the full set of
basis functions used in the KKT method is high, and thus it needs a high
$n_{\mathrm{sim}}$ for accurate coefficient estimates.

Figure 4: Obtained numerical results for 1D Logistic equation for increasing
number of samples in the domain. Top row shows mean squared errors and bottom
row shows maximum squared errors obtained by simple KKT (Algorithm 1), LAR-KKT
(Algorithm 2) and standard LAR method. Each column corresponds to a component
of total error measure $R^{2}$.

###  6.4 Wave equation: Time-dependent PDE

The previous 1D examples clearly show possibilities of PC2 for various types
of ODEs. Naturally, the PC2 can be used for higher number of input variables
as will be presented in the time-dependent wave equation:

$\displaystyle\frac{\partial^{2}u(x,t)}{\partial t^{2}}$
$\displaystyle=4\frac{\partial^{2}u(x,t)}{\partial x^{2}},\qquad
x\in[0,1],t\in[0,2]$ (28) $\displaystyle\frac{\partial u(x,0)}{\partial t}$
$\displaystyle=0,\qquad u(0,t)=u(1,t)=0,\qquad u(x,0)=\sin{(\pi x)}.$

Although the extension of PC2 for $n$-dimensional problems is straightforward
since it is based on standard PCE, it brings one significant difference in
comparison to 1D examples: boundary and initial conditions are now always
prescribed point-wise and thus the accuracy of PC2 is also dependent on number
and position of samples on the boundaries. In this paper, we use the DeepXDE
python package [34] for deterministic sampling on the boundaries and thus
their position does not influence the variance of the convergence plots.
Convergence plots for $n_{\mathrm{BC}}=10$ samples on the boundaries are shown
in Fig. 5. PC2 converges rapidly to very accurate approximation, although the
$R_{u}^{2}$ for very low $n_{\mathrm{sim}}$ is identical to LAR. Moreover, one
can see that the PDE error $R_{\mathcal{L}}^{2}$ is negligible for PC2, but
LAR does not respect the PDE, and thus it does not converge to the original
model. Note that the difference between PC2 and LAR is generally significantly
higher for increasing $n_{\mathrm{BC}}$ due to additional samples with Neumann
BC reflected only by PC2.

Figure 5: Obtained numerical results for wave equation for increasing number
of samples in the domain. Top row shows mean squared errors and bottom row
shows maximum squared errors obtained by simple KKT (Algorithm 1), LAR-KKT
(Algorithm 2) and standard LAR method. Each column corresponds to a component
of total error measure $R^{2}$.

A detailed comparison of the PC2 and standard LAR method for a selected
realization of the algorithm can be seen in Fig.6. The top row shows the
solution of the wave equation for a selected $t$ where the dashed line
represents the analytical reference solution and the colored line shows the
solution by LAR (left) and PC2 (right). The approximations are discretized to
200 points and their colors corresponds to the logarithm of the local squared
errors of the PDE $\log R^{2}_{\mathcal{L}}$. The bottom row shows
approximations of the wave equation over the entire input space together with
given the ED, $\bm{x}_{\mathrm{BC}}$ and $\bm{x}_{\mathrm{V}}$. It is clear
that standard LAR does not lead to sufficiently accurate approximation of
$u(x,t)$, which is also reflected by high $R^{2}_{\mathcal{L}}$. Furthermore,
even when LAR produces similar errors $R_{u}^{2}$, its PDE errors
$R_{\mathcal{L}}^{2}$ are much larger as illustrated in Fig. 7. In this case,
LAR leads to $R<10^{-3}$ though it is clear from $\log R^{2}_{\mathcal{L}}$
that the prescribed PDE is not respected. Note that PC2 leads to perfect
approximations in both selected realizations.

Figure 6: Selected realization of LAR (left) and PC2 (right) for the 2D wave
equation. The left pane shows significant error $R_{u}^{2}$ and
$R_{\mathcal{L}}^{2}$ for standard LAR. The right pane shows a very accurate
solution from PC2 using identical ED, $\bm{x}_{\mathrm{BC}}$ and
$\bm{x}_{\mathrm{V}}$. Figure 7: Selected realization of LAR (left) and PC2
(right) for the 2D wave equation. The left and right panes show similar errors
$R_{u}^{2}$ (bottom) but PC2 shows much lower error in PDE
$R_{\mathcal{L}}^{2}$ (top) using identical ED, $\bm{x}_{\mathrm{BC}}$ and
$\bm{x}_{\mathrm{V}}$.

###  6.5 Heat Equation: Uncertainty Quantification

The last numerical example presents the main advantage of PC2 – analytical UQ
of the original model containing deterministic and random variables as
proposed in section 4. This example is thus divided to two parts: (i) a
convergence study assuming all input variables to be deterministic, which
shows the efficiency of PC2 as a surrogate model. and (ii) UQ of the given
heat equation with uniformly distributed coefficient of thermal diffusivity
$\mathcal{D}$ in the following form:

$\displaystyle\frac{\partial u(x,t)}{\partial t}$
$\displaystyle=\mathcal{D}\frac{\partial^{2}u(x,t)}{\partial x^{2}},\qquad
x\in[0,1],\quad t\in[0,1],\quad\mathcal{D}\sim\pazocal{U}[0.2,0.8]$ (29)
$\displaystyle u(0,t)$ $\displaystyle=u(1,t)=0,\qquad u(x,0)=\sin(\pi x)$

#### Deterministic value of coefficient of thermal diffusivity

In the first part of this example, we assume a deterministic value of the
coefficient of thermal diffusivity $\mathcal{D}=0.4$. The results in Fig. 8
shows the convergence for $n_{\mathrm{BC}}=10$ and increasing number of
samples in the space-time domain $n_{\mathrm{sim}}$. The benefit of additional
information from the physical constraints can be clearly seen from the
comparison of the obtained convergence plots. The convergence rate of standard
LAR is significantly lower, while both PC2 algorithms achieve near-perfect
accuracy for $n_{\mathrm{sim}}=15$. Note that LAR-KKT has higher variance
since the selected sparse set of basis functions is dependent on the given ED
as expected, which is clearly visible especially in the second row showing
$\epsilon_{\mathrm{max}}$.

Figure 8: Obtained numerical results of heat equation for $n_{\mathrm{BC}=10}$
and increasing number of samples in the time-space domain. Top row shows mean
squared errors and bottom row shows maximum squared errors obtained by simple
KKT (Algorithm 1), LAR-KKT (Algorithm 2) and standard LAR method. Each column
corresponds to a component of total error measure $R^{2}$.

#### Uncertain coefficient of thermal diffusivity

In the present example – a heat equation with uncertain coefficient
$\mathcal{D}\sim\pazocal{U}[0.2,0.8]$ – there are 2 deterministic variables
(space coordinate $x$ and time $t$) and one random variable $\mathcal{D}$. To
efficiently quantify solution uncertainty, PC2 is based only on
$n_{\mathrm{BC}}=90$ samples and virtual points. We do not run the
deterministic simulation using samples of $\mathcal{D}$. This scenario is of
practical interest, since the traditional solution of PDEs with uncertain
inputs requires repetitive calculations with different realizations of the
random variables which could be costly. The number of samples on the
boundaries is selected relatively high to specifically present the advantages
of PC2 for UQ, since a general convergence for PC2 was shown in the previous
deterministic case. Note that the PC2 is based only on BC and virtual samples,
and thus the computational cost is not significantly affected by the
assumption of uncertain $\mathcal{D}$ since it affects only the virtual point.
An empty ED leads to very efficient UQ, since one can easily modify the
stochastic model of the input variables – number of random variables, their
distributions etc. without re-evaluation of the ED as necessary in standard
surrogate models. Once the PC2 is available, it is possible to use it for UQ –
estimation of local means and quantiles (e.g. $\pm\sigma$) as depicted in Fig.
9, which compares the mean solution and mean $\pm\sigma$ for PC2 with Monte
Carlo simulation using LHS with $n_{\mathrm{sim}=}10^{5}$ samples. Note that
there are negligible errors in the local mean values. The estimated local
variance can also be seen in a comparison to the LHS solution in Fig. 10.
Although the PC2 is based only on BC and PDE, it leads to a very accurate
approximation and it reflects the uncertainty of $\mathcal{D}$ very well both
in predictions (see the approximations and errors for selected values of
$\mathcal{D}$ in Fig. 11) and analytical estimation of local statistics.

Figure 9: Estimates of the local mean values at given coordinates
$\mathbb{E}\left[u\mid x,t\right]$ (middle), and local $\pm\sigma_{\left[u\mid
x,t\right]}$ quantiles (left and right). Comparison of analytical solution by
PC2 (top), original model (middle), and their absolute difference (bottom).
Figure 10: Estimates of the local variance at given coordinates
$\sigma^{2}_{\left[u\mid x,t\right]}$. Comparison of numerical estimation by
LHS (left), analytical solution by PC2 (middle), and their absolute difference
(right).

Figure 11: Approximation of heat equation by PC2 based on $n_{\mathrm{BC}}=90$
samples. Each row shows the original function, the PC2 approximation, and the
error corresponding to a selected realization of $\mathcal{D}$.

##  7 Conclusions & Further Work

A novel methodology for construction of physics-informed non-intrusive
regression-based PCE, referred to as PC2 was proposed in this paper. Given
physical constraints in the form of ordinary and partial differential
equations and their boundary conditions are imposed in a constrained least
squares solution to assemble a KKT system with Lagrange multipliers. The
proposed approach does not significantly increase the computational cost to
estimate the PCE coefficients, but improves accuracy considerably. The
proposed solution can be further employed with existing adaptive algorithms
instead of conventional ordinary least squares. An algorithm based on least
angle regression is then proposed to achieve a constrained sparse solution
with $p$-adaptivity. The presented approach was developed specifically for UQ
of costly mathematical models of physical systems using a small ED. PC2 allows
for analytical UQ similar to standard PCE by simple post-processing of the
coefficients, though the PC2 must first be reduced to exclude the influence of
deterministic variables. From the obtained numerical results, it is clear that
PC2 achieves significant improvement in the accuracy of PCE in terms of
surrogate model prediction error, but also with respect to error in the PDE
and boundary conditions for low-to-mid size ED with little additional
computational cost. We compared two algorithms for PC2 construction – sparse
solver LAR-KKT and a KKT-based PCE with a full set of basis functions, and it
can be concluded that while LAR-KKT does not always lead to superior results,
its benefits are more significant with increasing dimension of input random
vector. Note that both algorithms can be applied in tandem in practical
applications, since the PC2 using the full set of basis functions increases
the computational cost of LAR-KKT just by one additional iteration. Finally,
it was shown in the final numerical example that it is possible to create PC2
using only a set of virtual samples and boundary conditions, which leads to
extremely efficient UQ in comparison to other existing techniques.

This work opens the door for several open questions that should be
investigated in further research. First, the proposed approach is limited to
problems with linear and weakly nonlinear differential constraints. More
complicated non-linear PDEs will require advanced optimization techniques and
thus PC2 must be extended also for broader classes of applications and
suitable optimization techniques should be identified. Moreover, polynomial
basis functions restrict the application of PC2 to sufficiently smooth
operations, which could be alleviated by coupling the proposed method with
domain-decomposition techniques [21, 35, 20]. Further, PC2 will be combined
with active learning techniques such as $\Theta$ criterion[18] recently
proposed by authors of this paper, which can be used for construction of
optimal ED as well as optimal set of virtual samples.

## Acknowledgments

The first author acknowledges the financial support provided by the Czech
Science Foundation under project number 23-04712S. Additionally, a part of
this research was conducted during the research stay of the first author at
Johns Hopkins University supported by the stipend and partnership program of
Brno University of Technology for Excellence 2023. MDS and HS have been
supported by the Defense Threat Reduction Agency, Award HDTRA1202000.

## References

  * [1] G. E. Karniadakis, I. G. Kevrekidis, L. Lu, P. Perdikaris, S. Wang, L. Yang, Physics-informed machine learning, Nature Reviews Physics 3 (6) (2021) 422–440.
  * [2] M. Raissi, P. Perdikaris, 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 (2019) 686–707.
  * [3] S. Goswami, A. Bora, Y. Yu, G. E. Karniadakis, Physics-informed neural operators, arXiv preprint arXiv:2207.05748 (2022).
  * [4] N. B. Erichson, M. Muehlebach, M. W. Mahoney, Physics-informed autoencoders for lyapunov-stable fluid flow prediction, arXiv preprint arXiv:1905.10866 (2019).
  * [5] E. Kharazmi, Z. Zhang, G. E. Karniadakis, Variational physics-informed neural networks for solving partial differential equations, arXiv preprint arXiv:1912.00873 (2019).
  * [6] J. Sirignano, K. Spiliopoulos, Dgm: A deep learning algorithm for solving partial differential equations, Journal of computational physics 375 (2018) 1339–1364.
  * [7] L. P. Swiler, M. Gulian, A. L. Frankel, C. Safta, J. D. Jakeman, A survey of constrained gaussian process regression: Approaches and implementation challenges, Journal of Machine Learning for Modeling and Computing 1 (2) (2020).
  * [8] G. Pang, G. E. Karniadakis, Physics-informed learning machines for partial differential equations: Gaussian processes versus neural networks, Emerging Frontiers in Nonlinear Science (2020) 323–343.
  * [9] H. Sharma, J. A. Gaffney, D. Tsapetis, M. D. Shields, Learning thermodynamically constrained equations of state with uncertainty, arXiv preprint arXiv:2306.17004 (2023).
  * [10] N. Wiener, The homogeneous chaos, American Journal of Mathematics 60 (4) (1938) 897–936.  doi:10.2307/2371268.
  * [11] N. Lüthen, S. Marelli, B. Sudret, Sparse polynomial chaos expansions: Literature survey and benchmark, SIAM/ASA Journal on Uncertainty Quantification 9 (2) (2021) 593–649.  doi:10.1137/20M1315774.
  * [12] E. Torre, S. Marelli, P. Embrechts, B. Sudret, Data-driven polynomial chaos expansion for machine learning regression, Journal of Computational Physics 388 (2019) 601–623.
  * [13] B. Sudret, Global sensitivity analysis using polynomial chaos expansions, Reliability Engineering & System Safety 93 (7) (2008) 964–979.  doi:10.1016/j.ress.2007.04.002.
  * [14] T. Crestaux, O. L. Maître, J.-M. Martinez, Polynomial chaos expansion for sensitivity analysis, Reliability Engineering & System Safety 94 (7) (2009) 1161–1172.  doi:10.1016/j.ress.2008.10.008.
  * [15] N.-Z. Chen, C. Guedes Soares, Spectral stochastic finite element analysis for laminated composite plates, Computer methods in Applied Mechanics and Engineering 197 (51) (2008) 4830–4839.  doi:10.1016/j.cma.2008.07.003.
  * [16] N. Lüthen, S. Marelli, B. Sudret, Automatic selection of basis-adaptive sparse polynomial chaos expansions for engineering applications, International Journal for Uncertainty Quantification 12 (3) (2022) 49–74.
  * [17] N. Fajraoui, S. Marelli, B. Sudret, Sequential design of experiment for sparse polynomial chaos expansions, SIAM/ASA Journal on Uncertainty Quantification 5 (1) (2017) 1061–1085.  doi:10.1137/16m1103488.
  * [18] L. Novák, M. Vořechovský, V. Sadílek, M. D. Shields, Variance-based adaptive sequential sampling for polynomial chaos expansion, Computer Methods in Applied Mechanics and Engineering 386 (2021) 114105.  doi:https://doi.org/10.1016/j.cma.2021.114105.
  * [19] M. Thapa, S. B. Mulani, R. W. Walters, Adaptive weighted least-squares polynomial chaos expansion with basis adaptivity and sequential adaptive sampling, Computer methods in Applied Mechanics and Engineering 360 (2020) 112759.  doi:10.1016/j.cma.2019.112759.
  * [20] X. Wan, G. E. Karniadakis, An adaptive multi-element generalized polynomial chaos method for stochastic differential equations, Journal of Computational Physics 209 (2) (2005) 617–642.  doi:10.1016/j.jcp.2005.03.023.
  * [21] L. Novák, M. D. Shields, V. Sadílek, M. Vořechovský, Active learning-based domain adaptive localized polynomial chaos expansion (2023).  arXiv:2301.13635.
  * [22] S. Marelli, P.-R. Wagner, C. Lataniotis, B. Sudret, STOCHASTIC SPECTRAL EMBEDDING, International Journal for Uncertainty Quantification 11 (2) (2021) 25–47.  doi:10.1615/int.j.uncertaintyquantification.2020034395.
  * [23] S.-K. Choi, R. V. Grandhi, R. A. Canfield, C. L. Pettit, Polynomial chaos expansion with latin hypercube sampling for estimating response variability, AIAA Journal 42 (6) (2004) 1191–1198.  doi:10.2514/1.2220.
  * [24] G. Blatman, B. Sudret, Adaptive sparse polynomial chaos expansion based on least angle regression, Journal of Computational Physics 230 (6) (2011) 2345–2367.  doi:10.1016/j.jcp.2010.12.021.
  * [25] D. Xiu, G. E. Karniadakis, The Wiener–Askey polynomial chaos for stochastic differential equations, SIAM Journal on Scientific Computing 24 (2) (2002) 619–644.  doi:10.1137/s1064827501387826.
  * [26] R. G. Ghanem, P. D. Spanos, Stochastic Finite Elements: A Spectral Approach, Springer New York, 1991.  doi:10.1007/978-1-4612-3094-6.
  * [27] A. Cohen, G. Migliorati, Optimal weighted least-squares methods, The SMAI journal of computational mathematics 3 (2017) 181–203.  doi:10.5802/smai-jcm.24.
  * [28] A. C. Narayan, J. Jakeman, T. Zhou, A Christoffel function weighted least squares algorithm for collocation approximations, Math. Comput. 86 (2017) 1913–1947.  doi:10.1090/mcom/3192.
  * [29] B. Efron, T. Hastie, I. Johnstone, R. Tibshirani, Least angle regression, The Annals of Statistics 32 (2) (2004) 407–451.  doi:10.2307/3448465.
  * [30] G. Blatman, B. Sudret, An adaptive algorithm to build up sparse polynomial chaos expansions for stochastic finite element analysis, Probabilistic Engineering Mechanics 25 (2) (2010) 183–197.  doi:10.1016/j.probengmech.2009.10.003.
  * [31] L. Novák, On distribution-based global sensitivity analysis by polynomial chaos expansion, Computers & Structures 267 (2022) 106808.  doi:10.1016/j.compstruc.2022.106808.
  * [32] A. Olivier, D. Giovanis, B. Aakash, M. Chauhan, L. Vandanapu, M. D. Shields, UQpy: A general purpose python package and development environment for uncertainty quantification, Journal of Computational Science 47 (2020) 101204\.
  * [33] D. Tsapetis, M. D. Shields, D. G. Giovanis, A. Olivier, L. Novak, P. Chakroborty, H. Sharma, M. Chauhan, K. Kontolati, L. Vandanapu, D. Loukrezis, M. Gardner, Uqpy v4.1: Uncertainty quantification with python (2023).  arXiv:2305.09572.
  * [34] L. Lu, X. Meng, Z. Mao, G. E. Karniadakis, DeepXDE: A deep learning library for solving differential equations, SIAM Review 63 (1) (2021) 208–228.  doi:10.1137/19M1274067.
  * [35] A. Galetzka, D. Loukrezis, N. Georg, H. De Gersem, U. Römer, An hp-adaptive multi-element stochastic collocation method for surrogate modeling with information re-use, International Journal for Numerical Methods in Engineering 124 (12) (2023) 2902–2930.  doi:https://doi.org/10.1002/nme.7234.