# Can I Trust the Explanations? Investigating Explainable Machine Learning
Methods for Monotonic Models

Dangxing Chen  Kunshan Duke UniversityKunshanChina
<EMAIL_ADDRESS>
###### Abstract.

In recent years, explainable machine learning methods have been very
successful. Despite their success, most explainable machine learning methods
are applied to black-box models without any domain knowledge. By incorporating
domain knowledge, science-informed machine learning models have demonstrated
better generalization and interpretation. But do we obtain consistent
scientific explanations if we apply explainable machine learning methods to
science-informed machine learning models? This question is addressed in the
context of monotonic models that exhibit three different types of
monotonicity. To demonstrate monotonicity, we propose three axioms.
Accordingly, this study shows that when only individual monotonicity is
involved, the baseline Shapley value provides good explanations; however, when
strong pairwise monotonicity is involved, the Integrated gradients method
provides reasonable explanations on average.

explainable machine learning, monotonicity, fairness, neural networks

††doi:

##  1\. Introduction

In recent decades, machine learning (ML) models have achieved many successes.
In comparison with traditional methods, machine learning models are often
capable of increasing accuracy at the expense of black-box functionality. The
importance of model explanation is particularly important for highly regulated
industries such as the financial sector (OCC, 2021). As an example, the
Consumer Financial Protection Bureau (CFPB) confirmed that anti-discrimination
law requires companies to provide detailed explanations when denying an
application for credit when using machine learning methods
111https://www.consumerfinance.gov/about-us/newsroom/cfpb-acts-to-protect-the-
public-from-black-box-credit-models-using-complex-algorithms/. In response to
the growing regulatory requirements, researchers are investigating explainable
machine learning methods.

Explainable machine learning (XML) methods have been successfully used in the
past to achieve great success in machine learning. Popular used methods
include SHapley Additive exPlanations (SHAP) (Lundberg and Lee, 2017), Local
Interpretable Model-Agnostic Explanations (LIME) (Ribeiro et al., 2016),
Integrated Gradients (IG) (Sundararajan et al., 2017), Anchors (Ribeiro et
al., 2018), and Sensitivity-based methods (Horel and Giesecke, 2020). As a
result of these methods, we have gained a better understanding of how ML
models function (Boardman et al., 2022; Bussmann et al., 2021; Horel et al.,
2018).

In this paper, we address the attribution problem, which involves allocating
the prediction score of a model for a given input to its base features. The
attribution to a base feature can be understood as the importance of the
feature in the prediction. Credit scoring, for instance, can utilize
attribution to understand how each feature contributes to the credit score.
The SHAP and IG have been successfully applied to the attribution problem.
Further, their use has been demonstrated to comply with a number of
theoretical properties, as discussed in (Lundberg and Lee, 2017; Sundararajan
et al., 2017; Sundararajan and Najmi, 2020; Lundstrom et al., 2022). While
extensive analyses have been conducted, the majority of results have been
based on black-box machine learning models without domain knowledge.

Modeling and ensuring conceptual soundness require domain knowledge: the true
model should be consistent with the underlying theories. A number of studies
have demonstrated that physics-informed machine learning (PIML) (Karniadakis
et al., 2021; Greydanus et al., 2019) improved black-box machine learning
models in terms of interpretation and accuracy by enforcing conservation laws,
for example. Finance and other applications often require monotonicity. A
person’s credit score should be decreased when there is one more past due
balance on the account, for example. It is possible to achieve better
generalization and interpretation when monotonicity is successfully enforced
(Liu et al., 2020; Milani Fard et al., 2016; You et al., 2017; Repetto, 2022).
Such models can be categorized as finance-informed machine learning models
(FIML) or more generally science-informed machine learning models (SIML). In
addition, monotonicity is often associated with fairness. As an example, with
other factors being equal, a person with more past dues should have a lower
credit score. The violation of monotonicity may result in unfair consequences
that could cause damage to our society.

In this paper, we ask the following question. Can attribution methods deliver
consistent scientific explanations if SIML models contain certain scientific
knowledge? If so, to what extent? Specifically, do attribution methods
preserve monotonicity for monotonic ML models?

In the past, monotonicity has been considered for XML methods. In
(Sundararajan and Najmi, 2020), it is shown that among different SHAP methods,
the baseline Shapley value (BShap) method, which is a generalization of
Shapley-Shubik (Friedman and Moulin, 1999), preserves demand individual
monotonicity. However, as recently highlighted in (Chen and Ye, 2022, 2023;
Gupta et al., 2020; Cotter et al., 2019), individual monotonicity is not the
only problem to be addressed; pairwise monotonicity is just as important. The
concept of pairwise monotonicity refers to the comparison of features within a
given pair. For instance, a past due of more than 60 days should be considered
more serious than one of less than 60 days. Pairwise monotonicity is a
requirement of fairness that is informed by domain knowledge. It is
unfortunate that pairwise monotonicity has been neglected in the existing
literature.

This paper extends monotonicity results to a broader range of cases. As a
summary, we have made the following contributions:

  1. (1)

We propose three new axioms concerning the average preservation of three types
of monotonicity. Accordingly, a good attribution method should satisfy these
axioms if the underlying model is monotonic.

  2. (2)

As we show, IG preserves all average monotonicity axioms, but failing to
preserve demand individual monotonicity.

  3. (3)

While BShap preserves individual monotonicity and weak pairwise monotonicity,
it fails to preserve average strong pairwise monotonicity.

A number of examples are used to illustrate our results. BShap is appropriate
when only individual monotonicity is considered and is better suited if demand
individual monotonicity is required. On the other hand, when strong pairwise
monotonicity is involved, IG provides reasonable explanations on average.

The remainder of the paper is organized as follows. Section 2 introduces
attributions, Integrated gradients, (baseline) Shapley value, axioms, and
types of monotonicity. In Section 3, we propose three monotonic-related axioms
and analyze IG and BShap for these axioms. In Section 4, we present an
empirical example. Section 5 concludes.

##  2\. Preliminaries

For problem setup, assume we have $\mathcal{D}\times\mathcal{Y}$, where
$\mathcal{D}$ is the dataset with $n$ samples and $m$ features and
$\mathcal{Y}$ is the corresponding numerical values in regression and labels
in classification. For $\mathbf{a},\mathbf{b}\in\mathbb{R}^{m}$, define
$[\mathbf{a},\mathbf{b}]$ to be the hyperrectangle. We denote a class of
functions $f:[\mathbf{a},\mathbf{b}]\rightarrow\mathbb{R}$ by
$\mathcal{F}(\mathbf{a},\mathbf{b})$, or simply $\mathcal{F}$. For simplicity,
we assume $\mathbf{x}\in[\mathbf{a},\mathbf{b}]^{m}$ and $f$ is differentiable
almost everywhere.

###  2.1. Attribution

Following (Lundstrom et al., 2022), we call the point of interest $\mathbf{x}$
to explain as an explicand and $\mathbf{x}^{\prime}$ a baseline. The Baseline
Attribution Method that interprets features’ importance is defined as follows.

######  Definition 2.1 (Baseline Attribution Method (BAM)).

Given $\mathbf{x},\mathbf{x}^{\prime}\in[\mathbf{a},\mathbf{b}]$,
$f\in\mathcal{F}(\mathbf{a},\mathbf{b})$, a baseline attribution method is any
function of the form
$\mathbf{A}:[\mathbf{a},\mathbf{b}]\times[\mathbf{a},\mathbf{b}]\times\mathcal{F}(\mathbf{a},\mathbf{b})\rightarrow\mathbb{R}^{m}$
with $\mathbf{A}(\mathbf{x}^{\prime},\mathbf{x}^{\prime},f)=\mathbf{0}$. We
may also write $\mathbf{A}$ and denote $A_{i}$ as the $i$th attribution of
$\mathbf{A}$ for simplicity.

####  2.1.1. Integrated Gradients

Integrated Gradients (Sundararajan et al., 2017) is a simple yet useful tool
to measure features’ importance.

######  Definition 2.2 (Integrated Gradients (IG)).

Given $\mathbf{x},\mathbf{x}^{\prime}\in[a,b]$ and
$f\in\mathcal{F}(\mathbf{a},\mathbf{b})$, the integrated gradients attribution
of the $i$-th component of $\mathbf{x}$ is defined as

(1)
$\displaystyle\text{IG}_{i}(\mathbf{x},\mathbf{x}^{\prime},f)=(x_{i}-x_{i}^{\prime})\int_{0}^{1}\frac{\partial
f}{\partial
x_{i}}\left(\mathbf{x}^{\prime}+t(\mathbf{x}-\mathbf{x}^{\prime})\right)\ dt.$

In practice, this is usually calculated through numerical approximation

(2)
$\displaystyle\text{IG}_{i}(\mathbf{x},\mathbf{x}^{\prime},f)\approx(x_{i}-x_{i}^{\prime})\frac{1}{n}\sum_{k=1}^{n}\frac{\partial
f(\mathbf{x}^{\prime}+\frac{k}{n}(\mathbf{x}-\mathbf{x}^{\prime}))}{\partial
x_{i}},$

with 20 - 300 calls of gradient. For simplicity, we often use $\text{IG}_{i}$
for $\text{IG}_{i}(\mathbf{x},\mathbf{x}^{\prime},f)$.

####  2.1.2. (Baseline) Shapley Value (BShap)

The Shapley value (Lundberg and Lee, 2017) takes as input a set function
$v:2^{M}\rightarrow\mathbb{R}$. The Shapley value produces attributions
$s_{i}$ for each player $i\in M$ that add up to $v(M)$. The Shapley value of a
player $i$ is given by:

(3) $\displaystyle s_{i}=\sum_{S\subseteq M\backslash
i}\frac{|S|!(|M|-|S|-1)!}{M!}(v(S\cup i)-v(S)).$

Here, we focus on the Baseline Shapley (Shap), which calculates

(4) $\displaystyle v(S)=f(\mathbf{x}_{S};\mathbf{x}^{\prime}_{M\backslash
S}).$

That is, baseline values replace the feature’s absence. We denote Bshap
attribution by $\text{BS}_{i}(\mathbf{x},\mathbf{x}^{\prime},f)$ and
$\text{BS}_{i}$ sometimes.

###  2.2. Axioms

Many desirable characteristics of an attribution technique have been
identified in the literature, such as sensitivity, implementation variance,
and completeness. Interested readers are referred to (Lundstrom et al., 2022;
Sundararajan and Najmi, 2020; Sundararajan et al., 2017) for detailed
discussion. This paper focuses on the axiom of monotonicity. (Sundararajan and
Najmi, 2020) introduced demand individual monotonicity as one of the desired
axioms. The attributions for a feature increase with increasing feature values
if the model is monotonic with respect to that feature, all else being fixed,
including the baseline value. The paper (Sundararajan and Najmi, 2020) showed
that BShap is a powerful tool that preserves the monotonicity of demand
individual demands and other desirable characteristics.

###  2.3. Individual and pairwise monotonicity

Three types of monotonicity are discussed here. Without loss of generality, we
assume that all monotonic features are monotonically increasing throughout the
paper. Suppose $\boldsymbol{\alpha}$ is the set of all individual monotonic
features and $\neg\boldsymbol{\alpha}$ its complement, then the input
$\mathbf{x}$ can be partitioned into
$\mathbf{x}=(\mathbf{x}_{\boldsymbol{\alpha}},\mathbf{x}_{\neg\boldsymbol{\alpha}})$.
Individual monotonicity is defined as follows.

######  Definition 2.3 (Individual Monotonicity).

We say $f$ is individually monotonic with respect to
$\mathbf{x}_{\boldsymbol{\alpha}}$ if

$\displaystyle
f(\mathbf{x})=f(\mathbf{x}_{\boldsymbol{\alpha}},\mathbf{x}_{\neg\boldsymbol{\alpha}})\leq
f(\mathbf{x}^{*}_{\boldsymbol{\alpha}},\mathbf{x}_{\neg\boldsymbol{\alpha}})=f(\mathbf{x}^{*}),$
(5)
$\displaystyle\forall\mathbf{x},\mathbf{x}^{*}\in[\mathbf{a},\mathbf{b}]\text{
s.t.
}\mathbf{x}_{\boldsymbol{\alpha}}\leq\mathbf{x}^{*}_{\boldsymbol{\alpha}},$

where
$\mathbf{x}_{\boldsymbol{\alpha}}\leq\mathbf{x}_{\boldsymbol{\alpha}}^{\prime}$
means $x_{\alpha_{i}}\leq x_{\alpha_{i}}^{\prime},\forall{i}$.

######  Example 2.4.

In the case of credit scoring, given that $f$ represents the probability of
default and $x_{\alpha}$ represents the number of past dues, then $f$ should
increase monotonically with respect to $x_{\alpha}$.

In practice, certain features are intrinsically more important than others.
Analog to (2.3), we partition
$\mathbf{x}=(x_{\beta},x_{\gamma},\mathbf{x}_{\neg})$. Without sacrificing
generality, we assume that $x_{\beta}$ has greater significance than
$x_{\gamma}$. Lastly, we require that all features exhibiting pairwise
monotonicity also exhibit individual monotonicity. Pairwise monotonicity can
be categorized into two types: strong and weak. As a more general definition,
weak pairwise monotonicity is presented below.

######  Definition 2.5 (Weak Pairwise Monotonicity).

We say $f$ is weakly monotonic with respect to $x_{\beta}$ over $x_{\gamma}$
if

$\displaystyle f(\mathbf{x})=f(x_{\beta},x_{\gamma}+c,\mathbf{x}_{\neg})\leq
f(x_{\beta}+c,x_{\gamma},\mathbf{x}_{\neg})=f(\mathbf{x}^{*}),$ (6)
$\displaystyle\forall\mathbf{x},\mathbf{x}^{*}\in[\mathbf{a},\mathbf{b}]\text{
s.t. }x_{\beta}=x_{\gamma},c>0.$

Weak pairwise monotonicity compares the significance of $x_{\beta}$ and
$x_{\gamma}$ at the same magnitude. The following is an example.

######  Example 2.6.

An applicant who intends to study STEM in graduate school is required to take
the GRE general (170 in math and 170 in verbal) as one of the factors for
admission. Assume $f$ is the admitted probability, $x_{\beta}$ is the
student’s math score, and $x_{\gamma}$ is the student’s verbal score. Due to
the importance of math in STEM, $f$ should be weakly monotonic with respect to
$x_{\beta}$ over $x_{\gamma}$. In this case, it is necessary to fulfill the
condition $x_{\beta}=x_{\gamma}$. If the student has the same level of math
and verbal skills, it is more desirable to see improvement in math. It is a
different story when a student has a strong math score but a weak verbal
score. There is often a requirement for a minimum verbal score in schools in
order to ensure effective communication. A student who has a strong math score
but a very weak verbal score may be able to improve his/her chances of
admission more significantly if he/she improves his/her verbal score.

In addition, there is a stronger condition of pairwise monotonicity, known as
strong pairwise monotonicity, which is independent of the condition that
$x_{\beta}=x_{\gamma}$. Here is the definition.

######  Definition 2.7 (Strong Pairwise Monotonicity).

We say $f$ is strongly monotonic with respect to $x_{\beta}$ over $x_{\gamma}$
if

$\displaystyle f(\mathbf{x})=f(x_{\beta},x_{\gamma}+c,\mathbf{x}_{\neg})\leq
f(x_{\beta}+c,x_{\gamma},\mathbf{x}_{\neg})=f(\mathbf{x}^{*}),$ (7)
$\displaystyle\forall\mathbf{x},\mathbf{x}^{*}\in[\mathbf{a},\mathbf{b}]\text{
s.t. }c>0.$

Accordingly, $x_{\beta}$ is always more important than $x_{\gamma}$ regardless
of its value. Below is an example.

######  Example 2.8.

In credit scoring, suppose $f$ calculates the probability of default,
$x_{\beta}$ counts the number of past dues that have been outstanding for more
than two months, and $x_{\gamma}$ counts the number of past dues that are
outstanding between one and two months. As $x_{\beta}$ is always more
important than $x_{\gamma}$, $f$ should be strongly monotonic with respect to
$x_{\beta}$ over $x_{\gamma}$. The pairwise relationship is strong here since
whenever there is an additional past due, a longer past due is always more
serious.

##  3\. Monotonic Preserving

In related applications, monotonicity is an important component in domain
knowledge. This section discusses axioms derived from different types of
monotonicity and analyzes IG and BShap behaviors.

###  3.1. Demand Individual Monotonicity

Consider two explicands $\mathbf{x}=(x_{1},\dots,x_{m})$ and
$\mathbf{x}^{*}=(x_{1},\dots,x_{\alpha}+c,\dots,x_{m})$, where $c>0$. We
restate the axiom of demand individual monotonicity proposed in (Friedman and
Moulin, 1999).

######  Definition 3.1 (Demand Individual Monotonicity (DIM)).

Suppose $f$ is individually monotonic with respect to $x_{\alpha}$. We say a
BAM preserves demand individual monotonicity if

(8) $\displaystyle A_{\alpha}(\mathbf{x}^{*},\mathbf{x}^{\prime},f)\geq
A_{\alpha}(\mathbf{x},\mathbf{x}^{\prime},f),\forall\mathbf{x},\mathbf{x}^{*}\in[\mathbf{a},\mathbf{b}].$

######  Example 3.2.

Consider Example 2.4, whenever there is an additional past due, the
attribution of $x_{\alpha}$ should be increased since $x_{\alpha}$ caused the
change.

As discussed in (Sundararajan and Najmi, 2020), BShap preserves DIM. It is
presented here without the contents of the cost-sharing problem.

######  Theorem 3.3.

BShap preserves DIM.

###### Proof.

Without loss of generality, suppose $f$ is individually monotonic with respect
to $x_{1}$. Suppose $\mathbf{x}=(x_{1},\mathbf{x}_{\neg})$ and
$\mathbf{x}^{*}=(x_{1}+c,\mathbf{x}_{\neg})$ for $c>0$. Then

$\displaystyle\text{BS}_{1}(\mathbf{x}^{*},\mathbf{x}^{\prime},f)-\text{BS}_{1}(\mathbf{x},\mathbf{x}^{\prime},f)$
$\displaystyle=\sum_{S\subseteq M\backslash
1}\frac{|S|!(|M|-|S|-1)!}{M!}\Delta f_{S},$

where

$\displaystyle\Delta
f_{S}=f(x_{1}+c;\mathbf{x}_{S},\mathbf{x}_{M\backslash(S\cup
1)}^{\prime})-f(x_{1};\mathbf{x}_{S},\mathbf{x}_{M\backslash(S\cup
1)}^{\prime})\geq 0$

because of individual monotonicity. Thus, we conclude. ∎

IG has unfortunately produced a negative result.

######  Remark 3.4.

IG doesn’t preserve DIM.

The following example in (Friedman and Moulin, 1999) provides a counterexample
with a comparison between IG and BShap.

######  Example 3.5.

Consider the function $f(x_{1},x_{2})=\frac{x_{1}x_{2}}{x_{1}+x_{2}}$ with
baseline $\mathbf{x}^{\prime}=(0,0)$. For IG, it can be shown that
$\text{IG}_{1}(\mathbf{x},\mathbf{x}^{\prime},f)=\frac{x_{1}x_{2}^{2}}{(x_{1}+x_{2})^{2}}$.
Note $\frac{\partial}{\partial
x_{1}}\text{IG}_{1}=\frac{x_{2}^{2}(x_{2}-x_{1})}{(x_{1}+x_{2})^{3}}$,
therefore DIM is not preserved. It is important to note that the path for the
new explicand has been changed. As a result, the new path cannot guarantee
greater attributions for the explicand.

For BShap, we have
$\text{BS}_{1}(\mathbf{x},\mathbf{x}^{\prime},f)=\frac{x_{1}x_{2}}{2(x_{1}+x_{2})}$,
which preserves DIM. The main difference is that in the BShap calculation, the
new path covers the old path.

###  3.2. Average Individual Monotonicity

DIM concerns individual monotonicity for each change. The average result is
considered here as a weaker case. The following axiom is motivated by this
fact.

######  Definition 3.6 (Average Individual Monotonicity (AIM)).

Suppose $f$ is individually monotonic with respect to $x_{\alpha}$, then we
say a BAM preserves average individual monotonicity if

(9) $\displaystyle A_{\alpha}(\mathbf{x},\mathbf{x}^{\prime},f)\geq
0,\forall\mathbf{x}\in[\mathbf{a},\mathbf{b}]\text{ s.t.
}\mathbf{x}\geq\mathbf{x}^{\prime}.$

######  Example 3.7.

Consider Example 2.4, $x_{\alpha}$ should always provide positive
contributions.

Clearly, DIM is a stronger condition, as stated in the following.

######  Theorem 3.8.

If a BAM preserves DIM, then it preserves AIM.

###### Proof.

Suppose we want to explain $\mathbf{x}^{*}$, let
$\mathbf{x}=\mathbf{x}^{\prime}$ in Definition 3.1, then it follows from the
definition. ∎

As a result, we have the following corollary.

###### Corollary 3.9.

BShap preserves AIM.

In spite of the fact that the IG fails to preserve DIM, it preserves AIM.

######  Theorem 3.10.

IG preserves AIM.

###### Proof.

Suppose $f$ is individually monotonic with respect to $x_{\alpha}$, then

$\displaystyle\text{IG}_{\alpha}$
$\displaystyle=(x_{\alpha}-x_{\alpha}^{\prime})\int_{0}^{1}\frac{\partial
f}{\partial
x_{\alpha}}\left(\mathbf{x}^{\prime}+t(\mathbf{x}-\mathbf{x}^{\prime})\right)\
dt\geq 0,$

since $x_{\alpha}\geq x_{\alpha}^{\prime}$ and $\frac{\partial f}{\partial
x_{\alpha}}\geq 0$. ∎

As shown in Example 3.5, if we consider a new path, IG will still contribute
positively to the feature, but the contribution of the features may decrease.
To this end, we ask: does the failure of DIM pose a significant concern? DIM
is derived from a somewhat implicit baseline of $\mathbf{x}$. It is important
to recall that the increase is relative to $\mathbf{x}$. However, the
attribution we seek to calculate is essentially an average calculation.
Consider the following example.

######  Example 3.11.

Consider Example 2.8 with a concrete example, where
$f(x_{1},x_{2})=\min(0.2x_{1}+0.1x_{2},0.3)$ and $x_{1}$ is strongly monotonic
with respect to $x_{2}$. At $\mathbf{x}=(1,1)$ with baseline
$\mathbf{x}^{\prime}=(0,0)$, the corresponding attribution should be
$\mathbf{A}(\mathbf{x},\mathbf{x}^{\prime},f)=(0.2,0.1)$, due to linearity.
Now we consider $\mathbf{x}^{*}=(3,1)$, $f(\mathbf{x}^{*})=0.3$ and DIM
requires that $A_{1}(\mathbf{x}^{*},\mathbf{x}^{\prime},f)\geq 0.2$. In this
case, we wish to satisfy this requirement. However, if we consider
$\mathbf{x}^{!}=(0,4)$ with $f(x^{!})=0.3$. From the perspective of $x^{!}$,
$x_{1}$ does not incur an additional cost. Accordingly, a unit of $x_{1}$ and
$x_{2}$ should be attributed equally. This results in attributions at $(2,2)$
of $(0.15,0.15)$, which violates the principle of DIM.

######  Remark 3.12.

In light of Example 3.11, we argue that DIM may not be necessary in some
cases. Although DIM makes sense for some features, it can be achieved through
the use of a historical baseline, as in (Alam et al., 2022). Taking into
consideration the time at which there is an increase in monotonic features of
$\mathbf{x}_{t}$ and we have $\mathbf{x}_{t+1}$, then AIM provides positive
attributions for $A(\mathbf{x}_{t+1},\mathbf{x}_{t},f)$, which demand
individual monotonicity.

###  3.3. Average Weak Pairwise Monotonicity

Weak pairwise monotonicity compares the importance of features of the same
magnitude. Therefore, if two features have the same magnitude, we would expect
the more important feature to have more attributions. In light of this, we
propose the following axiom.

######  Definition 3.13 (Average Weak Pairwise Monotonicity (AWPM)).

Suppose $f$ is weakly monotonic with respect to $x_{\beta}$ over $x_{\gamma}$,
$x_{\beta}>x_{\beta}^{\prime}$ and $x_{\gamma}>x_{\gamma}^{\prime}$. Suppose
for an explicand $\mathbf{x}\in[\mathbf{a},\mathbf{b}]$, we have
$x_{\beta}=x_{\gamma}$. Then we say a BAM preserves average weak pairwise
monotonicity if

(10)
$\displaystyle\frac{1}{x_{\beta}-x_{\beta}^{\prime}}A_{\beta}(\mathbf{x},\mathbf{x}^{\prime},f)\geq\frac{1}{x_{\gamma}-x_{\gamma}^{\prime}}A_{\gamma}(\mathbf{x},\mathbf{x}^{\prime},f).$

######  Example 3.14.

Consider Example 2.6, math scores should contribute more attributions when a
student has the same verbal and math scores.

IG and BShap both preserve AWPM with appropriate baselines. For simplicity, we
only show results for $\mathbf{x}^{\prime}=\mathbf{0}$ here.

######  Theorem 3.15.

For $\mathbf{x}^{\prime}=\mathbf{0}$, IG preserves AWPM.

###### Proof.

Suppose $f$ is weakly monotonic with respect to $x_{\beta}$ over $x_{\gamma}$.
For $\mathbf{x}^{\prime}=\mathbf{0}$, we have
$\frac{1}{x_{\beta}-x_{\beta}^{\prime}}=\frac{1}{x_{\gamma}-x_{\gamma}^{\prime}}$.
Since $x_{\beta}=x_{\gamma}$ in $\mathbf{x}$ and $\mathbf{x}^{\prime}$,
$x_{\beta}=x_{\gamma}$ for
$\mathbf{x}^{\prime}+t(\mathbf{x}-\mathbf{x}^{\prime})$, $\forall t\in[0,1]$.
$f$ is weakly monotonic respect to $x_{\beta}$ over $x_{\gamma}$, therefore
$\frac{\partial f}{\partial x_{\beta}}(\mathbf{x})\geq\frac{\partial
f}{\partial x_{\gamma}}(\mathbf{x})$ if $x_{\beta}=x_{\gamma}$ in
$\mathbf{x}$. Hence,

$\displaystyle\text{IG}_{\beta}$
$\displaystyle=(x_{\beta}-x_{\beta}^{\prime})\int_{0}^{1}\frac{\partial
f}{\partial
x_{\beta}}\left(\mathbf{x}^{\prime}+t(\mathbf{x}-\mathbf{x}^{\prime})\right)\
dt$ $\displaystyle\geq(x_{\beta}-x_{\beta}^{\prime})\int_{0}^{1}\frac{\partial
f}{\partial
x_{\gamma}}\left(\mathbf{x}^{\prime}+t(\mathbf{x}-\mathbf{x}^{\prime})\right)\
dt$ $\displaystyle=(x_{\gamma}-x_{\gamma}^{\prime})\int_{0}^{1}\frac{\partial
f}{\partial
x_{\gamma}}\left(\mathbf{x}^{\prime}+t(\mathbf{x}-\mathbf{x}^{\prime})\right)\
dt$ $\displaystyle=\text{IG}_{\gamma}.$

∎

######  Theorem 3.16.

For $\mathbf{x}^{\prime}=\mathbf{0}$, BShap preserves AWPM.

###### Proof.

Without loss of generality, suppose $f$ is weakly monotonic with respect to
$x_{1}$ over $x_{2}$. For $\mathbf{x}^{\prime}=\mathbf{0}$, we have
$\frac{1}{x_{1}-x_{1}^{\prime}}=\frac{1}{x_{2}-x_{2}^{\prime}}$. Recall that

$\displaystyle\text{BS}_{1}=\sum_{S\subseteq M\backslash
1}\frac{|S|!(|M|-|S|-1)!}{M!}(v(S\cup 1)-v(S)).$

For $\text{BS}_{2}$, we use the symmetry, for each $S$ here, we consider
$S^{\prime}$ such that $x_{1}$ and $x_{2}$ are swapped within $S$, and
everything else is left unchanged. That is, if $2\notin S$, $S^{\prime}=S$; if
$2\in S$, then $S^{\prime}=(S\backslash 2)\cup 1$. For this arrangement, we
only need to show that $v(S\cup 1)-v(S)\geq v(S^{\prime}\cup 2)-v(S^{\prime})$
in the summation. If $2\notin S$, we have

$\displaystyle v(S\cup 1)-v(S)$
$\displaystyle=f(x_{1},0;\mathbf{x}_{S},\mathbf{x}_{M\backslash(S\cup\\{1,2\\})}’)$
$\displaystyle-f(0,0;\mathbf{x}_{S},\mathbf{x}_{M\backslash(S\cup\\{1,2\\}}’),$
$\displaystyle v(S^{\prime}\cup 2)-v(S^{\prime})$
$\displaystyle=f(0,x_{2};\mathbf{x}_{S^{\prime}},\mathbf{x}_{M\backslash(S^{\prime}\cup\\{1,2\\})}’)$
$\displaystyle-f(0,0;\mathbf{x}_{S^{\prime}},\mathbf{x}_{M\backslash(S^{\prime}\cup\\{1,2\\})}’).$

Since $S=S^{\prime}$ and $x_{1}=x_{2}$, by the definition of weak pairwise
monotonicity, we have

$\displaystyle
f(x_{1},0;\mathbf{x}_{S},\mathbf{x}_{M\backslash(S\cup\\{1,2\\})}’)\geq
f(0,x_{2};\mathbf{x}_{S^{\prime}},\mathbf{x}_{M\backslash(S^{\prime}\cup\\{1,2\\})}’),$

and therefore $v(S\cup 1)-v(S)\geq v(S^{\prime}\cup 2)-v(S^{\prime})$.

If $2\in S$, and we have

$\displaystyle v(S\cup 1)-v(S)$
$\displaystyle=f(x_{1},x_{2};\mathbf{x}_{S\backslash
2},\mathbf{x}_{M\backslash(S\cup 1)}’)$
$\displaystyle-f(0,x_{2};\mathbf{x}_{S\backslash
2},\mathbf{x}_{M\backslash(S\cup 1)}’),$ $\displaystyle v(S^{\prime}\cup
2)-v(S^{\prime})$
$\displaystyle=f(x_{1},x_{2};\mathbf{x}_{S^{\prime}\backslash
1},\mathbf{x}_{M\backslash(S^{\prime}\cup 2)}’)$
$\displaystyle-f(x_{1},0;\mathbf{x}_{S^{\prime}\backslash
1},\mathbf{x}_{M\backslash(S^{\prime}\cup 2)}’).$

Since $S\backslash 2=S^{\prime}\backslash 1$, $S\cup 1=S^{\prime}\cup 2$, and
$x_{1}=x_{2}$, by the definition of weak pairwise monotonicity, we have

$\displaystyle f(0,x_{2};\mathbf{x}_{S\backslash
2},\mathbf{x}_{M\backslash(S\cup 1)}’)\leq
f(x_{1},0;\mathbf{x}_{S^{\prime}\backslash
1},\mathbf{x}_{M\backslash(S^{\prime}\cup 2)}’),$

and therefore $v(S\cup 1)-v(S)\geq v(S^{\prime}\cup 2)-v(S^{\prime})$.

Since $v(S\cup 1)-v(S)\geq v(S^{\prime}\cup 2)-v(S^{\prime})$ for all $S$ with
corresponding $S^{\prime}$, we conclude that $\text{BS}_{1}\geq\text{BS}_{2}$.

∎

This result also suggests baseline points. We need
$x_{\beta}^{\prime}=x_{\gamma}^{\prime}$ in $\mathbf{x}^{\prime}$ in order to
preserve AWPM.

######  Example 3.17.

Consider $f(x_{1},x_{2})=4.5x_{1}-x_{1}^{2}+4x_{2}-x_{2}^{2}$ for
$x_{1},x_{2}$ in $[0,2]$. $f$ is weakly monotonic with respect to $x_{1}$ over
$x_{2}$. Consider the baseline $\mathbf{x}^{\prime}=(0,0)$ and the explicand
$\mathbf{x}=(2,2)$.

For IG, we have

$\displaystyle\textbf{IG}((2,2),(0,0),f)=\left[\begin{matrix}5\\\
4\end{matrix}\right].$

It can be seen that $\frac{1}{2}\text{IG}_{1}>\frac{1}{2}\text{IG}_{2}$, IG
preserves AWPM. Now suppose we consider $\mathbf{x}^{\prime}=(1,0)$, we have

$\displaystyle\textbf{IG}((2,2),(1,0),f)=\left[\begin{matrix}1.5\\\
4\end{matrix}\right].$

Since $\text{IG}_{1}<\frac{1}{2}\text{IG}_{2}$, IG fails to preserve AWPM. The
result of BShap is the same as that of $f$ being additively separable.

Consider Example 2.6. For $x_{\beta}^{\prime}=x_{\gamma}^{\prime}=0$, the weak
pairwise monotonicity in $f$ guarantees that math is more important than
verbal when they are equal when we determine the general importance of
features. If, on the other hand, we take the average historical statistics of
admitted students as our baseline, then it would be possible to see that
verbal skills are more important than mathematics skills as admitted students
are already competent in mathematics. It is important to keep in mind in this
case that we are comparing the cases of previously admitted students, and the
importance of features does not necessarily apply in general.

###  3.4. Average Strong Pairwise Monotonicity

When there is strong pairwise monotonicity, one feature is always more
important than another. Consequently, when considering attributions, the
average attribution of the more important ones should always be greater. In
light of this, we should consider the following axiom.

######  Definition 3.18 (Average Strong Pairwise Monotonicity (ASPM)).

Suppose $f$ is strongly monotonic with respect to $x_{\beta}$ over
$x_{\gamma}$, $x_{\beta}>x_{\beta}^{\prime}$,
$x_{\gamma}>x_{\gamma}^{\prime}$, and $\mathbf{x}\in[\mathbf{a},\mathbf{b}]$.
Then we say a BAM preserves average strong pairwise monotonicity if

(11)
$\displaystyle\frac{1}{x_{\beta}-x_{\beta}^{\prime}}A_{\beta}(\mathbf{x},\mathbf{x}^{\prime},f)\geq\frac{1}{x_{\gamma}-x_{\gamma}^{\prime}}A_{\gamma}(\mathbf{x},\mathbf{x}^{\prime},f).$

######  Example 3.19.

Consider Example 2.8. As a general rule, we should always expect $x_{\beta}$
to contribute more than $x_{\gamma}$ on average, since $x_{\beta}$ is of
greater importance. As an alternative, the interpretation may mislead people
to believe that payments between one and two months past due can be more
serious.

Clearly, AWPM preserving is a special case of ASPM for $x_{\beta}=x_{\gamma}$.
Thus, we have the following theorem.

######  Theorem 3.20.

If a BAM preserves ASPM, it also preserves AWPM for the same features.

As a result of IG, ASPM is preserved and we provide the results.

######  Theorem 3.21.

IG preserves ASPM.

###### Proof.

Suppose $f$ is strongly monotonic with respect to $x_{\beta}$ over
$x_{\gamma}$, then $\frac{\partial f}{\partial
x_{\gamma}}(\mathbf{x})\leq\frac{\partial f}{\partial x_{\beta}}(\mathbf{x})$,
$\forall\mathbf{x}\in[\mathbf{a},\mathbf{b}]$. Hence,

$\displaystyle\frac{1}{x_{\beta}-x_{\beta}^{\prime}}\text{IG}_{\beta}$
$\displaystyle=\int_{0}^{1}\frac{\partial f}{\partial
x_{\beta}}\left(\mathbf{x}^{\prime}+t(\mathbf{x}-\mathbf{x}^{\prime})\right)\
dt$ $\displaystyle\geq\int_{0}^{1}\frac{\partial f}{\partial
x_{\gamma}}\left(\mathbf{x}^{\prime}+t(\mathbf{x}-\mathbf{x}^{\prime})\right)\
dt$
$\displaystyle=\frac{1}{x_{\gamma}-x_{\gamma}^{\prime}}\text{IG}_{\gamma}.$

∎

BShap doesn’t preserve ASPM, as we documented here.

######  Remark 3.22.

BShap doesn’t preserve ASPM.

A counterexample is provided here, as well as comparisons with the IG.

######  Example 3.23.

Consider the function with two features
$f(x_{1},x_{2})=\log(1+10x_{1}+9x_{2})$, where $f$ is strongly monotonic with
respect to $x_{1}$ over $x_{2}$. This example mimics the diminishing marginal
effect, that is, the impact of $f$ decreases as $x_{1}+x_{2}$ increases.
Consider the explicand $\mathbf{x}=(4,1)$ and baseline
$\mathbf{x}^{\prime}=(0,0)$.

By calculations, we have $\textbf{BS}\approx\left[\begin{matrix}2.7\\\
1.3\end{matrix}\right]$. It is evident that ASPM is not preserved. BShap
requires a stronger condition to preserve the ASPM, for example, if the
averaged change of $x_{\beta}$ always has a greater impact on $f$. Such a
condition, however, is not realistic in applications such as credit scoring.
In addition, we apply the IG, and we obtain
$\textbf{IG}\approx\left[\begin{matrix}3.2\\\ 0.7\end{matrix}\right]$. ASPM is
successfully maintained.

Consider Example 2.8. Using BShap, at first glance, longer past dues appear to
be contributing more to the result. On second thought, it seems that averaged
longer past dues make a smaller contribution. Providing the customer with this
interpretation would cause the customer to incorrectly believe that longer
past dues have less impact, leading to a catastrophic outcome.

######  Remark 3.24.

As discussed in Lemma 4.8 in (Sundararajan and Najmi, 2020), BShap fits in the
SHAP framework. Therefore, SHAP doesn’t preserve ASPM as well.

##  4\. Empirical Examples

###  4.1. Data description

We use the Kaggle credit score dataset
222https://www.kaggle.com/c/GiveMeSomeCredit/overview. For simplicity, data
with missing variables are removed. Past dues greater than four times are
truncated. The dataset is randomly partitioned into $75\%$ training and $25\%$
test sets. The dataset contains 9 features as explanatory variables, whereas
the feature age is excluded from the original dataset to avoid potential
discrimination.

  * •

$x_{1}$: Number of times borrower has been 90 days or more past due (integer).

  * •

$x_{2}$: Number of times borrower has been 60-89 days past due but no worse in
the last 2 years (integer).

  * •

$x_{3}$: Number of times borrower has been 30-59 days past due but no worse in
the last 2 years (integer).

  * •

$x_{4}$: Total balance on credit cards and personal lines of credit except for
real estate and no installment debt such as car loans divided by the sum of
credit limits (percentage).

  * •

$x_{5}$: Monthly debt payments, alimony, and living costs divided by monthly
gross income (percentage).

  * •

$x_{6}$: Monthly income (real).

  * •

$x_{7}$: Number of open loans (installments such as car loan or mortgage) and
lines of credit (e.g., credit cards) (integer).

  * •

$x_{8}$: Number of mortgage and real estate loans including home equity lines
of credit (integer).

  * •

$x_{9}$: Number of dependents in the family, excluding themselves (spouse,
children, etc.) (integer).

  * •

$y$: Client’s behavior; 1 = Person experienced 90 days past due delinquency or
worse.

For simplicity, we consider the zero baseline as
$\mathbf{x}^{\prime}=\mathbf{0}$.

###  4.2. Model Performance

We apply black-box highly accurate fully-connected neural network (FCNN) and
transparent monotonic groves of neural additive model (MGNAM) (Chen and Ye,
2023). The FCNN is considered here due to its high degree of accuracy. The
MGNAM model was chosen because (i) it is transparent, and therefore easy to
analyze; (ii) monotonicity has been incorporated into the model. Accuracy is
measured by the area under the curve (AUC).

For FCNN, we use one hidden layer with ten neurons. For MGNAM, we consider
$g(\mathbb{E}[y|\mathbf{x}])=f(\mathbf{x})$ with the architecture

(12) $\displaystyle
f(\mathbf{x})=f_{1,2,3}(x_{1},x_{2},x_{3})+f_{4}(x_{4})+\dots+f_{9}(x_{9}).$

In other words, $x_{1}-x_{3}$ are grouped together, and the remaining features
are handled using 1-dimensional functions. As a convenience, we only impose
monotonicity on $x_{1}-x_{3}$. All functions are approximated by neural
networks with one hidden layer and two neurons. For $x_{1}-x_{3}$, we enforce
strong pairwise monotonicity. As a result of the domain knowledge, features
that are longer past due should have greater importance. In this paper, we
focus on simple architectures since there is no apparent improvement in
accuracy for more complicated models, and accuracy is not the primary concern.

Both models have an AUC of approximately $80\%$. This indicates that
transparent models are capable of modeling this dataset. We demonstrate the
application of IG and BShap to black-box FCNNs and then provide a more
detailed analysis of transparent MGNAMs regarding pairwise monotonicity.

###  4.3. FCNN Explanations

We demonstrate the application of IG and BShap to black-box FCNN. We
demonstrate how to interpret, debug, and avoid potential pitfalls with these
tools. An example of the application of IG and BShap is presented for

(13)
$\displaystyle\mathbf{x}=\left[\begin{matrix}3&1&1&1&0.11&6250&3&0&2\end{matrix}\right].$

Attributions for $\text{IG}(\mathbf{x},\mathbf{x}^{\prime},\text{FCNN})$ and
$\text{BShap}(\mathbf{x},\mathbf{x}^{\prime},\text{FCNN})$ are given in Figure
1 and Figure 2. In general, they produce similar results. The following is a
rough explanation of the result. Three times past due accounts with a lag time
of more than 90 days significantly increase the risk of default. Further, each
past due between 60-89 days and 30-59 days also increases the probability of
default, but at a smaller magnitude. Having a positive monthly income
decreases the probability and having two dependents increases it.

As a first step, we consider average monotonicity preservation. AIM is
preserved by both IG and BShap. In the case of pairwise monotonicity revealed
by IG, however, average attributions of $x_{1}$ are lower than those of
$x_{2}$. Based on Theorem 3.21, this indicates that strong pairwise
monotonicity is violated here. Therefore, this FCNN should not be used for
prediction in this scenario. Although BShap also observes a similar pattern,
it cannot guarantee a violation of strong pairwise monotonicity, as indicated
in Remark 3.22. In this regard, IG has an advantage over pairwise
monotonicity. If ASPM is violated in IG, then there must be a problem with the
model.

Figure 1. Instance explanation by IG for FCNN.  Figure 2. Instance explanation
by BShap for FCNN.

Next, we consider DIM preserving. We perturb $\mathbf{x}$ by increasing
$x_{3}$ by one unit,

(14)
$\displaystyle\mathbf{x}^{*}=\left[\begin{matrix}3&1&2&1&0.11&6250&3&0&2\end{matrix}\right].$

In principle, we would expect to observe an increase in the probability of
default. There is, however, a decrease of 0.015 in the probability of default.
As a result, the model does not adequately preserve individual monotonicity.
We then recompute IG and BShap in Figure 3 and Figure 4. IG and BShap continue
to provide positive attributions, suggesting that the model performs well on
average. Note that
$\text{IG}_{3}(\mathbf{x}^{*},\mathbf{x}^{\prime},\text{FCNN})$ decreases. In
spite of the fact that IG indicates that there may be a problem with the
model, we cannot conclude whether a problem exists, based on Remark 3.4.
Alternatively, $\text{BS}_{3}(\mathbf{x}^{*},\mathbf{x}^{\prime},\text{FCNN})$
increases, despite the decrease in the probability of default. BShap,
unfortunately, does not identify the model violation. As a result of these
results, although model explanations could help us debug, they cannot be
relied upon entirely, and it would be better to begin with the monotonic model
instead.

Figure 3. Instance explanation for the perturbed explicand (adding one
$x_{3}$) by IG for FCNN.  Figure 4. Instance explanation for the perturbed
explicand (adding one $x_{3}$) by BShap for FCNN.

###  4.4. MGNAM Explanations

We now check our results using MGNAM, where monotonicity is enforced.
Additionally, MGNAM is transparent, which allows us to examine explanations in
greater detail. Due to the fact that $x_{4}-x_{9}$ are modeled by one-
dimensional functions, attributions only represent values at the explicand, so
they are not very interesting. Because of this, we focus on the result for
$x_{1}-x_{3}$ and our focus is on $[0,2]$ for simplicity’s sake. Results by IG
and BShap for all cases are plotted in Figure 5 and Figure 6. Overall, they
perform similarly. All individual monotonicity is satisfied on average.
However, Bshap can be problematic when attempting to achieve pairwise
monotonicity. In the case of $(2,1,0)$,
$\frac{1}{2}\text{BS}_{1}<\text{BS}_{2}$, ASPM is violated by BShap.
Meanwhile, IG preserves ASPM. Furthermore, we would expect to see potential
more violations by BShap for larger values of $x_{1}-x_{3}$. Here, MGNAM has
already preserved monotonicity. IG provides appropriate explanations, whereas
BShap provides misleading explanations. Due to this, IG appears to be superior
in explaining pairwise monotonicity.

Figure 5. IG explanation for MGNAM for features $x_{1}-x_{3}$ in $[0,2]$.
X-axis are attributions and y-axis are features.  Figure 6. BShap explanation
for MGNAM for features $x_{1}-x_{3}$ in $[0,2]$. X-axis are attributions and
y-axis are features. At $(2,1,0)$, average strong pairwise monotonicity is
violated.

##  5\. Conclusion

We investigate attribution methods for monotonic models in this paper. New
axioms are proposed which must be satisfied by attribution methods. We analyze
IG and BShap methods. By analyzing and illustrating, we demonstrate that BShap
provides satisfactory explanations when only individual monotonicity is
involved, but IG provides reasonable explanations when strong pairwise
monotonicity is involved.

In the future, we would like to pursue several directions. (1) We have focused
on monotonicity in this paper, but it is necessary to investigate more domain
knowledge in order to gain a better understanding of the behavior of ML
models. (2) According to our analysis, neither IG nor BShap satisfies all
monotonic axioms. It is necessary to modify existing methods or propose new
methods in order to meet monotonicity requirements. (3) IG satisfies all
monotonicity axioms on average, however, it is necessary that the underlying
models are differentiable. Additionally, features such as the number of
defaults are naturally discrete. Although IG can still be applied in this case
using neural networks, it uses function values at points that are not useful
in practice. For this reason, it is desirable to have a modified Shapley value
method that is able to preserve strong pairwise monotonicity.

In general, we believe that domain knowledge, such as monotonicity, plays an
essential role in the interpretation of models. In sectors with high stakes,
such as the finance sector, we must take these factors into account when we
explain the models.

## References

  * (1)
  * OCC (2021) 2021\.  _Model risk management_.  https://www.occ.treas.gov/publications-and-resources/publications/comptrollers-handbook/files/model-risk-management/index-model-risk-management.html
  * Alam et al. (2022) Md Shafiul Alam, Jonathan Boardman, Xiao Huang, and Matthew Turner. 2022.  Applications of Integrated Gradients in Credit Risk Modeling. In _2022 IEEE International Conference on Big Data (Big Data)_. IEEE, 4997–5005.
  * Boardman et al. (2022) Jonathan Boardman, Md Shafiul Alam, Xiao Huang, and Ying Xie. 2022.  Integrated Gradients is a Nonlinear Generalization of the Industry Standard Approach to Variable Attribution for Credit Risk Models. In _2022 IEEE International Conference on Big Data (Big Data)_. IEEE, 5012–5023.
  * Bussmann et al. (2021) Niklas Bussmann, Paolo Giudici, Dimitri Marinelli, and Jochen Papenbrock. 2021.  Explainable machine learning in credit risk management.  _Computational Economics_ 57 (2021), 203–216.
  * Chen and Ye (2022) Dangxing Chen and Weicheng Ye. 2022.  Monotonic Neural Additive Models: Pursuing Regulated Machine Learning Models for Credit Scoring. In _Proceedings of the Third ACM International Conference on AI in Finance_. 70–78.
  * Chen and Ye (2023) Dangxing Chen and Weicheng Ye. 2023.  How to address monotonicity for model risk management?. In _Proceedings of the 40th International Conference on Machine Learning_ _(Proceedings of Machine Learning Research, Vol. 202)_. PMLR, 5282–5295.
  * Cotter et al. (2019) Andrew Cotter, Maya Gupta, Heinrich Jiang, Erez Louidor, James Muller, Tamann Narayan, Serena Wang, and Tao Zhu. 2019\.  Shape constraints for set functions. In _International conference on machine learning_. PMLR, 1388–1396.
  * Friedman and Moulin (1999) Eric Friedman and Herve Moulin. 1999.  Three methods to share joint costs or surplus.  _Journal of economic Theory_ 87, 2 (1999), 275–312.
  * Greydanus et al. (2019) Samuel Greydanus, Misko Dzamba, and Jason Yosinski. 2019\.  Hamiltonian neural networks.  _Advances in neural information processing systems_ 32 (2019).
  * Gupta et al. (2020) Maya Gupta, Erez Louidor, Oleksandr Mangylov, Nobu Morioka, Taman Narayan, and Sen Zhao. 2020\.  Multidimensional shape constraints. In _International Conference on Machine Learning_. PMLR, 3918–3928.
  * Horel and Giesecke (2020) Enguerrand Horel and Kay Giesecke. 2020.  Significance tests for neural networks.  _Journal of Machine Learning Research_ 21, 227 (2020), 1–29.
  * Horel et al. (2018) Enguerrand Horel, Virgile Mison, Tao Xiong, Kay Giesecke, and Lidia Mangu. 2018.  Sensitivity based Neural Networks Explanations.  _arXiv preprint arXiv:1812.01029_ (2018).
  * Karniadakis et al. (2021) George Em Karniadakis, Ioannis G Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang. 2021.  Physics-informed machine learning.  _Nature Reviews Physics_ 3, 6 (2021), 422–440.
  * Liu et al. (2020) Xingchao Liu, Xing Han, Na Zhang, and Qiang Liu. 2020\.  Certified monotonic neural networks.  _Advances in Neural Information Processing Systems_ 33 (2020), 15427–15438.
  * Lundberg and Lee (2017) Scott M Lundberg and Su-In Lee. 2017.  A unified approach to interpreting model predictions.  _Advances in neural information processing systems_ 30 (2017).
  * Lundstrom et al. (2022) Daniel D Lundstrom, Tianjian Huang, and Meisam Razaviyayn. 2022\.  A rigorous study of integrated gradients method and extensions to internal neuron attributions. In _International Conference on Machine Learning_. PMLR, 14485–14508.
  * Milani Fard et al. (2016) Mahdi Milani Fard, Kevin Canini, Andrew Cotter, Jan Pfeifer, and Maya Gupta. 2016.  Fast and flexible monotonic functions with ensembles of lattices.  _Advances in neural information processing systems_ 29 (2016).
  * Repetto (2022) Marco Repetto. 2022\.  Multicriteria interpretability driven deep learning.  _Annals of Operations Research_ (2022), 1–15.
  * Ribeiro et al. (2016) Marco Tulio Ribeiro, Sameer Singh, and Carlos Guestrin. 2016\.  ” Why should i trust you?” Explaining the predictions of any classifier. In _Proceedings of the 22nd ACM SIGKDD international conference on knowledge discovery and data mining_. 1135–1144.
  * Ribeiro et al. (2018) Marco Tulio Ribeiro, Sameer Singh, and Carlos Guestrin. 2018\.  Anchors: High-precision model-agnostic explanations. In _Proceedings of the AAAI conference on artificial intelligence_ , Vol. 32.
  * Sundararajan and Najmi (2020) Mukund Sundararajan and Amir Najmi. 2020.  The many Shapley values for model explanation. In _International conference on machine learning_. PMLR, 9269–9278.
  * Sundararajan et al. (2017) Mukund Sundararajan, Ankur Taly, and Qiqi Yan. 2017.  Axiomatic attribution for deep networks. In _International conference on machine learning_. PMLR, 3319–3328.
  * You et al. (2017) Seungil You, David Ding, Kevin Canini, Jan Pfeifer, and Maya Gupta. 2017.  Deep lattice networks and partial monotonic functions.  _Advances in neural information processing systems_ 30 (2017).


[2]The Grantham Institute, Imperial College London,
London, SW7 2AZ, UK
[3]Department of Computing, Imperial College London,
London, SW7 2AZ, UK
[4]Department of Mathematics, Imperial College London,
London, SW7 2AZ, UK

One approach to achieving correct finite element assembly is to ensure
that the local orientation of facets relative to each cell in the mesh
is consistent with the global orientation of that facet.
Rognes et al. have shown how to achieve this for any mesh
composed of simplex elements, and deal.II contains
a serial algorithm to construct a consistent orientation of any
quadrilateral mesh of an orientable manifold.

The core contribution of this paper is the extension of this algorithm
for distributed memory parallel computers, which facilitates its seamless
application as part of a parallel simulation system.

Furthermore, our analysis establishes a link between the
well-known Union-Find algorithm and the construction of a
consistent orientation of a quadrilateral mesh.
As a result, existing work on the parallelisation of the Union-Find
algorithm can be easily adapted to construct further parallel algorithms
for mesh orientations.

finite element assembly;
parallel algorithm;
quadrilateral mesh;

68W10; 68W15

M. HOMOLYA AND D. A. HAMPARALLEL QUADRILATERAL EDGE ORIENTATIONS

§ INTRODUCTION

A characteristic of the finite element method for the solution of partial
differential equations is that the representation of functions over the
domain can be chosen from a wide range of function spaces.  This choice is
achieved by selecting a particular local function space to be used for each
cell in the meshed domain. Information is communicated between adjacent
cells either by shared degrees of freedom on the cell boundaries, or by flux
integrals over the facets[By facet we mean a mesh entity of
codimension 1. In two dimensions these are the cell edges, in three
dimensions, the cell faces.] between cells. For every case except the
lowest degree continuous finite elements, these coupling terms require the
orientations of the two cells adjacent to each facet to be reconciled with
the orientation of that facet. Failure to do so will result in false
identification of degrees of freedom falling on the facets, or incorrect
accounting for the flux direction through facets.

Consequently, every finite element implementation which supports facet integrals, elements
of polynomial degree greater than one, or  or  conforming
elements must somehow ensure that adjacent cells agree on the orientation
of the intervening facet.  This can either be achieved by explicitly
recording orientation information and exploiting this in the local and/or
global assembly operations, or it can be achieved by ensuring, in a sense
which we will later make formal, that the local orientation of facets
relative to each cell in the mesh is consistent with the global orientation
of that facet.

The consistent numbering approach has the advantage that the integral
assembly code is simpler than that required by the explicit orientation
approach, at the cost of somehow establishing the consistent numbering and
with the limitation that a consistent numbering may not exist for all
meshes. Consistent numberings are simple to construct for any mesh composed
of simplex elements [21] but the same algorithm does not extend
to quadrilaterals.  Instead, the deal.II finite element library
[5] contains an algorithm to construct a consistent
orientation of any quadrilateral mesh of an orientable
manifold[A preprint documenting the deal.II algorithm has
recently appeared [1].].  The
algorithm employed in deal.II is serial. This is less than ideal on modern
supercomputers, for which it would be both more convenient and more
efficient to construct the orientation using a distributed parallel
algorithm. The core contribution of this paper is therefore to present a
distributed parallel extension of the algorithm in deal.II which creates a
global consistent numbering.  This algorithm is presented in
<ref>, with experimental evaluation in <ref>.
The formal definition of the problem is given in <ref>
and the proof of the existence of a solution in <ref>.
<ref> restates the deal.II serial algorithm in terms of
this proof.  The parallel part of this
algorithm is, in fact, an adaption of the well-known Union-Find
algorithm, so in <ref> we place this algorithm in the context
of that existing work.  We conclude the paper in <ref>.

Facet orientation has been addressed in different ways by various finite
element packages. Rognes et al. [21] give a more detailed
discussion of the features which make consistent orientations useful for
the efficient assembly of finite elements. They also provide an algorithm to
find consistent facet orientations for simplicial meshes, which is used in
FEniCS [13] and in Firedrake [20].  Other
finite element packages often use different approaches.  libMesh
[16] uses global vertex numbers to determine the direction of
edges for  conforming elements, however neither  conforming,
nor higher-order elements are supported.  FreeFem++ [12] does
not support quadrilaterals.  Nektar++ [7] stores and uses
explicit per-cell edge orientations, and also supports mixed cell meshes
(e.g. mixed triangles and quadrilaterals).  These orientations are
employed to match degrees of freedom on facets for higher-order elements,
however no sign flipping is necessary since  and  conforming
elements are not supported.  Raviart-Thomas elements in DUNE
[6] store explicit signs ($-1$ or $+1$) for each facet, which
are referenced from manually written formulæ.  deal.II [5]
employs the serial algorithm we describe in <ref>, as a
post-processing step after reading the mesh.

§ PROBLEM STATEMENT

This section aims to formally define the problem that we later
analyse and provide algorithmic solutions to.
Throughout this paper, we say that a mesh has consistent facet
orientations, if:

* Facet orientation is intrinsic to the facet, i.e. it
does not depend on the cell from which the facet is accessed.

* There is a fixed reference cell such that for each cell in the mesh,
there exists a mapping to the reference cell under which the orientation
of each facet incident to that cell matches the orientation of the
corresponding facet of the reference cell.

The orientation of facets on the reference cell can be arbitrary for
the purpose of the above definition,
although this choice can limit the scope of meshes for which a
consistent facet orientation exists.

In this paper, we tackle the problem of finding consistent
edge orientations for quadrilateral meshes on 2-dimensional
orientable manifolds.

§ PROOF OF THE EXISTENCE OF A SOLUTION

We prove that consistent edge orientations exist for any
quadrilateral mesh on an orientable 2-dimensional manifold.
Later we build on the analysis below to construct algorithms that
find consistent orientations.

As mentioned earlier, the choice of edge orientations on the
reference quadrilateral can affect the set of meshes for which
consistent orientations exist.
Up to  symmetry, there are four different edge orientations of a
quadrilateral. These are shown in Figure <ref>. The
symmetries arise since
we can arbitrarily choose the first vertex and the direction of
traversal when defining a mapping to the reference cell.
For each case, one can find a mesh on a 2-dimensional manifold which
does not have consistent edge orientations. [4] includes a
counter-example for (a). A structured grid with odd number of cells in
both direction, folded into a torus, is a counter-example for (b) and
(d). A structured grid bent into a Möbius strip is a counter-example
for (c). However, when we restrict our attention to meshes on
orientable manifolds (that is, the vast majority of domains actually
employed for finite element problems), then only (c) remains without a
counter-example. We here prove existence for this case:

Possible different quadrilateral edge orientations, after
considering symmetries.

Consistent edge orientations exist for any quadrilateral mesh on an
orientable 2-dimensional manifold, given the edge orientations on the
reference quadrilateral as in Figure <ref>.

Notice that on the reference quadrilateral, opposite edges always have the
same orientation. When constructing a global orientation of the cells in a
mesh, this implies that fixing the orientation of one edge
determines the orientation of the opposite edge. Conversely, the
orientation of an edge imposes no restriction on the orientation of the
two adjacent edges. Since each interior edge of the mesh is adjacent to
two cells, fixing the orientation of one edge actually fixes the
orientation of every edge reachable from the first edge by an unbroken
sequence of opposite cell edges. This defines a relation “the orientation
of edge $a$ determines the orientation of edge $b$”.  We will refer to
this as the orientation determination relation. It is easy to see
that this relation is reflexive, symmetric and transitive, and is
therefore an equivalence relation which defines equivalence classes on
the set of edges. The form of each equivalence class is a path or ribbon
(see Figure <ref>)
through the mesh connecting successive opposite faces. Within each class,
the orientation of any edge determines the orientation of all other
edges. Conversely, since the orientation of an edge only determines the
orientation of the opposite edge and not that of the adjacent edges,
separate classes can be oriented independently.

Unstructured quadrilateral mesh with consistent edge
orientations. Two “ribbons” (sets of edges where the orientation
of any edge determines the orientation of all edges) are
highlighted in dark red and light blue.

Note that if there exists a consistent edge orientation for a quadrilateral
mesh, one can invert the orientation of every edge in any equivalence class,
and still have a consistent edge orientation. Thus when assigning
orientations to a set of dependent edges, we can arbitrarily choose the
orientation of one edge, and that orientation propagates to all dependent
edges. Either both initial choices are fine, or a consistent orientation
does not exist.

This implies that a consistent orientation will always exist for a mesh,
unless there is a set of dependent edges for which the orientations necessarily
conflict. That is, there is an edge $u$, whose orientation forces a
different orientation on an edge $v$ through path $ p_1 $ than through path
$ p_2 $. From a slightly different point of view, there is a cycle $ p = u
\overset{p_1}{\leadsto} v \overset{p'_2}{\leadsto} u $ around $u$, where $
p'_2 $ is the reverse of $ p_2 $, such that the orientation of $u$ forces a
conflicting orientation on itself through path $p$. Without loss of
generality we can assume that $p$ is a simple cycle, i.e. $p$ contains
each of its edges only once (except for $u$, which appears as starting and
ending edge in $p$). In this case the quadrilateral cells in $p$ form a
Möbius strip (see Figure <ref>), which implies that the
mesh is not on an orientable manifold.

Edge orientations propagate in a conflicting manner on
a Möbius strip, and vice versa, if a cycle propagates
orientations in a conflicting way, then connecting the cells will
form a Möbius strip.

We have shown
that if the mesh does not admit consistent edge orientations, then its
manifold is non-orientable. It follows immediately that consistent orientations
exist for any quadrilateral mesh on an orientable manifold.

§ SERIAL ALGORITHM

Based on the above discussion, Alg. <ref> describes
the steps of consistently orienting a quadrilateral mesh. We assume
that initially $ \mathtt{visited}[e] = \mathtt{False} $ and
$ \mathtt{orientations}[e] $
is undefined for all $ e \in E $, where $ E $ denotes the set
of edges in the mesh. For generality, lines 13 and 14 lack the details
of mesh representation and orientation representation.

Serial algorithm

[1]
$ e \in E $
Orient$ e $, $ \uparrow $
Orient$ e $, $ o $
$ \mathtt{visited}[e] $
$ \mathtt{orientations}[e] \neq o $
Möbius strip found
$ \mathtt{visited}[e] \gets \mathtt{True} $
$ \mathtt{orientations}[e] \gets o $
$ c \in \text{cells incident to } e $
$ e' \gets \text{edge opposite to } e \text{ in } c $
$ o' \gets \text{orientation for } e' \text{ consistent with } e \text{ having orientation } o $
Orient$ e' $, $ o' $

Lines 10–16 run only $ |E| $ times, since each time they flip
a  to  in $ \mathtt{visited} $. Since the
quadrilateral mesh is on a 2-dimensional manifold, only 1 or 2 cells
can be incident to an edge (line 12), therefore lines 10–16 take
$ O(|E|) $ total runtime, not considering the time spent in recursive
Orient calls outside these lines. It also follows that
Orient is called $ O(|E|) $ times: $ |E| $ times in line 2,
and between $ |E| $ and $ 2 |E| $ times in line 15. Hence the serial
algorithm takes $ O(|E|) $ time.

The serial algorithm was first implemented in deal.II
[5], and briefly described in its
documentation [4].

§ PARALLEL ALGORITHM

This section describes the parallel extension of the above
algorithm for distributed memory systems, as it is
implemented in Firedrake.

The basic idea of Firedrake's parallel implementation is to first
solve edge orientation locally, then flip some local ribbon segments
until all local domains agree on the orientation of all shared edges
(see Figure <ref>).
To speed up flipping local ribbon segments, each parallel process
records connections between its shared edges. Interior edges are
re-aligned only when agreement on the orientations of all shared
edges is reached.

A “ribbon” in the global mesh, and its segments in the
local domains. It is assumed that the orientation of edges is
consistent inside these segments. However, the local domains
may (green, straight connection) or may not
(red, twisted connection) agree on the
orientation of shared edges.

Alg. <ref> describes how parallel processes negotiate the
orientation of shared edges. Before the negotiation starts, each
process aligns its edges locally, while populating the
and  arrays.  If $ u $
and $ v $ are edges shared with other processes, and connected by a
local ribbon segment, then  is populated with
entries $ \mathtt{affects\_edge}[u] = v $ and $
\mathtt{affects\_edge}[v] = u $, and  is
populated with entries $ \mathtt{affects\_orient}[u] =
\mathtt{affects\_orient}[v] ={} \uparrow $ when $ u $ and $ v $ have
the same orientation, and with entries $ \mathtt{affects\_orient}[u] =
\mathtt{affects\_orient}[v] ={} \downarrow $ otherwise. If $ u $ is a
shared edge, but it does not connect to any other shared edge, these
arrays are not updated. If we encounter a closed loop while locally
orienting edges, there is no need to negotiate the orientation its
edges with other parallel processes.

Another prerequisite is to initialise the  and
arrays with the weights and proposed orientations
of the shared edges of a parallel process. The weights shall be
globally unique for all local ribbon segments. We use $
\mathtt{our\_weight}[e] = \mathit{size} \cdot l(e) + \mathit{rank} $,
where $ l(e) $ is the local index of the ribbon segment which $ e $
belongs to, size is the number of parallel process, and
rank is the index of the current process. We need to ensure
that locally connected edges have the same weight and consistent
orientation initially. This invariant is maintained by
Alg. <ref>.

Parallel algorithm

[1]
$ \mathtt{conflict} \gets \mathtt{True} $

$ \mathtt{conflict} \gets \mathtt{False} $
$ e \in \text{shared edges} $
$ \mathtt{our\_weight}[e] = \mathtt{their\_weight}[e]$
$ \mathtt{our\_orient}[e] \neq \mathtt{their\_orient}[e] $
Möbius strip found
$ \mathtt{our\_weight}[e] < \mathtt{their\_weight}[e]$
$ \mathtt{our\_weight}[e] = \mathtt{their\_weight}[e] $
$ \mathtt{our\_orient}[e] = \mathtt{their\_orient}[e] $
$ e \in \mathtt{affects\_edge} $
$ e' \gets \mathtt{affects\_edge}[e] $
$ o' \gets \mathtt{affects\_orient}[e] \text{ XOR } \mathtt{our\_orient}[e] $
$ \mathtt{our\_orient}[e'] \neq o' $
$ \mathtt{conflict} \gets \mathtt{True} $
$ \mathtt{our\_weight}[e'] = \mathtt{our\_weight}[e] $
$ \mathtt{our\_orient}[e'] = o' $
$ \gets $ AllReduce, Or

Alg. <ref> consists of a main loop which terminates when there
are no more conflicts in the proposed orientations for shared edges. In
line 3–4 each parallel process exchanges its proposed orientations
and weights for its shared edges with its
neighbours.  and  have
type and size identical to  and
respectively, but they contain the values
proposed by the
neighbouring processes. The general rule is to adopt the orientation with
the larger weight.  Lines 11–23 handle the adoption of the remote weight and
orientation. Line 14 checks if such an update needs propagation via the
local ribbon segment. Lines 15–21 propagate the new information to the
other connected shared edge, and check if it causes that edge to flip.
If yes, then
another round of exchange is necessary, since the orientation for $ e' $
that this process forwarded in line 3, has now become outdated. When the
local weight is larger than the remote weight, there is nothing to do: the
other parallel process will adopt our orientation.  Since all weights are
initially unique, and then they are copied between neighbours, the local
weight being equal to the remote weight means that both parallel processes
have their local ribbon segment aligned to the same “source”. Therefore a
mismatch in the orientation is only possible if the global mesh contains a
Möbius strip, which is checked in lines 8–10.

This algorithm terminates in at most $ k $ communication rounds, where
$ k $ is the maximum number of local segments in any ribbon of the
global mesh. Since each ribbon is “untwisted” independently, further
discussion focuses on one arbitrary ribbon in the global mesh.
Any time there is a conflict between the orientation of two local
segments, the one with the higher weight “wins”, thus we can
conclude that the orientation of a ribbon in the global mesh is
ultimately determined by the local segment which has the highest
weight. Initially only one segment has this highest weight. In each
communication round one or two additional segments adopt this weight
and align their orientations, until the information has propagated
to all segments. At that point, there are either no conflicts along
the ribbon, or the algorithm will abort in the next communication
round due to finding a Möbius strip. Therefore the parallel
algorithm terminates in at most $ k $ communication rounds.

Theoretically, $ k $ can be large with respect to other parameters of the
mesh. For example, one can mesh an annulus domain with a very long,
spiraling ribbon. However, a far more typical scenario is for each ribbon to
either enter one side of the domain and exit the other, or to form a
loop. In either of these cases, assuming the domain has bounded aspect ratio
and that the parallel decomposition approximately minimises the surface area
of partitions, $ k $ will be $ O(\sqrt{P}) $, where $ P $ is the number of
parallel processes.

§ EXPERIMENTS

We have run several scaling experiments on the ARCHER UK National
Supercomputing Service. Our experimentation framework is available as
[14]. We used 6 different meshes to
evaluate the performance of the parallel algorithm described in

: structured grid on a square domain

: cubed sphere mesh (a quadrilateral
mesh of the surface of a sphere formed by deforming a cube to sphere shape
and refining)

: unstructured mesh on a square domain

: unstructured mesh on the surface of a sphere

, : unstructured meshes with
non-uniform resolution (see Figure <ref>)

The first two were generated using Firedrake utility functions, the other
four were generated with Gmsh [10], version
2.8.3.  and  are examples from the Gmsh tutorial. To
give an impression of the level of irregularity of these meshes,
low-resolution versions are presented in Figure <ref>.
We used meshes with cell count in the order of millions, unfortunately mesh
generation becomes very difficult for meshes bigger than that.

Low-resolution version of mesh samples  (left) and  (right).

We used Firedrake to construct consistent edge orientations for each
of these meshes on up to 24576 cores (1024 computing nodes, each with
24 CPU cores). We measured the number of required communication rounds
(Figure <ref>) and execution time
(Figure <ref>). Each experiment was repeated 4 times,
except for those involving the largest number of cores: for cost
efficiency reasons, they were carried out only once. The number of
communication rounds were the same during each instantiation of the
same experiment, however, the execution times varied greatly,
especially when they were under a second. Therefore reported execution
times should be taken with a grain of salt. Figure <ref>
shows the average of 4 measurements for each data point.

legend cell align=left,
legend pos=outer north east,
xlabel=# of parallel processes,
ylabel=# of communication rounds,
table [x=P, y=rounds, col sep=comma] rounds_s_square.csv;

table [x=P, y=rounds, col sep=comma] rounds_s_sphere.csv;

table [x=P, y=rounds, col sep=comma] rounds_u_square.csv;

table [x=P, y=rounds, col sep=comma] rounds_u_sphere.csv;

table [x=P, y=rounds, col sep=comma] rounds_t10.csv;

table [x=P, y=rounds, col sep=comma] rounds_t11.csv;

[black] table [x=P, y=rounds]
P     rounds
24    60
24576 1920
node[pos=0.4,pin=355:$\sqrt{P}$] ;

Number of required communication rounds for each mesh sample.

legend cell align=left,
legend pos=outer north east,
xlabel=# of parallel processes,
ylabel=Execution time [s],
table [x=P, y=time, col sep=comma] time_s_square.csv;

table [x=P, y=time, col sep=comma] time_s_sphere.csv;

table [x=P, y=time, col sep=comma] time_u_square.csv;

table [x=P, y=time, col sep=comma] time_u_sphere.csv;

table [x=P, y=time, col sep=comma] time_t10.csv;

table [x=P, y=time, col sep=comma] time_t11.csv;

Execution time to construct consistent edge orientations in
parallel for each mesh sample.

Unsurprisingly, the construction edge orientations for structured
meshes takes relatively few communication rounds.
Figure <ref> also shows that for a rectangular domain,
the number of required communication rounds stays the same order of
magnitude in case of non-uniform meshing as with uniform resolution.
However, the unstructured mesh on a spherical domain requires at least
one order of
magnitude more communication rounds. This is probably caused by the
fact that a sphere does not have a boundary, thus all ribbons must
form a closed loop, so traversing a ribbon must continue until the
starting edge is reached again. Regardless, our measurements confirm
our predicted $ O(\sqrt{P}) $ communication rounds for strong
scaling. Table <ref> collects the slopes of lines fitted to
the log-log data shown in Figure <ref>, all values
being close to the expected $ 0.5 $.

[tabular=l | l, table head=Test case     Slope
Slopes of lines fitted to the log-log data shown in Figure <ref>.

In most of the cases shown in Figure <ref>, the execution time
initially decreases as the number of parallel processes increases. It is
likely that this is due to the parallelisation of the local work each
process has to do. Eventually the increasing number of communication
rounds becomes dominant, and the execution times grow with the number
parallel processes. However, even the unstructured sphere mesh with
almost one million cells, running on about a quarter of ARCHER, takes only
105 seconds. For users running supercomputer-scale simulations, the
few minutes spent on edge orientations are not likely to become a

§ ALTERNATIVE APPROACHES TO PARALLELISM

The edge orientation problem can be reduced to the well-known
Union-Find algorithm, also known as the disjoint set data
structure [8].  We established in <ref>
that the orientation determination relation is an equivalence relation
which defines equivalence classes (the “ribbons”) on the set of edges.
To discover the ribbons, we can immediately use the Union-Find
algorithm. However, with an extension to it, we can also retrieve the
correct edge orientations.

First we briefly introduce the Union-Find algorithm. Let $ e_1,
e_2, \ldots, e_n $ be abstract elements, and $ S_1, S_2, \ldots, S_k $
disjoint non-empty sets over these elements, with $k \le n$. The sets are
initially singletons: $ S_i = \{ e_i \}, i = 1\ldots{}n $.
The Union operation destructively joins two sets.
The Find operation can be used to determine whether two elements,
$ e_i $ and $ e_j $, are in the same set.
Since the sets are disjoint, each element is in exactly one set.
Each set is identified by its representative element, which can
be any of its elements.
The Find operation returns, for any element $ e_i $, the
representative element of the set to which $ e_i $ belongs.
Thus $ \textsc{Find}(e_i) = \textsc{Find}(e_j) $ if $ e_i $
and $ e_j $ are in the same set.

To solve the edge orientation problem, the edges of the mesh become the
abstract elements for the Union-Find algorithm, starting with a singleton
set for each edge. Then, for each cell, we join both pairs of opposite
edges, more precisely we join the sets to which those opposite edges
belong. When done, each set corresponds to a “ribbon”, a set of edges
which determine each other's orientation.  However, we still need to
traverse the ribbon to determine the permitted relative orientations between
edges.  We will avoid this traversal by extending the Union-Find algorithm,
but let us first briefly present the usual representation of sets and the
algorithms for Find and Union.

Each element $ u $ links to its parent element $ p(u) $.
For representative elements, $ p(u) = u $.
For all other elements, $ p(u) $ is another element in the same set,
and following the parent links must eventually lead to the
representative element of the set.
Alg. <ref> and Alg. <ref> show the implementation of the
Find and Union operations, respectively.

Find operation for disjoint sets

[1]
Find$ u $
$ p(u) = u $
$ u $
Find$ p(u) $

Union operation for disjoint sets

[1]
Union$ u $, $ v $
$ r_u \gets $ Find$ u $
$ r_v \gets $ Find$ v $
$ r_u \neq r_v $
$ r_v $ becomes the representative element of the united set
$ p(r_u) \gets r_v $

Find simply traverses the parent links until it reaches the
representative element.  A variation of Find also reassigns the parent
links of each element on this path directly to the representative element –
this technique is called path compression.
Union first calls Find to look up the representative elements
of both sets. If the two sets are indeed different, then it assigns the
parent link of one of the representative elements to point to the other
representative element, thus effectively merging the sets.
Depending on the exact variant of the Union-Find algorithm [18],
there is generally a rule to determine which one of $ r_u $ and $ r_v $ remains
a representative element, but we ignore that detail for this discussion.
The textbook version of the Union-Find algorithm [22]
uses rank based merging in Union and path compression in
Find, and it is proven that $ m $ Find and $ n $
Union operations execute in $ O( (m + n) \alpha(m, n) ) $ time, where
$ \alpha(m, n) $, the functional inverse of Ackermann's function, is so
slow-growing that for all practical purposes it can be considered constant.

To address the edge orientation problem, we attach a binary value
$ \omega(u) $ to each parent link, which describes the
permitted relative orientation between the edge $ u $ and its parent
edge $ p(u) $.  Alg. <ref> and
Alg. <ref> describe Find and Union
with orientations, respectively.
The changes include the maintenance of orientation, and Möbius strip
detection in Union when one tries to connect a ribbon to
itself with the inconsistent orientation.

Find operation for disjoint sets with orientation

[1]
Find$ u $
$ p(u) = u $
$ u $, $ \uparrow $
$ r_u, o_p \gets $ Find$ p(u) $
$ r_u, \, \omega(u) \text{ XOR } o_p $

Union operation for disjoint sets with orientation

[1]
Union$ u $, $ v $, $ o $
$ r_u, o_u \gets $ Find$ u $
$ r_v, o_v \gets $ Find$ v $

$ r_u = r_v $
$ o_u \text{ XOR } o \text{ XOR } o_v ={} \downarrow $
Möbius strip found
$ r_v $ becomes the representative element of the united set
$ p(r_u) \gets (r_v, \, o_u \text{ XOR } o \text{ XOR } o_v ) $

This extension can be easily applied to most of the relevant variants
of the Union-Find algorithm, including parallel ones.
Alg. <ref> describes the algorithm to solve the edge
orientation problem using the Union-Find algorithm extended with
orientations. While Alg. <ref> itself does not feature any
parallel primitives, it will solve the edge orientation problem in
parallel if a parallel variant of the Union-Find algorithm
is employed. In that case each parallel process only iterates over its
local cells and edges.

Edge orientations using the extended Union-Find algorithm

[1]
$ c \in \text{\(local\) cells} $
$ e_1, e_2, e_1', e_2' \gets \text{edges of } c \text{, each consecutive pair sharing a vertex} $
$ o_1 \gets \text{permitted relative orientation between } e_1 \text{ and } e_1' $
$ o_2 \gets \text{permitted relative orientation between } e_2 \text{ and } e_2' $
Union$ e_1 $, $ e_1' $, $ o_1 $
Union$ e_2 $, $ e_2' $, $ o_2 $
$ e \in \text{\(local\) edges} $
$ r_e, o_e \gets $ Find$ e $
$ \mathtt{orientations}[e] \gets o_e $

The advantage of this approach is that we can easily re-use existing
work on parallelising the Union-Find algorithm to solve the edge
orientation problem.
Several attempts have been made to parallelise that algorithm both
on shared memory computers [9, 2, 3, 19], and on distributed memory systems [9, 17, 11, 15].
Since Firedrake is meant to scale on today's supercomputers, we are
only interested in distributed memory parallelism.
Cybenko et al. [9] attempt both shared memory and
distributed memory parallelisation, but their distributed memory
algorithm suffered from slowdown with increasing number of processors
and constant problem size.
Manne and Patwary [17] propose a distributed memory
parallel algorithm, which they demonstrate to scale up to 40
processors for certain graphs. This is, however, far from the scale
of supercomputers on which Firedrake is supposed to run.
Harrison et al. [11] achieve good scaling for their
purposes, which is connected component labeling for distributed
memory visualisation tools. However, the final step of their algorithm
uses all-to-all communication to merge the local components, which we
expect to become a serious bottleneck for our purposes, since
distributed quadrilateral meshes produce numerous local ribbon
segments and connections between them.
Iverson et al. [15] compare five parallel algorithms
for connected component labeling on distributed memory systems, both
theoretically and experimentally. The approach they call label
propagation is essentially equivalent to the parallelisation of
the edge orientation algorithm implemented in Firedrake.
Although it was not always the fastest in their evaluations,
it did not suffer from the explosive use of memory which affected some
other algorithms.

§ CONCLUSION & OUTLOOK

We have proposed a distributed mesh parallel algorithm, extending the
serial quadrilateral edge orientation algorithm in deal.II.
We have proven that there are consistent edge orientations for any
quadrilateral mesh of an orientable 2-dimensional manifold.
Our novel analysis establishes a link between the well-known
Union-Find algorithm and assigning orientations to the edges of
a quadrilateral mesh, thus existing parallelisation approaches for
the Union-Find algorithm can be easily adapted to construct further
parallel algorithms for mesh orientations.

[1]
Rainer Agelek, Michael Anderson, Wolfgang Bangerth, and William Barth,
On orienting edges of unstructured two- and three-dimensional meshes,
arXiv,  (2015).

[2]
Richard J. Anderson and Heather Woll, Wait-free parallel algorithms
for the union-find problem, in Proceedings of the Twenty-third Annual ACM
Symposium on Theory of computing - STOC '91, 1991, pp. 370–380.

[3]
David A. Bader and Guojing Cong, A fast, parallel spanning tree
algorithm for symmetric multiprocessors (SMPs), Journal of Parallel and
Distributed Computing, 65 (2005), pp. 994–1006.

[4]
Wolfgang Bangerth, Michael Anderson, and Ralf Hartmann, The deal.II
Library: GridReordering$<$dim, spacedim$>$ Class Template Reference.
Accessed: 2015-04-08.

[5]
W. Bangerth, R. Hartmann, and G. Kanschat, deal.II – A
general-purpose object-oriented finite element library, ACM Transactions on
Mathematical Software, 33 (2007), pp. 24/1–24/7.

[6]
Peter Bastian, Markus Blatt, Christian Engwer, Andreas Dedner, Robert
Klöfkorn, Sreejith P. Kuttanikkad, Mario Ohlberger, and Oliver Sander,
The Distributed and Unified Numerics Environment (DUNE), in
Proceedings of the 19th Symposium on Simulation Technique in Hannover, 2006.

[7]
C. D. Cantwell, D. Moxey, A. Comerford, A. Bolis, G. Rocco, G. Mengaldo,
D. de Grazia, S. Yakovlev, J.-E. Lombard, D. Ekelschot, B. Jordi, H. Xu,
Y. Mohamied, C. Eskilsson, B. Nelson, P. Vos, C. Biotto, R. M. Kirby, and
S. J. Sherwin, Nektar++: An open-source spectral/hp element
framework, Computer Physics Communications,  (2015).

[8]
Thomas H Cormen, Charles E Leiserson, Ronald L Rivest, Clifford Stein, and
Others, Introduction to algorithms, vol. 2, MIT Press Cambridge,

[9]
G. Cybenko, T. G. Allen, and J. E. Polito, Practical parallel
Union-Find algorithms for transitive closure and clustering, International
Journal of Parallel Programming, 17 (1988), pp. 403–423.

[10]
Christophe Geuzaine and Jean-François Remacle, Gmsh: a
three-dimensional finite element mesh generator with built-in pre- and
post-processing facilities, International Journal for Numerical Methods in
Engineering, 79 (2009), pp. 1309–1331.

[11]
Cyrus Harrison, Hank Childs, and Kelly P. Gaither, Data-parallel
mesh connected components labeling and analysis, in EG PGV'11 Proceedings
of the 11th Eurographics conference on Parallel Graphics and Visualization,
2011, pp. 131–140.

[12]
F Hecht, New development in FreeFem++, Journal of Numerical
Mathematics, 20 (2012), pp. 251–265.

[13]
J. Hoffman, J. Jansson, C. Johnson, M. G. Knepley, R. C. Kirby, A. Logg,
L. R. Scott, and G. N. Wells, FEniCS.
<http://fenicsproject.org/>, 2006.
Accessed: 2015-04-08.

[14]
Miklós Homolya, Experimentation framework for manuscript “A
parallel edge orientation algorithm for quadrilateral meshes”, version 2.
<http://dx.doi.org/10.5281/zenodo.31190>, Sept. 2015.

[15]
J. Iverson, C. Kamath, and G. Karypis, Evaluation of
connected-component labeling algorithms for distributed-memory systems,
Parallel Computing, 44 (2015), pp. 53 – 68.

[16]
Benjamin S. Kirk, John W. Peterson, Roy H. Stogner, and Graham F. Carey,
libMesh: a C++ library for parallel adaptive mesh refinement/coarsening
simulations, Engineering with Computers, 22 (2006), pp. 237–254.

[17]
Fredrik Manne and Md. Mostofa Ali Patwary, A scalable parallel
union-find algorithm for distributed memory computers, Lecture Notes in
Computer Science (including subseries Lecture Notes in Artificial
Intelligence and Lecture Notes in Bioinformatics), 6067 LNCS (2010),
pp. 186–195.

[18]
Md. Mostofa Ali Patwary, Jean Blair, and Fredrik Manne, Experiments
on union-find algorithms for the disjoint-set data structure, Lecture Notes
in Computer Science (including subseries Lecture Notes in Artificial
Intelligence and Lecture Notes in Bioinformatics), 6049 LNCS (2010),
pp. 411–423.

[19]
Md. Mostofa Ali Patwary, Peder Refsnes, and Fredrik Manne, Multi-core spanning forest algorithms using the disjoint-set data
structure, Proceedings of the 2012 IEEE 26th International Parallel and
Distributed Processing Symposium, IPDPS 2012,  (2012), pp. 827–835.

[20]
Florian Rathgeber, David A Ham, Lawrence Mitchell, Michael Lange, Fabio
Luporini, Andrew T T McRae, Gheorghe-Teodor Bercea, Graham R Markall, and
Paul H J Kelly, Firedrake: automating the finite element method by
composing abstractions, Submitted to ACM TOMS,  (2015).

[21]
Marie Rognes, Robert C Kirby, and Anders Logg, Efficient Assembly
of H(div) and H(curl) Conforming Finite Elements, SIAM Journal on
Scientific Computing, 31 (2012), pp. 4130–4151.

[22]
Robert Endre Tarjan, Efficiency of a good but not linear set union
algorithm, Journal of the ACM (JACM), 22 (1975), pp. 215–225.