=Paper=
{{Paper
|id=Vol-2964/article_67
|storemode=property
|title=Deep Autoencoders for Nonlinear Physics-Constrained Data-Driven Computational Framework with Application to Biological Tissue Modeling
|pdfUrl=https://ceur-ws.org/Vol-2964/article_67.pdf
|volume=Vol-2964
|authors=Xiaolong He,Qizhi He,Jiun-Shyan Chen
|dblpUrl=https://dblp.org/rec/conf/aaaiss/HeHC21
}}
==Deep Autoencoders for Nonlinear Physics-Constrained Data-Driven Computational Framework with Application to Biological Tissue Modeling==
https://ceur-ws.org/Vol-2964/article_67.pdf
Deep Autoencoders for Nonlinear Physics-Constrained Data-Driven
Computational Framework with Application to Biological Tissue Modeling
Xiaolong He1 , Qizhi He2 , Jiun-Shyan Chen1
1
Department of Structural Engineering, University of California, San Diego, La Jolla, CA 92093
xiaolong-he@ucsd.edu, js-chen@ucsd.edu
2
Physical and Computational Sciences Directorate, Pacific Northwest National Laboratory, Richland, WA 99354
qizhi.he@pnnl.gov
Abstract (2nd-PK) stress tensor. The superscript X denotes the refer-
ence (undeformed) configuration. F is the deformation gra-
Physics-constrained data-driven computing is an emerging
paradigm that directly integrates material database into physi-
dient related to u, defined as F(u) = ∂(X + u)/∂X, where
cal simulations of complex materials, bypassing the construc- X is the material coordinate. b, N, t, and g are the body
tion of classical constitutive models. However, most of the force, the surface normal on ΓX X
t , the traction on Γt , and the
developed data-driven computing approaches are based on prescribed displacement on ΓX u , respectively.
simplistic distance minimization and thus suffer from dealing To solve the boundary value problem (BVP) in Eq. (1),
with high-dimensional applications and lack generalization material laws describing the stress-strain relation are re-
ability. This study proposes a deep learning enhanced data- quired. However, it is difficult to construct phenomeno-
driven computing framework to address these fundamental logical material models for complex material systems.
issues for nonlinear materials modeling. To this end, an au- The physics-constrained data-driven computing framework
toencoder, a special multi-layer neural network architecture,
(Kirchdoerfer and Ortiz 2016; Ibanez et al. 2018; He and
is introduced to learn the underlying low-dimensional embed-
ding representation of the material database. Incorporating Chen 2020) offers an alternative to directly utilize material
the offline trained autoencoder and the discovered embedding data in physical simulations by formulating the BVP as an
space in online data-driven computation enables to search for optimization problem under physical constraints.
the optimal material state from database in low-dimensional In this framework, the material behavior is described by
embedding space, enhancing the robustness and predictabil- strain-stress pairs, ẑ = (Ê, Ŝ), defined as the material state
ity of data-driven computing on limited material data. To en- given by the material database, E = {ẑI }M I=1 ∈ E, where
hance stability and convergence of data-driven computing, a M is the number of material data points and E is an ad-
convexity-preserving interpolation scheme is introduced for
missible set of material database over the domain Ω. The
constructing the material state on the low-dimensional em-
bedding space given by autoencoders. The effectiveness and strain-stress pairs, z = (E, S) ∈ C, that satisfy the physical
enhanced generalization performance of the proposed ap- equations (Eq. (1)) are defined as the physical state, where
proach are examined by modeling biological tissue with ex- C = {z| Eq. (1)} is a physical admissible set. The data-
perimental data. driven solutions are obtained by fixed-point global-local it-
erations to minimize the distance between the material and
physical states. In the local step, the material data-driven
1 Nonlinear physics-constrained data-driven local solver searches for the optimal strain-stress pairs ẑ∗
modeling closest to a given physical state z by minimizing a distance
The physical equations governing the deformation of an functional defined as follows:
elastic solid in a domain ΩX bounded by a Neumann bound- Z
∗
ary ΓX X
t and a Dirichlet boundary Γu are given as F(z, ẑ ) = min d2z (z, ẑ)dΩ, (2)
ẑ∈E ΩX
DIV (F(u) · S) + b = 0, in ΩX ,
with
E = E(u) = (F F − I)/2, in ΩX ,
T
(1)
(F(u) · S) · N = t, on ΓX
t , d2z (z, ẑ) = d2E (E(u), Ê) + d2S (S, Ŝ), (3a)
on ΓX 1
u = g, u , d2E (E(u), Ê) = (E(u) − Ê) : Ĉ : (E(u) − Ê), (3b)
2
where u is the displacement vector, E is the Green La- 1
grangian strain tensor, and S is the second Piola-Kirchhoff dS (S, Ŝ) = (S − Ŝ) : Ĉ−1 : (S − Ŝ),
2
(3c)
2
Copyright © 2021 for this paper by its authors.Use permitted under
Creative Commons License Attribution 4.0 International (CC BY where Ĉ is a predefined symmetric and positive-definite ten-
4.0). sor used to properly regulate the distance between z and ẑ.
Given the optimal material data ẑ∗ = (Ê∗ , Ŝ∗ ) obtained
from the local step in Eq. (2), the global step of the data-
driven problem is to search for the closest physical state,
formulated as follows:
min F(E(u), S; Ê∗ , Ŝ∗ )
u,S
subject to: DIV (F(u) · S) + b = 0 in ΩX , (4)
(F(u) · S) · N = t on ΓX
t ,
which can be solved by the Lagrange multiplier method and
a nonlinear numerical solver. Note that the compatibility
constraint in Eq. (1b) is directly enforced by computing E
from u.
2 Autoencoders for nonlinear material
manifold learning
Autoencoders (DeMers and Cottrell 1993; Hinton and
Salakhutdinov 2006) aim to optimally copy its input to out- Figure 1: Schematic of an autoencoder consisting of an en-
put and retain the most representative features in a low- coder and a decoder, with the embedding dimension smaller
dimensional embedding layer. Hence, it allows for effective than the input dimension. The encoder learns an intrinsic
noise filtering, dimensionality reduction, and hidden struc- low-dimensional embedding of a high-dimensional input
ture discovery of data. As shown in Fig. 1, an autoencoder object, whereas the decoder optimally reconstructs the input
contains an encoder function henc (·; θ enc ) : Rd → Rp and object from the low-dimensional embedding.
a decoder function hdec (·; θ dec ) : Rp → Rd , such that the
autoencoder is
x̃ = h(x; θ enc , θ dec ) := hdec (·; θ dec ) ◦ henc (x; θ enc ), (5) {z0 ∈ Rp |z0 = henc (z; θ ∗enc ), ∀z ∈ Z}, in which the material
state is described by a lower-dimensional coordinate system
where d is the input dimension, p < d is the embedding di- z0 . Here, the prime symbol (·)0 is used to denote the quanti-
mension, θ enc and θ dec are the trainable parameters of the ties defined in the embedding space, and Z denotes the high-
encoder and the decoder, respectively. x̃ is the output of the dimensional phase space where the material states ẑ and the
autoencoder, a reconstruction of the original input x. With physical states z are defined. For example, the embedding
p < d, the encoder henc is trained to learn an intrinsic repre- 0
set of the given material data is E = {ẑ0I }M
0
I=1 ⊂ E , where
sentation of x ∈ Rd , denoted as the embedding x0 ∈ Rp , 0 ∗
ẑI = henc (ẑI ; θ enc ) for ẑI ∈ E.
whereas the decoder hdec is trained to reconstruct the in-
put data by mapping the embedding x0 back to the high- 3 Auto-embedding data-driven (AEDD)
dimensional space x̃ ∈ Rd .
In this study, autoencoders are employed to discover the solver
intrinsic low-dimensional material embedding of the given The nonlinear physics-constrained data-driven computing
material dataset E = {ẑI }M framework described in Section 1 is conducted in the high-
I=1 , where ẑI = (ÊI , ŜI ).
The optimal parameters (θ ∗enc , θ ∗dec ) of the autoencoder dimensional phase space Z (called data space), with the
h(·; θ enc , θ dec ) are obtained by minimizing the following physical state zα ∈ C, the material data ẑα ∈ E, and the
loss function: material dataset E defined in Z. The subscript ”α” is used to
M
denote the quantities at integration points. To enhance solu-
1 X tion accuracy and generalization ability of data-driven com-
(θ ∗enc , θ ∗dec ) = arg min ||h(ẑI ; θ enc , θ dec ) − ẑI ||2
θ enc ,θ dec M puting, deep manifold learning enabled by autoencoders is
I=1
introduced into the material data-driven local solver. The au-
L+1
X toencoders are trained offline and the trained encoder henc
+β ||W(l) ||2F , and decoder hdec functions are employed directly in the on-
l=1 line data-driven computation. As such, the encoder maps an
(6) arbitrary point from the data space to the embedding space,
where L is the number of hidden layers, W is the weight i.e. z0α = henc (zα ), whereas the decoder performs the re-
coefficient of the autoencoder, β is a regularization param- verse mapping, i.e. z̃α = hdec (z0α ).
eter, and || · ||F denotes the Frobenius norm. The first term In the proposed AEDD local solver, the local step defined
in the loss function is the reconstruction error of all training in Eq. (2) is reformulated by three steps:
data and the second term is a L2 -norm based weight regular-
ization term used to avoid over-fitting (Goodfellow, Bengio, Step 1 : z0α = henc (zα ), (7a)
and Courville 2016). ẑ0∗ 0 0
Step 2 : α = I {ΨI (zα ); ẑI }I∈Nk (z0α ) (7b)
With the trained autoencoder h(·; θ ∗enc , θ ∗dec ), a low-
dimensional embedding space can be defined, i.e., E 0 = Step 3 : ẑ∗α = hdec (ẑ0∗
α ), (7c)
0
for α = 1, ..., Nint , where ẑ0I ∈ E , Nint is the number of
integration points, and I is the convexity-preserving inter-
polation operator defined as
X
z0recon = I {ΨI (z0 ); ẑ0I }I∈Nk (z0 ) = ΨI (z0 )ẑ0I ,
I∈Nk (z0 )
(8)
where z0recon is the reconstruction of z0 , ẑ0I is the material
data embedding in E0 , Nk (z0 ) is the index set of the k near-
est neighbor points of z0 selected from E0 . The interpolation
functions are
φ(z0 − ẑ0I )
ΨI (z0 ) = P 0 0 , (9)
J=1 φ(z − ẑJ )
where φ(z0 − ẑ0I ) = 1/||z0 − ẑ0I ||2 is a positive kernel func-
tion representing the weight on the data set {ẑ0I }I∈Nk (z0 ) .
The positive interpolation functions satisfy the partition of Figure 2: Geometric schematic of the proposed auto-
unity, I∈Nk (z0 ) ΨI (z0 ) = 1, which ensures transforma-
P
embedding data-driven computational framework. The ma-
tion objectivity and convexity of the interpolation scheme, terial data points (the gray-filled circles), ẑI , in the phase
Eq. (8). space are related to the material embedding points (the
The schematic of data-driven computing with the AEDD white-filled circles) ẑ0I via the encoder function. The low-
local solver is illustrated in Fig. 2, where the integration dimensional embedding manifold is represented by the or-
point index α is dropped for brevity. For example, at the ange dash line.
v-th global-local iteration, after the physical state z(v) (the
blue-filled triangle) is obtained from the global step (Eq.
(4)), Step 1 of the local solver (Eq. (7a)) maps the sought (10 and 11). The normal components of the Green strain and
physical state from the data space to the embedding space the associated 2nd-PK stress tenors of all protocols are plot-
by the encoder, z0(v) = henc (z(v) ), depicted by the white- ted in Fig. 3(c-d). Considering the symmetry of specimen
filled triangle in Fig. 2. In Step 2, k nearest neighbors of geometry and loading conditions, a quarter model with sym-
metric boundary conditions is modelled, as shown in Fig. 3b.
z0(v) based on Euclidean distance are sought in the em-
The normalized root-mean-square deviation (NRMSD) with
bedding space and the optimal material embedding solution true
respect to the maximum experimental stress data Smax is
ẑ0∗(v) (the red square) is reconstructed by using the proposed employed to assess the prediction performance of the meth-
convexity-preserving interpolation (Eqs. (8-9)). Lastly, in ods, see Eq. (10),
Step 3, the optimal material embedding state ẑ0∗(v) is trans-
formed from the embedding space to the data space by the v
decoder, ẑ∗(v) = hdec (ẑ0∗(v) ) (the red star in Fig. 2). Sub-
uNeval
u X (S AEDD − S true )2
i i true
sequently, this material state ẑ∗(v) from the local solver in NRMSD = t /Smax , (10)
N eval
Eq. (7) is used in the next physical state update z(v+1) . This i
process completes one global-local iteration. The iterations where Neval = 200 is the number of evaluation points,
proceed until the distance between the physical and mate- SiAEDD and Sitrue are the predicted stress and the experi-
rial states is within a tolerance, yielding the data-driven so- mental stress data at the i-th evaluation point, respectively.
lution denoted by the green star in Fig. 2, which ideally is It is observed that autoencoders with an embedding di-
the intersection between the physical manifold and material mension p ≤ 2 could not capture a meaningful embedding
manifold in the data space. representation for the material dataset with two-dimensional
Here, the nearest neighbors searching and locally convex strain-stress pairs (d = 6). This is consistent with the ob-
reconstruction of the optimal material state are processed in servation in (He and Chen 2020) that the number of nearest
the filtered (noiseless) low-dimensional embedding space, neighbors for local convex reconstruction should be greater
resulting in the enhanced robustness against noise and accu- than the intrinsic dimensionality of data, which is 2 for
racy of the local solution. the two-dimensional strain-stress material dataset. Hence,
the autoencoder with an architecture of ”6-4-3-4-6” is em-
4 Data-driven biological tissue modeling ployed, where the first and the last values are the input and
The effectiveness of the proposed AEDD computational output dimensions, respectively, and the remaining values
framework is evaluated by modeling biological heart valve denote the number of neurons in hidden layers in sequence.
tissue with data from biaxial mechanical experiments on one Thus, the embedding dimension equals to 3. A hyperbolic
representative porcine mitral valve posterior leaflet (Jett et tangent function is adopted as the activation function for
al. 2018), see a schematic in Fig. 3a. There are eleven pro- all layers of autoencoders, except for the embedding layer
tocols, including nine biaxial tension protocols (1–9) with and the output layer, where a linear function is employed
various biaxial tension ratios and two pure shear protocols instead. The regularization parameter β in the loss function
(a) (b)
(a) AEDD (b) LCDD
Figure 4: (a) AEDD prediction on Protocol 5; (b) LCDD
prediction on Protocol 5. Protocols 1, 3, 4, 7, and 8 are used
to train the autoencoder applied in AEDD
References
DeMers, D., and Cottrell, G. W. 1993. Non-linear dimen-
(c) (d) sionality reduction. In Advances in neural information pro-
cessing systems, 580–587.
Figure 3: (a) Schematic of a mitral valve posterior leaflet Duchi, J.; Hazan, E.; and Singer, Y. 2011. Adaptive subgra-
specimen mounted on a biaxial testing system; (b) schematic dient methods for online learning and stochastic optimiza-
of the model of biaxial testing in data-driven computation; tion. Journal of machine learning research 12(Jul):2121–
(c) experimental Green strain data of all protocols; (d) ex- 2159.
perimental 2nd-PK stress data of all protocols Goodfellow, I.; Bengio, Y.; and Courville, A. 2016. Deep
learning. MIT press.
He, Q., and Chen, J.-S. 2020. A physics-constrained data-
driven approach based on locally convex reconstruction for
Eq. (6) is set as 10−5 . An adaptive gradient algorithm, Ada- noisy database. Computer Methods in Applied Mechanics
grad (Duchi, Hazan, and Singer 2011), is employed. The ini- and Engineering 363:112791.
tial learning rate is 0.1 and the number of training epochs is
He, Q.; Laurence, D. W.; Lee, C.-H.; and Chen, J.-S. 2020.
2000. The training dataset contains the strain-stress data of
Manifold learning based data-driven modeling for soft bio-
protocols 1, 3, 4, 7, and 8, see Fig. 3(c-d), which are stan-
logical tissues. Journal of Biomechanics 110124.
dardized to have zero mean and unit variance for accelerated
training process. Hinton, G. E., and Salakhutdinov, R. R. 2006. Reducing
the dimensionality of data with neural networks. science
The autoencoder is trained offline using the open-source 313(5786):504–507.
Pytorch library (Paszke et al. 2017) and then applied in the Ibanez, R.; Abisset-Chavanne, E.; Aguado, J. V.; Gonzalez,
AEDD solver during online computation of protocol 5. The D.; Cueto, E.; and Chinesta, F. 2018. A manifold learning
data-driven simulation predicts the stress responses of the approach to data-driven computational elasticity and inelas-
model (Fig. 3b) under the displacement-controlled loading ticity. Archives of Computational Methods in Engineering
prescribed by the deformation history of protocol 5 (Fig. 25(1):47–57.
3c). As shown in Fig. 4, AEDD achieves a higher predic-
tion accuracy than the local convexity data-driven (LCDD) Jett, S.; Laurence, D.; Kunkel, R.; Babu, A. R.; Kramer, K.;
computing approach (He and Chen 2020; He et al. 2020), Baumwart, R.; Towner, R.; Wu, Y.; and Lee, C.-H. 2018. An
(NRMSDAEDD = 0.061 < NRMSDLCDD = 0.158). The investigation of the anisotropic mechanical properties and
results demonstrate better extrapolative generalization abil- anatomical structure of porcine atrioventricular heart valves.
ity of AEDD, which is attributed to the underlying low- Journal of the mechanical behavior of biomedical materials
dimensional global material manifold learned by the autoen- 87:155–171.
coders. Specifically, AEDD performs local neighbor search- Kirchdoerfer, T., and Ortiz, M. 2016. Data-driven computa-
ing and locally convex reconstruction of optimal material tional mechanics. Computer Methods in Applied Mechanics
state based on geometric distance information in the low- and Engineering 304:81–101.
dimensional global embedding space, which contains the Paszke, A.; Gross, S.; Chintala, S.; Chanan, G.; Yang, E.;
underlying manifold structure of the material data and con- DeVito, Z.; Lin, Z.; Desmaison, A.; Antiga, L.; and Lerer,
tributes to a higher solution accuracy and better generaliza- A. 2017. Automatic differentiation in pytorch.
tion performance. In contrast, LCDD performs local neigh-
bor searching and locally convex reconstruction purely from
the existing material data points without any generalization,
leading to very limited extrapolative generalization ability.