# Data-driven discovery of interpretable causal relations for deep learning
material laws with uncertainty propagation

Xiao Sun 111Department of Applied Mathematics and Statistics, Johns Hopkins
University, Baltimore, MD. Bahador Bahmani 222Department of Civil Engineering
and Engineering Mechanics, Columbia University, New York, NY. Nikolaos N.
Vlassis 2 WaiChing Sun 2 Corresponding author<EMAIL_ADDRESS>Yanxun Xu 1
Corresponding author<EMAIL_ADDRESS>
###### Abstract

This paper presents a computational framework that generates ensemble
predictive mechanics models with uncertainty quantification (UQ). We first
develop a causal discovery algorithm to infer causal relations among time-
history data measured during each representative volume element (RVE)
simulation through a directed acyclic graph (DAG). With multiple plausible
sets of causal relationships estimated from multiple RVE simulations, the
predictions are propagated in the derived causal graph while using a deep
neural network equipped with dropout layers as a Bayesian approximation for
uncertainty quantification. We select two representative numerical examples
(traction-separation laws for frictional interfaces, elastoplasticity models
for granular assembles) to examine the accuracy and robustness of the proposed
causal discovery method for the common material law predictions in civil
engineering applications.

##  1 Introduction

Computer simulations for mechanics problems often require material
(constitutive) laws that replicate the local constitutive responses of the
materials. These material laws can be used to replicate the responses of an
interface (e.g., traction-separation laws or cohesive zone models) or bulk
materials (e.g., elastoplasticity models for solids, porosity-permeability
relationship and water retention curve). A computer model is then completed by
incorporating these local constitutive laws into a discretized form of balance
principles (balance of mass, linear momentum and energy) where discretized
numerical solutions can be sought by a proper solver.

Constitutive laws, such as stress-strain relationship for bulk materials,
traction-separation laws for interface, porosity-permeability for porous
media, are often derived following a set of axioms and rules
(truesdell2004non). In these hand-crafted models, phenomenological
observations are incorporated into constitutive laws (e.g., critical state
theory for soil mechanics (schofield1968critical; sun2013unified;
bryant2019micromorphically; na2019configurational), void-growth theory for
ductile damage (gurson1977continuum)). While the earlier simpler models are
often amended by newer and more comprehensive models (dafalias1984modelling)
in order to improve the performance (e.g. accuracy, more realistic
interpretation of mechanisms), these improvements are often a trade-off that
may unavoidably increase the number of parameters, leading to increasing
difficulty for the calibration, verification and validation processes, as well
as the uncertainty quantification (dafalias1984modelling; borja2007estimating;
borja2008coseismic; clement2013uncertainty; wang2019meta;
wang2016identifying).

The rise of big data and the great promises of machine learning have led to a
new generation of approaches that either bypass the usages of constitutive
laws via model-free data-driven methods (e.g., kirchdoerfer2016data;
he2020physics; bahmani2021kd; karapiperis2021data) or replace parts of the
modeling efforts/components with models generated from supervised learning
(e.g., furukawa1998implicit; lefik2003artificial; wang2018multiscale;
wang2019updated; zhang2020ai; vlassis2020geometric; vlassis2021sobolev;
logarzo2021smart; masi2021thermodynamics). However, one critical issue of
these machine learning and data-driven approaches is the lack of sufficient
interpretability of predictions. While there is no universally accepted
definition of interpretability, we will herein employ the definition used in
miller2019explanation which refers interpretability as the degree to which a
human can understand the cause of the prediction.

One possible way to boost the interpretability is to introduce a proper medium
to represent knowledge that can be understood by human (sussillo2013opening).
Causal graph, also known as causal Bayesian network, is one such medium in
which the causal relations among different entities are mathematically
represented by a directed graph. In the application of computational mechanics
wang2018multiscale; wang2019updated; heider2020so have derived a mathematical
framework to decompose a complex prediction task into multiple easier
predictions represented by subgraphs within graphs. More recent work such as
wang2019meta; wang2019cooperative introduce a deep reinforcement learning
approach that employs the Monte Carlo Tree Search (MCTS) to assemble a
directed graph that generates a sequence of interconnected predictions of
physical quantities to emulate a hand-crafted constitutive law. However, these
directed graphs are generated to optimize a given performance metric (e.g.
accuracy, calculation speed and forward prediction robustness), but not
necessarily reveal the underlying causal relations among physical quantities.

Discovering causal relations from observational data is an important problem
with applications in many fields of science, such as social science
(oktay2010causal), finance (gong2017causal), and biomedicine
(shen2020challenges). The standard way to discover causality is through
randomized controlled experiments. However, conducting such experiments can be
either impractical, unethical, and/or very expensive in many disciplines
(xu2012bayesian; xie2020bayesian). For mechanics problems, the major issues
include the time and labor cost for physical experiments, the lack of
facilities or equipment to complete the required tests and the difficulties to
obtain specimens (mitchell2005fundamentals; powrie2018soil; wood1990soil). As
a result, an alternative approach, which is adopted in this study, is to use
sub-scale simulations as the digital representation that generates auxiliary
data sets to build material laws or forecast engine for the macroscopic
material responses (liu2016determining; ma2020computational;
frankel2019predicting; wang2019updated). Classic methods for causal discovery
are based on probabilistic graphical modeling (pearl2000causality), the
structure of which is a directed acyclic graph (DAG) with nodes representing
random variables and edges representing conditional dependencies between
variables. Learning a DAG from observational data is highly challenging since
the number of possible DAGs is super-exponential to the number of nodes. There
are two main approaches for causal discovery: the constraint-based approach
and the score-based approach. The constraint-based approach aims to recover a
Markov equivalence class through inferring conditional independence
relationships among the variables, and the resulting Markov equivalence class
may contain multiple DAGs that indicate the same conditional independence
relationships (le2016fast; cui2016copula). On the other hand, the score-based
approach uses a scoring function, such as the Bayesian Information Criterion
(BIC), to search for the DAG that best fits the data (heckerman1995learning;
huang2018generalized).

In this paper, we aim to discover causal relations that can explain the
underlying mechanism of a history-dependent macroscopic constitutive law
upscaled from direct numerical simulations at the meso-scale. The most common
method for constructing the causal relations from time-series data is Granger
causality (granger1969), which assumes a number of lagged effects and analyzes
the data in a unit no more than that number of lags. See runge2018causal for a
review of causal discovery methods on time-series data. However, most of these
causal discovery methods assume that the data are generated from a stationary
process, meaning that the data are generated by a distribution that does not
change with time. Such an assumption does not hold in many physical processes,
in which the mechanisms or parameters in the causal model may change over
time. Several methods have been proposed recently to tackle time-varying
causal relations in non-stationary processes (ghassami2018multi;
huang2019causal; huang2020causal). However, they either assume linear causal
models, or does not offer the flexibility of incorporating known physical
knowledge, limiting their applicability to the nonlinear path-dependent
relations in learning material constitutive laws.

In this work, we offer two major innovations. First, we introduce a new
decoupled discovery/training approach where the discovery of causal relations
represented by a DAG is enabled by a causal discovery algorithm that deduces
plausible causal relations from non-stationary time series data and
incorporates known physical knowledge. Second, we leverage the obtained causal
graph as the representation of mechanics knowledge and adopt a Bayesian
approximation using the dropout layer technique first introduced by
gal2016dropout to propagate epistemic uncertainty in the causal graph and
generate quantitative predictions with uncertainty quantification.

The rest of the paper is organized as follows. Section 2 first introduces the
two data sets (learning traction-separation law and hypo-plasticity of
granular materials) used for our numerical experiments. This is followed by
the description of theory and implementation of the proposed causal discovery
algorithm used to deduce the causal relations from non-stationary time series
data (Section 3). The setup of the deep neural network model for the
prediction tasks and the uncertainty propagation are included in Sections 4
and 5 respectively. The proposed framework is then tested against two
numerical experiments (Section 6), which is then followed by the conclusion.

##  2 Causal relations and constitutive laws

As demonstrated in previous studies such as wang2018multiscale; wang2019meta;
wang2019updated; wang2019cooperative; heider2020so; vlassis2020geometric, the
relationships in a constitutive model can be represented by a network of
unidirectional information flow, i.e., a DAG $G=({\bm{V}},\bm{E})$ where
${\bm{V}}$ represents a vertex set and $\bm{E}$ denotes an edge set. With
appropriate assumptions that will be discussed later, the DAG can be
identified as a causal graph (pearl2000causality). The causal relations are
not only useful to explain the underlying mechanism of a process but also
provide us a basis to formulate multi-step transfer learning to predict
constitutive responses. This strategy can be beneficial because one can
leverage more data gathered from physical numerical experiments to train the
prediction model. For instance, while a black-box prediction of stress-strain
curves only leverages the stress-strain pair for supervised learning, the
introduction of knowledge graphs may introduce multiple supervised learning
tasks where measurements of porosity, fabric tensors or any other physical and
geometrical attributes measured during the experiments can be leveraged to
improve the training. For completeness, we briefly describe the procedure to
consider the data set as vertex sets in graphs and the causal discovery
process used to create the directed edge set in a knowledge graph through two
examples.

Note that many of the physical quantities that become the vertices in the
knowledge graphs are graph metrics obtained from analyzing the connectivity
topology of the granular system. For brevity, we will not provide a review on
the applications of graph theory for granular matter here. Interested readers
may refer to Appendix LABEL:sec:graph_metrics for the definitions of the graph
metrics and satake1978constitution; bagi1996stress; walker2010topological;
tordesillas2010force; o2011particulate; kuhn2015stress and
vlassis2020geometric for reviews on the graph theory applied to particulate
and granular systems.

###  2.1 Dataset for traction-separation law

In the first example, our goal is to conduct a numerical experiment to verify
whether the causal discovery algorithm is able to re-discover the well-known
causal relation that links the plastic dilatancy and contraction to the
frictional behaviors (scholz1998earthquakes; popov2010contact) with a small
data set.

Following wang2018multiscale; wang2019meta; wang2019updated;
wang2019cooperative, we consider the vertex set consists of five elements, the
displacement jump/separation $\bm{U}$, the traction $\bm{T}$, and three
geometric measures, i.e.,

  1. 1.

Displacement jump $\bm{U}$, the relative displacement of an interface of two
in-contact bodies.

  2. 2.

Porosity $\phi$, the ratio between the volume of the void and the total volume
of RVE.

  3. 3.

Coordination number (averaged) $CN=N_{\text{contact}}/N_{\text{particle}}$
where $N_{\text{contact}}$ is the number of particle contacts and
$N_{particle}$ is the number of particles in the RVE.

  4. 4.

Fabric tensor
$\mathbf{A}_{f}={\frac{1}{N_{\text{contact}}}}\sum_{c=1}^{N_{\text{contact}}}\mathbf{n}^{c}\otimes\mathbf{n}^{c}$
where $\mathbf{n}^{c}$ is the normal vector of a particle contact $c$ in the
RVE. The symbol ‘$\otimes$’ denotes a juxtaposition of two vectors (e.g.,
$\bm{a}\otimes\bm{b}=a_{i}b_{j}$) or two symmetric second order tensors [e.g.,
$(\bm{\alpha}\otimes\bm{\beta})_{ijkl}=\alpha_{ij}\beta_{kl}$].

  5. 5.

Traction $\bm{T}$, the traction vector acts on the interface.

To generate a machine learning based traction-separation law, we identify the
displacement jump as the root and the traction as the leaf of the causal
graph. The causal graph is a DAG $G=({\bm{V}},\bm{E})$ where ${\bm{V}}$ is the
set consisting of the physical quantities $\bm{U},\phi,CN,\bm{A}_{f}$ and
$\bm{T}$. Meanwhile, $\bm{E}\subseteq{\bm{V}}\times{\bm{V}}$ is a set of
directed edges that connect any two elements from ${\bm{V}}$, and $\bm{E}$ is
determined from the causal discovery algorithm outlined in Section 3.

The dataset is generated using an open-source code YADE. In total, there are
100 traction-separation law simulations run with different loading paths
performed on the same RVE. This RVE consists of spherical particles with radii
between $1\pm 0.3$ mm with uniform distribution. The RVE has a height of 20 mm
in the normal direction of the frictional surface and is initially
consolidated to an isotropic pressure of 10 MPa. The inter-particle
interaction is controlled by Cundall’s elastic-frictional contact model
(cundall1979discrete) with an inter-particle elastic modulus of
$E_{eq}=1\,GPa$, a ratio between shear and normal stiffness of
$k_{s}/k_{n}=0.3$, a frictional angle of $\phi=30^{\circ}$, a density
$\rho=2600\,kg/m^{3}$, and a Cundall damping coefficient
$\alpha_{\text{damp}}=0.2$. For brevity, the generation and setup of the
simulations are not included in this paper. Interested readers please refer to
wang2019meta for more information. The data required to replicate the results
of this paper and for 3rd-party validation can be found in the Mendeley Data
repository (wang2019discrete).

###  2.2 Dataset for hypo-plasticity of granular materials

While the first data set is used to determine whether the causal graph
algorithm may re-discover known physical relations in the literature, the
second problem is designed to test whether the causal graph algorithm may
successfully investigate new plausible causal relations not known a prior in
the literature.

For this purpose, we run 60 discrete element simulations and use 30 of them
for calibrations and 30 for blind forward predictions. In addition to the
conventional microstructural attributes (e.g., porosity and fabric tensor)
typically used for hand-crafted constitutive laws (manzari1997critical;
dafalias2004simple; sun2013unified; bryant2019micromorphically;
na2019configurational), we have also recorded the evolution of the particle
contact pairs in each incremental time step of the discrete element
simulations. The particle contact connectivity is itself an undirected graph
$G_{\text{contact}}=({\bm{V}}_{\text{particle}},\bm{E}_{\text{contact}})$
where ${\bm{V}}_{\text{particle}}$ is the set of particles and
$\bm{E}_{\text{contact}}$ is the set of particle contacts, one for each
contact between two contacting convex particle represented. They are
undirected edges. To facilitate new discovery, we compute 15 different graph
metrics of $G_{\text{contact}}$ (see Appendix LABEL:sec:graph_metrics for
definition) that have not been used for composing constitutive laws and see if
(1) whether the causal discovery algorithm may discover causal relations among
these new physical quantities and (2) whether the new discovery helps improve
the accuracy, robustness and consistency of the forward predictions enabled by
neural networks trained according to the discovered causal relations.

In total, there are 11 types of time-history data in which 3 of them are
second-order tensors (strain, stress and the strong fabric tensor), and the
rest are scalar (porosity, coordination number graph density, graph local
efficiency, graph average clustering, graph degree assortativity coefficient,
graph transitivity and graph clique number). As such there are 11 elements in
the vertex set and the goal of the causal discovery is to establish the edge
set to complete the causal graph. A sequence of supervised learning is then
used to generate predictions via deep learning.

##  3 Causal discovery and knowledge graph constructions

###  3.1 Notations and assumptions

Let $G=({\bm{V}},\bm{E})$ be a DAG containing only directed edges and has no
directed cycles. For each $V_{i}\in{\bm{V}}$, let $\mathrm{PA}^{i}$ denotes
the set of parents of $V_{i}$ in $G$. Since our data are time-history
dependent, we assume that the joint probability distribution of ${\bm{V}}$ at
each time point according to $G$ can factorize as
$p({\bm{V}})=\prod_{i=1}^{m}p(V_{i}\mid\mathrm{PA}^{i}),$ where $m$ is the
number of vertices in $G$. Here $p(V_{i}\mid\mathrm{PA}^{i})$ can be regarded
as “causal mechanism.” For non-stationary time series data, the causal
mechanism $p(V_{i}\mid\mathrm{PA}^{i})$ can change over time, and the changes
may be due to the involved functional models or the causal strengths.

Throughout this section, we use the example of traction-separation law to
illustrate the proposed causal discovery method without loss of generality.
Therefore, ${\bm{V}}$ is the set consisting of displacement jump $\bm{U}$,
porosity $\phi$, coordination number $CN$, fabric tensor $\bm{A}_{f}$, and
traction $\bm{T}$. Zhangkun2019 developed a constraint-based causal discovery
algorithm for non-stationary time series data to identify changing causal
modules and recover the causal structure. In this paper, we extend the
algorithm in Zhangkun2019 such that the proposed causal discovery algorithm
not only handles non-stationary time-history data but also incorporates
certain physical constraints. For example, in constructing the traction-
separation law, we have the prior knowledge that the dynamic changes in
$\bm{U}$ can cause changes in other variables, not vice versa. Therefore, if
there exists a directed edge between the displacement jump $\bm{U}$ and any
other variable $V_{i}$, then $\bm{U}\rightarrow{V_{i}}$.

Denote ${\bm{V}}_{-\bm{U}}$ to include all other variables in ${\bm{V}}$
excluding $\bm{U}$ (e.g., porosity, fabric tensor). Since the causal mechanism
can change over time, we assume that the changes can be explained by certain
time-varying confounders, which can be written as functions of time. As we
have the prior knowledge that $\bm{U}$ itself is a time-dependent variable and
could affect all other variables, we regard $\bm{U}$ as such a confounder and
assume that the causal relation for each $V_{i}\in{\bm{V}}_{-\bm{U}}$ can be
represented by the following structural equation model:

$\displaystyle V_{i}$ $\displaystyle=$ $\displaystyle
g_{i}(\mathrm{PA}^{i},\theta_{i}(\bm{U}),\epsilon_{i}),$ (1)

where $\mathrm{PA}^{i}$ includes $\bm{U}$ if the changes in $\bm{U}$ can
affect the changes in $V_{i}$, $\theta_{i}(\bm{U})$ denotes a function of
$\bm{U}$ that influences $V_{i}$ as effective parameters, $\epsilon_{i}$ is a
noise term that is independent of $\bm{U}$ and $\mathrm{PA}^{i}$. The
$\epsilon_{i}$’s are assumed to be independent. As we treat $\bm{U}$ as a
random variable, there is a joint distribution over
${\bm{V}}\cup\\{\theta_{i}(\bm{U})\\}_{i:V_{i}\in{\bm{V}}_{-\bm{U}}}$. Denote
$G^{aug}$ to be the graph by adding
$\\{\theta_{i}(\bm{U})\\}_{i:V_{i}\in{\bm{V}}_{-\bm{U}}}$ to $G$, and for each
$i$, adding an arrow from $\theta_{i}(\bm{U})$ to $V_{i}$. Note that $G$ is
the induced subgraph of $G^{aug}$ over ${\bm{V}}$. Denote the joint
distribution of $G^{aug}$ to be $p^{aug}$.

In order to apply any conditional independence test on the variable set
${\bm{V}}$ for recovering causal structure, we set the following assumptions
(spirtes2000causation).

Assumption 1. (Causal Markov condition) $G^{aug}$ and the joint distribution
$p^{aug}$ on
${\bm{V}}\cup\\{\theta_{i}(\bm{U})\\}_{i:V_{i}\in{\bm{V}}_{-\bm{U}}}$ satisfy
the causal Markov condition if and only if a vertex of $G^{aug}$ is
probabilistically independent of all its non-descendants in $G^{aug}$ given
the set of all its parents.

Assumption 2. (Faithfulness) $G^{aug}$ and the joint distribution $p^{aug}$
satisfy the faithfulness condition if and only if no conditional independence
holds unless entailed by the causal Markov condition.

Assumption 3. (Causal sufficiency) The common causes of all variables in
${\bm{V}}\cup\\{\theta_{i}(\bm{U})\\}_{i:V_{i}\in{\bm{V}}_{-\bm{U}}}$ are
measured.

###  3.2 Recovery of the causal skeleton

In this section, we propose a constraint-based method building upon the PC
algorithm (spirtes2000causation) to first identify the skeleton of $G$,
defined as the obtained undirected graph if we ignore the directions of edges
in a DAG $G$. We prove that given Assumptions 1-3, we can apply conditional
independence tests to ${\bm{V}}$ to recover the skeleton of $G$. Algorithm 1
describes the proposed method, which is supported by Theorem 1. The proof is
provided in Appendix LABEL:sec:proof following Zhangkun2019.

1:Object: To obtain the undirected skeleton of $G$

2:Build a complete undirected graph $U_{G}$ with variables ${\bm{V}}$

3:for each node $V_{i}\in{\bm{V}}_{-\bm{U}}$ do

4: if $V_{i}$ and $\bm{U}$ are independent given a subset of
$\\{V_{k}|V_{k}\in{\bm{V}}_{-\bm{U}},k\neq i\\}$ then

5: Remove the edge between $V_{i}$ and $\bm{U}$

6:for every $V_{i},V_{j}\in{\bm{V}}_{-\bm{U}}$ do

7: if $V_{i}$ and $V_{j}$ are independent given a subset of $\\{V_{k}\mid
V_{k}\in{\bm{V}}_{-\bm{U}},k\neq i,k\neq j\\}\cup{\bm{U}}$ then

8: Remove the edge between $V_{i}$ and $V_{j}$

9:return $U_{G}$

Algorithm 1 Obtain the undirected skeleton of $G$

######  Theorem 1.

Given Assumptions 1-3, for every $V_{i},V_{j}\in{\bm{V}}_{-\bm{U}}$, $V_{i}$
and $V_{j}$ are not adjacent in $G$ if and only if they are independent
conditional on some subset of $\\{V_{k}\mid V_{k}\in{\bm{V}}_{-\bm{U}},k\neq
i,k\neq j\\}\cup\\{\bm{U}\\}$.

In lines 3-7 of Algorithm 1, we determine whether the changes in $\bm{U}$
cause changes in $V_{i}$. If not, $\bm{U}$ is not in the parent set of $V_{i}$
and there is no edge between $\bm{U}$ and $V_{i}$ in $G$. The lines 8-12 of
Algorithm 1 aims to identify the causal skeleton between variables in
${\bm{V}}$ except $\bm{U}$. Since how other variables change with $\bm{U}$ and
the relations between these variables are usually unknown and potentially very
complex, we use a nonparametric conditional independence test, kernel-based
condition independence (KCI) test developed by zhang2012kernel, to determine
the dependence between variables throughout this paper. This nonparametric
approach can not only capture the linear/nonlinear correlations between
variables by testing for zero Hilbert-Schmidt norm of the partial cross-
covariance operator, but also handle multidimensional data that are common in
mechanics problems.

###  3.3 Determination of causal directions

After obtaining the skeleton $U_{G}$, we need to determine the causal
directions of edges. meek1995causal provided a set of orientation rules to
determine the directions of undirected edges in a graph based on conditional
independence tests. However, the Meek rule (meek1995causal) is only applicable
to edges that satisfy its conditions. In this section, we first introduce the
Meek rule (meek1995causal), then propose an algorithm to orient the edges that
are not covered by the Meek rule after incorporating known physical knowledge.

Denote $\leftrightarrow$ to be an undirected edge. The Meek rule has the
following principles:

  1. 1.

For all triples $V_{i}\leftrightarrow{V_{j}}\leftrightarrow{V_{k}}$, if
$V_{i}$ and $V_{k}$ are marginally independent but conditionally dependent
given $V_{j}$, then $V_{i}\rightarrow V_{j}\leftarrow V_{k}$;

  2. 2.

If $V_{i}\rightarrow{V_{j}}\leftrightarrow{V_{k}}$ and there is no edge
between $V_{i}$ and $V_{k}$, then orient $V_{j}\rightarrow{V_{k}}$;

  3. 3.

If $V_{i}\rightarrow{V_{j}}\leftrightarrow{V_{k}}$ and there is an edge
between $V_{i}$ and $V_{k}$, then orient $V_{i}\rightarrow{V_{k}}$;

  4. 4.

If $V_{i}\rightarrow V_{j}\leftarrow V_{k}$, $V_{i}\leftrightarrow
V_{k}\leftrightarrow V_{k}$, and $V_{k}\leftrightarrow V_{j}$, then
$V_{k}\rightarrow V_{j}$.

Now we describe our algorithm on how to determine the edge directions in the
obtained skeleton $U_{G}$. Firstly, for any node $V_{i}$ adjacent to $\bm{U}$,
we orient $\bm{U}\rightarrow V_{i}$ due to the prior physical knowledge that
only $\bm{U}$ affects other variables, not vice versa. Then we apply the Meek
rule to the obtained graph after orienting the edges from $\bm{U}$ to its
neighbours. For instance, suppose $\bm{U}\rightarrow V_{i}\leftrightarrow
V_{j}$, if $V_{j}$ and $\bm{U}$ are independent given a set of variables
including $V_{i}$, then we orient $V_{i}\rightarrow V_{j}$; if $V_{j}$ and
$\bm{U}$ are independent given a set of variables excluding $V_{i}$, then we
have $V_{j}\rightarrow V_{i}$.

Next, we discuss how to determine the edge direction between two adjacent
variables if they are both adjacent to $\bm{U}$, i.e., $V_{i}\leftrightarrow
V_{j}$, $\bm{U}\rightarrow V_{i}$, and $\bm{U}\rightarrow V_{j}$, since such a
scenario is not covered by the Meek rule. The modularity property of causal
systems (pearl2000causality) demonstrated that if there are no confounders for
cause and effect, then $p(\mathrm{\it cause})$ and $p(\mathrm{\it
effect}\mid\mathrm{\it cause})$ are either fixed or change independently.
Based on this principle, since both $V_{i}$ and $V_{j}$ change with $\bm{U}$,
we can test the conditional independence between
$p(V_{i}\mid\theta_{i}(\bm{U}))$ and $p(V_{j}\mid V_{i},\theta_{j}(\bm{U}))$,
as well as between $p(V_{j}\mid\theta_{j}(\bm{U}))$ and $p(V_{i}\mid
V_{j},\theta_{i}(\bm{U}))$ to determine the direction between $V_{i}$ and
$V_{j}$. That says, if $p(V_{i}\mid\theta_{i}(\bm{U}))$ and $p(V_{j}\mid
V_{i},\theta_{j}(\bm{U}))$ are conditionally independent but
$p(V_{j}\mid\theta_{j}(\bm{U}))$ and $p(V_{i}\mid V_{j},\theta_{i}(\bm{U}))$
are not, then $V_{i}\rightarrow V_{j}$. Zhangkun2019 developed a kernel
embedding of non-stationary conditional distributions and extended the Hilbert
Schmidt Independence Criterion (HSIC, gretton2008kernel) to measure the
dependence between distributions, based on which the causal directions can be
determined. For example, if we have two random variables $V_{1}$ and $V_{2}$,
we can compute the dependence between $p(V_{1})$ and $p(V_{2}\mid V_{1})$
using the normalized HSIC, denoted by $\hat{\Delta}_{V_{1}\rightarrow V_{2}}$.
By the same token, we can compute the dependence between $p(V_{2})$ and
$p(V_{1}\mid V_{2})$ using the normalized HSIC, denoted by
$\hat{\Delta}_{V_{2}\rightarrow V_{1}}$. If $\hat{\Delta}_{V_{1}\rightarrow
V_{2}}<\hat{\Delta}_{V_{2}\rightarrow V_{1}}$, we orient $V_{1}\rightarrow
V_{2}$; otherwise $V_{2}\rightarrow V_{1}$. After orienting all possible
edges, we can get the Markov equivalent class of the DAG $G$. Algorithm 2
summarizes the proposed method on how to determine causal directions.

1:Object: To orient the directions in the causal skeleton $U_{G}$

2:Input: The undirected skeleton output from Algorithm 1

3:for any node $V_{i}$ adjacent to $\bm{U}$ do

4: Orient $\bm{U}\rightarrow{V_{i}}$

5:for all other undirected edges do

6: Apply the Meek rule

7:Find all nodes that are adjacent to $\bm{U}$ and have undirected edges with
other nodes, denoted by ${\cal S}$

8:if $\cal S$ is empty then

9: return $G$

10:else

11: repeat

12: for each node $V\in{\cal S}$ do

13: Consider the set $\cal Z$ of nodes that either are directed parents of $V$
or have undirected edges to $V$

14: Calculate the normalized HSIC using the node $V$ as the effect and the set
of nodes $\cal Z$ as the cause

15: Pick the node $V$ with the smallest normalized HSIC

16: Orient all edge directions from the nodes in $\cal Z$ to node $V$

17: Remove the node $V$ from ${\cal S}$

18: until $\cal S$ is empty

19:return $G$

Algorithm 2 Obtain the Markov equivalence class of the DAG $G$

######  Remark 1.

In the numerical examples, we setup a threshold inclusion probability ($20\%$)
below which the causality relation is not included in the hierarchical neural
network models. This treatment allows us to ensure that the causalities with
sufficient likelihoods are included but the less prominent relationship is
omitted to improve the efficiency and simplicity of the resultant model. This
threshold can be viewed as a hyperparameter. A highly threshold may yield a
DAG with less vertices and therefore reduce the total number of required
supervised training at the expense of being less precise on the causality
relations among the data.

##  4 Supervised learning for path-dependent material laws

Once causal relations are identified, a directed graph $G=({\bm{V}},\bm{E})$
can be established where there is an edge $e_{ij}\in\bm{E}$ from the node
$V_{i}\in{\bm{V}}$ to $V_{j}\in{\bm{V}}$ if $V_{i}$ is a direct cause of
$V_{j}$. Denote the leaf node to be the vertex that is not the cause to any
other vertices, the root node to be the vertex that is not the target of any
other vertices. Figure 1 demonstrates a directed graph indicating an
information flow how leaf node(s) is related to root node(s) via some
intermediate nodes, e.g., in Figure 1 $\\{V_{1},V_{2}\\}$, $\\{V_{6}\\}$, and
$\\{V_{3},V_{4},V_{5}\\}$ are sets of root, leaf, and intermediate nodes.

Figure 1: (a) shows a direct graph partitioned into three sub-graphs. (b)
indicates how each sub-graph is used in a separate supervised machine learning
task to predict the downstream node(s) from the upstream node(s).

Along with the similar idea introduced in wang2019meta, we aim to discover all
sub-graphs that sequentially pass the information from the root to leaf, see
Algorithm 3. Each of these subgraphs will contain leaf and root but without
any intermediate nodes. As such, supervised learning can help us train neural
network that predict the root of each subgraph with the corresponding leaf (or
leaves) as input(s). To identify all sub-graphs, we traverse the graph
backward from leaf to root nodes. The leaf with its immediate predecessors
formed a directed tree which will be added into the list of potential sub-
graphs. Then, we remove edges of the founded tree in the graph $G$. In the
next step, we select a new leaf node from the updated $G$ and do the process
just described until there is not any edges in $G$. For the case shown in
Figure 1 we have the following potential sub-graphs:

$\displaystyle
G_{a}=(\\{V_{2},V_{3},V_{5},V_{6}\\},\\{e_{26},e_{36},e_{56}\\}),$ (2)
$\displaystyle G_{b}=(\\{V_{1},V_{4},V_{5}\\},\\{e_{14},e_{15}\\}),$ (3)
$\displaystyle G_{c}=(\\{V_{2},V_{4},V_{5}\\},\\{e_{24},e_{25}\\}),$ (4)
$\displaystyle G_{d}=(\\{V_{1},V_{3}\\},\\{e_{13}\\}).$ (5)

Those sub-graphs (directed trees) that share common upstream nodes will be
merged into a bigger sub-graph. For the graph in Figure 1 final sub-graphs
are:

$\displaystyle
G_{1}=G_{a}=(\\{V_{2},V_{3},V_{5},V_{6}\\},\\{e_{26},e_{36},e_{56}\\}),$ (6)
$\displaystyle G_{2}=G_{b}\cup
G_{c}=(\\{V_{1},V_{2},V_{4},V_{5}\\},\\{e_{14},e_{15},e_{24},e_{25}\\}),$ (7)
$\displaystyle G_{3}=G_{d}=(\\{V_{1},V_{3}\\},\\{e_{13}\\}).$ (8)

For each sub-graph, we have a separate supervised machine learning (ML) task.
In each ML task, inputs and outputs features are upstream and downstream nodes
in each directed sub-graph, as shown in Figure 1(b).

1:Input: Directed graph $G=({\bm{V}},\bm{E})$ $\triangleright$ result of
causal graph

2:$\mathcal{S}\leftarrow\emptyset$ $\triangleright$ a set of 2-tuples (input
nodes, output nodes) for ML tasks

3:while $\bm{E}\neq\emptyset$ do

4: $V_{l}\leftarrow$ get a leaf node of $G$ $\triangleright$ if $G$ has
multiple leaf nodes returns a random leaf

5: $\bar{{\bm{V}}}\leftarrow\\{V_{l}\\}$

6: $\hat{{\bm{V}}}\leftarrow\\{V_{i}|e_{il}\in\bm{E}\\}$

7: $\mathcal{S}\leftarrow\mathcal{S}\cup\\{(\hat{{\bm{V}}},\bar{{\bm{V}}})\\}$
$\triangleright$ $\hat{{\bm{V}}}$ is input node(s), and $\bar{{\bm{V}}}$ is
the output node for a potential ML task

8: $\bm{E}\leftarrow\bm{E}\setminus\\{e_{il}|V_{i}\in\hat{{\bm{V}}}\\}$

9:Modification: Elements $i$ and $j$ of $\mathcal{S}$ known as
$(\hat{{\bm{V}}}_{i},\bar{{\bm{V}}}_{i})$ and
$(\hat{{\bm{V}}}_{j},\bar{{\bm{V}}}_{j})$, respectively, are merged into one
element as $(\hat{{\bm{V}}}_{i},\bar{{\bm{V}}}_{j}\cup\bar{{\bm{V}}}_{j})$ if
$\hat{{\bm{V}}}_{i}\equiv\hat{{\bm{V}}}_{j}$.

10:return $\mathcal{S}$

Algorithm 3 Obtain all supervised learning input-output pairs

In the training step, we use the same architecture for all ML tasks consisting
of five layers GRU-Dropout-GRU-Dropout-Dense. Each GRU layer has 32 neurons,
and the linear activation function is used for the Dense layer. Gated
Recurrent Unit (GRU) is one type of Recurrent Neural Networks (RNN) to model
history dependence (chung2014empirical). We use the GRU to strengthen the
robustness of the ML black-box in dealing with any possible path-dependence
(wang2019cooperative). The dropout rate in the GRU is controlled to quantify
uncertainty in the model prediction, which will be detailed in Section 5. All
the input and output feature columns for each supervised learning task are
normalized to zero mean and unity standard deviation. The loss function is
defined as the mean squared root error between output prediction and ground
truth. This loss is minimized by the Adam optimizer inside Keras library with
Tensorflow 2.0 as its backend. Optimization process is done by the mini-batch
stochastic gradient descent algorithm with a batch size 256 during 1000
epochs. The described neural network architecture and hyperparameters are
chosen to be as close as possible to the ones used in wang2019meta.

Note that, the tuning of the hyperparameters (e.g. number of neurons, number
of layers, type of activation) may have a significant effect on the
performance of the neural network models. The best combination of
hyperparameters can be estimated via a variety of approaches such as the
greedy search, random search (bergstra2012random), random forest
(probst2019hyperparameters), Bayesian optimization (klein2017fast), meta-
gradient iteration or deep reinforcement learning (wang2020non; Fuchs:2021aa).
In this work, we adopt the random search approach in bergstra2012random to
fine-tune the hyperparameters (cf. Sec. 6.1) A rigorous hyperparameter study
that compares different hyperparameter tuning for neural networks that
generates constitutive laws may provide further insights on the optimal setup
of the hyperparameters but is out of the scope of this study.

For the blind prediction, after training, we start from the root (e.g.,
$\bm{U}$ in the traction-separation law) and sequentially predict intermediate
nodes via their corresponding sub-graph trained neural networks (NN) until
reaching the leaf node. For example in the case shown in Figure 1: NN 3
predicts $V_{3}$ from input $V_{1}$; $V_{4}$ and $V_{5}$ are predicted by NN 2
from inputs $V_{1}$ and $V_{2}$; and finally NN 1 is used to predict the
target variable $V_{6}$ from input $V_{2}$ and the obtained intermediate nodes
$V_{3}$ and $V_{5}$.

##  5 Uncertainty propagation in causal graph with dropout layers

As described in Section 4, we use the deep learning method, GRU in the
training step and prediction to handle path-dependent predictions. However,
the GRU itself is not designed to capture prediction uncertainty, which is of
crucial importance in learning material law. In machine learning and
statistics, Bayesian methods as probabilistic models provide us a natural way
to quantify the model uncertainty through computing the posterior distribution
of unknown parameters (xie2020bayesian). However, these methods often suffer
from a prohibitively high computational cost. In this paper, we show that the
dropout technique (srivastava2014dropout) used in the GRU can quantify
uncertainty in prediction as a Bayesian approximation.

Dropout, a regularization method that randomly masks or ignores neurons during
training, has been widely used in many deep learning models to avoid over-
fitting and improve prediction (hinton2012improving; li2016improved;
boluki2020learnable). gal2016dropout firstly prove the link between dropout
and a well known probabilistic model, the Gaussian process
(rasmussen2003gaussian), and show that the use of dropout in the feed forward
neural networks can be interpreted as a Bayesian approximation of Gaussian
processes. In the context of RNNs, gal2016theoretically treated RNNs as
probabilistic models by assuming network weights as random variables with a
Gaussian mixture prior (with one component fixed at zero with a small
variance). Such a technique is similar to the spike-and-slab prior in Bayesian
statistics for variable selection (ishwaran2005spike). Then
gal2016theoretically show that optimizing the objective in the variational
inference (blei2017variational) for approximating the posterior distribution
over the weights is equivalent to conducting dropout in the respective RNNs,
and demonstrate the implementation in one commonly-used RNN model, the long
short-term memory (hochreiter1997long). In this section, we propose to extend
the technique developed in gal2016theoretically in the context of GRUs for
uncertainty quantification (UQ) in the prediction.

Given training inputs ${\bm{X}}$ and the corresponding output ${\bm{Y}}$,
suppose that we aim to predict an output ${\bm{y}}^{*}$ for a new input
${\bm{x}}^{*}$. From the Bayesian point of view, the prediction uncertainty
can be characterized by the posterior predictive distribution of
${\bm{y}}^{*}$ as follows:

$p\left({\bm{y}}^{*}\mid{\bm{x}}^{*},{\bm{X}},{\bm{Y}}\right)=\int
p\left({\bm{y}}^{*}\mid{\bm{x}}^{*},{\bm{\omega}}\right)p({\bm{\omega}}\mid{\bm{X}},{\bm{Y}})\mathrm{d}{\bm{\omega}},$
(9)

where ${\bm{\omega}}$ includes all unknown model parameters,
$p({\bm{\omega}}\mid{\bm{X}},{\bm{Y}})$ is the posterior distribution of
${\bm{\omega}}$. In the GRU, all unknown weights can be viewed as
${\bm{\omega}}$. As the posterior distribution
$p({\bm{\omega}}\mid{\bm{X}},{\bm{Y}})$ is generally intractable, the
variational inference method approximates it by proposing a variational
distribution $q({\bm{\omega}})$ and then finding the optimal parameters in the
variational distribution through minimizing the Kullback-Leibler (KL)
divergence between the approximating distribution and the full posterior
distribution:

$\mathrm{KL}(q(\bm{\omega})\|p(\bm{\omega}\mid{\bm{X}},{\bm{Y}}))\propto-\int
q(\bm{\omega})\log
p\left({\bm{Y}}\mid{\bm{X}},{\bm{\omega}}\right)\mathrm{d}{\bm{\omega}}+\mathrm{KL}(q(\bm{\omega})\|p(\bm{\omega})),$
(10)

where $p(\bm{\omega})$ is the prior distribution of ${\bm{\omega}}$.

Given an input sequence ${\bm{X}}=[{\bm{x}}_{1},\dots,{\bm{x}}_{T}]$ of length
$T$, the hidden state ${\bm{h}}_{t}$ at time step $t$ in the GRU neural
network can be generated as follows:

$\displaystyle{\bm{z}}_{t}$
$\displaystyle=\sigma\left({\bm{W}}_{{\bm{z}}}{\bm{x}}_{t}+{\bm{U}}_{{\bm{z}}}{\bm{h}}_{t-1}+{\bm{b}}_{{\bm{z}}}\right),$
(11) $\displaystyle{\bm{r}}_{t}$
$\displaystyle=\sigma\left({\bm{W}}_{{\bm{r}}}{\bm{x}}_{t}+{\bm{U}}_{{\bm{r}}}{\bm{h}}_{t-1}+{\bm{b}}_{r}\right),$
$\displaystyle{\bm{i}}_{t}$
$\displaystyle=\tanh\left({\bm{W}}_{{\bm{h}}}{\bm{x}}_{t}+{\bm{U}}_{{\bm{h}}}\left({\bm{r}}_{t}\odot{\bm{h}}_{t-1}\right)+{\bm{b}}_{{\bm{h}}}\right),$
$\displaystyle{\bm{h}}_{t}$
$\displaystyle={\bm{z}}_{t}\odot{\bm{h}}_{t-1}+\left(1-{\bm{z}}_{t}\right)\odot{\bm{i}}_{t},$

where $\sigma$ denotes the sigmoid function and $\odot$ denotes the element-
wise product. Also, we assume that the model output at time step $t$ can be
written as
$f_{{\bm{Y}}}({\bm{h}}_{t})={\bm{h}}_{t}{\bm{W}}_{{\bm{Y}}}+{\bm{b}}_{{\bm{Y}}}$.
Then the unknown parameters in the GRU are
$\bm{\omega}=\left\\{{\bm{W}}_{{\bm{Y}}},{\bm{W}}_{{\bm{z}}},{\bm{U}}_{{\bm{z}}},{\bm{W}}_{{\bm{r}}},{\bm{U}}_{{\bm{r}}},{\bm{W}}_{{\bm{h}}},{\bm{U}}_{{\bm{h}}},{\bm{b}}_{{\bm{Y}}},{\bm{b}}_{{\bm{z}}},{\bm{b}}_{{\bm{r}}},{\bm{b}}_{{\bm{h}}}\right\\}$.
We write
${\bm{h}}_{t}=f^{\bm{\omega}}_{{\bm{h}}}\left({\bm{x}}_{t},{\bm{h}}_{t-1}\right)$
and $f^{{\bm{\omega}}}_{{\bm{Y}}}$ for the output in order to make the
dependence on ${\bm{\omega}}$ clear.

Then the right hand of (10) can be written as follows:

$\displaystyle-\int q(\bm{\omega})\log p\left({\bm{Y}}\mid
f_{{\bm{Y}}}^{{\bm{\omega}}}\left({\bm{x}}_{1},\dots,{\bm{x}}_{T},f_{{\bm{h}}}^{{\bm{\omega}}}\left({\bm{x}}_{T},f_{{\bm{h}}}^{{\bm{\omega}}}(\dots
f_{{\bm{h}}}^{{\bm{\omega}}}({\bm{x}}_{1},{\bm{h}}_{0})\dots)\right)\right)\right)\mathrm{d}{\bm{\omega}}+\mathrm{KL}(q(\bm{\omega})\|p(\bm{\omega})),$
(12)

which can be approximated by Monte Carlo integration with the generated
samples $\hat{{\bm{\omega}}}^{b}\sim q({\bm{\omega}})$ and plug in the sampled
$\hat{{\bm{\omega}}}^{b}$’s to (12).

Following gal2016theoretically, we use a mixture of Gaussian distributions as
the variational distribution for every weight matrix row ${\bm{\omega}}_{k}$:

$q({\bm{\omega}})=\prod_{k=1}^{K}q({\bm{\omega}}_{k}),\ \
q\left({\bm{\omega}}_{k}\right)=\pi\mathcal{N}\left({\bm{\omega}}_{k};\mathbf{0},\tau^{2}I\right)+(1-\pi)\mathcal{N}\left({{\bm{\omega}}}_{k};{\bm{m}}_{k},\tau^{2}I\right),$
(13)

where $\pi$ is the dropout probability, ${\bm{m}}_{k}$ is the variational
parameter (row vector), and $\tau^{2}$ is a small variance. We optimize over
${\bm{m}}_{k}$ by minimizing the KL divergence in (12). Sampling each row of
$\hat{{\bm{\omega}}}^{b}$ is equivalent to randomly mask rows in each weight
matrix, i.e., conducting dropout. Then the predictive posterior distribution
can be approximated by

$p\left({\bm{y}}^{*}\mid{\bm{x}}^{*},{\bm{X}},{\bm{Y}}\right)=\int
p\left({\bm{y}}^{*}\mid{\bm{x}}^{*},{\bm{\omega}}\right)p({\bm{\omega}}\mid{\bm{X}},{\bm{Y}})\mathrm{d}{\bm{\omega}}\approx\frac{1}{B}\sum_{b=1}^{B}p({\bm{y}}^{*}\mid{\bm{x}}^{*},\hat{{\bm{\omega}}}^{b}),$
(14)

where $\hat{{\bm{\omega}}}^{b}\sim q({\bm{\omega}})$ and $B$ is the total
number of generated samples.

To implement the dropout in the GRU, we re-parametrize (11) as follows:

$\left(\begin{array}[]{c}{\bm{z}}_{t}\\\ {\bm{r}}_{t}\\\
{\bm{i}}_{t}\end{array}\right)=\left(\begin{array}[]{c}\sigma\\\ \sigma\\\
\tanh\end{array}\right)\left(\left(\begin{array}[]{c}{\bm{x}}_{t}\circ{\bm{m}}_{{\bm{x}}}\\\
{\bm{h}}_{t-1}\circ{\bm{m}}_{{\bm{h}}}\end{array}\right)\cdot{\bm{\omega}}\right),$
(15)

where ${\bm{m}}_{{\bm{x}}}$ and ${\bm{m}}_{{\bm{h}}}$ denote randomly masks
repeated at all time steps.

We take the prediction tasks in Fig. 1 for example. The posterior predictive
distribution of $V_{3}$ given $V_{1}$ is straightforward by (14) since only
root and leaf nodes are involved. When involving the intermediate nodes, e.g.,
the subgroup $G_{1}$, the posterior predictive distribution of $V_{6}$ given
the root nodes ($V_{1}$ and $V_{2})$ can be computed by

$\displaystyle p\left(V_{6}\mid V_{1},V_{2}\right)=\int p\left(V_{6}\mid
V_{2},V_{3},V_{5},{\bm{\omega}}\right)p(V_{3}\mid
V_{1},{\bm{\omega}})p(V_{5}\mid V_{1},V_{2},{\bm{\omega}})p({\bm{\omega}}\mid
V_{1},V_{2})\mathrm{d}V_{3}\mathrm{d}V_{5}\mathrm{d}{\bm{\omega}}.$ (16)

Specifically, when we generate $B$ samples from the posterior distribution of
${\bm{\omega}}$, we also generate $B$ samples for the intermediate nodes.
Those $B$ samples of ${\bm{\omega}}$ and the corresponding intermediate nodes
can be then used in Monte Carlo integration to calculate (16). Throughout this
paper, we set $B$=200.

##  6 Numerical Examples

In this section, we conduct two numerical experiments to test our proposed
framework that combines causal discovery with deep learning to build
constitutive laws for granular materials. In Section 6.1, the causal discovery
is conducted to determine the constitutive relationships for an RVE interface
composed of spherical grains. The subsequent supervised machine learning then
leverages the causal relations learned from the causal discovery algorithm to
establish a serial of supervised learning that constitutes a forecast engine
for traction. The propagation of uncertainty is enabled by the dropout
technique that approximates Monte Carlo simulations to determine the
confidence intervals for a given dropout rate. In Section
LABEL:sec:bulkplasticity, the same exercise is repeated for another data set
to generate hypoplasticity surrogate model for a discrete element assembly
where new topological measures are computed and incorporated into the proposed
framework to (1) discover new physical mechanisms and (2) determine the
benefit of the new discovery on the accuracy, robustness, and consistency of
the forward predictions on unseen events.

###  6.1 Numerical Example 1: Machine Learning traction-separation law

Traction-separation laws are known as one of the main ingredients of cohesive
fracture models used for brittle materials (pandolfi2000three;
park2011cohesive). Generally, a traction-separation law constitutes a relation
between traction and displacement jump fields over the fracture surface. There
exist many hand-crafted models developed by experts for different applications
(park2011cohesive), while no unified framework had been developed until
recently by wang2019meta who suggest an approach based on reinforcement
learning. Also, in some applications such as granular materials more
descriptors, e.g., porosity or fabric tensor, should be considered in these
constitutive laws (sun2013multiscale) to derive more predictive models. Lack
of robustness in adding more descriptors is another weakness of classical
models.

Our first test for the data-driven causal discovery model with dropout UQ is
on the traction-separation law data publicly available in the repository
Mendeley data (cf. wang2019discrete). This dataset has also been used in
wang2019meta where the traction-separation law is determined from
reinforcement learning. Our major point of departure is three-fold. Firstly,
we develop a causal discovery algorithm to identify causal relations among
history-dependent physical quantities in RVE simulations. Secondly, we
decouple the causal discovery from the training of the neural network such
that we now first discover causal relations, then utilize the discovered
relationships to generate quantitative predictions using the method detailed
in section 4. Thirdly, we introduce the Bayesian approximation using the
dropout technique to propagate the uncertainty in the causal graph, building
upon the theoretical framework established in gal2016theoretically.

The database includes 100 DEM experiments. In each DEM experiment, the time
history of all the variables included in the DAG are recorded. Each experiment
is conducted by a different ratio of normal to tangential loading rate and
loading-unloading cycles on the same representative volume element of granular
materials. As such, the total number of time-history data points in these
experiments may vary from 51 to 111.

The interested reader is referred to the appendix in wang2019meta for more
information. In our study, feature space consists of displacement jump vector,
traction vector, coordination number, symmetric part of fabric tensor, and
porosity. We use half of the experiments for causal discovery and training
artificial neural networks, and the rest is used for test and validation. In
the causal discovery step, as different experimental setups may lead to
different causal relations among variables, we apply the proposed causal
discovery algorithm (Algorithms 1 and 2 in Section 3) to each of the training
experiment, and then report the final causal graph by calculating the
inclusion probabilities of directed edges appearing in all training
experiments. The inclusion probability of one edge is defined as the
proportion of causal graphs containing this edge. The directed edges with
inclusion probabilities being larger than a pre-defined threshold (20% in our
paper) are kept in the final causal graph. When both edge directions between
two variables appear with positive inclusion probabilities (e.g.,
$V_{i}\rightarrow V_{j}$ and $V_{j}\rightarrow V_{i}$ both exist), we keep the
edge direction that has a higher inclusion probability. The goal of the
resultant model is to predict the same granular assembly responds to a
different cyclic loading path unseen in the training. As such, the focus of
this model is to generate a surrogate for one representative element volume.

Fig. 2 plots the final causal graph on the training data sets with edge
inclusion probabilities. The strong confidence ($96\%$) in the edge starting
from the displacement jump vector to porosity is consistent with the common
field knowledge, i.e., the immediate consequence of displacement jump is the
volume change.

The displacement jump vector, as the only control variable, affects all the
intermediate physical quantities and traction vector. This observation may
seem to be trivial, but it is not always the case which will be shown in the
next example. The causal effects of fabric and coordination number on traction
is aligned with the modern Critical State Theory (li2012anisotropic) which is
obtained without expert interpretation by the causal algorithm. Note that
fabric encodes microstructural information in more detail such as directional
dependence due to its tensorial nature, rather than porosity which smears out
information into one scalar quantity. Therefore, it is reasonable to see that
fabric has a considerable contribution in describing material behavior with a
complex arrangement of force chains at the microstructural level.

Figure 2: Final causal graph for the traction-separation law deduced from
time-history of displacement, traction, porosity, coordination number, and
fabric tensor. The number on each edge represents the edge inclusion
probabilities among all possible causal relations from the training data sets.