Ballester-Ripoll and Leonelli

Global Sensitivity Analysis in Probabilistic Graphical Models

Global Sensitivity Analysis in Probabilistic Graphical Models

Rafael Ballester-Ripoll School of Human Sciences & Technology, IE University,
Madrid<EMAIL_ADDRESS>Manuele Leonelli School of Human
Sciences & Technology, IE University, Madrid, Spain<EMAIL_ADDRESS>
We show how to apply Sobol’s method of global sensitivity analysis to measure
the influence exerted by a set of nodes’ evidence on a quantity of interest
expressed by a Bayesian network. Our method exploits the network structure so
as to transform the problem of Sobol index estimation into that of
marginalization inference. This way, we can efficiently compute indices for
networks where brute-force or Monte Carlo based estimators for variance-based
sensitivity analysis would require millions of costly samples. Moreover, our
method gives exact results when exact inference is used, and also supports the
case of correlated inputs. The proposed algorithm is inspired by the field of
tensor networks, and generalizes earlier tensor sensitivity techniques from
the acyclic to the cyclic case. We demonstrate the method on three medium to
large Bayesian networks that cover the areas of project risk management and
reliability engineering.

Bayesian networks, Global sensitivity analysis, Sobol indices, Tensor
networks, Uncertainty quantification

##  1 Introduction

The assessment of the effect of the choice of inputs on the outputs of a
mathematical model is a fundamental but often overlooked step of any real-
world analysis. The set of techniques to perform such a critical task is
usually referred to under the broad name of _sensitivity analysis_. There is
an increasing recognition of the importance of performing sensitivity analysis
after the definition of a mathematical model and there are now a wide array of
recent reviews of the methods available (e.g. Borgonovo and Plischke 2016,
Borgonovo 2017, Razavi et al. 2021, Saltelli et al. 2019).

Probabilistic graphical models and specifically Bayesian networks (BNs) (e.g.
Darwiche 2009, Koller et al. 2009) are a class of models that are widely used
for risk assessment of complex operational systems in a variety of domains
(e.g. Bielza and Larrañaga 2014, Cai et al. 2018, Drury et al. 2017, McLachlan
et al. 2020). They provide an efficient and intuitive framework to represent
the relationship between a vector of variables of interest by simply using a
graph. Sensitivity analysis is increasingly becoming part of any real-world BN
modeling (e.g. Hänninen and Kujala 2012, Kleemann et al. 2017, Makaba et al.
2020) and general-purpose software that implements it is now available
(Leonelli et al. 2021, Darwiche 2010).

Sensitivity analysis in BNs can be broken down into two main steps. First,
some parameters of the model are varied and the effect of these variations on
output probabilities of interest are investigated. For this purpose, a simple
mathematical function, usually termed _sensitivity function_ , describes an
output probability of interest as a function of the BN parameters (Castillo et
al. 1997, Coupé and Van Der Gaag 2002). Second, once parameter variations are
identified, their effect is summarized by a distance or divergence measure
between the original and the varied distributions underlying the BN, most
commonly the Chan-Darwiche distance (Chan and Darwiche 2005) or the well-known
Kullback-Leibler divergence.

As noticed by Rohmer (2020) all the above-mentioned methods are local since
they most often measure the effect of a parameter variation on an output of
interest, while the other parameters are kept fixed. Multi-way sensitivity
analyses, where more than one parameter is varied contemporaneously, have been
proposed (e.g. Bolt and Renooij 2014, Chan and Darwiche 2004, Kjærulff and van
der Gaag 2000, Leonelli et al. 2017, Leonelli and Riccomagno 2019). However,
these quickly become computationally intractable.

Outside of the realm of BNs, the most common approach is the so-called global
sensitivity analysis (Saltelli et al. 2000, 2008) which provides a “global”,
instead of local, representation of how different factors jointly influence
some function of the model’s output. General GSA is computationally intensive
and remains largely unused for probabilistic graphical models. Here we bridge
this gap and contribute a practical computation of Sobol indices (Sobol 1990)
for variables associated to vertices of a network. To our knowledge, the only
other attempt at introducing global sensitivity methods for BNs is due to Li
and Mahadevan (2017). However, their approach suffers of a number of
drawbacks: i. it is computationally very expensive and requires a complex
Monte Carlo simulation; ii. its implementation is not freely available; iii.
it can only compute Sobol indices for one vertex of the network.

Our method takes advantage of the similarities between tensor networks (TNs;
see e.g. Ye and Lim 2018) and probabilistic graphical models and is part of an
ongoing effort to foster interactions between these two communities (Robeva
and Seigal 2018). Like probabilistic networks, TNs have long had an important
role in knowledge discovery and mathematical modeling, and are experiencing a
surge in popularity along with the latest trends in machine learning
(Stoudenmire and Schwab 2016). Although topologies admitting cycles are much
more expressive, much of the extant TN literature focuses on learning acyclic
structures from observable multidimensional arrays (tensors). Further, the
Sobol indices of an acyclic tensor networks can be computed efficiently,
especially in the particular case of the _linear tensor network_ or _tensor
train_ (Rai 2014, Zhang et al. 2015, Ballester-Ripoll et al. 2018, 2019). Our
contribution is _to make exact Sobol sensitivity analysis possible in the case
where the function of interest is given by a Bayesian network_ ; note that
such networks can only be cast as tensor networks by introducing cycles (via
moralization).

Importantly, by taking advantage of the underlying structure of the BN, the
methods proposed here can formally account for correlated inputs, although
some restrictions need to be imposed. Global sensitivity methods usually
ignore input correlations which can greatly falsify the outcome of the
sensitivity analysis (Do and Razavi 2020), although some recent proposals
attempt to embed correlated inputs (e.g. Chastaing et al. 2012, Iooss and
Prieur 2019, Sheikholeslami et al. 2020).

The implementation of our methods is freely available in Python (available at
https://github.com/rballester/bnsobol) contributing to the recent effort of
promoting global sensitivity analysis (Douglas-Smith et al. 2020).

###  1.1 Notation

We use the following throughout the paper:

  * •

The capital $Y$ represents random variables, $y$ their realizations and
$\mathbb{Y}$ their sample spaces.

  * •

Bold letters for vectors, e.g. $\boldsymbol{y}$, $\boldsymbol{Y}$, and
blackboard letters for sets, e.g. $\mathbb{A}$, $\mathbb{B}$;
$\mathcal{P}(\mathbb{A})$ is the power set of $\mathbb{A}$.

  * •

$[n]$ denotes $\\{1,\dots,n\\}$.

  * •

For a subset $\mathbb{A}\subset[n]$, we denote
$\boldsymbol{Y}_{\mathbb{A}}=(Y_{i})_{i\in\mathbb{A}}$ and
$\mathbb{Y}_{\mathbb{A}}=\times_{i\in\mathbb{A}}\mathbb{Y}_{i}$.

  * •

Calligraphic letters are used for graphical models: $\mathcal{B},\mathcal{M}$
and $\mathcal{T}$ for Bayesian networks, Markov random fields and tensor
networks, respectively.

  * •

If $\mathcal{T}$ is a graphical model and $\mathbb{A}\subset[n]$,
$\mathcal{T}(\boldsymbol{y}_{\mathbb{A}})$ is the model that results after
setting $\boldsymbol{Y}_{\mathbb{A}}=\boldsymbol{y}_{\mathbb{A}}$.

  * •

A subscript, as in $\mathcal{T}_{\mathbb{A}}$, denotes variable elimination,
i.e. removing $\boldsymbol{Y}_{\mathbb{A}}$ from a model $\mathcal{T}$.

  * •

The backslash ($\setminus$) denotes the absence of an element in a set. We use
it whenever the set in question is clear from the context; for example, for a
network $\mathcal{T}$ on variables $[n]$, $\mathcal{T}_{\setminus i}$ denotes
removing $\boldsymbol{Y}_{[n]\setminus\\{i\\}}$.

##  2 Bayesian Networks and Graphical Models

BNs are one of the most common statistical methods to investigate the
dependence structure for a random vector of interest and to efficiently answer
inferential queries. Next we review the theory of BNs and discuss their
connection to more general probabilistic graphical models.

###  2.1 Bayesian Networks

Consider a random vector $\boldsymbol{Y}=(Y_{1},\dots,Y_{n})$ of interest,
where for the purposes of this paper we assume $Y_{i}$ to be taking values in
a discrete space $\mathbb{Y}_{i}$, $i\in[n]$. For three random vectors
$\boldsymbol{Y}_{\mathbb{A}}$, $\boldsymbol{Y}_{\mathbb{B}}$ and
$\boldsymbol{Y}_{\mathbb{C}}$, where
$\mathbb{A},\mathbb{B},\mathbb{C}\subset[n]$, we say that
$\boldsymbol{Y}_{\mathbb{A}}$ is conditionally independent of
$\boldsymbol{Y}_{\mathbb{B}}$ given $\boldsymbol{Y}_{\mathbb{C}}$ and write
$\boldsymbol{Y}_{\mathbb{A}}\perp\\!\\!\\!\perp\boldsymbol{Y}_{\mathbb{B}}|\boldsymbol{Y}_{\mathbb{C}}$
if

$\textnormal{Pr}(\boldsymbol{Y}_{\mathbb{A}}=\boldsymbol{y}_{\mathbb{A}}|\boldsymbol{Y}_{\mathbb{B}}=\boldsymbol{y}_{\mathbb{B}},\boldsymbol{Y}_{\mathbb{C}}=\boldsymbol{y}_{\mathbb{C}})=\textnormal{Pr}(\boldsymbol{Y}_{\mathbb{A}}=\boldsymbol{y}_{\mathbb{A}}|\boldsymbol{Y}_{\mathbb{C}}=\boldsymbol{y}_{\mathbb{C}}),$

for all $\boldsymbol{y}_{\mathbb{A}}\in\mathbb{Y}_{\mathbb{A}}$,
$\boldsymbol{y}_{\mathbb{B}}\in\mathbb{Y}_{\mathbb{B}}$ and
$\boldsymbol{y}_{\mathbb{C}}\in\mathbb{Y}_{\mathbb{C}}$. In the following, as
a shorthand, we write
$\textnormal{Pr}(\boldsymbol{Y}_{\mathbb{A}}=\boldsymbol{y}_{\mathbb{A}}|\boldsymbol{Y}_{\mathbb{B}}=\boldsymbol{y}_{\mathbb{B}})$
as $\textnormal{Pr}(\boldsymbol{y}_{\mathbb{A}}|\boldsymbol{y}_{\mathbb{B}})$,
for any sets $\mathbb{A},\mathbb{B}\subseteq[n]$.

A BN gives a visual representation of conditional independence by means of a
directed acyclic graph (DAG). For a review of DAGs refer to Appendix 1. The
DAG associated to the BN has vertex set $[n]$, i.e. a vertex is associated to
each variable of $\boldsymbol{Y}$ and edges represent dependence as formalized
next.

###### Definition 2.1

A BN $\mathcal{B}$ for a discrete random vector $\boldsymbol{Y}$ consists of:

  * •

a DAG $G$ with vertex set $[n]$ and edge set
$\\{(j,i):i\in[n],j\in\mathbb{P}_{i}\\}$, encoding $m-1$ conditional
independences $Y_{i}\perp\\!\\!\\!\perp
Y_{\mathbb{ND}_{i}}\;|\,Y_{\mathbb{P}_{i}}$, where $\mathbb{P}_{i}$ and
$\mathbb{ND}_{i}$ are the parent and non-descendant sets of $i$ in $G$;

  * •

conditional probabilities
$\textnormal{Pr}(y_{i}|\boldsymbol{y}_{\mathbb{P}_{i}})$, for
$y_{i}\in\mathbb{Y}_{i}$ and
$\boldsymbol{y}_{\mathbb{P}_{i}}\in\mathbb{Y}_{\mathbb{P}_{i}}$.

We let $\mathbb{V}(\mathcal{B})$ and $\mathbb{E}(\mathcal{B})$ denote the
vertex and edge sets, respectively, of the DAG of the BN $\mathcal{B}$.

Figure 2.1 gives an illustration of a DAG over five discrete variables
$(Y_{1},\dots,Y_{5})$. This DAG embeds the conditional independences
$Y_{2}\perp\\!\\!\\!\perp Y_{1}$, $Y_{3}\perp\\!\\!\\!\perp Y_{2}|Y_{1}$,
$Y_{4}\perp\\!\\!\\!\perp Y_{3}|Y_{1},Y_{2}$ and $Y_{5}\perp\\!\\!\\!\perp
Y_{1},Y_{2}|Y_{3},Y_{4}$. The definition of the BN is completed by a numerical
specification of the probabilities of each variable conditional on the
possible values of the parents.

${1}$${2}$${3}$${4}$${5}$

A simple DAG of a BN with vertex set $\\{1,2,3,4,5\\}$ and edge set
$\\{(1,3),(1,4),(2,4),(3,5),(4,5)\\}$.

Importantly, a BN specifies a factorization of the joint probability of
$\boldsymbol{Y}$ in terms of simpler conditional distributions of individual
variables $Y_{i}$ conditional on their parents
$\boldsymbol{Y}_{\mathbb{P}_{i}}$, $i\in[n]$, i.e.

$\textnormal{Pr}(\boldsymbol{y})=\prod_{i=2}^{n}\textnormal{Pr}(y_{i}|\boldsymbol{y}_{\mathbb{P}_{i}})\textnormal{Pr}(y_{1}),$
(1)

for $\boldsymbol{y}=(y_{1},\dots,y_{n})\in\mathbb{Y}$. For the BN in Figure
2.1, it holds that

$\textnormal{Pr}(\boldsymbol{y})=\textnormal{Pr}(y_{5}|y_{4},y_{3})\textnormal{Pr}(y_{4}|y_{2},y_{1})\textnormal{Pr}(y_{3}|y_{1})\textnormal{Pr}(y_{2})\textnormal{Pr}(y_{1})$

###  2.2 Using Bayesian Networks

The DAG associated to a BN provides an intuitive overview of the relationships
existing between variables of interest. However, it does also provide a
framework to assess if any generic conditional independence holds for a
specific subset of the variables via the so-called _d-separation_ criterion
(see e.g. Darwiche 2009). Furthermore, the DAG provides a framework for fast
exact propagation of probabilities and evidence for the computation of any
(conditional) probability involving a specific subset of the variables.

Both of the above tasks usually require a manipulation of the original DAG,
called _moralization_. In a nutshell, moralization adds an undirected edge
between any two parents of a vertex which are not joined by an edge (hence the
word moralize). This is illustrated for the DAG in Figure 2.1. The vertices 3
and 4 are parents of 5, but there is no edge between them. The same argument
holds for the vertices 1 and 2, which are parents of 4. Thus the moralization
process adds an undirected edge between 1 and 2 and between 3 and 4. The
moralized DAG is reported in Figure 2.2 (left). The _skeleton_ of a DAG
consists of the moralized DAG where all directed edges are replaced by
undirected ones. The skeleton of the DAG in Figure 2.1 is given in Figure 2.2
(right).

The skeleton is an undirected graph, with the extra property of being
decomposable (see Appendix 1), giving the same probability to events
$\boldsymbol{y}\in\mathbb{Y}$ as the original BN. Such graphs enable for fast
probability computations taking advantage of an alternative factorization of
$\textnormal{Pr}(\boldsymbol{Y}=\boldsymbol{y})$ to the one in Equation (1).
Let $\mathbb{C}$ be the set of cliques of the skeleton and $\mathbb{S}$ the
set of separators of the skeleton. Then

$\textnormal{Pr}(\boldsymbol{y})=\frac{\prod_{C\in\mathbb{C}}\textnormal{Pr}(\boldsymbol{y}_{C})}{\prod_{S\in\mathbb{S}}\textnormal{Pr}(\boldsymbol{y}_{S})}$

For example, the skeleton in Figure 2.2 has cliques $\\{1,2,4\\}$,
$\\{1,3,4\\}$ and $\\{3,4,5\\}$, and separators $\\{1,4\\}$ and $\\{3,4\\}$.
Therefore, its probabilities factorize as

$\textnormal{Pr}(\boldsymbol{y})=\frac{\textnormal{Pr}(y_{1},y_{2},y_{4})\textnormal{Pr}(y_{1},y_{3},y_{4})\textnormal{Pr}(y_{3},y_{4},y_{5})}{\textnormal{Pr}(y_{1},y_{4})\textnormal{Pr}(y_{3},y_{4})}$

where the probabilities on the right-hand side are easily computed from the
conditional probabilities defining the BN: for instance,
$\textnormal{Pr}(y_{1},y_{2},y_{4})=\textnormal{Pr}(y_{4}|y_{1},y_{2})\textnormal{Pr}(y_{2})\textnormal{Pr}(y_{1})$.

${1}$${2}$${3}$${4}$${5}$

${1}$${2}$${3}$${4}$${5}$

Moralization process for the BN in Figure 2.1: moralized DAG (left) and
skeleton (right).

###  2.3 Markov Random Fields

Starting from a BN, the moralization process returns an undirected graph which
provides an efficient framework to carry out a variety of tasks including
d-separation, probability propagation, etc. However, in many applications it
is customary to define a probabilistic graphical model based on undirected
graphs directly (see e.g. Guo et al. 2015, Ma et al. 2019, L. Salemi et al.
2019). Such models are usually called _Markov Random Fields_ (MRFs). Here we
follow Robeva and Seigal (2018) and define them via hypergraphs (see the
Appendix for details).

###### Definition 2.2

An MRF $\mathcal{M}$ with respect to an hypergraph with vertex set $[n]$ and
hyperedges $\mathbb{C}_{\mathcal{M}}$ is defined by

$\displaystyle\phi_{\mathcal{M}}:\mathbb{Y}\rightarrow$
$\displaystyle\mathbb{R}_{+}$
$\displaystyle\phi_{\mathcal{M}}(\boldsymbol{y})\mapsto$
$\displaystyle\frac{1}{z}\prod_{C\in\mathcal{P}([n])}\phi_{\mathcal{M}}(\boldsymbol{y}_{C}),$

where $\phi_{\mathcal{M}}(\boldsymbol{y}_{C})=1$ if
$C\not\in\mathbb{C}_{\mathcal{M}}$ and otherwise is a generic function from
$\mathbb{Y}_{\mathbb{C}}$ to $\mathbb{R}_{+}$. The functions
$\phi_{\mathcal{M}}(\boldsymbol{y}_{\mathbb{C}})$ are called _potentials_ and
$z$ is a normalizing constant to impose that
$\sum_{\boldsymbol{y}\in\mathbb{Y}}\phi_{\mathcal{M}}(\boldsymbol{y})=1$. We
write $\mathbb{V}(\mathcal{M})=[n]$ and
$\mathbb{E}(\mathcal{M})=\mathbb{C}_{\mathcal{M}}$.

The previous definition is very general but is needed when working with TNs.
More often, MRFs are defined via a simple undirected graph, and the function
$\phi$ factorizes over the cliques of the graph (for this reason we denoted
the hyperedges as $\mathbb{C}_{\mathcal{M}}$). Conversely, Definition 2.2
allows the factorization to work over any subsets of the vertex set. Notice
that $\phi_{\mathcal{M}}(\boldsymbol{y})$ is actually a probability
distribution since it is non-negative and sums to one.

As an illustration, consider the undirected graph in Figure 2.2, whose cliques
are $\\{1,2,4\\}$, $\\{1,3,4\\}$, $\\{3,4,5\\}$. An MRF over this undirected
graph can be defined by constructing an hypergraph with the same vertex set
and as hyperedges its cliques. Its probability distribution can be written as

$\textnormal{Pr}(\boldsymbol{y})=\phi_{\mathcal{M}}(\boldsymbol{y})=\frac{1}{z}\phi_{\mathcal{M}}(y_{1},y_{2},y_{4})\phi_{\mathcal{M}}(y_{1},y_{3},y_{4})\phi_{\mathcal{M}}(y_{3},y_{4},y_{5}).$

Every BN can be cast as an MRF, as formalized next.

###### Proposition 2.3

A BN $\mathcal{B}$ with $\mathbb{V}(\mathcal{B})=[n]$ is equivalent to an MRF
$\mathcal{M}$ such that $\mathbb{V}(\mathcal{M})=[n]$ and
$\mathbb{E}(\mathcal{M})$ is the set of cliques of the skeleton of the DAG of
$\mathcal{B}$.

Our new global sensitivity analysis methods take advantage of the
transformation of a BN to its skeleton to perform efficient probability
computations. For this reason, although our methods will be illustrated for
BNs only, they also apply for MRFs as well as chain graphs (which again are
usually transformed to an undirected graph).

In particular, an operation required for probability propagation and that our
methods will heavily use is that of _marginalization_. Given an MRF
$\mathcal{M}$ over the vector $\boldsymbol{Y}$ and a set
$\mathbb{A}\subset[n]$ we want to compute the probability
$\textnormal{Pr}(\boldsymbol{y}_{\mathbb{A}})$ which can be derived as

$\textnormal{Pr}(\boldsymbol{y}_{\mathbb{A}})=\sum_{\boldsymbol{y}_{\setminus\mathbb{A}}\in\mathbb{Y}_{\setminus\mathbb{A}}}\textnormal{Pr}(\boldsymbol{y}_{\mathbb{A}},\boldsymbol{y}_{\setminus\mathbb{A}})=\frac{1}{z}\sum_{\boldsymbol{y}_{\setminus\mathbb{A}}\in\mathbb{Y}_{\setminus\mathbb{A}}}\prod_{C\in\mathbb{C}_{\mathcal{M}}}\phi_{\mathcal{M}}(\boldsymbol{y}_{C}).$

Since the normalizing constant $z$ is often hard to compute, we will simply
work with potentials and un-normalized probabilities and call:

$\mathcal{M}_{\setminus\mathbb{A}}(\boldsymbol{y}_{\mathbb{A}})=\sum_{\boldsymbol{y}_{\setminus\mathbb{A}}\in\mathbb{Y}_{\setminus\mathbb{A}}}\prod_{C\in\mathbb{C}_{\mathcal{M}}}\phi_{\mathcal{M}}(\boldsymbol{y}_{C}),$

so that
$\textnormal{Pr}(\boldsymbol{y}_{\mathbb{A}})=\frac{1}{z}\mathcal{M}_{\setminus\mathbb{A}}(\boldsymbol{y}_{\mathbb{A}})$.

###  2.4 Tensor Networks

Since MRFs encode probabilities in a non-normalized way, they are still
subject to the non-negativity constrain. _Tensor networks_ drop this condition
and thus become a graphical representation for general functions Robeva and
Seigal (2018). Therefore, a TN $\mathcal{T}$ is defined as an MRF
$\mathcal{M}$, but whose potentials $\phi_{\mathcal{T}}$ are generic functions
to $\mathbb{R}$. For instance, the problem of tensor decomposition seeks to
find a low-parametric network that encodes, either exactly or approximately, a
given real or complex-valued input tensor.

In this paper we do not aim to learn a TN out of data. Instead, we turn the
input BN $\mathcal{B}$, which may have been learned from data, into its
equivalent MRF as in Proposition 2.3, and then derive a TN $\mathcal{T}$ from
$\mathcal{M}$ by simply letting
$\mathbb{V}(\mathcal{T})=\mathbb{V}(\mathcal{M})$,
$\mathbb{E}(\mathcal{T})=\mathbb{C}_{\mathcal{M}}$ and
$\phi_{\mathcal{T}}(\boldsymbol{y}_{C})=\phi_{\mathcal{M}}(\boldsymbol{y}_{C})$
if $C\in\mathbb{C}_{\mathcal{M}}$ and 1 otherwise. Such a TN is said to be
_derived_ from $\mathcal{M}$. Given this TN $\mathcal{T}$ we then manipulate
it in order to extract relevant sensitivity indices for the original BN
$\mathcal{B}$: see Section 4 below.

##  3 Sobol Indices

The method of Sobol (1990) is one of the most powerful frameworks for global
sensitivity analysis of a general function
$f:\Omega\subseteq\mathbb{R}^{k}\to\mathbb{R}$, whose domain is usually
restricted to a $k$-dimensional rectangle $\Omega=\times_{i\in[k]}\Omega_{i}$.
In our case, as formalized in Section 2, $\Omega$ will be a discrete Cartesian
product $\mathbb{Y}_{\mathbb{E}}$ for a set $\mathbb{E}$ of variables in a BN.

There are multiple kinds of Sobol indices, each with different nuances of
interpretation. Next, we review the two most important ones.

###  3.1 Variance Components

The _variance component_ associated to variable $Y_{i}$ is defined as

$S_{i}:=\frac{\mathrm{Var}_{i}\left[\textnormal{E}_{\setminus
i}[f]\right]}{\mathrm{Var}[f]},$ (2)

where $\textnormal{E}_{\setminus i}$ is the expectation with respect to
$\boldsymbol{Y}_{\setminus i}$ and $\mathrm{Var}_{i}$ is the variance with
respect to $Y_{i}$. This index measures the _additive effect_ that is due to
$Y_{i}$ or, more plainly, how the average model output changes as we vary
$y_{i}$.

###  3.2 Total Indices

The _total index_ of the variable $Y_{i}$ swaps the order of expectation and
variance in Equation (2):

$S^{T}_{i}:=\frac{\textnormal{E}_{\setminus
i}\left[\mathrm{Var}_{i}[f]\right]}{\mathrm{Var}[f]}.$ (3)

The $S^{T}_{i}$ measures the overall effect due to $Y_{i}$ or, in plain words,
the average variability that results from changing $y_{i}$ while leaving all
other parameters fixed. This effect generalizes the additive effect. When
$f$’s inputs are all mutually independent, this implies $S_{i}\leq S^{T}_{i}$;
we have $S_{i}=S^{T}_{i}$ if and only if $Y_{i}$ can be separated from $f$ as
an additive component, that is, if the original function can be decomposed as
$f(\boldsymbol{y})=g(y_{i})+h(\boldsymbol{y}_{\setminus i})$.

Among other uses, total indices are a popular tool for model simplification: a
small $S^{T}_{i}$ guarantees that variable $Y_{i}$ can be safely _frozen_ to
any value of its domain and removed from the model.

###  3.3 Algorithmic Challenges

Being able to capture global effects, Sobol indices are highly expressive and
have become a gold-standard in sensitivity analysis. However, since both the
expectation and the variance in Equations (2) and (3) are multidimensional
integrals (or multi-index summations, in the discrete case), they are
notoriously difficult to compute. Estimation for black-box functions is often
undertaken via Monte Carlo sampling routines which, due to the Central Limit
Theorem, have a slow convergence rate of $O(\sqrt{n})$ in general where $n$ is
the number of samples taken. Even when a Quasi-Monte Carlo sampling plan is
used, for example via Sobol or Saltelli sequences (Saltelli et al. 2008),
large sampling budgets (say, with $n$ between $10^{4}$ and $10^{8}$) are often
needed in practice for each and every Sobol index. Alternatively, more
efficient algorithms exist when $f$ is represented as a _surrogate model_ ,
which is the traditional sensitivity analysis approach to approximate
expensive function evaluations. This includes _polynomial chaos expansions_
(Sudret 2008) and _Gaussian processes_ (Marrel et al. 2009), among others.
Unfortunately, such surrogate models are often poorly suited to approximate
Bayesian networks, which are the target of this paper.

##  4 Proposed Method for Global Sensitivity Analysis in BNs

As customary in sensitivity analysis for BNs (see e.g. Gómez-Villegas et al.
2013), we assume a BN $\mathcal{M}$ and a partition of its vertices into an
output node $\mathbb{O}$, evidential nodes $\mathbb{E}$ and chance nodes
$\mathbb{U}$, i.e. $[n]=\mathbb{O}\cup\mathbb{E}\cup\mathbb{U}$. In more
details:

  * •

$Y_{\mathbb{O}}$ is an _output variable_ whose expected value is a quantity of
interest.

  * •

$\boldsymbol{Y}_{\mathbb{E}}$ is a vector of _evidential variables_ whose
influence on $Y_{\mathbb{O}}$ we wish to quantify.

  * •

$\boldsymbol{Y}_{\mathbb{U}}$ is a vector of _chance variables_ which are
unobserved and cannot be fixed.

Formally, given such a network, we wish to study the following function of
interest:

$\displaystyle\begin{split}f:\mathbb{Y}_{\mathbb{E}}&\to\mathbb{R}\\\
\boldsymbol{y}_{\mathbb{E}}&\mapsto\textnormal{E}_{\mathbb{U}}[Y_{\mathbb{O}}|\boldsymbol{Y}_{\mathbb{E}}=\boldsymbol{y}_{\mathbb{E}}],\end{split}$
(4)

where $\textnormal{E}_{\mathbb{U}}$ refers to the expectation with respect to
$\boldsymbol{Y}_{\mathbb{U}}$. Without loss of generality we assume there is a
single output variable $Y_{\mathbb{O}}$ (see the remark below).

Specifically, our goal is to assess the _global effect_ of observed evidence
$\boldsymbol{Y}_{\mathbb{E}}=\boldsymbol{y}_{\mathbb{E}}$ on the quantity of
interest $f(\boldsymbol{y}_{\mathbb{E}})$ from Equation (4). In Sobol’s
variance-based sense, the effect of an evidential variable $Y_{i}$ is captured
by two indices $S_{i}$ and $S^{T}_{i}$ for every $i\in\mathbb{E}$. As argued
earlier, such effects are key to important tasks including model
simplification, elicitation, interpretation, factor prioritization, risk
analysis, etc. (Saltelli et al. 2008). Ours is the first method to address the
case in which $f$ is defined by a probabilistic graphical model as a function
of evidence in several of its nodes. Note that explicitly evaluating
$f(\boldsymbol{y}_{\mathbb{E}})$ requires inference in a BN. This operation,
which bears significant computational cost for one sample
$\boldsymbol{y}_{\mathbb{E}}$, becomes quickly impractical for thousands or
millions of samples, as is typically required by Monte Carlo based methods.
Building a surrogate model is rarely viable and accurate enough for general
BNs, and the case of correlated input variables is even more challenging. Our
algorithms avoid the explicit evaluation of $f$ and, instead, exploit the
network’s structure to efficiently obtain all variances required in Equations
(2) and (3).

One must also account for the way evidential variables are distributed: they
follow a joint probability mass function
$\textnormal{Pr}(\boldsymbol{y}_{\mathbb{E}})=\mathcal{M}_{\setminus\mathbb{E}}(\boldsymbol{y}_{\mathbb{E}})$
for $\boldsymbol{y}_{\mathbb{E}}\in\mathbb{Y}_{\mathbb{E}}$. If all vertices
$\mathbb{E}$ are root nodes (i.e. they have no parents), Pr is separable,
which is an important particular case in the sensitivity analysis literature;
in fact, many applications of the method of Sobol assume uncorrelated inputs.
For example, when input values can be set at will (as in an experimental
design), it is customary to assume independent, uniformly distributed inputs
in a rectangle $\Omega\subset\mathbb{R}^{|\mathbb{E}|}$. The method we propose
supports correlated inputs as well by means of deriving the network
$\mathcal{M}_{\setminus\mathbb{E}}$, which we do via marginalization. While
marginalizing the offspring of $\mathbb{E}$ is straightforward (no operations
are required), marginalizing its ancestors is more computationally intensive.
Therefore, our method is most useful when the evidential nodes are located in
a relatively upstream position in the graph.

#### Remark.

Note that our formulation is flexible enough to allow for arbitrary output
functions to be considered, including those that represent the _utility_ of a
set of nodes. For example, suppose we initially have a network with nodes
$\mathbb{E}\cup\mathbb{U}$ and a utility function
$g:\mathbb{Y}_{\mathbb{D}}\to\mathbb{Y}_{\mathbb{O}}$ where
$\mathbb{D}\subseteq\mathbb{E}\cup\mathbb{U}$ is a subset of the network’s
variables, and assume we want to study $g$’s output as a function of variables
$\mathbb{E}$. To reduce this task to Equation (4), we encode $g$ into a new
node $\mathbb{O}$ that we add to the graph as in Luque and Díez (2010). The
domain of $Y_{\mathbb{O}}$ is $\mathbb{Y}_{\mathbb{O}}$, its parents are
$\boldsymbol{Y}_{\mathbb{D}}$, and its conditional probability table Pr is
built so as to reflect $g$:

$\textnormal{Pr}(y_{\mathbb{O}}|\boldsymbol{y}_{\mathbb{D}})=\begin{cases}1\mbox{
if }y_{\mathbb{O}}=g(\boldsymbol{y}_{\mathbb{D}})\\\ 0\mbox{
otherwise}.\end{cases}$

This function encoding strategy is demonstrated in one of our experiments in
Section 5.

###  4.1 Network Arithmetics

The proposed algorithm relies on a few building blocks for TN manipulation,
which are introduced next.

###### Definition 4.1

Given a TN $\mathcal{T}$ with vertex set $[n]$ and $\mathbb{A}\subseteq[n]$,
the square of $\mathcal{T}$ with respect to $\mathbb{A}$ is a TN
$\mathcal{T}^{2}$ such that
$(\mathcal{T}^{2})_{\setminus\mathbb{A}}(\boldsymbol{y}_{\mathbb{A}})=(\mathcal{T}_{\setminus\mathbb{A}}(\boldsymbol{y}_{\mathbb{A}}))^{2}$
for any $\boldsymbol{y}_{\mathbb{A}}\in\mathbb{Y}_{\mathbb{A}}$. Since the
squaring commutes with the marginalization, we may simply write
$\mathcal{T}^{2}_{\setminus\mathbb{A}}$.

###### Proposition 4.2

Let $\mathbb{B}=[n]\setminus\mathbb{A}$ and $\mathcal{T}$ be a TN. The square
$\mathcal{T}^{2}$ w.r.t. $\mathbb{A}$ can be constructed as follows:

  * •

$\mathbb{V}(\mathcal{T}^{2})=\mathbb{A}\cup\mathbb{B}\cup\widehat{\mathbb{B}}$
where $\mathbb{A}\cap\widehat{\mathbb{B}}=\emptyset$,
$\mathbb{B}\cap\widehat{\mathbb{B}}=\emptyset$ and there exists a bijection
$\eta:\widehat{\mathbb{B}}\rightarrow\mathbb{B}$;

  * •

$\mathbb{E}(\mathcal{T}^{2})=\mathbb{E}(\mathcal{T})\cup\mathbb{E}^{\prime}$
where $\mathbb{E}^{\prime}$ contains the hyperedge
$\\{\mathbb{A}\cap\mathbb{D}\\}\cup\eta(\mathbb{D}\setminus\mathbb{A})$ for
every $\mathbb{D}\in\mathbb{C}$;

  * •

if $\mathbb{D}\subseteq\mathbb{A}\cup\mathbb{B}$ then
$\phi_{\mathcal{T}^{2}}(\boldsymbol{y}_{\mathbb{D}})=\phi_{\mathcal{T}}(\boldsymbol{y}_{\mathbb{D}})$,
otherwise $\mathbb{D}\subseteq\mathbb{A}\cup\widehat{\mathbb{B}}$ and
$\phi_{\mathcal{T}^{2}}(\boldsymbol{y}_{\mathbb{D}})=\phi_{\mathcal{T}}(\boldsymbol{y}_{C\cap\eta(\widehat{\mathbb{B}})},\boldsymbol{y}_{\mathbb{D}\cap\mathbb{A}})$.

###### Proof 4.3

Proof. For any $\boldsymbol{y}_{\mathbb{A}}\in\mathbb{Y}_{\mathbb{A}}$ we have

$\displaystyle(\mathcal{T}_{\setminus\mathbb{A}}(\boldsymbol{y}_{\mathbb{A}}))^{2}$
$\displaystyle=\left(\sum_{\boldsymbol{y}_{\mathbb{B}}\in\mathbb{Y}_{\mathbb{B}}}\prod_{\mathbb{D}\in\mathcal{P}(\mathbb{V}(\mathcal{T}))}\phi_{\mathcal{T}}(\boldsymbol{y}_{\mathbb{A}\cap\mathbb{D}},\boldsymbol{y}_{\mathbb{B}\cap\mathbb{D}})\right)^{2}$
$\displaystyle=\sum_{\boldsymbol{y}_{\mathbb{B}},\boldsymbol{y}_{\mathbb{B}}^{\prime}\in\mathbb{Y}_{\mathbb{B}}}\prod_{\mathbb{D}\in\mathcal{P}(\mathbb{V}(\mathcal{T}))}\phi_{\mathcal{T}}(\boldsymbol{y}_{\mathbb{A}\cap\mathbb{D}},\boldsymbol{y}_{\mathbb{B}\cap\mathbb{D}})\phi_{\mathcal{T}}(\boldsymbol{y}_{\mathbb{A}\cap\mathbb{D}},\boldsymbol{y}_{\mathbb{B}\cap\mathbb{D}}^{\prime})$
$\displaystyle=\sum_{\boldsymbol{y}_{\mathbb{B}},\boldsymbol{y}_{\mathbb{B}}^{\prime}\in\mathbb{Y}_{\mathbb{B}}}\prod_{\mathbb{D}\in\mathcal{P}(\mathbb{V}(\mathcal{T}))}\phi_{\mathcal{T}^{2}}(\boldsymbol{y}_{\mathbb{A}\cap\mathbb{D}},\boldsymbol{y}_{\mathbb{B}\cap\mathbb{D}})\phi_{\mathcal{T}^{2}}(\boldsymbol{y}_{\mathbb{A}\cap\mathbb{D}},\boldsymbol{y}_{\mathbb{B}\cap\eta(\widehat{\mathbb{D}})}^{\prime})$
$\displaystyle=\sum_{\boldsymbol{y}_{\mathbb{B}\cup\widehat{\mathbb{B}}}\in\mathbb{Y}_{\mathbb{B}\cup\widehat{\mathbb{B}}}}\prod_{\mathbb{D}\in\mathcal{P}(\mathbb{V}(\mathcal{T}^{2}))}\phi_{\mathcal{T}^{2}}(\boldsymbol{y}_{\mathbb{A}\cap\mathbb{D}},\boldsymbol{y}_{(\mathbb{B}\cup\widehat{\mathbb{B}})\cap\mathbb{D}})$
(5)

where the last equality holds since $\mathbb{E}(\mathcal{T}^{2})$ does not
include any element with vertices in both $\mathbb{B}$ and
$\widehat{\mathbb{B}}$. The result easily follows by noticing that the r.h.s.
of Equation (5) is exactly
$\mathcal{T}^{2}_{\setminus\mathbb{A}}(\boldsymbol{y}_{\mathbb{A}})$.

Figure 4.1 illustrates the squaring procedure of Proposition 4.2. A TN
$\mathcal{T}$ has hypergraph with vertex set $[5]$ and hyperedges
$\\{1,2,4\\}$, $\\{1,3,4\\}$ and $\\{3,4,5\\}$. This can be represented by the
undirected graph in Figure 4.1 (left). Suppose that we are interested in the
squared TN $\mathcal{T}^{2}$ w.r.t. $\\{1,2\\}$, which, using Proposition 4.2,
can be constructed as the one in Figure 4.1 (right). The idea is related to a
known algorithm to compute the dot product between two tensor networks
$\mathcal{T}$ and $\mathcal{T}^{\prime}$ by contracting all nodes to yield a
scalar $\left<\mathcal{T},\mathcal{T}^{\prime}\right>$ (Lee and Cichocki
2017). Our version sets $\mathcal{T}=\mathcal{T}^{\prime}$ and skips the node
contraction part.

${1}$${2}$${3}$${4}$${5}$

${1}$${2}$${3}$${4}$${5}$${\widehat{3}}$${\widehat{4}}$${\widehat{5}}$

Squaring process for the tensor network $\mathcal{T}$ (left) with respect to
$\\{1,2\\}$ (right).

###### Definition 4.4

Let $\mathcal{T}$ and $\mathcal{T}^{\prime}$ be two TNs for $\boldsymbol{Y}$
whose hypergraphs have the same vertex set $[n]$. The quotient of
$\mathcal{T}$ over $\mathcal{T}^{\prime}$ denoted as
$\mathcal{T}/\mathcal{T}^{\prime}$ is such that
$(\mathcal{T}/\mathcal{T}^{\prime})(\boldsymbol{y})=\mathcal{T}(\boldsymbol{y})/\mathcal{T^{\prime}}(\boldsymbol{y})$
for any $\boldsymbol{y}\in\mathbb{Y}$.

###### Proposition 4.5

The quotient $\mathcal{T}/\mathcal{T}^{\prime}$ is a TN whose hypergraph has
vertex set $[n]$ and hyperedges
$\mathbb{E}(\mathcal{T})\cup\mathbb{E}(\mathcal{T}^{\prime})$, and potentials
$\phi_{\mathcal{T}/\mathcal{T}^{\prime}}(\boldsymbol{y})=\phi_{\mathcal{T}}(\boldsymbol{y})/\phi_{\mathcal{T}^{\prime}}(\boldsymbol{y})$
for any $\boldsymbol{y}\in\mathbb{Y}$.

###### Proof 4.6

Proof. For any $\boldsymbol{y}\in\mathbb{Y}$ we have

$\frac{\mathcal{T}(\boldsymbol{y})}{\mathcal{T}^{\prime}(\boldsymbol{y})}=\frac{\prod_{\mathbb{A}\in\mathcal{P}(\mathbb{V}(\mathcal{T}))}\phi_{\mathcal{T}}(\boldsymbol{y}_{\mathbb{A}})}{\prod_{\mathbb{A}\in\mathcal{P}(\mathbb{V}(\mathcal{T^{\prime}}))}\phi_{\mathcal{T}^{\prime}}(\boldsymbol{y}_{\mathbb{A}})}=\prod_{\mathbb{A}\in\mathcal{P}(\mathbb{V}(\mathcal{T}))}\frac{\phi_{\mathcal{T}}(\boldsymbol{y}_{\mathbb{A}})}{\phi_{\mathcal{T}^{\prime}}(\boldsymbol{y}_{\mathbb{A}})}=\prod_{{\mathbb{A}\in\mathcal{P}(\mathbb{V}(\mathcal{T}/\mathcal{T}^{\prime}))}}\phi_{\mathcal{T}/\mathcal{T}^{\prime}}(\boldsymbol{y}_{\mathbb{A}})=(\mathcal{T}/\mathcal{T}^{\prime})(\boldsymbol{y})\hfill\Halmos$

###### Proposition 4.7

Let $\mathbb{A}\subset[n]$. Then,
$(\mathcal{T}/\mathcal{T}_{\mathbb{A}}^{\prime})_{\mathbb{A}}=\mathcal{T}_{\mathbb{A}}/\mathcal{T}_{\mathbb{A}}^{\prime}$.
In addition,
$(\mathcal{T}/\mathcal{T}_{\mathbb{A}}^{\prime})^{2}_{\mathbb{A}}=(\mathcal{T}_{\mathbb{A}})^{2}/(\mathcal{T}_{\mathbb{A}}^{\prime})^{2}$,
where the square is w.r.t. $\setminus\mathbb{A}$.

###### Proof 4.8

Proof. For the first statement and any
$\boldsymbol{y}_{\setminus\mathbb{A}}=\mathbb{Y}_{\setminus\mathbb{A}}$,

$\displaystyle(\mathcal{T}/\mathcal{T}_{\mathbb{A}}^{\prime})_{\mathbb{A}}(\boldsymbol{y}_{\setminus\mathbb{A}})$
$\displaystyle=\sum_{\boldsymbol{y}_{\mathbb{A}}\in\mathbb{Y}_{\mathbb{A}}}(\mathcal{T}(\boldsymbol{y}_{\mathbb{A}},\boldsymbol{y}_{\setminus\mathbb{A}})/\mathcal{T}_{\mathbb{A}}^{\prime}(\boldsymbol{y}_{\setminus\mathbb{A}}))$
$\displaystyle=\left(\sum_{\boldsymbol{y}_{\mathbb{A}}\in\mathbb{Y}_{\mathbb{A}}}\mathcal{T}(\boldsymbol{y}_{\mathbb{A}},\boldsymbol{y}_{\setminus\mathbb{A}})\right)/\mathcal{T}_{\mathbb{A}}^{\prime}(\boldsymbol{y}_{\setminus\mathbb{A}})=\mathcal{T}_{\mathbb{A}}(\boldsymbol{y}_{\setminus\mathbb{A}})/\mathcal{T}_{\mathbb{A}}^{\prime}(\boldsymbol{y}_{\setminus\mathbb{A}}).$

The second statement follows immediately from the above and Definition 4.1.

###  4.2 Representing the Target Function

The quotient operation is also useful to build a TN representation of the
function of interest $f$.

###### Definition 4.9

Let $\mathcal{B}$ be our starting point BN, i.e. with a vertex set that
decomposes into evidential, chance and output nodes:
$[n]=\mathbb{E}\cup\mathbb{U}\cup\mathbb{O}$. Let $\mathcal{M}$ be the MRF
derived from $\mathcal{B}$ (Prop. 2.3). The _function TN_ of $\mathcal{M}$ is
the tensor network that results by i. dropping the non-negativity condition on
$\mathcal{M}$; ii. adding the potential $\phi(y_{\mathbb{O}})=y_{\mathbb{O}}$.

###### Proposition 4.10

Let $\mathcal{T}$ be the function TN of the MRF $\mathcal{M}$ (Def. 4.9).
Then, the function of interest $f$ (Equation 4) equals
$(\mathcal{T}/\mathcal{M}_{\mathbb{U}\cup\mathbb{O}})_{\mathbb{U}\cup\mathbb{O}}=\mathcal{T}_{\mathbb{U}\cup\mathbb{O}}/\mathcal{M}_{\mathbb{U}\cup\mathbb{O}}$.

###### Proof 4.11

Proof. For any $\boldsymbol{y}_{\mathbb{E}}\in\mathbb{Y}_{\mathbb{E}}$ we have

$\displaystyle f(\boldsymbol{y}_{\mathbb{E}})$
$\displaystyle=\textnormal{E}_{\mathbb{U}\cup\mathbb{O}}[Y_{\mathbb{O}}|\boldsymbol{Y}_{\mathbb{E}}=\boldsymbol{y}_{\mathbb{E}}]=\sum_{\boldsymbol{y}_{\mathbb{O}\cup\mathbb{U}}\in\mathbb{Y}_{\mathbb{O}\cup\mathbb{U}}}y_{\mathbb{O}}\cdot\textnormal{Pr}(\boldsymbol{y}_{\mathbb{E}},\boldsymbol{y}_{\mathbb{U}\cup\mathbb{O}})/\textnormal{Pr}(\boldsymbol{y}_{\mathbb{E}})$
$\displaystyle=\sum_{\boldsymbol{y}_{\mathbb{O}\cup\mathbb{U}}\in\mathbb{Y}_{\mathbb{O}\cup\mathbb{U}}}y_{\mathbb{O}}\cdot\mathcal{M}(\boldsymbol{y}_{\mathbb{E}},\boldsymbol{y}_{\mathbb{U}\cup\mathbb{O}})/\mathcal{M}_{\mathbb{U}\cup\mathbb{O}}(\boldsymbol{y}_{\mathbb{E}})=\sum_{\boldsymbol{y}_{\mathbb{O}\cup\mathbb{U}}\in\mathbb{Y}_{\mathbb{O}\cup\mathbb{U}}}\mathcal{T}(\boldsymbol{y}_{\mathbb{E}},\boldsymbol{y}_{\mathbb{U}\cup\mathbb{O}})/\mathcal{M}_{\mathbb{U}\cup\mathbb{O}}(\boldsymbol{y}_{\mathbb{E}})=(\mathcal{T}/\mathcal{M}_{\mathbb{U}\cup\mathbb{O}})_{\mathbb{U}\cup\mathbb{O}}(\boldsymbol{y}_{\mathbb{E}})\hfill\Halmos$

###### Proposition 4.12

Let $\mathcal{T}$ be the TN derived from a MRF $\mathcal{M}$. For any
partition $\\{\mathbb{A},\mathbb{B},\mathbb{D}\\}$ of $[n]$ we have
$\textnormal{E}_{\mathbb{A}}[\mathcal{T}_{\mathbb{B}}]=(\mathcal{T}_{\mathbb{B}}\cdot\mathcal{M}_{\mathbb{B}}/\mathcal{M}_{\mathbb{A}\cup\mathbb{B}})_{\mathbb{A}}$.

###### Proof 4.13

Proof. Note that
$\textnormal{E}_{\mathbb{A}}(\mathcal{T}(\boldsymbol{y}_{\setminus\mathbb{A}}))=\sum_{\boldsymbol{y}_{\mathbb{A}\in\mathbb{Y}_{\mathbb{A}}}}\mathcal{T}(\boldsymbol{y}_{\mathbb{A}},\boldsymbol{y}_{\setminus\mathbb{A}})\cdot\textnormal{Pr}(\boldsymbol{y}_{\mathbb{A}}|\boldsymbol{Y}_{\setminus\mathbb{A}}=\boldsymbol{y}_{\setminus\mathbb{A}})$
and
$\textnormal{Pr}(\boldsymbol{y}_{\mathbb{D}})=\mathcal{M}_{\mathbb{A}\cup\mathbb{B}}(\boldsymbol{y}_{\mathbb{D}})$
for any $\boldsymbol{y}_{\mathbb{D}}\in\mathbb{Y}_{\mathbb{D}}$. Therefore

$\displaystyle\textnormal{E}_{\mathbb{A}}(\mathcal{T}_{\mathbb{B}}(\boldsymbol{y}_{\mathbb{D}}))$
$\displaystyle=\sum_{\boldsymbol{y}_{\mathbb{A}}\in\mathbb{Y}_{\mathbb{A}}}\mathcal{T}_{\mathbb{B}}(\boldsymbol{y}_{\mathbb{A}},\boldsymbol{y}_{\mathbb{D}})\cdot\textnormal{Pr}(\boldsymbol{y}_{\mathbb{A}}|\boldsymbol{Y}_{\mathbb{D}}=\boldsymbol{y}_{\mathbb{D}})=\sum_{\boldsymbol{y}_{\mathbb{A}}\in\mathbb{Y}_{\mathbb{A}}}\mathcal{T}_{\mathbb{B}}(\boldsymbol{y}_{\mathbb{A}},\boldsymbol{y}_{\mathbb{D}})\cdot\textnormal{Pr}(\boldsymbol{y}_{\mathbb{A}},\boldsymbol{y}_{\mathbb{D}})/\textnormal{Pr}(\boldsymbol{y}_{\mathbb{D}})$
$\displaystyle=\sum_{\boldsymbol{y}_{\mathbb{A}}\in\mathbb{Y}_{\mathbb{A}}}(\mathcal{T}_{\mathbb{B}}(\boldsymbol{y}_{\mathbb{A}},\boldsymbol{y}_{\mathbb{D}})\cdot\mathcal{M}_{\mathbb{B}}(\boldsymbol{y}_{\mathbb{A}},\boldsymbol{y}_{\mathbb{D}})/\mathcal{M}_{\mathbb{A}\cup\mathbb{B}}(\boldsymbol{y}_{\mathbb{D}}))=(\mathcal{T}_{\mathbb{B}}\cdot\mathcal{M}_{\mathbb{B}}/\mathcal{M}_{\mathbb{A}\cup\mathbb{B}})_{\mathbb{A}}(\boldsymbol{y}_{\mathbb{D}}).\hfill\Halmos$

###### Proposition 4.14

$\textnormal{E}[f]=\mathcal{T}_{[n]}$.

###### Proof 4.15

Proof. Combining Propositions. 4.7 to 4.12 and using the fact that
$\mathcal{M}_{[n]}=1$, we have

$\textnormal{E}[f]=\textnormal{E}_{\mathbb{E}}[(\mathcal{T}/\mathcal{M}_{\setminus\mathbb{E}})_{\setminus\mathbb{E}}]=((\mathcal{T}/\mathcal{M}_{\setminus\mathbb{E}})_{\setminus\mathbb{E}}\cdot\mathcal{M}_{\setminus\mathbb{E}}/\mathcal{M}_{\setminus\mathbb{E}\cup\mathbb{E}})_{\mathbb{E}}=(\mathcal{T}_{\setminus\mathbb{E}}/\mathcal{M}_{\setminus\mathbb{E}}\cdot\mathcal{M}_{\setminus\mathbb{E}})_{\mathbb{E}}=(\mathcal{T}_{\setminus\mathbb{E}})_{\mathbb{E}}=\mathcal{T}_{[n]}.\hfill\Halmos$

Once equipped with the tools described above, we are ready to compute the
Sobol indices as defined in Section 3.

###  4.3 Computing the Global Variance of $f$

Both types of Sobol indices (Equations 2 and 3) require computing
$\mathrm{Var}[f]$ in the denominator. We do so via the formula
$\mathrm{Var}[f]=\textnormal{E}[f^{2}]-\textnormal{E}[f]^{2}$ and deriving
each of these terms as follows:

  * •

First term:

$\displaystyle\textnormal{E}[f^{2}]$
$\displaystyle=\textnormal{E}_{\mathbb{E}}[(\mathcal{T}/\mathcal{M}_{\setminus\mathbb{E}})^{2}_{\setminus\mathbb{E}}]\mbox{
(using Prop.~{}\ref{prop:f})}$
$\displaystyle=((\mathcal{T}/\mathcal{M}_{\setminus\mathbb{E}})^{2}_{\setminus\mathbb{E}}\cdot\mathcal{M}_{\setminus\mathbb{E}}/\mathcal{M}_{\setminus\mathbb{E}\cup\mathbb{E}})_{\mathbb{E}}\mbox{
(using Prop.~{}\ref{prop:expected_value})}$
$\displaystyle=(\mathcal{T}^{2}_{\setminus\mathbb{E}}/\mathcal{M}_{\setminus\mathbb{E}}/\mathcal{M}_{\setminus\mathbb{E}\cup\mathbb{E}})_{\mathbb{E}}\mbox{
(using Prop.~{}\ref{prop:quotient_projection})}$
$\displaystyle=(\mathcal{T}^{2}_{\setminus\mathbb{E}}/\mathcal{M}_{\setminus\mathbb{E}})_{\mathbb{E}}$
(6)
$\displaystyle=(\mathcal{T}^{2}/\mathcal{M}_{\setminus\mathbb{E}})_{\mathbb{E}\cup\setminus\mathbb{E}\cup\widehat{\setminus\mathbb{E}}}\mbox{,
(where the square of }\mathcal{T}\mbox{ is taken w.r.t. }\mathbb{E}\mbox{)}.$
(7)

Note the computational advantage of Equation (7) over Equation (6): the latter
requires marginalization over variables $\setminus\mathbb{E}$, then
$\mathbb{E}$, while the former marginalizes once over all variables, which
gives marginalization heuristics full freedom to find the best order.

  * •

Second term: $\textnormal{E}^{2}[f]=(\mathcal{T}_{[n]})^{2}$ (using Prop.
4.14).

###  4.4 Computing the Variance Components

We obtain index $S_{i}$ via Equation (2), where the numerator is found via the
identity

$\mathrm{Var}_{i}[\textnormal{E}_{\setminus
i}[f]]=\textnormal{E}_{i}[E^{2}_{\setminus
i}[f]]-\textnormal{E}^{2}_{i}[\textnormal{E}_{\setminus i}[f]]$ (8)

and by observing that:

  * •

The first term is equal to

$\displaystyle\textnormal{E}_{i}[\textnormal{E}^{2}_{\setminus i}[f]]$
$\displaystyle=\textnormal{E}_{i}[\textnormal{E}^{2}_{\setminus
i}[(\mathcal{T}/\mathcal{M}_{\setminus\mathbb{E}})_{\setminus\mathbb{E}}]]\mbox{
(using Prop.~{}\ref{prop:f})}$
$\displaystyle=\textnormal{E}_{i}[((\mathcal{T}/\mathcal{M}_{\setminus\mathbb{E}})_{\setminus\mathbb{E}}\cdot\mathcal{M}_{\setminus\mathbb{E}}/\mathcal{M}_{\setminus\mathbb{E}\cup\setminus
i})_{\setminus i}^{2}]\mbox{ (using Prop.~{}\ref{prop:expected_value})}$
$\displaystyle=\textnormal{E}_{i}[(\mathcal{T}_{\setminus\mathbb{E}}/\mathcal{M}_{\setminus\mathbb{E}}\cdot\mathcal{M}_{\setminus\mathbb{E}}/\mathcal{M}_{\setminus\mathbb{E}\cup\setminus
i})_{\setminus i}^{2}]$
$\displaystyle=\textnormal{E}_{i}[(\mathcal{T}_{\setminus\mathbb{E}}/\mathcal{M}_{\setminus
i})^{2}_{\setminus i}]\mbox{ (since }\setminus\mathbb{E}\cup\setminus
i=\setminus i)$
$\displaystyle=((\mathcal{T}_{\setminus\mathbb{E}}/\mathcal{M}_{\setminus
i})^{2}_{\setminus i}\cdot\mathcal{M}_{\setminus i}/\mathcal{M}_{\setminus
i\cup\\{i\\}})_{i}\mbox{ (using Prop.~{}\ref{prop:expected_value})}$
$\displaystyle=(\mathcal{T}^{2}_{\setminus i}/\mathcal{M}_{\setminus i})_{i}.$

Note that, since $\mathcal{T}_{\setminus i}$ is a one-variable network,
$\mathcal{T}^{2}_{\setminus i}$ is most efficiently computed as
$(\mathcal{T}_{\setminus i})^{2}$.

  * •

For the second term, $\textnormal{E}^{2}_{i}[\textnormal{E}_{\setminus
i}[f]]=(\textnormal{E}_{i}[\textnormal{E}_{\setminus
i}[f]])^{2}=\textnormal{E}^{2}[f]=(\mathcal{T}_{[n]})^{2}$ (using Prop. 4.14).

The denominator is simply the model’s global variance $\mathrm{Var}[f]$
(Section 4.3).

###  4.5 Computing the Total Indices

To obtain index $S^{T}_{i}$, we use the following alternative definition:

$S^{T}_{i}=1-\frac{\mathrm{Var}_{\setminus
i}\left[\textnormal{E}_{i}[f]\right]}{\mathrm{Var}[f]},$ (9)

which is equivalent (Saltelli et al. 2008) to the one given in Equation (3).
To obtain the numerator, we use the identity in Equation (8) and observe that:

  * •

For the first term, analogously to Section 4.4, substituting $i$ for
$\setminus i$ and vice versa we arrive at $\textnormal{E}_{\setminus
i}[\textnormal{E}^{2}_{i}[f]]=(\mathcal{T}^{2}_{i}/\mathcal{M}_{\setminus\mathbb{E}\cup\\{i\\}})_{\setminus
i}$, where the square is taken w.r.t. $\mathbb{E}\setminus i$.

  * •

For the second term, $\textnormal{E}^{2}_{\setminus
i}[\textnormal{E}_{i}[f]]=(\textnormal{E}_{\setminus
i}[\textnormal{E}_{i}[f]])^{2}=\textnormal{E}^{2}[f]=(\mathcal{T}_{[n]})^{2}$
(using Prop. 4.14).

###  4.6 Summarized Algorithm

We put together all steps in Algorithm 1 so as to minimize the cost of
intermediate calculations. Note that our algorithm can also obtain Sobol
indices for any vector of variables $\boldsymbol{\alpha}$ (the so-called
higher-order indices), since all mathematical arguments above hold true when
replacing a single variable $i$ by a vector of them $\boldsymbol{\alpha}$.

Input: Bayesian network $\mathcal{B}$ with variables
$\mathbb{O}\cup\mathbb{U}\cup\mathbb{E}$, variable of interest
$i\in\mathbb{E}$, domain $\mathbb{Y}_{\mathbb{O}}$ of the output variable
$Y_{\mathbb{O}}$

Result: Variance component $S_{i}$ and total index $S^{T}_{i}$ for variable
$Y_{i}$

$\blacktriangleright$ Derive the MRF $\mathcal{M}$ that is equivalent to
$\mathcal{B}$ (Prop. 2.3)

$\blacktriangleright$ Construct the function TN $\mathcal{T}$ of $\mathcal{M}$
(Def. 4.9)

$\blacktriangleright$ Find $\mathcal{J}:=\mathcal{M}_{\setminus\mathbb{E}}$
via marginalization

$\blacktriangleright$ Compute
$\mathcal{T}_{i},\mathcal{J}_{i},\mathcal{T}_{\setminus
i},\mathcal{J}_{\setminus i},\mathcal{T}_{[n]}$ via marginalization

$\blacktriangleright$ Compute $\mathcal{T}^{2},\mathcal{T}^{2}_{i}$, and
$\mathcal{T}^{2}_{\setminus i}$ via network squaring (Prop. 4.2)

$\blacktriangleright$ Compute
$\mathcal{T}^{2}/\mathcal{J},\mathcal{T}^{2}_{i}/\mathcal{J}_{i}$, and
$\mathcal{T}^{2}_{\setminus i}/\mathcal{J}_{\setminus i}$ via network quotient
(Prop. 4.5)

$\blacktriangleright$
$V:=(\mathcal{T}^{2}/\mathcal{J})_{\mathbb{E}\cup\setminus\mathbb{E}\cup\widehat{\setminus\mathbb{E}}}-(\mathcal{T}_{[n]})^{2}$

$\blacktriangleright$ $S_{i}:=((\mathcal{T}^{2}_{\setminus
i}/\mathcal{J}_{\setminus i})_{i}-(\mathcal{T}_{[n]})^{2})/V$

$\blacktriangleright$
$S^{T}_{i}:=1-\left((\mathcal{T}^{2}_{i}/\mathcal{J}_{i})_{\setminus
i}-(\mathcal{T}_{[n]})^{2}\right)/V$

return $S_{i},S^{T}_{i}$

Algorithm 1 Compute first-order Sobol indices for a given BN

##  5 Illustrations

In order to assess the validity of our algorithms, we now consider two
networks based on simulated data and one real-world network. In all cases, we
use exact marginalization via the _minimal weight_ heuristic for variable
ordering (Kjærulff 1990): at each step, we marginalize the node such that the
product of its neighbors’ cardinalities is minimal.

###  5.1 Concrete Resistance

As a first example, we consider a 24-node Bayesian network due to Castillo and
Kjærulff (2003) designed to assess the damage of reinforced concrete
structures of buildings. There is one output variable (the damage of a
reinforced concrete beam), 16 observable variables influencing the damage of
reinforced concrete structures and seven intermediate unobservable variables
that define some partial states of the structure. The list of variables and
their definition is reported in Table 5.1 and their cause-effect relationships
are summarized by the BN in Figure 5.1. The original variables in Castillo and
Kjærulff (2003) are continuous and their model is defined as a Gaussian BN
(e.g. Darwiche 2009). From their definition we derived a discrete BN as
follows: (i) we simulate 1000000 observations from their Gaussian BN; (ii)
each variable is independently discretized into three categories (low, medium,
high) using the equal frequency method (Nojavan et al. 2017); (iii) retaining
the same DAG as in Castillo and Kjærulff (2003) the conditional probabilities
for the discrete data are learned. All these steps are carried out with the
bnlearn package (Scutari 2010).

Variable labels and names for the concrete resistance network.  Variable
Definition $O$ Damage assessment $U_{7}$ Cracking state $U_{6}$ Cracking state
in shear domain $U_{5}$ Steel corrosion $U_{4}$ Cracking state in flexure
domain $U_{3}$ Shrinkage cracking $U_{2}$ Worst cracking in flexure domain
$U_{1}$ Corrosion state $E_{16}$ Weakness of the beam $E_{15}$ Deflection of
the beam $E_{14}$ Position of the worst shear crack $E_{13}$ Breadth of the
worst shear crack $E_{12}$ Position of the worst flexure crack $E_{11}$
Breadth of the worst flexure crack $E_{10}$ Length of the worst flexure cracks
$E_{9}$ Cover $E_{8}$ Structure age $E_{7}$ Humidity $E_{6}$ pH value in the
air $E_{5}$ Content of chlorine in the air $E_{4}$ Number of shear cracks
$E_{3}$ Number of flexure cracks $E_{2}$ Shrinkage $E_{1}$ Corrosion

The 16 evidential variables are the roots of the DAG (that is, without
parents), for a total of $3^{16}\approx 4.31\cdot 10^{7}$ grid points in the
domain of the function of interest $f$. The resulting concrete resistances are
given by output variable $O$, whose values we map to $\mathbb{R}$ as follows:
$\mathrm{low}\equiv 0,\mathrm{medium}\equiv 1$, and $\mathrm{high}\equiv 2$.
Table 5.1 and Figure 5.1 summarize the global sensitivity indices. Overall,
computing all 16 variance components and total indices took $1.12$s and
$1.79$s, respectively. Each index requires four marginalizations in our
implementation, for a total of $16\cdot 2\cdot 4=128$ such operations that
took 2.24s globally. Based on the average marginalization time, a brute-force
computation evaluating all entries of $f$ would have taken well above $10^{5}$
seconds, and even a Monte Carlo approximation using $10^{6}$ samples would
take ca. 18’000s.

The results suggest that the most critical variables for the output are
$E_{12}$ (position of the worst flexure crack), $E_{13}$ (breadth of the worst
shear crack) and $E_{14}$ (position of the worst shear crack). Notice that the
significant differences between variance components and total indices indicate
that much of those effects are due to variable interactions, not to single
variables alone, therefore highlighting the need of global sensitivity
analyses to uncover the actual relationships between input and output
variables.

Sobol indices and computation times for the concrete resistance network.  Var.
$i$ $S_{i}$ Time ($S_{i}$) $S^{T}_{i}$ Time ($S^{T}_{i}$) E1 0.05464 0.06325
0.06335 0.09594 E2 0.03081 0.05967 0.03569 0.10103 E3 0.06339 0.06109 0.07750
0.09831 E4 0.04156 0.05868 0.05470 0.09775 E5 0.04837 0.06028 0.06217 0.09880
E6 0.06281 0.06016 0.07860 0.09927 E7 0.01783 0.05837 0.02503 0.09913 E8
0.02344 0.05937 0.03393 0.09868 E9 0.02357 0.05922 0.03412 0.09612 E10 0.05574
0.06039 0.07461 0.09657 E11 0.05560 0.05826 0.07452 0.09817 E12 0.12943
0.05956 0.15069 0.09695 E13 0.14802 0.05955 0.17370 0.09784 E14 0.08480
0.05865 0.10512 0.10081 E15 0.07363 0.05868 0.07535 0.10183 E16 0.01338
0.05963 0.01388 0.09872

Total indices for the concrete resistance network: chance nodes are shown in
gray, output node in orange, while evidential nodes are color-coded according
to their respective total indices.

###  5.2 Synthetic Network with Dependent Inputs

As a second simulated example, we created a 40-node Gaussian Bayesian network
with three output nodes, 15 evidential nodes, and 22 chance nodes. Both the
structure of the DAG (reported in Figure 5.2) and the conditional
distributions of the nodes where chosen a priori so as to best illustrate the
methodology. In this case, however, the evidential variables exhibit some
dependence. A discrete BN is obtained using the same routine as in the
previous example. The output is a new additional node
$O(O_{1},O_{2},O_{3}):=O_{1}+O_{2}+O_{3}$, where $\mathrm{low}\equiv
0,\mathrm{medium}\equiv 1$ and $\mathrm{high}\equiv 2$. Since there are 15
evidential variables, the domain of the function of interest $f$ contains
$3^{15}\approx 1.43\cdot 10^{7}$ grid points. Table 5.2 (right) and Figure 5.2
report the global sensitivity indices. Note that $S_{i}>S^{T}_{i}$ for several
variables; this is a well-known consequence of the non-independence of the
function’s inputs. Overall, computing all 15 variance components and total
indices took $5.88$s and $11.48$s, respectively. Although, the computation
time is slightly larger than in the previous example, it is still feasible to
compute the global sensitivity indices in the much more challenging case of
dependent inputs.

To illustrate the consequences of overlooking dependence in the inputs, we
created the same discrete BN with 15 inputs, but where the dependence between
the evidential variables is deleted. The global sensitivity indices for the
independent network are reported on the left side of Table 5.2. We can notice
two things. First, $S_{i}^{T}\geq S_{i}$ for all variables, due to missed
dependence between inputs. Second, evidential variables that are important for
the output through the dependence with other inputs (for instance $E_{2}$ or
$E_{5}$) have now a much smaller effect on the output and their relevance
would be overlooked by considering the inputs independent.

Sobol indices and computation times for the synthetic network (with both
dependent and independent inputs).  Independent inputs Dependent inputs Var.
$i$ $S_{i}$ Time ($S_{i}$) $S^{T}_{i}$ Time ($S^{T}_{i}$) $S_{i}$ Time
($S_{i}$) $S^{T}_{i}$ Time ($S^{T}_{i}$) E1 0.00812 0.38456 0.00832 0.66675
0.05757 0.36064 0.00437 0.63704 E10 0.00617 0.35793 0.00658 0.67984 0.00277
0.35715 0.00311 0.65014 E11 0.03099 0.35240 0.03118 0.65987 0.01198 0.35403
0.01237 0.64012 E12 0.06793 0.34724 0.06883 0.74234 0.04408 0.35532 0.02202
0.71773 E13 0.00471 0.35308 0.00497 0.71207 0.02735 0.35807 0.00482 0.72696
E14 0.00336 0.35799 0.00337 0.69925 0.00121 0.36493 0.00124 0.65699 E15
0.03169 0.36250 0.03193 0.60682 0.01366 0.35283 0.01395 0.61534 E2 0.04169
0.35959 0.04200 0.64520 0.15364 0.35414 0.01426 0.64304 E3 0.02196 0.35772
0.02223 0.64319 0.10733 0.35941 0.01452 0.63806 E4 0.05410 0.35593 0.05445
0.64833 0.13931 0.35352 0.03504 0.65013 E5 0.08062 0.35327 0.08140 0.58844
0.27517 0.35414 0.03071 0.60475 E6 0.34252 0.35654 0.34322 0.59323 0.37965
0.35600 0.17480 0.61408 E7 0.11156 0.35384 0.11220 0.64536 0.16543 0.41148
0.05473 0.64713 E8 0.08040 0.35663 0.08076 0.59313 0.13329 0.35635 0.02463
0.61735 E9 0.11108 0.35002 0.11224 0.64339 0.17694 0.35846 0.06924 0.64184

Total indices for the synthetic networks with dependent inputs: chance nodes
are shown in gray, output node in orange, while evidential nodes are color-
coded according to their respective total indices.

###  5.3 Risk Assessment in Project Complexity

Our last example is a BN due to Qazi et al. (2016) representing critical risks
specific to a construction project and adapted from an existing model proposed
by Eybpoosh et al. (2011). Qazi et al. (2016) considers eight project
complexity elements which can be chosen at the starting of the project as well
as four project objectives: timeliness; cost; quality; and market share. All
these objectives are henceforth presented as negative counterparts in order to
align these to the notion of risks and are assumed to be equally important.
There are 14 additional variables that are not observed at the start of the
project (chance variables). Table 5.3 reports a full list of the variables,
which are all defined as binary. Although the eight project complexity
elements are considered as decision variables in Qazi et al. (2016), giving an
influence diagram representation of the problem (e.g. Darwiche 2009), here
they are considered as evidential random variables and are assigned a priori
an equal probability to their outcomes.

Variable labels and names for the project complexity network.  Variable
Definition $E_{1}$ Lack of experience with the involved team $E_{2}$ Use of
innovative technology $E_{3}$ Lack of experience with technology $E_{4}$
Strict quality requirements $E_{5}$ Multiple contracts $E_{6}$ Multiple
stakeholder and variety of perspectives $E_{7}$ Political instability $E_{8}$
Susceptibility to natural disasters $U_{1}$ Contractors’ lack of experience
$U_{2}$ Suppliers default $U_{3}$ Delays in design and regulatory approval
$U_{4}$ Contract related problems $U_{5}$ Economic issues in country $U_{6}$
Major design changes $U_{7}$ Delays in obtaining raw material $U_{8}$ Non-
availability of local resources $U_{9}$ Unexpected events $U_{10}$ Increase in
raw material prices $U_{11}$ Change in project specifications $U_{12}$
Conflicts with project stakeholders $U_{13}$ Decrease in productivity $U_{14}$
Delays/interruptions $O_{1}$ Decrease in quality of work $O_{2}$ Low market
share/reputational issues $O_{3}$ Time overrun $O_{4}$ Cost overrrun

Table 5.3 and Figure 5.3 report a summary of the results. Overall, computing
all 8 variance components and total indices took $1.06$s and $2.00$s,
respectively, thus highlighting the computational efficiency of our
algorithms. Clearly, the lack of experience with the involved team ($E_{1}$)
has the biggest impact on the possible difficulties involved with the project,
with Sobol indices more than twice larger than for all other variables. Other
complexity elements that are relevant for the output are (in order of
importance) strict quality requirements ($E_{4}$), lack of experience with
technology ($E_{3}$) and political instability ($E_{7}$).

Sobol indices and computation times for the project complexity network.  Var.
$i$ $S_{i}$ Time ($S_{i}$) $S^{T}_{i}$ Time ($S^{T}_{i}$) E1 0.48218 0.12326
0.48513 0.22200 E2 0.00502 0.12252 0.02151 0.21638 E3 0.12845 0.12056 0.13723
0.21051 E4 0.16822 0.12094 0.18604 0.21015 E5 0.00193 0.11998 0.00466 0.21473
E6 0.00093 0.12192 0.00224 0.22305 E7 0.09408 0.12048 0.10671 0.22147 E8
0.08309 0.12162 0.09678 0.22114

Total indices for the synthetic networks with dependent inputs: chance nodes
and initial output nodes are shown in gray, final output node in orange, while
evidential nodes are color-coded according to their respective total indices.

##  6 Conclusions

We have addressed efficient Sobol sensitivity analysis of Bayesian networks,
more specifically when the output is the expected value of one of its nodes as
a function of a set of other nodes. Our algorithms compute the two most
important Sobol indices (variance components and total indices) by first
transforming the input BN into a TN that captures the target function.
Afterwards, via various manipulations of the TN, we compute each required
variance using only a few inference queries (marginalization) over the
modified TN. We showed this procedure to scale well: it is able to handle
networks where both brute-force and Monte Carlo estimation would be
impossible, as they would require millions of inference queries over the BN.
In addition, our method is exact, provided that exact inference is affordable
on the modified TN.

Our method is a generalization of previous SA algorithms for acyclic tensor
networks such as the Tucker or tensor train decompositions. Our version works
on cyclic networks, which arise naturally (due to moralization) when casting a
general Bayesian network into a TN. The key enabling concepts are those of
network squaring and quotient, thanks to which we can manipulate the function
of interest so as to produce networks whose marginalization yields the desired
statistics.

As a venue of future research, we will explore the feasibility of computing
more sophisticated variance-based sensitivity indices from BN-defined
functions. These potentially include the mean dimension, the Shapley values,
the dimension distribution, and other more recent metrics that are based on
the method of Sobol.

##  7 Elements of Graph Theory

###  7.1 Graphs

A graph $G=(\mathbb{V}(G),\mathbb{E}(G))$ is defined by two sets
$\mathbb{V}(G)$ and $\mathbb{E}(G)$, where $\mathbb{V}(G)$ is a finite set of
vertices, in this paper $[n]$, and $\mathbb{E}(G)$ is a set of edges
consisting of a subset of all possible ordered pairs of distinct vertices. The
edges of a graph can be directed on undirected, depending on whether or not
the order of the involved vertices matters.

Let $G=(\mathbb{V}(G),\mathbb{E}(G))$ be a graph. When $(i,j)\in\mathbb{E}(G)$
and $(j,i)\not\in\mathbb{E}(G)$, we say there is a directed edge from $i$ to
$j$. Conversely, if $(i,j),(j,i)\in\mathbb{E}(G)$, we say there is an
undirected edge between $i$ and $j$. A graph which all edges are directed is
called a directed graph, a graph in which all edges are undirected is called
undirected graph, and if both types of edges are present the graph is called
mixed.

A path from vertex $i$ to vertex $j$ is an ordered sequence of vertices
$(i_{1},\dots,i_{r})$ starting in $i=i_{1}$ and ending in $j=i_{r}$ such that
there is an edge $(i_{k},i_{k+1})$ in $G$, $k=1,\dots,r-1$. The length of the
path is $r-1$. A path is said to be closed if it has the same starting and
ending vertices, that is $i=j$.

###  7.2 Directed Graphs

If $(i,j)\in\mathbb{V}(G)$, $i$ is said to be a parent of $j$ and $j$ is said
to be a child of $i$. The set of all parents of a node $i$ is denoted by
$\mathbb{P}_{i}$. The descendant set of a node $i$ is made of all vertices $j$
such that there is a path from $i$ to $j$, and denoted by $\mathbb{D}_{i}$.
The set of non-descendants of a vertex $i$ is
$\mathbb{ND}_{i}=\mathbb{V}(G)\setminus\mathbb{D}_{i}$.

A cycle is a closed path in a directed path. A directed acyclic graph is a
directed graph with no cycles. The mixed graph obtained by first joining every
pair of vertices with a common child in a directed graph is called a moral
graph. The undirected graph which substitutes every directed edge in a moral
graph with an undirected edge is called a skeleton.

###  7.3 Undirected Graphs

A set of vertices $\mathbb{A}\subseteq\mathbb{V}(G)$ is said to be complete if
there are edges between every pair of nodes in $\mathbb{A}$. A complete set of
nodes $\mathbb{C}$ is called a clique if it is maximal, that is, it is not a
proper subset of another complete set. A loop is a closed path in a directed
graph. A chord is and edge between two vertices in a loop that is not
contained in that loop. An undirected graph is said to be triangulated if
every loop of length four or more has at least one chord. The skeleton of a
directed acyclic graph is triangulated.

Let $(\mathbb{C}_{1},\dots,\mathbb{C}_{k})$ be an ordered sequence of all
cliques of an undirected graph and $S_{i}=C_{i}\cap\\{C_{1}\cup\cdots\cup
C_{i-1}\\}$ be the separator. The sequence
$(\mathbb{C}_{1},\dots,\mathbb{C}_{k})$ is said to respect the running
intersection property if $S_{i}$ is a subset of at least one element in
$\\{C_{1},\dots,C_{i-1}\\}$, for all $i=1,\dots,k$. An undirected graph has at
least one ordered sequence of cliques respecting the running intersection
property if and only if it is triangulated. Triangulated graphs, which respect
the running intersection property, have great computational advantages, and
this is the reason why a directed acyclic graph is usually transformed into
its skeleton. We will not give details here about this (which can be found in
Darwiche 2009).

###  7.4 Hypergraphs

A hypergraph is a generalization of a graph in which an edge can join any
number of vertices. Conversely, in ordinary graphs an edge can only join two
vertices. Formally, an hypergraph $H=(\mathbb{V}(H),\mathbb{E}(H))$ where
$\mathbb{V}(H)$ is a set of vertices and $\mathbb{E}(H)$ is a set of non-empty
subsets of $\mathbb{V}(H)$ called hyperedges. Therefore,
$\mathbb{E}(H)\subset\mathcal{P}(\mathbb{V}(H))\setminus\emptyset$.

## References

  * Ballester-Ripoll et al. (2018) Ballester-Ripoll R, Paredes EG, Pajarola R (2018) Tensor algorithms for advanced sensitivity metrics. _SIAM/ASA Journal on Uncertainty Quantification_ 6(3):1172–1197.
  * Ballester-Ripoll et al. (2019) Ballester-Ripoll R, Paredes EG, Pajarola R (2019) Sobol tensor trains for global sensitivity analysis. _Reliability Engineering & System Safety_ 183:311–322.
  * Bielza and Larrañaga (2014) Bielza C, Larrañaga P (2014) Bayesian networks in neuroscience: a survey. _Frontiers in Computational Neuroscience_ 8:131.
  * Bolt and Renooij (2014) Bolt JH, Renooij S (2014) Local sensitivity of Bayesian networks to multiple simultaneous parameter shifts. _European Workshop on Probabilistic Graphical Models_ , 65–80 (Springer).
  * Borgonovo (2017) Borgonovo E (2017) _Sensitivity analysis: an introduction for the management scientist_ (Springer).
  * Borgonovo and Plischke (2016) Borgonovo E, Plischke E (2016) Sensitivity analysis: a review of recent advances. _European Journal of Operational Research_ 248(3):869–887.
  * Cai et al. (2018) Cai B, Kong X, Liu Y, Lin J, Yuan X, Xu H, Ji R (2018) Application of Bayesian networks in reliability evaluation. _IEEE Transactions on Industrial Informatics_ 15(4):2146–2157.
  * Castillo et al. (1997) Castillo E, Gutiérrez JM, Hadi AS (1997) Sensitivity analysis in discrete Bayesian networks. _IEEE Transactions on Systems, Man, and Cybernetics-Part A: Systems and Humans_ 27(4):412–423.
  * Castillo and Kjærulff (2003) Castillo E, Kjærulff U (2003) Sensitivity analysis in Gaussian Bayesian networks using a symbolic-numerical technique. _Reliability Engineering & System Safety_ 79(2):139–148.
  * Chan and Darwiche (2004) Chan H, Darwiche A (2004) Sensitivity analysis in Bayesian networks: from single to multiple parameters. _Proocedings of the 16th Conference on Uncertainty in Artificial Intelligence_ , 317–325.
  * Chan and Darwiche (2005) Chan H, Darwiche A (2005) A distance measure for bounding probabilistic belief change. _International Journal of Approximate Reasoning_ 38:149–174.
  * Chastaing et al. (2012) Chastaing G, Gamboa F, Prieur C (2012) Generalized Hoeffding-Sobol decomposition for dependent variables-application to sensitivity analysis. _Electronic Journal of Statistics_ 6:2420–2448.
  * Coupé and Van Der Gaag (2002) Coupé VMH, Van Der Gaag LC (2002) Properties of sensitivity analysis of Bayesian belief networks. _Annals of Mathematics and Artificial Intelligence_ 36(4):323–356.
  * Darwiche (2009) Darwiche A (2009) _Modeling and reasoning with Bayesian networks_ (Cambridge University Press).
  * Darwiche (2010) Darwiche A (2010) _samIam: Sensitivity analysis, modeling, inference and more_. URL http://reasoning.cs.ucla.edu/samiam/.
  * Do and Razavi (2020) Do NC, Razavi S (2020) Correlation effects? A major but often neglected component in sensitivity and uncertainty analysis. _Water Resources Research_ 56(3):e2019WR025436.
  * Douglas-Smith et al. (2020) Douglas-Smith D, Iwanaga T, Croke BF, Jakeman AJ (2020) Certain trends in uncertainty and sensitivity analysis: an overview of software tools and techniques. _Environmental Modelling & Software_ 124:104588.
  * Drury et al. (2017) Drury B, Valverde-Rebaza J, Moura MF, de Andrade Lopes A (2017) A survey of the applications of Bayesian networks in agriculture. _Engineering Applications of Artificial Intelligence_ 65:29–42.
  * Eybpoosh et al. (2011) Eybpoosh M, Dikmen I, Talat Birgonul M (2011) Identification of risk paths in international construction projects using structural equation modeling. _Journal of Construction Engineering and Management_ 137(12):1164–1175.
  * Gómez-Villegas et al. (2013) Gómez-Villegas MA, Main P, Susi R (2013) he effect of block parameter perturbations in Gaussian Bayesian networks: Sensitivity and robustness. _Information Sciences_ 222:439–458.
  * Guo et al. (2015) Guo J, Cheng J, Levina E, Michailidis G, Zhu J (2015) Estimating heterogeneous graphical models for discrete data with an application to roll call voting. _The Annals of Applied Statistics_ 9(2):821.
  * Hänninen and Kujala (2012) Hänninen M, Kujala P (2012) Influences of variables on ship collision probability in a Bayesian belief network model. _Reliability Engineering & System Safety_ 102:27–40.
  * Iooss and Prieur (2019) Iooss B, Prieur C (2019) Shapley effects for sensitivity analysis with correlated inputs: comparisons with Sobol’ indices, numerical estimation and applications. _International Journal for Uncertainty Quantification_ 9(5).
  * Kjærulff (1990) Kjærulff U (1990) Triangulation of graphs - algorithms giving small total state space. Technical Report R 90-09, Department of Mathematics and Computer Science, Strandvejen, DK 9000 Aalborg, Denmark.
  * Kjærulff and van der Gaag (2000) Kjærulff U, van der Gaag LC (2000) Making sensitivity analysis computationally efficient. _Proceedings of the Sixteenth Conference on Uncertainty in Artificial Intelligence_ , 317–325.
  * Kleemann et al. (2017) Kleemann J, Celio E, Fürst C (2017) Validation approaches of an expert-based Bayesian belief network in northern Ghana, West Africa. _Ecological modelling_ 365:10–29.
  * Koller et al. (2009) Koller D, Friedman N, Bach F (2009) _Probabilistic graphical models: principles and techniques_ (MIT press).
  * L. Salemi et al. (2019) L Salemi P, Song E, Nelson BL, Staum J (2019) Gaussian Markov random fields for discrete optimization via simulation: framework and algorithms. _Operations Research_ 67(1):250–266.
  * Lee and Cichocki (2017) Lee N, Cichocki A (2017) Fundamental tensor operations for large-scale data analysis using tensor network formats. _Multidimensional Systems and Signal Processing_ 1–40.
  * Leonelli et al. (2017) Leonelli M, Görgen C, Smith JQ (2017) Sensitivity analysis in multilinear probabilistic models. _Information Sciences_ 411:84–97.
  * Leonelli et al. (2021) Leonelli M, Ramanathan R, Wilkerson RL (2021) Sensitivity and robustness analysis in Bayesian networks with the bnmonitor R package. _arXiv:2107.11785_ .
  * Leonelli and Riccomagno (2019) Leonelli M, Riccomagno E (2019) A geometric characterisation of sensitivity analysis in monomial models. _arXiv:1901.02058_ .
  * Li and Mahadevan (2017) Li C, Mahadevan S (2017) Sensitivity analysis of a Bayesian network. _ASCE-ASME Journal of Risk and Uncertainty in Engineering Systems, Part J: Mechanical Engineering_ 4(1):011003.
  * Luque and Díez (2010) Luque M, Díez FJ (2010) Variable elimination for influence diagrams with super value nodes. _International Journal of Approximate Reasoning_ 51(6):615–631.
  * Ma et al. (2019) Ma YP, Sudakov I, Strong C, Golden KM (2019) Ising model for melt ponds on arctic sea ice. _New Journal of Physics_ 21(6):063029.
  * Makaba et al. (2020) Makaba T, Doorsamy W, Paul BS (2020) Bayesian network-based framework for cost-implication assessment of road traffic collisions. _International Journal of Intelligent Transportation Systems Research_ 19:240–253.
  * Marrel et al. (2009) Marrel A, Iooss B, Laurent B, Roustant O (2009) Calculations of Sobol indices for the Gaussian process metamodel. _Reliability Engineering & System Safety_ 94:742–751.
  * McLachlan et al. (2020) McLachlan S, Dube K, Hitman GA, Fenton N, Kyrimi E (2020) Bayesian networks in healthcare: distribution by medical condition. _Artificial Intelligence in Medicine_ 101912\.
  * Nojavan et al. (2017) Nojavan F, Qian SS, Stow CA (2017) Comparative analysis of discretization methods in Bayesian networks. _Environmental Modelling & Software_ 87:64–71.
  * Qazi et al. (2016) Qazi A, Quigley J, Dickson A, Kirytopoulos K (2016) Project complexity and risk management (ProCRiM): Towards modelling project complexity driven risk paths in construction projects. _International Journal of Project Management_ 34(7):1183–1198.
  * Rai (2014) Rai P (2014) _Sparse Low Rank Approximation of Multivariate Functions – Applications in Uncertainty Quantification_. Doctoral thesis, École Centrale Nantes, URL https://tel.archives-ouvertes.fr/tel-01143694.
  * Razavi et al. (2021) Razavi S, Jakeman A, Saltelli A, Prieur C, Iooss B, Borgonovo E, Plischke E, Lo Piano S, Iwanaga T, Becker W, et al. (2021) The future of sensitivity analysis: an essential discipline for systems modeling and policy support. _Environmental Modelling & Software_ 137:104954.
  * Robeva and Seigal (2018) Robeva E, Seigal A (2018) Duality of graphical models and tensor networks. _Information and Inference: A Journal of the IMA_ 8(2):273–288.
  * Rohmer (2020) Rohmer J (2020) Uncertainties in conditional probability tables of discrete Bayesian belief networks: a comprehensive review. _Engineering Applications of Artificial Intelligence_ 88:103384.
  * Saltelli et al. (2019) Saltelli A, Aleksankina K, Becker W, Fennell P, Ferretti F, Holst N, Li S, Wu Q (2019) Why so many published sensitivity analyses are false: a systematic review of sensitivity analysis practices. _Environmental Modelling & Software_ 114:29–39.
  * Saltelli et al. (2008) Saltelli A, Ratto M, Andres T, Campolongo F, Cariboni J, Gatelli D, Saisana M, Tarantola S (2008) _Global Sensitivity Analysis: The Primer_ (John Wiley & Sons, Ltd.).
  * Saltelli et al. (2000) Saltelli A, Tarantola S, Campolongo F (2000) Sensitivity anaysis as an ingredient of modeling. _Statistical Science_ 15(4):377–395.
  * Scutari (2010) Scutari M (2010) Learning Bayesian networks with the bnlearn R package. _Journal of Statistical Software_ 35(3):1–22.
  * Sheikholeslami et al. (2020) Sheikholeslami R, Gharari S, Papalexiou SM, Clark MP (2020) Viscous: A variance-based sensitivity analysis using copulas for efficient identification of dominant hydrological processes. _Earth and Space Science Open Archive_ 37\.
  * Sobol (1990) Sobol IM (1990) Sensitivity estimates for nonlinear mathematical models (in Russian). _Mathematical Models_ 2:112–118.
  * Stoudenmire and Schwab (2016) Stoudenmire E, Schwab DJ (2016) Supervised learning with tensor networks. Lee D, Sugiyama M, Luxburg U, Guyon I, Garnett R, eds., _Advances in Neural Information Processing Systems_ , volume 29.
  * Sudret (2008) Sudret B (2008) Global sensitivity analysis using polynomial chaos expansions. _Reliability Engineering & System Safety_ 93(7):964 – 979.
  * Ye and Lim (2018) Ye K, Lim LH (2018) Tensor network ranks. _arXiv:1801.02662_ .
  * Zhang et al. (2015) Zhang Z, Yang X, Oseledets IV, Karniadakis GE, Daniel L (2015) Enabling high-dimensional hierarchical uncertainty quantification by ANOVA and tensor-train decomposition. _IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems_ 34(1):63–76.