∎

11institutetext: 1 HEDPS and LTCS, College of Engineering, Peking University,
Beijing 100871, China. 22institutetext: 2 Department of Mechanical
Engineering, Northwestern University, Evanston, USA. 3 Visiting student at
Department of Mechanical Engineering, Northwestern University 33institutetext:
* Corresponding Author<EMAIL_ADDRESS><EMAIL_ADDRESS>
# HiDeNN-PGD: reduced-order hierarchical deep learning neural networks

Lei Zhang1,3 Ye Lu2 Shaoqiang Tang1* Wing Kam Liu2*

(Received: date / Accepted: date)

###### Abstract

This paper presents a proper generalized decomposition (PGD) based reduced-
order model of hierarchical deep-learning neural networks (HiDeNN). The
proposed HiDeNN-PGD method keeps both advantages of HiDeNN and PGD methods.
The automatic mesh adaptivity makes the HiDeNN-PGD more accurate than the
finite element method (FEM) and conventional PGD, using a fraction of the FEM
degrees of freedom. The accuracy and convergence of the method have been
studied theoretically and numerically, with a comparison to different methods,
including FEM, PGD, HiDeNN and Deep Neural Networks. In addition, we
theoretically showed that the PGD converges to FEM at increasing modes, and
the PGD error is a direct sum of the FEM error and the mode reduction error.
The proposed HiDeNN-PGD performs high accuracy with orders of magnitude fewer
degrees of freedom, which shows a high potential to achieve fast computations
with a high level of accuracy for large-size engineering problems.

###### Keywords:

hierarchical deep-learning neural networks, proper generalized decomposition,
canonical tensor decomposition, finite element method, reduced order methods,
convergence study, error bound

##  1 Introduction

Despite the constantly increasing computer power, numerical simulations of
physical systems with numerous degrees of freedom remain computationally
prohibitive. These kinds of problems arise usually in simulation-based
engineering applications and the repetitive manipulation (or modification) of
the mesh system has been identified as a key time-costly issue in standard
finite element method (FEM) belytschko2013nonlinear . This has been a
motivation for developing the isogeometric approaches hughes2005isogeometric .

In recent years, the deep neural network (DNN) has shown some interesting
features in handling the solution of physics constrained systems. The
universal approximation theorem hornik1989multilayer ;
cybenko1989approximation and the natural scalability of DNN have been the
foundation of its superior performance for large systems. This has thus
motivated the use of DNN to approximate the solution of partial differential
equations (PDEs) weinan2017deep ; raissi2019physics ; li2019clustering . A
recently developed Hierarchical Deep-learning Neural Network (HiDeNN) method
zhang2021hierarchical ; saha2021hierarchical falls within this perspective.
The so-called HiDeNN is developed by constraining the weights and biases of
DNN to mesh coordinates to build multiple dimensions finite element, meshfree,
isogeometric, B-spline, and NURBs interpolation functions. HiDeNN allowed the
automatic mesh adaptivity and showed a good potential to prevent large mesh
systems and the standard time-consuming mesh refinement procedure. In order to
further enhance the efficiency of HiDeNN, this work proposed HiDeNN-PGD, a
reduced-order model of HiDeNN using the proper generalized decomposition
(PGD).

The PGD-based model reduction methods rely on the idea of separation of
variables, and are usually written in the format of canonical decomposition.
This kind of method was originally proposed in an a priori setting
ammar2006new ; chinesta2011short ; bhattacharyya2018multi ; modesto2015proper
, in which the separated functions are computed on-the-fly by solving the
PDEs. It has gained increased popularity in recent years. For overcoming the
intrusiveness and extending the applicability of the method, the a posteriori
data-driven PGD lu2018adaptive ; lu2019datadriven ; lu2018multi ;
diez2018algebraic has also been developed more recently. In contrast to a
priori PGD, the a posteriori method uses a database to learn the separated
functions and thus can be used as regression for constructing reduced order
surrogate models.

In our work, we adopted the same idea of separation of variables, in
particular, the separation of space variables for solving PDEs in the HiDeNN
framework, leading to the so-called HiDeNN-PGD method, which is expected to
have reduced degrees of freedom with a high accuracy. Indeed, the space
separated PGD, leading to lower dimensional space functions, is usually
considered for reducing the computational complexity of 3D separable domains
(see e.g. bognet2014separated ). However, the convergence aspect with respect
to the mesh refinement and number of modes has been less studied.

We investigated the convergence aspect of the PGD approach in this paper.
Based on the approximation function spaces, we analyzed the numerical error
and convergence associated with different approaches and compare their error
bounds. It can be shown that the HiDeNN-PGD is more accurate than both FEM and
conventional PGD, thanks to the adaptivity achieved by HiDeNN. For further
enhancing the optimality of PGD modes, we suggested fixing the number of modes
firstly and solve them together. Hence, HiDeNN-PGD can require fewer modes
than PGD. This is advantageous for high-dimensional problems where the
optimality of modes is crucial. The numerical examples have confirmed our
theoretical analysis. In addition, we numerically investigated the
relationship between the approximation error and the modes and proposed a
strategy to select the prescribed mode number in HiDeNN-PGD. The proposed
HiDeNN-PGD has shown a high potential to achieve high performance computing
with high accuracy.

The paper is organized as follows. Section 2 gives a brief overview of
different numerical methods for partial differential equations (PDEs), in
which the approximation function spaces are described. The error analysis
based on a class of PDEs is given in Section 3. Section 4 presents the
proposed HiDeNN-PGD method. Section 5 provides some numerical examples and
discussions. Finally, the paper closes with some concluding remarks.

##  2  DNN, HiDeNN, FEM and canonical decomposition based function
approximation

Function approximation is a key component in numerical solutions of partial
differential equations (PDEs). In this section, we briefly review how such an
approximation can be performed in terms of FEM, DNN, HiDeNN and Canonical
tensor Decomposition (CD). We present their approximation function sets, which
will be used in the theoretical analysis in Subsection 2.2. We restrict the
discussion to a scalar-value function
${u}(\bm{x}):\mathbb{R}^{3}\rightarrow\mathbb{R}$. The conclusions should be
straightforwardly extended to vector functions.

(a) DNN-based interpolation function

(b) HiDeNN interpolation function

(c) FEM interpolation function

(d) CD interpolation function

Figure 1: Illustration for the interpolation functions in the form of DNN with
$\bm{x}=(x,y,z)$ as input and $u^{h}$ as output. Weights and biases inside
dashed line-box are constrained, and those inside the red solid line-box are
fixed. $\mathcal{A}_{0}$ is the identity activation function defined by
$\mathcal{A}_{0}(x)=x$.

###  2.1 Overview of the approximation function sets

Deep neural network-based method

According to the universal approximation theorem, a deep neural network (DNN)
can be designed to approximate any given continuous function to desired
accuracy hornik1989multilayer ; cybenko1989approximation ; tang2021neural .
Thus it can be a candidate to approximate solutions for solving PDEs
wu2016physics ; xiao2016quantifying ; weinan2017deep ; berg2018unified ;
raissi2017physics ; raissi2019physics ; sirignano2018dgm ; weinan2018deep ,
i.e.,

$u^{h}(\bm{x})=\mathcal{F}^{NN}(\bm{x}),$ (1)

where $\mathcal{F}^{NN}$ represents the neural network with $\bm{x}$ as input
and $u^{h}$ as output. Note that $u^{h}$ can be a multidimensional vector. For
instance, in a classical feedforward neural network (FFNN) Goodfellow-et-
al-2016 ; haykin1994neural ; oishi2017computational with $N_{L}$ layers,
recursive relations among neurons are as follows

$\displaystyle a^{l}_{j=1}=x,a^{l}_{j=2}=y,a^{l}_{j=3}=z,\text{ if }l=1\text{
(input layer)};$ (2) $\displaystyle
a^{l}_{j}=\mathcal{A}(\sum_{i=1}^{N_{N}(l-1)}{W^{l}_{ij}a^{l-1}_{i}+b^{l}_{j}}),\text{
if }l\in\\{2,...,N_{L}-1\\}\text{ (hidden layer)}.$ (3)

Hence, the output layer can be defined as

$\displaystyle\mathcal{F}^{NN}_{j}$ $\displaystyle=$ $\displaystyle
a^{N_{L}}_{j}=\sum_{i=1}^{N_{N}(N_{L}-1)}{W^{N_{L}}_{ij}a^{N_{L}-1}_{i}+b^{N_{L}}_{j}},\text{
if }l=N_{L}\text{ (output layer)},$ (4)

with the detailed definition of the notations in Table 1. Therefore, once the
weights $\bm{W}$, biases $\bm{b}$ and activation functions $\bm{\mathcal{A}}$
have been chosen, $\mathcal{F}^{NN}$ can serve as an approximation function
with the input variable as $\bm{x}=(x,y,z)$.

Table 1: Notation table of variables used in the feed forward neural network $\bm{x}=(x,y,z)$ | Space coordinates
---|---
$l$ | Counting index for number of layers
$i$ |  Counting index for neurons in layer $l-1$
$j$ |  Counting index for neurons in layer $l$
$N_{L}$ | Number of layers in the neural network
$N_{N}(l)$ |  Number of neurons in layer $l$
$W_{ij}^{l}$ |  Weight connecting the $i^{\text{th}}$ neuron in layer $l-1$ to the $j^{\text{th}}$ in layer $l$
$b_{j}^{l}$ |  Bias of the $j^{\text{th}}$ neuron in layer $l$
$a_{j}^{l}$ |  Neuron value for $j^{th}$ neuron in $l^{th}$ layer
$\mathcal{A}$ | Activation function
$\mathcal{F}^{NN}$ | Feedforward neural network function

The approximation function set formed by a general NN is

$\mathcal{N}^{h}=\left\\{u^{h}(\bm{x})\big{|}u^{h}=\mathcal{F}^{NN}(\bm{x};\bm{W},\bm{b},\bm{\mathcal{A}}),W^{l}_{ij}\in\mathbb{R},b^{l}_{j}\in\mathbb{R}\right\\},$
(5)

where $\mathcal{F}^{NN}(\bm{x};\bm{W},\bm{b},\bm{\mathcal{A}})$ denotes an NN
with the input $\bm{x}$ and depends on weights $\bm{W}$, biases $\bm{b}$ and
activation functions $\bm{\mathcal{A}}$.

For interpretation, DNN can be somehow rewritten in the form of shape
functions associated to nodal values of $u^{h}$. In this way, the function
approximation reads

$u^{h}=\sum_{\mathcal{I}=1}^{np}\mathcal{F}^{NN}_{\mathcal{I}}(\bm{x};\bm{W},\bm{b},\bm{\mathcal{A}})\mathscr{u}_{\mathcal{I}},$
(6)

as illustrated in Fig. 1(a). $\mathcal{F}^{NN}_{\mathcal{I}}$ represents the
value of the $\mathcal{I}$-th neuron in the last hidden layer, i.e., the
output of the previous hidden layers. $\mathscr{u}_{\mathcal{I}}$ is the
corresponding weight connecting the output layer with the last hidden layer.
This interpolation form may provide a more interpretable structure for DNN.

In multidimensional cases, such as for 3D mechanical problems, the above
equation can be straightforwardly applied to each component of displacement
fied as follows

$u_{x}^{h}=\sum_{\mathcal{I}=1}^{np}\mathcal{F}^{NN}_{\mathcal{I}}(\bm{x};\bm{W},\bm{b},\bm{\mathcal{A}})\mathscr{u}_{x\mathcal{I}},\\\
$ (7)
$u_{y}^{h}=\sum_{\mathcal{I}=1}^{np}\mathcal{F}^{NN}_{\mathcal{I}}(\bm{x};\bm{W},\bm{b},\bm{\mathcal{A}})\mathscr{u}_{y\mathcal{I}},$
(8)
$u_{z}^{h}=\sum_{\mathcal{I}=1}^{np}\mathcal{F}^{NN}_{\mathcal{I}}(\bm{x};\bm{W},\bm{b},\bm{\mathcal{A}})\mathscr{u}_{z\mathcal{I}},$
(9)

HiDeNN

The recently developed HiDeNN method zhang2021hierarchical uses a similar DNN
structure of (6) with additional constraints to build a family of function
approximations. Similar to the FEM, the continuous domain $\Omega$ is
discretized by a mesh with $np$ nodes
$\bm{x}_{1},\bm{x}_{2},\cdots,\bm{x}_{np}$. Then the finite element shape
functions $N_{\mathcal{I}}$ can be constructed by the neural network block,
namely,

$\mathcal{N}_{\mathcal{I}}(\bm{x};\bm{W},\bm{b},\bm{\mathcal{A}})$ (10)

where $\mathcal{I}=1,2,\cdots,np$. Different from
$\mathcal{F}_{\mathcal{I}}^{NN}(\bm{x};\bm{W},\bm{b},\bm{\mathcal{A}})$ in
(6), $\mathcal{N}_{\mathcal{I}}(\bm{x};\bm{W},\bm{b},\bm{\mathcal{A}})$
precisely equals to the finite element shape function
$N_{\mathcal{I}}(\bm{x})$ with inputs $\bm{x}$ and an output
$N_{\mathcal{I}}$, satisfying the following constraints for shape functions
automatically, i.e.

$\displaystyle\sum_{\mathcal{I}=1}^{np}\mathcal{N}_{\mathcal{I}}(\bm{x};\bm{W},\bm{b},\bm{\mathcal{A}})=1,\mathcal{N}_{\mathcal{I}}(\bm{x}_{\mathcal{J}};\bm{W},\bm{b},\bm{\mathcal{A}})=\delta_{\mathcal{I}\mathcal{J}}.$
(11)

With Kronecker Delta constraints, we can apply Dirichlet boundary conditions
directly similar to that of the finite element method, so that all the weights
$\bm{W}$, and biases $\bm{b}$ are functions of nodal coordinates $\bm{x}_{I}$.
Thus we can rewrite the shape function explicitly in terms of $\bm{x}_{I}$ as

$\mathcal{N}_{\mathcal{I}}(\bm{x};\bm{x}_{\mathcal{I}}^{*},\bm{\mathcal{A}})=\mathcal{N}_{\mathcal{I}}(\bm{x};\bm{W},\bm{b},\bm{\mathcal{A}}),$
(12)

where $\bm{x}_{\mathcal{I}}^{*}$ denotes the support of
$N_{\mathcal{I}}(\bm{x})$, e.g. in linear 1D cases
$\bm{x}_{\mathcal{I}}^{*}=[x_{{\mathcal{I}}-1},x_{\mathcal{I}},x_{{\mathcal{I}}+1}]$.

Combining such neural network blocks for the entire mesh gives the final form
of HiDeNN, as shown in Fig. 8(b). This results in the approximation function
set

$\mathcal{H}^{h}=\left\\{u^{h}(\bm{x})\bigg{|}u^{h}=\sum_{\mathcal{I}=1}^{np}\mathcal{N}_{\mathcal{I}}(\bm{x};\bm{x}_{\mathcal{I}}^{*},\bm{\mathcal{A}})\mathscr{u}_{\mathcal{I}},\mathscr{u}_{\mathcal{I}}\in\mathbb{R}\right\\}.$
(13)

The parametric expression with nodal positions $\bm{x}_{\mathcal{I}}^{*}$
allows automatic r-adaptivity, and accordingly improves the local and global
accuracy of the interpolant.

Finite element method

The approximation function set $\mathcal{H}^{h}$ degenerates to the FE
approximation function set $\mathcal{V}^{h}$ belytschko2013nonlinear when the
nodal position is fixed, which reads

$\mathcal{V}^{h}=\left\\{u^{h}(\bm{x})\bigg{|}u^{h}=\sum_{\mathcal{I}=1}^{np}N_{\mathcal{I}}(\bm{x})u_{\mathcal{I}},u_{\mathcal{I}}\in\mathbb{R}\right\\}.$
(14)

From the DNN viewpoint, this corresponds to fixing the weights and biases that
are functions of $\bm{x}_{\mathcal{I}}^{*}$, as shown in Fig. 1(b).

Canonical tensor decomposition

Under the assumption of separation of variables, the function $\bm{u}$ may be
approximated by the sum of the products of multiple 1D functions, i.e.,

$u^{h}(\bm{x})=u^{h}(x,y,z)=\sum_{q=1}^{Q}X^{(q)}(x)Y^{(q)}(y)Z^{(q)}(z),$
(15)

where $Q$ is the number of modes, and the product of $X^{(q)},Y^{(q)},Z^{(q)}$
provides a mode for the interpolation function. This form or concept is known
as CD kolda2009tensor . The so-called PGD method has adopted this concept for
solving PDEs ammar2006new ; gonzalez2010recent and for data learning
lu2018adaptive ; lu2019datadriven ; blal2019non .

Thanks to the separation of variables, only shape functions of reduced
dimension are needed. In (15), 1D FE shape functions can be used for a 3D
problem, namely,

$\displaystyle X^{(q)}(x)$ $\displaystyle=$
$\displaystyle\sum_{I=1}^{n_{1}}N_{I}(x)\beta_{I}^{(q)},$ (16) $\displaystyle
Y^{(q)}(y)$ $\displaystyle=$
$\displaystyle\sum_{J=1}^{n_{2}}N_{J}(y)\gamma_{J}^{(q)},$ (17) $\displaystyle
Z^{(q)}(z)$ $\displaystyle=$
$\displaystyle\sum_{K=1}^{n_{3}}N_{K}(z)\theta_{K}^{(q)}.$ (18)

Here, $n_{1},n_{2},n_{3}$ are the number of nodes in $x,y,z$ directions,
respectively. Thus the corresponding approximation function set for a given
$Q$ modes is

$\displaystyle\mathcal{M}_{Q}^{h}=$ $\displaystyle\Bigg{\\{}$ $\displaystyle
u^{h}(\bm{x})\bigg{|}u^{h}=\sum_{q=1}^{Q}\left(\sum_{I=1}^{n_{1}}N_{I}(x)\beta_{I}^{(q)}\right)\left(\sum_{J=1}^{n_{2}}N_{J}(y)\gamma_{J}^{(q)}\right)\left(\sum_{K=1}^{n_{3}}N_{K}(z)\theta_{K}^{(q)}\right),$
$\displaystyle\beta_{I}^{(q)},\gamma_{J}^{(q)},\theta_{K}^{(q)}\in\mathbb{R}\Bigg{\\}}.$

Note that the interpolation function in (15) can be rearranged as

$u^{h}(x,y,z)=\sum_{I=1}^{n_{1}}\sum_{J=1}^{n_{2}}\sum_{K=1}^{n_{3}}N_{I}(x)N_{J}(y)N_{K}(z)\left(\sum_{q=1}^{Q}\beta^{(q)}_{I}\gamma^{(q)}_{J}\theta^{(q)}_{K}\right)$
(20)

which is regarded as a finite element interpolation function with
$N_{I}(x)N_{J}(y)N_{K}(z)$ as shape functions and
$\left(\sum_{q=1}^{Q}\beta^{(q)}_{I}\gamma^{(q)}_{J}\theta^{(q)}_{K}\right)$
as coefficients. Fig. 1(d) illustrates a DNN format of (20), which will be the
basis for the proposed HiDeNN-PGD method.

Multidimensional shape functions of CD, i.e., the product of 1D shape
functions, are fixed and determined by nodal positions along each direction.
In addition, the nodal values in the last layer are constraint in the form of
the tensor product, i.e.,

$u_{(I,J,K)}=\sum_{q=1}^{Q}\beta^{(q)}_{I}\gamma^{(q)}_{J}\theta^{(q)}_{K}.$
(21)

When a few modes may represent the function $u^{h}(x,y,z)$, this method is
advantageous in terms of lower integration complexity and less degrees of
freedom (DoFs) ammar2006new ; chinesta2013proper . The DoFs in (20) is of the
order $O((n_{1}+n_{2}+n_{3})Q)$, which is linear with the spatial dimension
and far smaller than traditional methods (e.g., FEM).

###  2.2 HiDeNN-PGD: Reduced order HiDeNN via PGD

HiDeNN gets better accuracy compared with classical FEM due to the adaptivity
of nodal positions. More DoFs might result in more cost. On the other hand,
representation of separated variables provides a reduced order model to
improve the efficiency but might lose accuracy. Here, we propose HiDeNN-PGD, a
reduced-order model of HiDeNN via PGD, which seeks to accomplish an optimized
balance between accuracy and computational cost.

In HiDeNN-PGD, the shape functions in each direction are written in the DNN
format, namely, (16-18) are replaced by 1D HiDeNN interpolants (refer to
Appendix C),

$\displaystyle X^{(q)}(x)$ $\displaystyle=$
$\displaystyle\sum_{I=1}^{n_{1}}\mathcal{N}_{I}(x;\bm{x}_{I}^{*},\bm{\mathcal{A}})\beta_{I}^{(q)},$
(22) $\displaystyle Y^{(q)}(y)$ $\displaystyle=$
$\displaystyle\sum_{J=1}^{n_{2}}\mathcal{N}_{J}(y;\bm{y}_{J}^{*},\bm{\mathcal{A}})\gamma_{J}^{(q)},$
(23) $\displaystyle Z^{(q)}(z)$ $\displaystyle=$
$\displaystyle\sum_{K=1}^{n_{3}}\mathcal{N}_{K}(z;\bm{z}_{K}^{*},\bm{\mathcal{A}})\theta_{K}^{(q)},$
(24)

where
$\mathcal{N}_{I}(x;\bm{x}_{I}^{*},\bm{\mathcal{A}}),\mathcal{N}_{J}(y;\bm{y}_{J}^{*},\bm{\mathcal{A}}),\mathcal{N}_{K}(z;\bm{z}_{K}^{*},\bm{\mathcal{A}})$
are the 1D HiDeNN shape functions in $x,y,z$ directions, respectively. Thus
the interpolation function set is defined by

$\displaystyle\mathcal{G}_{Q}^{h}=$ (25)
$\displaystyle\Bigg{\\{}u^{h}\bigg{|}u^{h}=\sum_{q=1}^{Q}\left(\sum_{I=1}^{n_{1}}\mathcal{N}_{I}(x;\bm{x}_{I}^{*},\bm{\mathcal{A}})\beta_{I}^{(q)}\right)\left(\sum_{J=1}^{n_{2}}\mathcal{N}_{J}(y;\bm{y}_{J}^{*},\bm{\mathcal{A}})\gamma_{J}^{(q)}\right)\left(\sum_{K=1}^{n_{3}}\mathcal{N}_{K}(z;\bm{z}_{K}^{*},\bm{\mathcal{A}})\theta_{K}^{(q)}\right)\Bigg{\\}}.$

Since adaptivity occurs only in each direction, the mesh is always regular.

###  2.3 Relationship among the NN, HiDeNN, FEM, CD and HiDeNN-PGD
approximation function sets

In this subsection, we explore the relationship among FEM, NN, HiDeNN, CD and
HiDeNN-PGD.

Assume that CD, FEM are based on the same regular mesh with
$n_{1},n_{2},n_{3}$ nodes along $x,y,z$ directions. This mesh also serves as
an initial guess of $\bm{x}_{I}^{*}$ in HiDeNN and HiDeNN-PGD. The shape
functions of FEM is the product of 1D shape functions, i.e., the shape
function associated with the node $(x_{i},y_{j},z_{k})$ is
$N_{(I,J,K)}(x,y,z)=N_{I}(x)N_{J}(y)N_{K}(z)$. NN has a more general structure
than HiDeNN, and might be fully-connected.

By definition, we have the following relationship among approximation function
sets of NN, HiDeNN, FEM and CD:

$\mathcal{M}_{Q}^{h}\subset\mathcal{V}^{h}\subset\mathcal{HI}^{h}\subset\mathcal{N}^{h}.$
(26)

Especially when $Q$ is big enough ($Q\geq\min\\{n_{1},n_{2}\\}$ for 2D and
$Q\geq\min\\{n_{1}n_{2},n_{2}n_{3},n_{3}n_{1}\\}$ for 3D), we have

$\mathcal{M}^{h}_{Q}=\mathcal{V}^{h}\subset\mathcal{G}^{h}_{Q}\subset\mathcal{H}^{h}\subset\mathcal{N}^{h},$
(27)

as illustrated in Fig. 2.

Figure 2: Illustration for relationship among interpolation function sets of
CD, FEM, HiDeNN-PGD, HiDeNN and NN. Especially, when $Q$ tends to infinity,
$\mathcal{M}_{Q}^{h}$ approaches to $\mathcal{V}^{h}$. When $Q$ is big enough,
FEM interpolation set $\mathcal{V}^{h}$ is the subset of HiDeNN-PGD one
$\mathcal{G}^{h}_{Q}$.

The above conclusions (26)-(27) are based on the following observations:

$\bullet$

According to (20), the CD interpolation functions can be regarded as finite
element interpolation functions, which belong to $\mathcal{V}^{h}$, so
$\mathcal{M}_{Q}^{h}\subset\mathcal{V}^{h}$. Especially, when $Q$ is big
enough ($Q\geq\min\\{n_{1},n_{2}\\}$ for 2D and
$Q\geq\min\\{n_{1}n_{2},n_{2}n_{3},n_{3}n_{1}\\}$ for 3D),
$\mathcal{M}_{Q}^{h}$ approaches to $\mathcal{V}^{h}$, i.e.,
$\mathcal{M}^{h}_{Q}=\mathcal{V}^{h}$. Detailed results will be shown in
Section 4. Proofs can be found in Appendix B.

$\bullet$

In HiDeNN, an optimization of the nodal positions is performed. Thus FEM may
be regarded as a specific case in the HiDeNN with nodal coordinates fixed .

$\bullet$

HiDeNN is a class of structured NN with weights and biases as functions of
nodal values and nodal positions.

$\bullet$

HiDeNN-PGD requires the mesh regular, while HiDeNN optimizes nodal positions
freely, so $\mathcal{G}^{h}_{Q}\subset\mathcal{H}^{h},\forall Q\in\mathbb{N}$.
On the other hand, HiDeNN-PGD has more DoFs than CD under the same number of
modes, so $\mathcal{M}^{h}_{Q}\subset\mathcal{G}^{h}_{Q},\forall
Q\in\mathbb{N}$. When $Q$ is small, $\mathcal{M}^{h}_{Q}$ is the subset of the
intersection of $\mathcal{G}^{h}_{Q}$ and $\mathcal{V}^{h}$.

We can also summarize the DoFs for different methods in Table 2. It is shown
that HiDeNN-PGD and PGD have only a linear growth in terms of DoFs, whereas
the DoFs of FEM and HiDeNN may grow in a polynomial manner.

Table 2: Comparison of DoFs for different methods on a 3D mesh | FEM | PGD/CD | HiDeNN-PGD | HiDeNN
---|---|---|---|---
DoFs | $n_{1}\times n_{2}\times n_{3}$ | $(n_{1}+n_{2}+n_{3})\times Q$ | $(n_{1}+n_{2}+n_{3})\times Q+n_{1}+n_{2}+n_{3}$ | $n_{1}\times n_{2}\times n_{3}+n_{1}\times n_{2}\times n_{3}$

##  3 Error analysis of FEM, CD, HiDeNN, NN-based solutions and HiDeNN-PGD for
PDEs

We consider a partial differential equation with homogeneous boundary
conditions

$\left\\{\begin{array}[]{l}\mathcal{L}\bm{u}(\bm{x})+\bm{b}(\bm{x})=0\text{ in
}\Omega\subset\mathbb{R}^{d},\\\
\bm{u}|_{\partial\Omega}=\bm{0},\end{array}\right.$ (28)

where $\bm{u}\in\mathbb{R}^{m}$ denotes an $m$-dimensional vector-valued
function in a certain Hilbert space $H(\mathbb{R}^{d},\mathbb{R}^{m})$,
$\bm{b}$ the source term, $\bm{x}\in\mathbb{R}^{d}$ the $d$-dimensional space
coordinates, and $\mathcal{L}$ a second-order differential operator.

We assume that an energy potential $\Pi(u)$ exists, formulated in the
following form

$\Pi(\bm{u})=\dfrac{1}{2}a(\bm{u},\bm{u})-(\bm{b},\bm{u}),$ (29)

where $a(\cdot,\cdot)$ is the symmetric and bilinear form corresponding to the
second-order differential operator $\mathcal{L}$, and
$(\bm{f},\bm{g})=\int_{\Omega}\bm{f}\cdot\bm{g}\mathrm{d}\bm{x}$ denotes the
inner product. For example, let
$a(\bm{f},\bm{g})=\int_{\Omega}\nabla\bm{f}\cdot\nabla\bm{g}\mathrm{d}\bm{x}$
for Poisson equation. The minimization of $\Pi(u)$ gives the solution to (28).
This leads to a weak form of (28), which reads

$\bm{u}=\operatorname*{arg\,min}_{\bm{u}^{*}\in
H(\mathbb{R}^{m},\mathbb{R}^{n})}\Pi[\bm{u^{*}}].$ (30)

Such a weak form is commonly adopted in interpolation theory based numerical
approaches, such as the methods shown in Section 2. Denoted by
$\mathcal{S}^{h}\subset H(\mathbb{R}^{d},\mathbb{R}^{m})$ the discretized
approximation solution set with a characteristic mesh size $h$, the
approximate solution based on this given interpolation function set is then

$\bm{u}^{h}=\operatorname*{arg\,min}_{\bm{u}^{h*}\in\mathcal{S}^{h}}\Pi[\bm{u^{h*}}].$
(31)

In the following, we shall take $\mathcal{S}^{h}$ to be
$\mathcal{M}_{0,Q}^{h},\mathcal{V}_{0}^{h},\mathcal{H}_{0}^{h},\mathcal{N}_{0}^{h},\mathcal{G}_{0,Q}^{h}$,
which are the subsets of
$\mathcal{M}_{Q}^{h},\mathcal{V}^{h},\mathcal{H}^{h},\mathcal{N}^{h},\mathcal{G}_{Q}^{h}$
under homogeneous boundary conditions, respectively.

###  3.1 Error analysis

We assert the following relations among error bounds for FEM, NN, HiDeNN and
CD

$\left\|u^{CD}-u^{exact}\right\|_{E}\geq\left\|u^{FEM}-u^{exact}\right\|_{E}\geq\left\|u^{HiDeNN}-u^{exact}\right\|_{E}\geq\left\|u^{NN}-u^{exact}\right\|_{E},$
(32)

Here, $\left\|\cdot\right\|_{E}=\sqrt{a(\cdot,\cdot)}$ is called as the energy
norm, and $u^{exact}$ is the real solution of the problem, i.e.,

$u^{exact}=\arg\min_{u\in H^{1}_{0}}\Pi[u].$ (33)

Especially when $Q$ is big enough ($Q\geq\min\\{n_{1},n_{2}\\}$ for 2D and
$Q\geq\min\\{n_{1}n_{2},n_{2}n_{3},n_{3}n_{1}\\}$ for 3D), the error bounds
for five methods become

$\displaystyle\left\|u^{CD}-u^{exact}\right\|_{E}\geq\left\|u^{FEM}-u^{exact}\right\|_{E}\geq\left\|u^{HiDeNN-
PGD}-u^{exact}\right\|_{E}$ (34)
$\displaystyle\geq\left\|u^{HiDeNN}-u^{exact}\right\|_{E}\geq\left\|u^{NN}-u^{exact}\right\|_{E},$

The relationship (26) is inherited by
$\mathcal{M}_{0,Q}^{h},\mathcal{V}_{0}^{h},\mathcal{H}_{0}^{h},\mathcal{N}_{0}^{h}$,
i.e.,

$\mathcal{M}_{0,Q}^{h}\subset\mathcal{V}_{0}^{h}\subset\mathcal{H}_{0}^{h}\subset\mathcal{N}_{0}^{h}.$
(35)

These four methods are all based on the minimal energy principle,

$\bm{u}=\arg\min_{\bm{u}^{h}\in\mathcal{S}^{h}}\Pi[\bm{u}^{h}],$ (36)

where $\mathcal{S}^{h}$ is selected as
$\mathcal{M}_{0,Q}^{h},\mathcal{V}_{0}^{h},\mathcal{H}_{0}^{h},\mathcal{N}_{0}^{h}$
for CD, FEM, HiDeNN, and NN, respectively. Thus due to the relationship (35)
among them, we have

$\Pi[u^{CD}]\geq\Pi[u^{FEM}]\geq\Pi[u^{HiDeNN}]\geq\Pi[u^{NN}].$ (37)

This leads to the error bounds (32). In the same manner, the relationship (27)
leads to the error bounds (3.1). (3.1) shows that HiDeNN-PGD might reach
better accuracy than FEM with regular mesh at increasing number of modes.

We remark that the above theoretical analysis does not account for numerical
aspects, such as the difficulties to find out the global minimizer of the
energy potential.

###  3.2 Proof of the mesh in-dependency of mode reduction error in the PGD
method

The PGD based model reduction induces two kinds of errors: mesh discretization
error and mode reduction error.

###### Theorem 3.1

Let $u^{exact},u^{PGD}$,$u^{FEM}$ be the exact solution, the numerical
solution of PGD and FEM, respectively. PGD and FEM take the same regular mesh,
and FEM takes the shape functions as the product of 1D shape functions of PGD
in each dimension. Then the following error decomposition holds:

$\|u^{PGD}-u^{exact}\|_{E}^{2}=\|u^{FEM}-u^{exact}\|_{E}^{2}+\|u^{PGD}-u^{FEM}\|_{E}^{2}.$
(38)

###### Proof

We calculate

$\displaystyle\|u^{PGD}-u^{exact}\|_{E}^{2}$ (39) $\displaystyle=$
$\displaystyle\int_{\Omega}(\nabla u^{PGD}-\nabla u^{FEM}+\nabla
u^{FEM}-\nabla u^{exact})^{2}\mathrm{d}\bm{x}$ $\displaystyle=$
$\displaystyle\|u^{FEM}-u^{exact}\|_{E}^{2}+\|u^{PGD}-u^{FEM}\|_{E}^{2}+\int_{\Omega}2(\nabla
u^{PGD}-\nabla u^{FEM})\cdot(\nabla u^{FEM}-\nabla
u^{exact})\mathrm{d}\bm{x}.$

By the Gauss theorem, we obtain

$\displaystyle\int_{\Omega}(\nabla u^{PGD}-\nabla u^{FEM})\cdot(\nabla
u^{FEM}-\nabla u^{exact})\mathrm{d}\bm{x}$ $\displaystyle=$
$\displaystyle\int_{\Omega}(\nabla u^{PGD}-\nabla u^{FEM})\cdot\nabla
u^{FEM}+(u^{PGD}-u^{FEM})\nabla^{2}u^{exact}\mathrm{d}\bm{x}$ $\displaystyle=$
$\displaystyle\int_{\Omega}(\nabla u^{PGD}-\nabla u^{FEM})\cdot\nabla
u^{FEM}-(u^{PGD}-u^{FEM})b(\bm{x})\mathrm{d}\bm{x}.$

Since PGD and FEM share the same mesh and shape functions, $v=u^{PGD}-u^{FEM}$
belongs to the test function space of FEM. By the weak form of the FEM
problem, we have

$\int_{\Omega}\nabla v\cdot\nabla
u^{FEM}-v(\bm{x})b(\bm{x})\mathrm{d}\bm{x}=0.$ (42)

That is to say, (3.2) vanishes and hence (38) holds. $\hfill\square$.

This theorem asserts that the PGD error is a direct sum of the FEM error and
the mode reduction error (the difference between FEM and PGD).

##  4 The formulation of HiDeNN-PGD: the 2D Poisson problem as illustration

In subsection 2.2, we defined HiDeNN-PGD in terms of the approximation
function space. Here, we give the detailed formulation of the method.

For the sake of simplicity and without loss of generality, we consider 2D
Poisson problem,

$\left\\{\begin{array}[]{l}\Delta u(x,y)+b(x,y)=0\text{ in
}\Omega_{(x,y)}\subset\mathbb{R}^{2},\\\
u|_{\partial\Omega}=0,\end{array}\right.$ (43)

(43) is solved in the regular domain $\Omega_{(x,y)}=[a,b]\times[c,d]$ with
homogeneous boundary conditions. Note that if inhomogeneous boundary
conditions are under consideration, we can separate the solution into two
parts,

$u=u^{0}+\tilde{u},$ (44)

where $u^{0}$ is an arbitrary function satisfying boundary conditions, and
$\tilde{u}$ is the solution to the new equations with homogeneous boundary
condition. The variational formula of (43) is

$\Pi(u)=\dfrac{1}{2}\int_{\Omega_{(x,y)}}\left|\nabla
u\right|^{2}\mathrm{d}x\mathrm{d}y-\int_{\Omega_{(x,y)}}u(x,y)b(x,y)\mathrm{d}x\mathrm{d}y$
(45)

Substituting HiDeNN-PGD interpolation function into (45), we obtain

$\displaystyle\Pi(u^{h})$ $\displaystyle=$
$\displaystyle\dfrac{1}{2}\sum_{p=1}^{Q}\sum_{q=1}^{Q}\left(\int_{x_{1}}^{x_{n_{1}}}\dfrac{d}{dx}X^{(p)}(x)\dfrac{d}{dx}X^{(q)}(x)\mathrm{d}x\right)\left(\int_{y_{1}}^{y_{n_{2}}}Y^{(p)}(y)Y^{(q)}(y)\mathrm{d}y\right)$
$\displaystyle+$
$\displaystyle\dfrac{1}{2}\sum_{p=1}^{Q}\sum_{q=1}^{Q}\left(\int_{x_{1}}^{x_{n_{1}}}X^{(p)}(x)X^{(q)}(x)\mathrm{d}x\right)\left(\int_{y_{1}}^{y_{n_{2}}}\dfrac{d}{dy}Y^{(p)}(y)\dfrac{d}{dy}Y^{(q)}(y)\mathrm{d}y\right)$
$\displaystyle-$
$\displaystyle\sum_{q=1}^{Q}\int_{\Omega_{(x,y)}}X^{(q)}(x)Y^{(q)}(y)b(x,y)\mathrm{d}x\mathrm{d}y,$

with the discrete mesh
$[x_{1}=a,x_{2},\cdots,x_{n_{1}}=b]\times[y_{1}=c,y_{2},\cdots,y_{n_{2}}=d]$.
Notice that there exist cross terms in Eq. (4). For convenience and
considering the difficulties to do exact integration between different
discrete meshes, all the modes share the same mesh
$[x_{1},x_{2},\cdots,x_{n_{1}}]\times[y_{1},y_{2},\cdots,y_{n_{2}}]$ and the
same shape functions. We use Gauss quadrature for the source term.

Once the interpolation function set (25) is obtained, minimal variational
principle gives the approximated solution. The process in HiDeNN-PGD is
formulated as

find
$\displaystyle\beta_{I}^{(1)},\gamma_{J}^{(1)},\cdots,\beta_{I}^{(Q)},\gamma_{J}^{(Q)},x_{I},y_{J},I=1,\cdots,n_{1},J=1,\cdots,n_{2}$
(47) min $\displaystyle\dfrac{1}{2}\int_{\Omega_{(x,y)}}\left|\nabla
u\right|^{2}\mathrm{d}x\mathrm{d}y-\int_{\Omega_{(x,y)}}u(x,y)b(x,y)\mathrm{d}x\mathrm{d}y$
$\displaystyle
u^{h}=\sum_{q=1}^{Q}\left(\sum_{I=1}^{n_{1}}\mathcal{N}_{I}(x;\bm{x}_{I}^{*},\bm{\mathcal{A}})\beta_{I}^{(q)}\right)\left(\sum_{J=1}^{n_{2}}\mathcal{N}_{J}(y;\bm{y}_{J}^{*},\bm{\mathcal{A}})\gamma_{J}^{(q)}\right)$
and
$\displaystyle\sum^{n_{1}}_{I=1}\mathcal{N}(\bm{x}^{*}_{I},x,\mathcal{A})=1,\sum^{n_{2}}_{J=1}\mathcal{N}(\bm{y}^{*}_{J},y,\mathcal{A})=1.$

The gradient descent method is applied to iteratively minimize
$\Pi(\mathscr{u}^{h})$ and solve for all parameters together. In the following
numerical examples, we choose Adam algorithm kingma2014adam , i.e.,

  1. 1.

Initialization: Set number of modes $Q$, initial nodal positions
$x_{I},y_{J},I=1,2,\cdots,n_{1},J=1,2,\cdots,n_{2}$, initial coefficients
$\beta^{(q)}_{I},\gamma^{(q)}_{J},q=1,2,\cdots,Q$ and maximal iteration step
$M$

  2. 2.

Algorithm:

While $k\leq M$ do

    1. (a)

Compute gradient $\dfrac{\partial}{\partial x_{I}},\dfrac{\partial}{\partial
y_{J}},\dfrac{\partial}{\partial\beta_{I}^{(q)}},\dfrac{\partial}{\partial\gamma_{J}^{(q)}},I=1,\cdots,n_{1},J=1,\cdots,n_{2},q=1,\cdots,Q$;

    2. (b)

Update
$x_{I},y_{J},\beta_{I}^{(q)},\gamma_{J}^{(q)},I=1,\cdots,n_{1},J=1,\cdots,n_{2},q=1,\cdots,Q$
by using Adam algorithm;

End while.

##  5 Numerical examples

In this section, we study the performance of HiDeNN-PGD with comparison to
FEM, HiDeNN and PGD methods. This study mainly focuses on the accuracy
comparison and the convergence behavior, as the computational cost can be
strongly affected by the under-optimized implementation. The computational
efficiency of the proposed HiDeNN-PGD method will be investigated in our
future work on the basis of GPU-computing.

###  5.1 2D case

The HiDeNN-PGD is applied to the Poisson problem with a concentrated load. To
demonstrate the capability of the method, the domain analyzed by HiDeNN-PGD is
initialized with a uniform mesh of 40 by 40 elements. As shown in Figure 3,
the final solution agrees with the reference FE solution. This reference
solution is obtained with a very fine mesh containing $4,000\times 4,000$
elements.

Table 3 illustrates the evolution of accuracy with an increasing number of
modes. For comparison purposes, we also applied the PGD, CD, FEM and HiDeNN
for the same problem on the uniform mesh of 40 by 40. The errors of different
methods are computed based on the energy norm with comparison to the reference
solution. As expected, HiDeNN is the most accurate one but leads to a
significantly larger number of degrees of freedom (DoFs). The proposed HiDeNN-
PGD method can have the same level of accuracy and only requires a small
number of modes. Compared to PGD and CD, the HiDeNN-PGD is much more accurate
at a limited number of modes. Taking a closer look at PGD and CD, these two
methods converge overall to the FEM on the coarse mesh. HiDeNN-PGD and HiDeNN
can overcome this limitation imposed by the mesh size with the adaptivity.
This observation is consistent with our theoretical analysis. Moreover, it
should be noticed that, unlike HiDeNN, the HiDeNN-PGD only increases slightly
the DoFs when compared with PGD and CD. This attributes to the separation of
variables. Indeed, the mesh adaptation is performed only in the two separated
axes, as shown in Figure 3(d). For comparison purposes, the final optimized
mesh of HiDeNN is illustrated in Figure 4. It is shown that the HiDeNN enables
a full adaptivity for the entire mesh, which leads to a significantly
different and more accurate final result. Nevertheless, the HiDeNN-PGD shows
attractive advantages in terms of DoFs.

Figure 5 illustrates the first four modes of HiDeNN-PGD, PGD, and CD. The PGD
modes remain similar to CD in this example. However, the modes of HiDeNN-PGD
seem to be more concentrated in the region of interest. This difference mainly
comes from the mesh adaptivity.

To further confirm the performance of HiDeNN-PGD, we have compared the
accuracy of different methods on different meshes. In Table 4, the PGD, CD and
HiDeNN-PGD results are obtained from the final converged mode. It is shown
that the HiDeNN-PGD always gives more accurate results with fewer degrees of
freedom. This confirms the previous observation.

(a) Reference FE solution obtained with a extremely fine mesh

(b) HiDeNN-PGD solution

(c) Initial uniform mesh for HiDeNN-PGD

(d) Final optimized mesh for HiDeNN-PGD

Figure 3: FE solution versus HiDeNN-PGD solution.

(a) Inital mesh for HiDeNN

(b) Final optimized mesh for HiDeNN

Figure 4: Mesh optimization in HiDeNN.

(a) First four modes in HiDeNN-PGD

(b) First four modes in PGD

(c) First four modes in CD

Figure 5: Mode comparison for HiDeNN-PGD, PGD and CD Table 3: Accuracy comparison for different methods on the coarse 40 by 40 mesh | PGD | CD | FEM | HiDeNN-PGD | HiDeNN
---|---|---|---|---|---
Mode number | DoFs | Err | DoFs | Err | DoFs | Err | DoFs | Err | DoFs | Err
1 | 78 | $38.167\%$ | 78 | $38.167\%$ | 1521 | $\bm{11.659\%}$ | 156 | $37.357\%$ | 4719 | $\bm{2.102\%}$
2 | 156 | $16.500\%$ | 156 | $14.422\%$ |  | - | 234 | $\bm{9.293\%}$ |  | -
3 | 234 | $13.188\%$ | 234 | $11.789\%$ |  | - | 312 | $3.674\%$ |  | -
4 | 312 | $11.811\%$ | 312 | $11.664\%$ |  | - | 390 | $3.662\%$ |  | -
5 | 390 | $11.685\%$ | 390 | $\bm{11.659\%}$ |  | - | 468 | $3.661\%$ |  | -
6 | 468 | $11.666\%$ | 468 | $11.659\%$ |  | - | 546 | $3.661\%$ |  | -
8 | 624 | $\bm{11.659\%}$ |  | - |  | - |  | - |  | -
20 | 1560 | $11.659\%$ |  | - |  | - |  | - |  | -
Table 4: Accuracy comparison for different methods with different meshes

| PGD | CD | FEM | HiDeNN-PGD | HiDeNN
---|---|---|---|---|---
Mesh | DoFs | Err | DoFs | Err | DoFs | Err | DoFs | Err | DoFs | Err
$40\times 40$ | 624(8) | $11.659\%$ | 390(5) | $11.659\%$ | 1,521 | $\bm{11.659\%}$ |  468(5) | $\bm{3.661\%}$ | 4719 | $\bm{2.102\%}$
$80\times 80$ | 1,422(9) | $5.887\%$ | 790(5) | $5.887\%$ | 6,241 | $5.887\%$ | 1,106(6) | $1.851\%$ | 19,039 | $1.406\%$
$160\times 160$ | 2,862(9) | $2.948\%$ | 1,908(6) | $2.948\%$ | 25,281 | $2.948\%$ | 2,226(6) | $1.174\%$ | 76,479 | $1.081\%$
$320\times 320$ | 7,018(11) | $1.469\%$ | 5,104(8) | $1.469\%$ | 101,761 | $1.469\%$ | 3,828(5) | $0.896\%$ | 306,559 | $0.889\%$
$640\times 640$ | 11,502(9) | $0.724\%$ | 7,668(6) | $0.724\%$ | 408,321 | $0.724\%$ | 10,224(7) | $0.606\%$ | 1,227,519 | $0.597\%$

Note that the number in the parentheses indicates the number of modes.

####  5.1.1 Convergence studies

In the HiDeNN-PGD method, the mode number has to be prescribed. In general,
this is unknown for a given problem and can vary significantly from one to
another. Thus, we want to study the convergence property of this method and
the PGD to get a general idea about how to choose the mode number.

As a first attempt, we restrict ourselves to a one-mode solution problem. This
eliminates the effect of the number of modes, and allows us to study the
convergence rate of the methods with respect to mesh refinement. This kind of
error is usually known as discretization error in FEM. To do so, the body
force term is manufactured so that the final solution is analytically known in
a separated form. As shown in Figure 6, the PGD and HiDeNN-PGD results
converge with respect to the element size at a rate similar to FEM. However,
if we consider the DoFs, this error converges much faster for PGD and HiDeNN-
PGD. This confirms that doing separation of variables for PGD and HiDeNN-PGD
does not degrade the convergence rate in terms of mesh refinement.

The PGD based model reduction induces two kinds of errors: mesh discretization
error and mode reduction error. In particular, we have theoretically shown
that the latter one should be independent of the mesh. In order to show this
numerically, we use a manufactured load to compute the multi-modes PGD
solution in different meshes. The results are shown in Figure 7. As expected,
the convergence rates on the number of modes remain similar regardless of the
mesh size. It seems that the log error is linearly proportional to the mode
number. The decreasing slop remained unchanged from the very coarse mesh to
the fine one. This implies the coarse mesh has the same mode reduction error
as a fine mesh. Hence, it can be used to choose the mode number.

From the above observation, we may consider using a coarse mesh PGD, which is
very cheap, to study the mode reduction error and choose an appropriate mode
number for HiDeNN-PGD. Since the the HiDeNN-PGD is always more accurate than
the usual PGD, the selected prescribed number should be large enough.

(a) Convergence with the element size

(b) Convergence with DoFs

Figure 6: Convergence rate of FEM, PGE, HiDeNN-PGD with respect to mesh
refinement Figure 7: Convergence of PGD with respect to the increasing number
of modes for different meshes

###  5.2 3D case

The proposed HiDeNN-PGD method has been tested in three dimensional cases.
Similar to the previous two-dimensional example, Table 5 reports the evolution
of error at an increasing mode number for different methods on a coarse mesh.
Again, the HiDeNN-PGD outperforms the other methods in terms of accuracy and
DoFs. The same conclusion can be drawn on finer meshes, as reported in Table
6.

Table 5: Accuracy comparison for different methods on the coarse $40\times 40\times 40$ mesh | PGD | CD | FEM | HiDeNN-PGD | HiDeNN
---|---|---|---|---|---
Mode number | DoFs | Err | DoFs | Err | DoFs | Err | DoFs | Err | DoFs | Err
1 | 117 | $72.474\%$ | 117 | $72.474\%$ | 59,319 | $\bm{27.870\%}$ | 234 | $72.349\%$ | 255,996 | $\bm{10.780\%}$
2 | 234 | $29.883\%$ | 234 | $27.927\%$ |  | - | 351 | $\bm{10.797\%}$ |  | -
3 | 351 | $28.026\%$ | 351 | $27.923\%$ |  | - | 468 | $10.797\%$ |  | -
4 | 468 | $27.931\%$ | 468 | $27.920\%$ |  | - | 585 | $10.797\%$ |  | -
5 | 585 | $27.895\%$ | 585 | $27.872\%$ |  | - | 702 | $10.796\%$ |  | -
6 | 702 | $27.881\%$ | 702 | $27.871\%$ |  | - | - | - |  | -
7 | 819 | $27.877\%$ | 819 | $27.871\%$ |  | - | - | - |  | -
8 | 936 | $27.874\%$ | 936 | $\bm{27.870\%}$ |  | - | - | - |  | -
17 | 1,989 | $\bm{27.870\%}$ |  | - |  | - |  | - |  | -
20 | 2,340 | $27.870\%$ |  | - |  | - | - | - |  | -
Table 6: Accuracy comparison for different methods with different meshes for 3D problem  | PGD | CD | FEM | HiDeNN-PGD | HiDeNN
---|---|---|---|---|---
Mesh | DoFs | Err | DoFs | Err | DoFs | Err | DoFs | Err | DoFs | Err
$40\times 40\times 40$ | 1,989(17) | $27.870\%$ | 936(8) | $27.870\%$ | 59,319 | $\bm{27.870\%}$ |  702(5) | $\bm{10.796\%}$ | 255,996 | $\bm{10.780\%}$
$80\times 80\times 80$ | 2,607(11) | $14.416\%$ | 1185(5) | $14.416\%$ | 493,039 | $14.416\%$ | 1185(4) | $6.771\%$ | 2,047,996 |
$160\times 160\times 160$ | 4,770(10) | $7.247\%$ | 2,385(5) | $7.247\%$ | 4,019,679 | $7.247\%$ | 3,339(6) | $4.036\%$ | 16,383,996 |
$320\times 320\times 320$ | 9,570(10) | $3.628\%$ | 5,742(6) | $3.628\%$ | 32,461,759 | $3.628\%$ | 5,742(5) | $1.831\%$ | 131,071,996 |

##  6 Conclusion

A reduced-order hierarchical deep learning network has been proposed. The so-
called HiDeNN-PGD is a combination of HiDeNN and PGD with separated spatial
variables. This combined method presents several advantages over HiDeNN and
PGD methods. First, it allows to leverage the automatic mesh adaptivity of the
HiDeNN method for reducing the mode number in PGD approximation. Second,
combining PGD with HiDeNN reduces significantly the number of degrees of
freedom for HiDeNN and potentially leads to a high computational efficiency.
Furthermore, we have demonstrated that both HiDeNN and HiDeNN-PGD can provide
more accurate solutions than FEM and PGD (or CD), through an error analysis
with the help of analyzing the approximation function spaces.

The numerical results have confirmed the mathematical analysis. These examples
have been performed based on 2D and 3D Poisson problems. It is shown that the
proposed HiDeNN-PGD method can provide accurate solutions with the least
degrees of freedom. In order to have an idea for the prescribed number of
modes in HiDeNN-PGD, we have studied numerically the convergence rate on PGD
approximation. It has been found that the convergence rate on the mode number
is insensitive to the mesh size. Therefore, we can expect to use a coarse mesh
PGD to compute a roughly estimated mode number for HiDeNN-PGD. This finding is
interesting and provides a useful guideline on the choice of the number of
modes for HiDeNN-PGD or other PGD-based methods that may require a better
optimality in terms of basis.

In the future, the computational efficiency of HiDeNN-PGD will be further
explored based on realistic problems, such as thermo-mechanical analysis in
additive manufacturing or multi-scale composite simulations. In terms of
convergence studies, theoretical results need to be derived through a rigorous
mathematical analysis. The numerical results provided in this paper can serve
as the first evidence for demonstrating the capabilities of the method.

## Acknowledgement

L. Zhang, and S. Tang are supported by National Natural Science Foundation of
China grant number 11890681, 11832001, 11521202 and 11988102. W. K. Liu and Y.
Lu are supported by National Science Foundation grant number CMMI-1934367 and
CMMI-1762035.

##  Appendix A Convergence for general canonical decomposition method at
increasing number of modes

In this section, we discuss the convergence of the canonical decomposition
method at increasing $Q$ (number of modes). In Sect. 2, we have shown that
$\mathcal{M}_{Q}^{h}\subset\mathcal{V}^{h}$, provided that their
interpolations are based on the same basis functions. Here, we make some
further discussions.

For 2D case, we compare

$\mathcal{M}_{Q}^{h}=\left\\{u^{h}\bigg{|}u^{h}=\sum_{q=1}^{Q}\left(\sum_{I=1}^{n_{1}}N_{I}(x)\beta_{I}^{(q)}\right)\left(\sum_{J=1}^{n_{2}}N_{J}(y)\gamma_{J}^{(q)}\right),\beta_{I}^{(q)},\gamma_{J}^{(q)}\in\mathbb{R},I=1,\cdots,n_{1},J=1,\cdots,n_{2}\right\\}$
(48)

with

$\mathcal{V}^{h}=\left\\{u^{h}\bigg{|}u^{h}=\sum_{I=1}^{n_{1}}\sum_{J=1}^{n_{2}}N_{I}(x)N_{J}(y)u_{(I,J)},u_{(I,J)}\in\mathbb{R},I=1,\cdots,n_{1},J=1,\cdots,n_{2}\right\\}.$
(49)

As the same basis functions are used, it follows that:

1.

$\forall Q\in\mathbb{N},\mathcal{M}_{Q}^{h}\subset\mathcal{V}^{h}$;

2.

If $Q_{1}\leq Q_{2}$, $Q_{1},Q_{2}\in\mathbb{N}$,
$\mathcal{M}_{Q_{1}}^{h}\subset\mathcal{M}_{Q_{2}}^{h}$;

3.

$\forall
Q\geq\min\\{n_{1},n_{2}\\},Q\in\mathbb{N},\mathcal{M}_{Q}^{h}=\mathcal{V}^{h}$.

The first two statements are straightforward. The last property is proved in
Appendix A.

We then conclude the following relationship:

$\mathcal{M}^{h}_{Q=1}\subset\mathcal{M}^{h}_{Q=2}\subset\cdots\subset\mathcal{M}^{h}_{Q=\min\\{n_{1},n_{2}\\}}=\mathcal{M}^{h}_{Q=\min\\{n_{1},n_{2}\\}+1}=\cdots=\mathcal{V}^{h}.$
(50)

In other words, $\mathcal{M}^{h}_{Q}$ is always a subset of $\mathcal{V}^{h}$,
and it becomes the same as $\mathcal{V}^{h}$ when $Q$ increases. Consequently,
when enough number of modes are taken, the canonical decomposition result
reaches the same accuracy as the FEM solution. We remark that
$\min\\{n_{1},n_{2}\\}$ is precisely the minimal number of modes to ensure
$\mathcal{M}^{h}_{Q}=\mathcal{V}^{h}$.

The above discussions extend to 3D readily.

1.

$\forall Q\in\mathbb{N},\mathcal{M}_{Q}^{h}\subset\mathcal{V}^{h}$;

2.

If $Q_{1}\leq
Q_{2},Q_{1},Q_{2}\in\mathbb{N},\mathcal{M}_{Q_{1}}^{h}\subset\mathcal{M}_{Q_{2}}^{h}$;

3.

$\forall
Q\geq\min\\{n_{1}n_{2},n_{2}n_{3},n_{1}n_{3}\\},Q\in\mathbb{N},\mathcal{M}_{Q}^{h}=\mathcal{V}^{h}$.

We note that the minimal number of modes to ensure
$\mathcal{M}_{Q}^{h}=\mathcal{V}^{h}$ in 3D is essentially to find a best
rank-$r$ approximation to order-3 tensor, which is an open mathematical
problem. An upper bound is given in property 3.

##  Appendix B Convergence for 2D canonical decomposition method at increasing
number of modes

###### Theorem B.1

For 2D case, $\mathcal{M}_{Q}^{h}$ and $\mathcal{V}^{h}$ are defined in
Eq.(48) and Eq.(49), respectively. We have

$\forall
Q\geq\min(n_{1},n_{2}),Q\in\mathbb{N},\mathcal{M}_{Q}^{h}=\mathcal{V}^{h}.$
(51)

###### Proof

For any interpolation function

$u^{h,FEM}=\sum_{i=1}^{n_{1}}\sum_{j=1}^{n_{2}}N_{i}(x)N_{j}(y)u_{i,j}$ (52)

in the set $\mathcal{V}^{h}$, we write the nodal values in the form of matrix,

$\bm{U}=\left[\begin{array}[]{cccc}u_{1,1}&u_{1,2}&\cdots&u_{1,n_{2}}\\\
u_{2,1}&u_{2,2}&\cdots&u_{2,n_{2}}\\\ \vdots&\vdots&&\vdots\\\
u_{n_{1},1}&u_{n_{1},2}&\cdots&u_{n_{1},n_{2}}\end{array}\right].$ (53)

According to the singular value decomposition (SVD), $\bm{U}$ is represented
by

$\bm{U}=\sum_{q=1}^{rank(\bm{U})}\sigma^{(q)}\bm{w}^{(q)}\otimes\bm{v}^{(q)},\sigma^{(1)}\geq\sigma^{(2)}\geq\cdots\geq\sigma^{(rank(\bm{U}))}>0,$
(54)

where $\bm{w}^{(q)}$ is the $n_{1}$-dimensional vector, and $\bm{v}^{(q)}$ is
the $n_{2}$-dimensional vector. Thus $u^{h,FEM}$ is rewritten in the form of
the separation of variables, i.e.,

$u^{h,FEM}=\sum_{i=1}^{n_{1}}\sum_{j=1}^{n_{2}}N_{i}(x)N_{j}(y)\left(\sum_{q=1}^{rank(\bm{U})}\sigma^{(q)}w_{i}^{(q)}v_{j}^{(q)}\right).$
(55)

So we have $u^{h,FEM}\in\mathcal{M}_{Q}^{h}$, if
$Q\geq\min\\{n_{1},n_{2}\\}\geq rank(\bm{U})$. Combining with
$\mathcal{M}_{Q}^{h}\subset\mathcal{V}^{h}$, we obtain Eq. (51).

We remark that SVD tells us the minimal number of modes to reproduce the FE
solution, i.e. $\min\\{n_{1},n_{2}\\}$.

##  Appendix C 1D HiDeNN Formulation

In standard 1D FEM, the computational domain $\Omega$ is discretized by a grid
with $n$ nodes and the shape function associated with an internal node $x_{I}$
is

$N_{I}(x)=\left\\{\begin{array}[]{cc}\dfrac{x-x_{I-1}}{x_{I}-x_{I-1}},&x_{I-1}\leq
x\leq x_{I},\\\ \dfrac{x_{I+1}-x}{x_{I+1}-x_{I}},&x_{I}\leq x\leq x_{I+1},\\\
0,&elsewhere,\end{array}\right.$ (56)

where $x_{I-1}$ and $x_{I+1}$ are the two neighbor points of the node
$\>x_{I}$ from the left side and right side, respectively.

We rewrite $N_{I}(x)$ in a DNN format consists of weights, biases, and
activation functions. Considering the shape function is a piecewise linear
function, the activation function is selected as ReLU function, i.e.,
$\mathcal{A}_{1}=\max(0,x)$. Fig. 8(a) shows the DNN representation of the
linear shape function. The corresponding formula is

$\displaystyle\mathscr{N}_{I}(x;\bm{W},\bm{b},\>\bm{\mathcal{A}})$
$\displaystyle=$ $\displaystyle
W_{11}^{l=4}\mathcal{A}_{1}\left(W_{11}^{l=3}\mathcal{A}_{1}\left(W_{11}^{l=2}x+b_{1}^{l=2}\right)+b_{1}^{l=3}\right)$
$\displaystyle+W_{21}^{l=4}\mathcal{A}_{1}\left(W_{22}^{l=3}\mathcal{A}_{1}\left(W_{12}^{l=2}x+b_{2}^{l=2}\right)+b_{2}^{l=3}\right)+b_{1}^{l=4}$
$\displaystyle=$
$\displaystyle\mathcal{A}_{1}\left(\dfrac{-1}{x_{I}-x_{I-1}}\mathcal{A}_{1}\left(-x+x_{I}\right)+1\right)+\mathcal{A}_{1}\left(\dfrac{-1}{x_{I+1}-x_{I}}\mathcal{A}_{1}\left(x-x_{I}\right)+1\right)-1,$

where
$\bm{W}=[W_{11}^{l=2},W_{12}^{l=2},W_{11}^{l=3},W_{22}^{l=3},W_{11}^{l=4},W_{21}^{l=4}]$,
and $\bm{b}=[b_{1}^{l=2},b_{2}^{l=2},b_{1}^{l=3},b_{2}^{l=3},b_{1}^{l=4}]$ are
the weights and biases of the connected neurons. Note that all the weights and
biases are functions of nodal coordinates. The formula can be rewritten as the
form of

$\mathcal{N}_{I}(\bm{x};\bm{x}_{I}^{*},\bm{\mathcal{A}}),$ (58)

where $\bm{x}_{I}^{*}$ denotes the vector that represents the neighbor nodes
of node $\bm{x}_{I}$ involved in $N_{I}(\bm{x})$. For 1D linear shape
function, it should be $\bm{x}_{I}^{*}=[x_{I-1},\>x_{I},\>x_{I+1}]$. For the
sake of clarity, one more layer is added to introduce the nodal value $u_{I}$,
i.e., the formula becomes

$\displaystyle\mathscr{u}_{I}^{h}$ $\displaystyle=$
$\displaystyle\mathscr{N}_{I}(x;\>\bm{W},\bm{b},\>\bm{\mathcal{A}})\mathscr{u}_{I}=\mathscr{N}_{I}(x;\>\bm{x}_{I}^{*},\>\bm{\mathcal{A}})\mathscr{u}_{I};\mbox{
no summation on }{I}$ $\displaystyle=$
$\displaystyle\mathcal{A}_{0}\left(\mathcal{A}_{1}\left(\dfrac{-1}{x_{I}-x_{I-1}}\mathcal{A}_{1}\left(-x+x_{I}\right)+1\right)-0.5\right)\mathscr{u}_{I}$
$\displaystyle+\mathcal{A}_{0}\left(\mathcal{A}_{1}\left(\dfrac{-1}{x_{I+1}-x_{I}}\mathcal{A}_{1}\left(x-x_{I}\right)+1\right)-0.5\right)\mathscr{u}_{I},$

where $\mathscr{u}_{I}^{h}$ and $\mathscr{u}_{I}$ are the interpolated
displacement and nodal displacement at node $x_{I}$,
$\bm{\mathcal{A}}=[\mathcal{A}_{0},\>\mathcal{A}_{1}]$ are the activation
functions used for the construction of the DNN approximation.
$\mathcal{A}_{0}(x)=x$ is an identical function. Fig. 8(b) gives the DNN
representation of the interpolation of the nodal displacement at node $x_{I}$.

(a) DNN-based 1D shape function

(b) DNN-based 1D interpolation function

Figure 8: Deep neural network (DNN) representation of the 1D global shape
function and interpolation function.

Once the shape function with nodal value for an arbitrary node $x_{I}$ is
constructed, the interpolation is obtained by assembling all DNNs, i.e.,

$\mathscr{u}^{h}(x)=\sum_{I=1}^{n}\mathscr{N}_{I}(x;\>\bm{x}_{I}^{*},\>\bm{\mathcal{A}})\mathscr{u}_{I}.$
(60)

Compared with classical FEM, nodal positions are introduced as additional DoFs
in the optimization for HiDeNN, which increases both the local and global
accuracy of the interpolants.

Reference zhang2021hierarchical also presented the DNN representation of
various rational functions including Lagrange polynomials, B-spline,
Reproducing Kernel Particle Method (RKPM), NURBS, Isogeometric analysis (IGA),
etc., and multidimensional shape functions.

##  Appendix D Space separated PGD

When the domain is not intrinsically separable, fully separated representation
can not be applied directly. gonzalez2010recent immersed the non-separable
domain onto a fully separable one. ghnatios2019advanced used geometrical
mapping to deal with layered domain, where interfaces are not planar. Now we
combine representation of separated variables with FE mapping to deal with
more complex geometrical cases. In section 4.1, we define the mapping to a
regular parameter space by means of FE geometrical mapping. Then section 4.2
addresses several solution schemes. Finally, we take a 2D Poisson problem as
illustration for the whole solution procedure in section 4.3.

###  D.1 Mesh mapping and recovering for irregular domains

Figure 9: A quarter of ring is transformed in to a rectangular by the polar
coordinates transformation.

The main idea is to map original irregular domain $\Omega_{(\bm{x})}$ to a
regular one $\tilde{\Omega}_{(\tilde{\bm{x}})}$, and then apply the separated
representation. Fig, 9 illustrates a simple example. A quarter of ring becomes
a rectangular through polar transformation. Then the new representation of
separated variables is shown as below,

$u^{h}=\sum_{q=1}^{Q}\tilde{X}^{(q)}(\theta)\tilde{Y}^{(q)}(r),$ (61)

where $\theta$ and $r$ are the functions of space coordinates $x,y$.

Figure 10: Illustration for the geometrical mapping. The irregular domain with
irregular mesh is related to a regular domain with regular mesh by a 2-step
mapping.

By virtue of parametric transformation in FEM, we propose a general way to
define this mapping as illustrated in Fig 10. We present a FE mesh over the 2D
irregular computational domain $\Omega_{\bm{x}}$ first with nodes
$\bm{x}_{i,j},i=1,2,\cdots,n_{1},j=1,2,\cdots,n_{2}$. Then we define a mapping
to its corresponding lattice $(i,j)$. $\tilde{\bm{x}}$ is the coordinates of
the transformed domain $\tilde{\Omega}_{\tilde{\bm{x}}}$. The mapping consist
of two steps:

1\. Mapping each element to a square or cubic

The first mapping is the classical parametric mapping in FEM. We make a change
of coordinates which maps the 4-node element into a square $[-1,1]^{2}$ for 2D
or maps the 8-node element into a cubic $[-1,1]^{3}$ for 3D. The coordinates
of a point $\bm{\xi}$ in the square is related to the physical coordinates of
a point $\bm{x}$ in the element by mappings of the form

$\bm{x}=\sum_{a=1}^{n_{e}}N_{a}^{e}(\bm{\xi})\bm{x}_{a}^{e}$ (62)

where $n_{e}$ is the number of nodes of the element ($n_{e}=4$ for 2D and
$n_{e}=8$ for 3D), $\bm{x}_{a}^{e}$ is the coordinates of the $a$-th node of
the element, and $N_{a}^{e}(\bm{\xi})$ is the corresponding shape function.
$\bm{\xi}$ is called as natural coordinates.

2\. Mapping the square to a lattice

For the sake of separated representation, we define the second mapping to
assemble the square into a lattice. The transformed formula is

$\tilde{\bm{x}}=\sum_{a=1}^{n_{e}}N_{a}^{e}(\bm{\xi})\tilde{\bm{x}}_{a}^{e}$
(63)

and the inverse transformation is

$\displaystyle\xi$ $\displaystyle=$
$\displaystyle\dfrac{2\tilde{x}-\tilde{x}^{e}_{1}-\tilde{x}^{e}_{2}}{\tilde{x}^{e}_{2}-\tilde{x}^{e}_{1}}$
(64) $\displaystyle\eta$ $\displaystyle=$
$\displaystyle\dfrac{2\tilde{y}-\tilde{y}^{e}_{1}-\tilde{y}^{e}_{2}}{\tilde{y}^{e}_{2}-\tilde{y}^{e}_{1}}$
(65) $\displaystyle\zeta$ $\displaystyle=$
$\displaystyle\dfrac{2\tilde{z}-\tilde{z}^{e}_{1}-\tilde{z}^{e}_{2}}{\tilde{z}^{e}_{2}-\tilde{z}^{e}_{1}}.$
(66)

The final transformed domain $\tilde{\Omega}_{\tilde{\bm{x}}}$ is called as
reference domain, which is a regular domain with a regular mesh
$[\tilde{x}_{1},\tilde{x}_{2},\cdots,\tilde{x}_{n_{1}}]\times[\tilde{y}_{1},\tilde{y}_{2},\cdots,\tilde{y}_{n_{2}}]\times[\tilde{z}_{1},\tilde{z}_{2},\cdots,\tilde{z}_{n_{3}}]$.
$\tilde{x}^{e}_{1},\tilde{x}^{e}_{2},\tilde{y}^{e}_{1},\tilde{y}^{e}_{2},\tilde{z}^{e}_{1},\tilde{z}^{e}_{2}$
are the coordinates of the element
$[\tilde{x}^{e}_{1},\tilde{x}^{e}_{2}]\times[\tilde{y}^{e}_{1},\tilde{y}^{e}_{2}]\times[\tilde{z}^{e}_{1},\tilde{z}^{e}_{2}]$
in the reference domain. For convenience, we might take the mesh in the
reference domain as a lattice corresponding to the index of the nodes in the
physical domain, i.e.,
$\tilde{x}_{i}=i,i=1,2,\cdots,n_{1},\tilde{y}_{j}=j,j=1,2,\cdots,n_{2},\tilde{z}_{k}=k,k=1,2,\cdots,n_{3}$.

The whole mapping is defined as below,

$\bm{x}=\sum_{a=1}^{n_{e}}N_{a}^{e}(\bm{\xi}(\tilde{\bm{x}}))\bm{x}_{a}^{e}.$
(67)

Then separation of spatial variables is applicable to the reference domain,
i.e., the interpolation function set is

$\tilde{\mathcal{M}}^{h}_{Q}=\left\\{u^{h}\Bigg{|}u^{h}=\sum_{q=1}^{Q}\tilde{X}(\tilde{\bm{x}})\tilde{Y}(\tilde{y})\tilde{Z}(\tilde{z})=\sum_{q=1}^{Q}\left(\sum_{i=1}^{n_{1}}N_{i}(\tilde{x})\beta_{i}^{(q)}\right)\left(\sum_{j=1}^{n_{2}}N_{j}(\tilde{y})\gamma_{j}^{(q)}\right)\left(\sum_{k=1}^{n_{3}}N_{k}(\tilde{z})\theta_{k}^{(q)}\right)\right\\},$
(68)

where $N_{i}(\tilde{x}),N_{j}(\tilde{y}),N_{i}(\tilde{z})$ are shape
functions, and $\beta_{i}^{(q)},\gamma_{j}^{(q)},\theta_{k}^{(q)}$ are the
corresponding coefficients of the $q$-th mode. Note that
$\tilde{x},\tilde{y},\tilde{z}$ are the functions of physical coordinates
$x,y,z$.

###  D.2 Solution schemes

There are different solution schemes to solve the problem. A straight way is
to find the solution by minimizing the variational formula with the given
number of modes $Q^{*}$ directly, i.e.,

$u^{h}=\operatorname*{arg\,min}_{u^{h}\in\tilde{\mathcal{M}}^{h}_{Q=Q^{*}}}\Pi(u^{h}(\bm{x};\bm{\beta}^{(q)},\bm{\gamma}^{(q)},\bm{\theta}^{(q)})).$
(69)

All the parameters $\beta$ are solved at the same time.

Yet this global optimization might be expensive, so we might borrow the idea
from PGD ammar2006new ; chinesta2013proper , i.e., incremental solution
scheme. More precisely, the solution scheme is

$\begin{array}[]{lll}\text{The first mode:}&u^{PGD,(1)}=\arg\min_{\Delta
u\in\tilde{\mathcal{M}}^{h}_{Q=1}}\Pi[u^{0}+\Delta u];\\\ \text{For $m>1$, the
$m$-th mode:}&\Delta u^{(m)}=\arg\min_{\Delta
u\in\tilde{\mathcal{M}}^{h}_{Q=1}}\Pi[u^{PGD,(m-1)}+\Delta u],\\\ \text{The
PGD solution with $m$ modes:}&u^{PGD,(m)}=u^{PGD,(m-1)}+\Delta
u^{(m)},\end{array}$ (70)

We remarked it is also possible to solve several modes simultaneously in one
incremental step.

In general, the initial guess $u_{0}$ is set to be zero. When dealing with the
boundary conditions, $u^{0}$ can be arbitrary continuous functions satisfying
boundary conditions. We can also set an appropriate initial guess to improve
the efficiency.

These two solution schemes have their advantages and disadvantages. With the
same number of modes $Q^{*}$, we solve the modes one by one using the latter
solution scheme while optimize all the modes together using the previous one.
Thus, we have $\Pi(u^{CD})<\Pi(u^{PGD})$ and then obtain

$\|u^{CD}-u^{exact}\|_{E}\leq\|u^{PGD}-u^{exact}\|_{E}.$ (71)

This indicates the solution of the previous solution scheme might be better
than that of the incremental way. However, the previous one might cost more.

###  D.3 Illustrating the solution procedure: the 2D Poisson problem

For the sake of simplicity and without loss of generality, we consider 2D
Poisson problem with incremental solution scheme for illustration,

$\left\\{\begin{array}[]{l}\nabla^{2}\bm{u}(x,y)+b(x,y)=0\text{ in
}\Omega_{(x,y)}\subset\mathbb{R}^{2},\\\
\bm{u}|_{\partial\Omega}=\bm{0}.\end{array}\right.$ (72)

(72) is solved in the irregular domain $\Omega_{(x,y)}$ with homogeneous
boundary conditions.

The solution is assumed in the form of (48). Then we solve it with incremental
solution scheme. Let previous $q-1$ modes solved. The $q$-th mode is obtained
by

$\Delta u^{(q)}=\operatorname*{arg\,min}_{\Delta
u\in\tilde{\mathcal{M}}^{h}_{Q=1}}\Pi[u^{PGD,(q-1)}+\Delta u^{(q)}],$ (73)

where $u^{PGD,(q-1)}$ is the sum of the previous $q-1$ modes. We rewrite the
interpolation function in the following matrix form

$u^{PGD,(q)}=u^{PGD,(q-1)}+\Delta u^{(q)},\Delta
u^{(q)}=\left((\bm{\beta}^{(q)})^{T}\bm{N}^{\beta}(\tilde{x})\right)\left((\bm{\gamma}^{(q)})^{T}\bm{N}^{\gamma}(\tilde{y})\right),$
(74)

where $\bm{\beta}^{(q)},\bm{\gamma}^{(q)}$ are the coefficient vector, and
$\bm{N}^{\beta}(\tilde{x}),\bm{N}^{\gamma}(\tilde{y})$ denotes the vector
containing shape functions.

Substituting (74) into the variational formula (45), we have

$\displaystyle\Pi(u^{PGD,(q)})$ $\displaystyle=$
$\displaystyle\dfrac{1}{2}\int_{\Omega_{(x,y)}}\left(\nabla(\Delta
u^{(q)})\right)^{2}\mathrm{d}x\mathrm{d}y+\int_{\Omega_{(x,y)}}\left(\nabla(\Delta
u^{(q)})\right)\left(\nabla(u^{PGD,(q-1)})\right)\mathrm{d}x\mathrm{d}y$
$\displaystyle-\int_{\Omega_{(x,y)}}u(x,y)b(x,y)\mathrm{d}x\mathrm{d}y+\Pi(u^{PGD,(q-1)}).$

The quadratic term with respect to $\bm{\beta}^{(q)},\bm{\gamma}^{(q)}$ in the
variational formula is given by

$\displaystyle\int_{\Omega_{(x,y)}}\left(\nabla(\Delta
u^{(q)})\right)^{2}\mathrm{d}x\mathrm{d}y$ $\displaystyle=$
$\displaystyle\int_{\tilde{\Omega}_{(\tilde{x},\tilde{y})}}\left(\nabla_{(\tilde{x},\tilde{y})}(\Delta
u^{(q)})\right)^{T}\bm{J}^{-T}\bm{J}^{-1}\nabla_{(\tilde{x},\tilde{y})}(\Delta
u^{(q)})\mathrm{det}(\bm{J})\mathrm{d}\tilde{x}\mathrm{d}\tilde{y}.$ (76)

with the Jacobi matrix

$\bm{J}=\dfrac{\partial(x,y)}{\partial(\tilde{x},\tilde{y})}.$ (77)

The gradient of $\Delta u^{(q)}$ is

$\nabla_{(\tilde{x},\tilde{y})}\Delta
u^{(q)}=\left[\begin{array}[]{c}\dfrac{\partial}{\partial\tilde{x}}\\\
\dfrac{\partial}{\partial\tilde{y}}\end{array}\right]\Delta
u^{(q)}=\left[\begin{array}[]{c}(\bm{\beta}^{(q)})^{T}\dfrac{d\bm{N}^{\beta}(\tilde{x})}{d\tilde{x}}(\bm{N}^{\gamma}(\tilde{y}))^{T}\gamma^{(q)}\\\
(\bm{\beta}^{(q)})^{T}\bm{N}^{\beta}(\tilde{x})(\dfrac{d\bm{N}^{\gamma}(\tilde{y})}{d\tilde{y}})^{T}\gamma^{(q)}\end{array}\right].$
(78)

Note that the expression is a nonlinear algebraic system. It is hard to solve
it directly, so we use the alternating direction strategy as below:

In each iteration step,

1\. Fix $\gamma$ and solve $\beta$

The quadratic term becomes

$\int_{\Omega_{(x,y)}}\left(\nabla(\Delta
u^{(q)})\right)^{2}\mathrm{d}x\mathrm{d}y=\int_{\tilde{\Omega}_{(\tilde{x},\tilde{y})}}\beta^{(q),T}\bm{B}^{\beta,T}(\tilde{x},\tilde{y})\bm{B}^{\beta}(\tilde{x},\tilde{y})\beta^{(q)}\mathrm{det}(\bm{J})\mathrm{d}\tilde{x}\mathrm{d}\tilde{y}$
(79)

with

$\bm{B}^{\beta}(\tilde{x},\tilde{y})=\bm{J}^{-1}\left[\begin{array}[]{c}(\bm{\gamma}^{(q)})^{T}\bm{N}^{\gamma}(\tilde{y})(\dfrac{d\bm{N}^{\beta}(\tilde{x})}{d\tilde{x}})^{T}\\\
(\bm{\gamma}^{(q)})^{T}\dfrac{d\bm{N}^{\gamma}(\tilde{y})}{d\tilde{y}}(\bm{N}^{\beta}(\tilde{x}))^{T}\end{array}\right].$
(80)

The stiffness matrix for $\beta^{(q)}$ is

$\bm{K}^{\beta}=\int_{\tilde{\Omega}_{(\tilde{x},\tilde{y})}}\bm{B}^{\beta,T}(\tilde{x},\tilde{y})\bm{B}^{\beta}(\tilde{x},\tilde{y})\mathrm{det}(\bm{J})\mathrm{d}\tilde{x}\mathrm{d}\tilde{y}.$
(81)

2\. Fix $\beta$ and solve $\gamma$

The quadratic term becomes

$\int_{\Omega_{(x,y)}}\left(\nabla(\Delta
u^{(q)})\right)^{2}\mathrm{d}x\mathrm{d}y=\int_{\tilde{\Omega}_{(\tilde{x},\tilde{y})}}\gamma^{(q),T}\bm{B}^{\gamma,T}(\tilde{x},\tilde{y})\bm{B}^{\gamma}(\tilde{x},\tilde{y})\gamma^{(q)}\mathrm{det}(\bm{J})\mathrm{d}\tilde{x}\mathrm{d}\tilde{y}$
(82)

with

$\bm{B}^{\gamma}(\tilde{x},\tilde{y})=\bm{J}^{-1}\left[\begin{array}[]{c}(\bm{\beta}^{(q)})^{T}\dfrac{d\bm{N}^{\beta}(\tilde{x})}{d\tilde{x}}(\bm{N}^{\gamma}(\tilde{y}))^{T}\\\
(\bm{\beta}^{(q)})^{T}\bm{N}^{\beta}(\tilde{x})(\dfrac{d\bm{N}^{\gamma}(\tilde{y})}{d\tilde{y}})^{T}\end{array}\right].$
(83)

The stiffness matrix for $\gamma^{(q)}$ is

$\bm{K}^{\gamma}=\int_{\tilde{\Omega}_{(\tilde{x},\tilde{y})}}\bm{B}^{\gamma,T}(\tilde{x},\tilde{y})\bm{B}^{\gamma}(\tilde{x},\tilde{y})\mathrm{det}(\bm{J})\mathrm{d}\tilde{x}\mathrm{d}\tilde{y}.$
(84)

If we present a regular mesh over a regular domain $\Omega_{(x,y)}$, the
mapping is a linear transformation for coordinates. Let one element of the
regular mesh $[x^{e}_{1},x^{e}_{2}]\times[y^{e}_{1},y^{e}_{2}]$. The mapping
(67) reduces to

$x=\dfrac{x^{e}_{2}-x^{e}_{1}}{\tilde{x}^{e}_{2}-\tilde{x}^{e}_{1}}(\tilde{x}-\tilde{x}^{e}_{1})+x^{e}_{1},y=\dfrac{y^{e}_{2}-y^{e}_{1}}{\tilde{y}^{e}_{2}-\tilde{y}^{e}_{1}}(\tilde{y}-\tilde{y}^{e}_{1})+y^{e}_{1}.$
(85)

$\bm{J}$ reduces to a diagonal matrix
$\mathrm{diag}(\dfrac{x^{e}_{2}-x^{e}_{1}}{\tilde{x}^{e}_{2}-\tilde{x}^{e}_{1}},\dfrac{y^{e}_{2}-y^{e}_{1}}{\tilde{y}^{e}_{2}-\tilde{y}^{e}_{1}})$.
Thus we have

$\bm{J}^{-T}\bm{J}^{-1}\mathrm{det}(J)=\dfrac{x^{e}_{2}-x^{e}_{1}}{\tilde{x}^{e}_{2}-\tilde{x}^{e}_{1}}\dfrac{y^{e}_{2}-y^{e}_{1}}{\tilde{y}^{e}_{2}-\tilde{y}^{e}_{1}},$
(86)

which is constant in each element and separated representation. Thus this
method degenerates to the classical PGD. For irregular mesh,
$\bm{J}^{-T}\bm{J}^{-1}\mathrm{det}(J)$ in (76) is the function of $\bm{x}$
and mostly non-separated representation. In the numerical implementation, we
usually approximate it by a separated form using SVD technique, i.e.,

$\bm{J}^{-T}\bm{J}^{-1}\mathrm{det}(J)\approx\left[\begin{array}[]{cc}\sum_{a}\phi_{11}^{(a)}(\tilde{x})\psi_{11}^{(a)}(\tilde{y})&\sum_{a}\phi_{12}^{(a)}(\tilde{x})\psi_{12}^{(a)}(\tilde{y})\\\
\sum_{a}\phi_{12}^{(a)}(\tilde{x})\psi_{12}^{(a)}(\tilde{y})&\sum_{a}\phi_{22}^{(a)}(\tilde{x})\psi_{22}^{(a)}(\tilde{y})\end{array}\right].$
(87)

This converts the 2D integration (76) to the product of 1D integration along
different directions ($\tilde{x}$ and $\tilde{y}$ directions in the reference
domain), which might reduce the computational cost for integration.

###  D.4 Numerical exmples

(a) Physical domain

(b) FE mesh

(c) Reference solution

Figure 11: Geometry and mesh for the plate. (a) Model of the plate. (b) $100$
elements mesh as an illustration. (c) Reference solution with $4.90\times
10^{7}$ elements.

(a) Physical domain

(b) FE mesh

(c) Reference mesh

(d) Reference solution

Figure 12: Geometry and mesh for the plate with one hole. (a) Model of the
plate with one hole. (b) Discretization by a $84$ elements mesh. (c) Reference
mesh corresponding to $84$ elements.

## References

  * (1) T. Belytschko, W. K. Liu, B. Moran, K. Elkhodary, Nonlinear Finite Elements for Continua and Structures, John wiley & sons, 2013.
  * (2) T. J. Hughes, J. A. Cottrell, Y. Bazilevs, Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement, Computer methods in applied mechanics and engineering 194 (39-41) (2005) 4135–4195.
  * (3) K. Hornik, M. Stinchcombe, H. White, et al., Multilayer feedforward networks are universal approximators, Neural Networks 2 (5) (1989) 359–366.
  * (4) G. Cybenko, Approximation by superpositions of a sigmoidal function, Mathematics of Control, Signals and Systems 2 (4) (1989) 303–314.
  * (5) W. E, J. Han, A. Jentzen, Deep learning-based numerical methods for high-dimensional parabolic partial differential equations and backward stochastic differential equations, Communications in Mathematics and Statistics 5 (4) (2017) 349–380.
  * (6) 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.
  * (7) H. Li, O. L. Kafka, J. Gao, C. Yu, Y. Nie, L. Zhang, M. Tajdari, S. Tang, X. Guo, G. Li, et al., Clustering discretization methods for generation of material performance databases in machine learning and design optimization, Computational Mechanics 64 (2) (2019) 281–305.
  * (8) L. Zhang, L. Cheng, H. Li, J. Gao, C. Yu, R. Domel, Y. Yang, S. Tang, W. K. Liu, Hierarchical deep-learning neural networks: finite elements and beyond, Computational Mechanics 67 (1) (2021) 207–230.
  * (9) S. Saha, Z. Gan, L. Cheng, J. Gao, O. L. Kafka, X. Xie, H. Li, M. Tajdari, H. A. Kim, W. K. Liu, Hierarchical deep learning neural network (HiDeNN): An artificial intelligence (AI) framework for computational science and engineering, Computer Methods in Applied Mechanics and Engineering 373 (2021) 113452.
  * (10) A. Ammar, B. Mokdad, F. Chinesta, R. Keunings, A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids, Journal of Non-Newtonian Fluid Mechanics 139 (3) (2006) 153–176.
  * (11) F. Chinesta, P. Ladeveze, E. Cueto, A short review on model order reduction based on proper generalized decomposition, Archives of Computational Methods in Engineering 18 (4) (2011) 395–404.
  * (12) M. Bhattacharyya, A. Fau, U. Nackenhorst, D. Néron, P. Ladevèze, A multi-temporal scale model reduction approach for the computation of fatigue damage, Computer Methods in Applied Mechanics and Engineering 340 (2018) 630–656.
  * (13) D. Modesto, S. Zlotnik, A. Huerta, Proper generalized decomposition for parameterized helmholtz problems in heterogeneous and unbounded domains: application to harbor agitation, Computer Methods in Applied Mechanics and Engineering 295 (2015) 127–149.
  * (14) Y. Lu, N. Blal, A. Gravouil, Adaptive sparse grid based hopgd: Toward a nonintrusive strategy for constructing space-time welding computational vademecum, International Journal for Numerical Methods in Engineering 114 (13) (2018) 1438–1461.
  * (15) Y. Lu, N. Blal, A. Gravouil, Datadriven hopgd based computational vademecum for welding parameter identification, Computational Mechanics 64 (1) (2019) 47–62.
  * (16) Y. Lu, N. Blal, A. Gravouil, Multi-parametric space-time computational vademecum for parametric studies: Application to real time welding simulations, Finite Elements in Analysis and Design 139 (2018) 62–72.
  * (17) P. Díez, S. Zlotnik, A. García-González, A. Huerta, Algebraic pgd for tensor separation and compression: an algorithmic approach, Comptes Rendus Mécanique 346 (7) (2018) 501–514.
  * (18) B. Bognet, A. Leygue, F. Chinesta, Separated representations of 3d elastic solutions in shell geometries, Advanced Modeling and Simulation in Engineering Sciences 1 (1) (2014) 1–34.
  * (19) S. Tang, Y. Yang, Why neural networks apply to scientific computing?, Theoretical and Applied Mechanics Letters (2021) 100242.
  * (20) J. L. Wu, J. X. Wang, H. Xiao, J. Ling, Physics-informed machine learning for predictive turbulence modeling: A priori assessment of prediction confidence (2016).
  * (21) H. Xiao, J. L. Wu, J. X. Wang, R. Sun, C. Roy, Quantifying and reducing model-form uncertainties in reynolds-averaged navier–stokes simulations: A data-driven, physics-informed bayesian approach, Journal of Computational Physics 324 (2016) 115–136.
  * (22) J. Berg, K. Nyström, A unified deep artificial neural network approach to partial differential equations in complex geometries, Neurocomputing 317 (2018) 28–41.
  * (23) M. Raissi, P. Perdikaris, G. E. Karniadakis, Physics informed deep learning (part i): Data-driven solutions of nonlinear partial differential equations, arXiv preprint arXiv:1711.10561.
  * (24) J. Sirignano, K. Spiliopoulos, DGM: A deep learning algorithm for solving partial differential equations, Journal of Computational Physics 375 (2018) 1339–1364.
  * (25) W. E, B. Yu, The deep Ritz method: a deep learning-based numerical algorithm for solving variational problems, Communications in Mathematics and Statistics 6 (1) (2018) 1–12.
  * (26) I. Goodfellow, Y. Bengio, A. Courville, Deep Learning, MIT Press, 2016, http://www.deeplearningbook.org.
  * (27) S. Haykin, Neural Networks: a Comprehensive Foundation, Prentice Hall PTR, 1994\.
  * (28) A. Oishi, G. Yagawa, Computational mechanics enhanced by deep learning, Computer Methods in Applied Mechanics and Engineering 327 (2017) 327–351.
  * (29) T. G. Kolda, B. W. Bader, Tensor decompositions and applications, SIAM Review 51 (3) (2009) 455–500.
  * (30) D. González, A. Ammar, F. Chinesta, E. Cueto, Recent advances on the use of separated representations, International Journal for Numerical Methods in Engineering 81 (5) (2010) 637–659.
  * (31) N. Blal, A. Gravouil, Non-intrusive data learning based computational homogenization of materials with uncertainties, Computational Mechanics 64 (3) (2019) 807–828.
  * (32) F. Chinesta, R. Keunings, A. Leygue, The Proper Generalized Decomposition for Advanced Numerical Simulations: a Primer, Springer Science & Business Media, 2013\.
  * (33) D. P. Kingma, J. Ba, Adam: A method for stochastic optimization.
  * (34) C. Ghnatios, E. Abisset, A. Ammar, E. Cueto, J.-L. Duval, F. Chinesta, Advanced separated spatial representations for hardly separable domains, Computer Methods in Applied Mechanics and Engineering 354 (2019) 802–819.