# Globally Interpretable Classifiers via Boolean Formulas with Dynamic
Propositions

Reijo Jaakkola  Tampere University  Tomi Janhunen  Tampere University  Antti
Kuusisto  Tampere University  Masood Feyzbakhsh Rankooh  Tampere University
Miikka Vilander111The authors are given in the alphabetical order. Tampere
University

###### Abstract

Interpretability and explainability are among the most important challenges of
modern artificial intelligence, being mentioned even in various legislative
sources. In this article, we develop a method for extracting immediately human
interpretable classifiers from tabular data. The classifiers are given in the
form of short Boolean formulas built with propositions that can either be
directly extracted from categorical attributes or dynamically computed from
numeric ones. Our method is implemented using Answer Set Programming. We
investigate seven datasets and compare our results to ones obtainable by
state-of-the-art classifiers for tabular data, namely, XGBoost and random
forests. Over all datasets, the accuracies obtainable by our method are
similar to the reference methods. The advantage of our classifiers in all
cases is that they are very short and immediately human intelligible as
opposed to the black-box nature of the reference methods.

##  1 Introduction

Human interpretability is one of the main challenges of modern artificial
intelligence (AI). This has lead to an increasing interest in explainable AI
(or XAI for short). The right for receiving explanations is already mentioned
even in various legislative sources such as the European Data Protection
Regulation [13] and the California Consumer Privacy Act [10].

Interpretability is typically divided into _local_ and _global_ variants. The
former relates to providing reasons why a classifier produced a particular
decision with a particular input, while global interpretability means that the
entire behaviour of the classifier—on any input—can be understood due to a
single and sufficiently clear explication.

The current paper introduces and investigates a method for producing globally
interpretable classifiers for binary classification tasks on tabular data. The
method is implemented declaratively by deploying the paradigm of Answer Set
Programming (ASP; see [9] for an overview) and, more specifically its theory-
based extension ASP _modulo difference logic_ , abbreviated ASP(DL) [19].
Using simple difference constraints, the discretization of numeric values
takes place dynamically during the search for actual (optimal) explanations.
As shown in this work, Boolean statements resulting from such dynamic
discretization offer more fine-grained primitives for classifying data points
and give rise to even shorter classifying formulas compared to an earlier
static approach based on medians [18].

The global nature of the interpretability stems from the fact that the
formulas acting as classifiers are very short Boolean assertions about the
attributes of the original datasets. Thus, the formulas read as very simple
statements in natural language about the data studied. E.g., the formula

$\neg(p_{\geq 6}\lor q_{\geq 7}\lor r_{\geq 5})$

is an explanation that we obtain from our experiments on a dataset concerning
breast cancer. The predicate $p_{\geq 6}$ states that the value of $p$, a
measure relating to bare nuclei, is at least $6$. The interpretations are
similar for the other predicates, with $q$ and $r$ denoting measures of clump
thickness and uniformity of cell size, respectively. The formula is clearly
short and it is immediately clear (at least to an expert from the respective
medical field) what the general meaning of the formula is. The accuracy of
this formula on the holdout data is $97.1$ percent.

While this paper is focused on obtaining explanations from data, we note that
our method could also be used to explain black-box classifiers. In this case,
one would first use the black-box classifier to generate data and then explain
that data as described in this paper.

###  1.1 The formula-size method

We next give an overview of our method; for the full details, please see
Section 3. Let $M$ denote a tabular dataset with the set
$\\{Y_{1},\dots,Y_{s}\\}\cup\\{Z_{1},\dots,Z_{k}\\}$ of attributes, each
$Y_{i}$ being numeric and $Z_{i}$ categorical. To enable Boolean
representations of classifiers, we model each $Y_{i}$ and $Z_{j}$ by Boolean
propositions defined as follows.

  1. 1.

For each number $r$ in the range of the attribute $Y_{i}$ in the dataset $M$,
we define a corresponding proposition $p_{i}^{\geq r}$.

  2. 2.

Applying a natural one-hot encoding, we define, for each category $C$ of
$Z_{j}$, a corresponding proposition $p_{j}^{C}$.

These propositions are interpreted in $M$ in a natural way, meaning that a
data point $w$ satisfies $p_{i}^{\geq r}$ (respectively, $p_{j}^{C}$) if $w$
has the value of $Y_{i}$ at least $r$ (respectively, belongs to the category
$C$ of $Z_{j}$).

Given this typically huge set $\Pi$ of proposition symbols at our disposal, we
do the following. We randomly split the data into two parts: a 70 percent part
$W_{70}$ for training and a 30 percent part $W_{30}$ for validation. We
consider the accuracies of formulas over the training and validation sets. By
accuracy we mean the percentage of points on which the formula agrees with the
target attribute of the binary classification task. In the training part, we
do the following for increasing _length bounds_ $\ell$, starting with
$\ell=1$.

  1. 1.

Among all Boolean formulas of length at most $\ell$ and using the propositions
in $\Pi$, we find a formula $\varphi_{\ell}$ with maximum accuracy in
$W_{70}$.

  2. 2.

We record the accuracy $\mathit{acc}_{30}(\varphi_{\ell})$ of $\varphi_{\ell}$
in the validation set $W_{30}$. We then repeat the above step $1$ with
$\ell:=\ell+1$.

Using early stopping, we scan for longer and longer formulas $\varphi_{\ell}$
until overfitting. Our early stopping trigger checks when
$\mathit{acc}_{30}(\varphi_{\ell})$ is smaller than the best accuracy obtained
so far, twice in a row. Supposing $\ell$ is the greatest value we reach before
overfitting, we then compute the final classifier formula $\psi$ by using the
the formula bound $\ell$ and the union $W$ of the training and valuation data.

We also test two other variants of the formula-size method. A more restricted
version we introduced in [18] uses a system based on propositions $p_{\geq
r^{\prime}}$ where $r^{\prime}$ is always the median of the corresponding
numeric attribute. We discuss this further below in the section devoted for
related work. Another more general version has proposition symbols that
specify an interval $[r,r^{\prime}]$ instead of a single pivot value $r$. All
other details are the same for our experiments on these variants of the method
and we compare all three in terms of formula size as well as accuracy.

To test our method, we use a ten-fold cross validation procedure. This means
that each dataset $W$ in the actual experiments is 90 percent of the very
original dataset studied, with 10 percent reserved as holdout data.

Our method is implemented in ASP [9] which brings along the strengths of
declarative programming: the logic programs can be understood as _logical
specifications_ of the concepts involved such as formulas forming the
hypothesis space, the objective function used to measure their quality, and
the schemes used to discretize the data provided as input. Since logic
programs are highly elaboration tolerant, it is easy to formalize these
concepts one-by-one in an incremental way. Moreover, we advocate the overall
modularity of the implementation by the three-layer architecture illustrated
in Figure 1. As a consequence, it is straightforward to experiment with
different encodings as well as alternative concepts.

It is worth noting that the number of Boolean formulas grows asymptotically
extremely fast, with the number of non-equivalent formulas of size at most $k$
being exponential in relation to both $k$ and the number of propositions used.
While symmetry reductions and other techniques can be deployed to reduce the
hypothesis space significantly, it is still natural to question how the method
described above could be feasible in the first place. The answer relates to
overfitting. It is clear that if too long formulas are allowed, they will
overfit. Therefore, we need to use early stopping to avoid overfitting to the
training data. In a large majority of our experiments on real-life datasets,
signs of overfitting emerge very fast, triggering early stopping and leading
to short formulas and feasible runtimes. Furthermore, such short formulas
yield accuracies similar to reference methods. The fact that short formulas
suffice indicates that, in a sense, asymptotic complexity is not the right
measure here. This is important to our method, since too long formulas would
not be inherently interpretable anyway.

###  1.2 Results

We run our method on seven datasets originating from the UCI machine learning
repository. As tabular data is still challenging to deep learning approaches
[15, 7], we use XGBoost [11] and Sklearns implementation of random forests [8]
as our comparison methods. These methods are esteemed state-of-the-art
approaches for tabular data [24, 15] and thus suitable for comparisons.

As mentioned above, we compare three different variants of our method: the
static median-based formula-size method of [18], a dynamic discretization
based on single pivot points $r$ and a version with intervals
$[r,r^{\prime}]$. Here the focus is seeing how much improvement the different
approaches to dynamic discretization can provide. From Figure 2 and Tables 1,
2 and 3, one can see that for several datasets, the dynamic approaches either
improve the accuracy or decrease the formula size needed to obtain a similar
accuracy.

The average accuracy and standard deviation of each method over ten-fold
cross-validation are reported in Figure 2. For all datasets, we obtain
immediately interpretable formulas with accuracy similar to the reference
methods. Formulas obtained as classifiers are reported in Table 1 for the
interval method, in Table 2 for the pivot method and in Table 3 for the median
method. The tables also list average formula lengths over all runs. The
example formulas in the tables were taken from the first round of the ten fold
cross validation in each case.

As an example result, consider the classifier

$p_{[25,60]}\land q_{[128,196]}$

obtained with the interval method from the Diabetes dataset. Here the
attribute $p$ is the age of the person and the attribute $q$ is their measured
glucose value. The formula offers high glucose and age between 25 and 60 as an
explanation for diabetes. Another example is the formula

$p\lor q_{\geq 767}$

obtained with the pivot method from four different splits of the BankMarketing
dataset (with slight variations on the bound of $q$). The classification task
concerns whether a marketing call resulted in subscription for a term deposit.
Here the attribute $p$ indicates whether the customer has previously
subscribed for a deposit and the attribute $q$ gives the duration of the call
in seconds. So a customer is likely to subscribe to a deposit if they have
done so previously or the telemarketer can keep them on the phone for a long
duration. The accuracies of this formula on the holdout data ranged from
$88.1$ percent to $90.3$ percent for different splits.

We conclude that in all cases, our classifiers are immediately globally
interpretable and the accuracy is similar to the ones obtained by XGBoost and
random forests.

The interpretability is a real gain, as the classifiers produced by XGBoost
and random forests are of a black-box nature. Indeed, the classifiers returned
by these methods are ensembles of decision trees, with at least tens or even
thousands of members.

A shortcoming of our method is that it is computationally demanding and thus
not feasible in scenarios where classifiers must be obtained fast. In the
worst case, the experiments with our method took days in a high-power
computation environment, while the comparison methods could be completed much
faster with an ordinary laptop. The runtimes are discussed in Section 5.3.

###  1.3 Related work

Explainable AI is a large and active field, with logic-based explainability
being a growing subfield. See [20] for a thorough survey on logic-based
explainability. Much of the recent work on this topic has been initiated in
[23] and [17], where natural minimality notions such as prime implicants are
utilized as explanations.

Much of the work on logic-based explainability focuses on local aspects of
explainability. The current paper, on the other hand, focuses on global
explainability using short formulas of Boolean logic. It is generally agreed
that small classifiers can be regarded as global explanations [2, 5]. In this
work we aim for very small formulas, following a “simplicity first”
philosophy. Relating to this ideology, already [16] reported that in many of
the 16 datasets studied there, even a single attribute was sometimes enough to
get an accuracy not too drastically smaller than one obtainable by the then
state-of-the-art decision trees.

Other approaches to global explainability using small classifiers include the
articles [5, 2]. These studies investigate the use of small rule lists and
decision trees, which are optimal with respect to size and training error. The
empirical results reported in these papers also demonstrate the surprising
effectiveness of interpretable models on real-world tabular data. For further
related work on global explanations and interpretable AI, see the survey [22].

The quest for short explanations is compatible with the principle of Occam’s
razor [6] and the early search-based approaches, such as current best
hypothesis search [21], already implement this idea. However, the performance
of computers has improved drastically since then and the solver technology
used to implement the search for hypotheses in the present work has emerged by
the side.

In this work, we continue the line of research initiated in [18] where numeric
attributes are still made Boolean statically, with the intuitive reading “the
value of the attribute is above median” for Boolean variables. That paper
presents three experiments with single split cross validation, obtaining
related Boolean classifiers. In the current paper, we redo the experiments of
[18] with the seven datasets studied, two of which were used in [18]. As
mentioned above, the related results are in Figure 2 and Tables 1, 2 and 3.
While reasonable classifiers are obtained, the ones using dynamic predicates
$p_{[r,r^{\prime}]}$ and $p_{\geq r}$ lead to formulas that are shorter or
more accurate, while keeping the running times feasible.

##  2 Preliminaries

A finite set $\tau=\\{p_{1},\dots,p_{k}\\}$ of proposition symbols is called a
Boolean vocabulary. The syntax of propositional logic over the vocabulary
$\tau$, denoted $\mathrm{PL}[\tau]$, is given by the following grammar:

$\varphi::=p\mid\neg\varphi\mid\varphi\lor\varphi\mid\varphi\land\varphi$

where $p\in\tau$. We also define the equivalence connective as the
abbreviation
$\varphi\leftrightarrow\psi:=(\varphi\land\psi)\lor(\neg\varphi\land\neg\psi)$.

A Boolean $\tau$-model is a structure $M=(W,V)$, where the set
$W\neq\emptyset$ is the domain of $M$ and the function
$V:\tau\to\mathcal{P}(W)$ is a $\tau$-valuation. For $p\in\tau$, the set
$V(p)\subseteq W$ is the set of points, where the proposition $p$ is
considered to be true. Given a $\tau$-valuation $V$, we extend it in the
standard way to a valuation $V:\mathrm{PL}[\tau]\to\mathcal{P}(W)$ concerning
all $\mathrm{PL}[\tau]$-formulas.

The size of a formula $\varphi\in\mathrm{PL}[\tau]$ is defined as the number
of occurrences of proposition symbols and connectives in the formula. For
example the size of $\neg(p\land r)$ is 4.

Let $M=(W,V)$ be a $\tau\cup\\{q\\}$-model and let $\varphi$ be a
$\tau$-formula. The accuracy of $\varphi$ over $M$ is

$\mathrm{acc}(\varphi,M):=\frac{|V(\varphi\leftrightarrow q)|}{|W|}.$ (1)

The accuracy of $\varphi$ is the percentage of points where $\varphi$ and the
target attribute $q$ agree.

We move on from propositional structures to structures with both Boolean and
numeric attributes. The datasets we consider also have categorical attributes
with more than two categories, but we use one-hot encoding to Booleanize the
categorical attributes as a preprocessing step. See Subsection 5.1 for more
discussion on this.

A Boolean attribute over a set of data points $W$ is simply a subset of $W$.
As seen above, the interpretations of proposition symbols $p$ in propositional
logic are Boolean attributes $V(p)\subseteq W$. A numeric attribute over $W$
and with a finite codomain $R$ is a function $f:W\rightarrow R$. We denote the
set of functions $f:W\to R$ by $R^{W}$. In this paper, $R$ is a set of
floating point numbers and thus always finite. We assume implicitly that the
codomain $R$ has a linear order $\leq$.

A general vocabulary $\mathcal{X}=\mathcal{B}\uplus\mathcal{N}$ is a disjoint
union of a Boolean vocabulary $\mathcal{B}=\\{p_{1},\dots,p_{\ell}\\}$ and a
numeric vocabulary $\mathcal{N}=\\{s_{1},\dots,s_{m}\\}$. A general valuation
$V$ over $\mathcal{X}$ and $R$ is a function
$V:\mathcal{X}\rightarrow\mathcal{P}(W)\cup R^{W}$ which splits naturally into
two functions $V_{B}:\mathcal{B}\rightarrow\mathcal{P}(W)$ and
$V_{N}:\mathcal{N}\rightarrow R^{W}$. For convenience, we write $s_{i}^{V}$
instead of $V(s_{i})$. A general structure is a tuple $(W,V)$, where $W$ is a
set of data points and $V$ is a general valuation.

We define two different ways to discretize general vocabularies and structures
based on a single pivot value $r\in R$ and an interval
$[r,r^{\prime}]\subseteq R$, respectively. We first define the pivot
discretizations.

Let $\mathcal{X}=\mathcal{B}\uplus\mathcal{N}$ be a general vocabulary. The
complete pivot vocabulary over $\mathcal{N}=\\{s_{1},\dots,s_{m}\\}$ and $R$
is the (finite) set of proposition symbols $p_{s_{i}}^{\geq r}$ where
$s_{i}\in\mathcal{N}$ and $r\in R$. Intuitively, the complete dynamic
vocabulary contains all possible propositions $p_{s_{i}}^{\geq r}$ one might
need when Booleanizing the numeric attributes of $\mathcal{N}$ using a single
pivot point.

A Boolean vocabulary $\eta$ is a $R$-pivot discretization of a numeric
vocabulary $\mathcal{N}$ if $\eta$ is the range of an injection $f$ from
$\mathcal{N}$ into the complete dynamic vocabulary over $\mathcal{N}$ and $R$
such that for all $s\in\mathcal{N}$, $f(s)=p_{s}^{\geq r}$ for some $r\in R$.
Following the intuition of discretization, this turns every numeric attribute
symbol $s$ into a Boolean symbol $p_{s}^{\geq r}$.

Let $M=(W,V)$ be a general structure over
$\mathcal{X}=\mathcal{B}\uplus\mathcal{N}$. A Boolean structure
$M^{\prime}=(W,V^{\prime})$ is a pivot discretization of $M$ if the following
conditions hold.

  1. 1.

$V^{\prime}$ is a function
$V^{\prime}:\mathcal{B}\uplus\eta\to\mathcal{P}(W)$, where $\eta$ is a
$R$-discretization of $\mathcal{N}$.

  2. 2.

$V^{\prime}(p)=V(p)$ for every $p\in\mathcal{B}$.

  3. 3.

$V^{\prime}(p_{s}^{\geq r})=\\{w\in W\mid s^{V}(w)\geq r\\}$ for all
$p_{s}^{\geq r}\in\eta$. Recall that $s^{V}$ denotes the function
$V(s):W\rightarrow R$ that we obtain from $V$ when given the input $s$.

The pivot class over $M$ is the set of all pivot discretizations of $M$,
denoted by $\mathrm{PD}(M)$.

We similarly define the notions of complete interval vocabulary, $R$-interval
discretization, interval discretization of $M$ and interval class
$\mathrm{ID}(M)$ over $M$ by simply replacing the condition $\geq r$ with
inclusion in an interval $[r,r^{\prime}]$.

Let $\mathcal{F}$ be a set of Boolean formulas and $\mathcal{M}$ a set of
Boolean models (possibly over different vocabularies). The maximum accuracy
over $\mathcal{F}$ and $\mathcal{M}$ is the maximum number
$\mathrm{acc}(\varphi,M)$ such that

  1. 1.

$\varphi\in\mathcal{F}$ and $M\in\mathcal{M}$,

  2. 2.

only propositions from the vocabulary of $M$ occur in $\varphi$.

The Boolean formula $\varphi$ is then said to realize the maximum accuracy
(over the model $M$).

Let $M$ be a general structure over $\mathcal{X}=\mathcal{B}\uplus\mathcal{N}$
and $R$. A Boolean formula $\varphi$ is a $M$-pivot-formula (resp.
$M$-interval formula), if $\varphi$ is a Boolean formula over a Boolean
vocabulary $\mathcal{B}\uplus\eta$, where $\eta$ is a $R$-pivot discretization
(resp. $R$-interval discretization) of $\mathcal{N}$. We denote the set of all
$M$-pivot-formulas by $\mathrm{PF}(M)$ and the set of all $M$-interval-
formulas by $\mathrm{IF}(M)$.

For $\ell\in\mathbb{N}$, we denote the set of $M$-pivot-formulas with size at
most $\ell$ by $\mathrm{PF}_{\ell}(M)$. We define $\mathrm{IF}_{\ell}(M)$ in
the same way.

When discussing real-life datasets in the other sections below, we will
suppress the model notation $M=(W,V)$ and only discuss the set $W$ of data
points as well as some subsets of $W$, to be interpreted as submodels of $M$.

##  3 The dynamic formula length method

We next give a step-by-step description of our method. We describe the method
for pivot discretization; the case of interval discretization is identical
with pivot-formulas switched for interval-formulas. The median method uses a
static discretization based on the median; the other details are the same as
with the other variants.

Note that this procedure assumes we have already separated ten percent of the
original dataset to be used as a final testing data. Thus the input to the
process is a dataset $W$ consisting of the remaining 90 percent.

In the procedure, we scan through all formulas up to increasing formula length
bounds. We shall describe in Section 4 how this can be done reasonably
efficiently, avoiding a naive brute force search. We use early stopping to
avoid overfitting and formulas too long to be interpretable.

The method consists of the following stages.

  1. 1.

An input to the method is a tabular dataset $W$ with binary and numeric
attributes $X_{1},\dots,X_{m},q$, where $q$ is a binary target attribute. From
this data $W$, we randomly separate a 30 percent set $W_{v}\subseteq W$ to be
used as _validation data_. The remaining 70 percent $W_{\mathit{tr}}\subseteq
W$ is called the _training data_.

  2. 2.

For a length bound $\ell$ starting with 1, we scan through the set
$\mathrm{PF}_{\ell}(W_{\mathit{tr}})$ of all $W_{\mathit{tr}}$-pivot-formulas
$\varphi$ with maximum length $\ell$, choosing the formula $\varphi_{\ell}$
that realizes the maximum accuracy over the training set $W_{\mathit{tr}}$.

  3. 3.

For each chosen formula $\varphi_{\ell}$, we record its accuracy over the
validation set $W_{v}$. Then, we update the best validation accuracy so far,
denoted $\Delta$. Now, if the accuracies of both $\varphi_{\ell-1}$ and
$\varphi_{\ell}$ are smaller than $\Delta$, we trigger early stopping and move
on to step 4. Otherwise, we return to step 2 with $\ell$ increased by one.

  4. 4.

We scan through the formulas

$\varphi_{1},\dots,\varphi_{\ell-1}$

again and find the smallest number $L\in\\{1,\dots,\ell-1\\}$ such that the
accuracy of $\varphi_{L}$ over $W_{v}$ is the same as that of
$\varphi_{\ell-1}$. This step ensures that if many formulas of different
lengths obtain the same accuracy, we choose the smallest length, leading to
shorter final formulas.

  5. 5.

Finally, we use the data $W=W_{tr}\cup W_{v}$ and the length bound $L$ to
compute the final output formula by scanning through the set
$\mathrm{PF}_{L}(W)$ of all $W$-pivot-formulas $\psi$ with maximum length $L$,
and then choosing the formula with the maximum accuracy over $W$.

We emphasize that the search conducted in steps 2 and 5 for the optimal pivot-
formula is not carried out as a brute force search. Even though the search
remains a challenging problem, we use symmetry reductions and other techniques
to reduce the size of the hypothesis space to make it practically feasible to
conduct the search on real-life datasets. See Section 4.1 for more details.
Our results in Section 5.2, along with the training times we report in Section
5.3, show that the method is useful in practice.

In step 3, we stop once the accuracy of the formulas on the 30-percent
validation data is worse than the best accuracy so far, twice in a row. This
early stopping is done to avoid overfitting. In general, as the length bound
of the formulas increases, the accuracy on both the training data
$W_{\mathit{tr}}$ and the validation data $W_{v}$ improves at the start. At
some point, however, the validation accuracy will start to decline while the
training accuracy keeps on improving. This is when the formulas start to
overfit to the training data. Our simple early stopping trigger aims to avoid
overfitting and also obtain short, interpretable formulas.

Finally in step 5 we use the formula length found in step 4 to compute the
final classifier using all of the available data. This is to ensure that our
early stopping procedure does not lead to throwing away 30 percent of the
data.

Note that even though our method gives a way to select a unique classifier, it
can also be worthwhile to consider the full sequence of formulas with
increasing lengths. This sequence gives an evolving picture of explanations
that can be found from the data.

##  4 Implementation

In this section, we present an implementation of dynamic propositions in ASP
and, in particular, its theory-style extension by simple linear constraints
known as _difference constraints_. The respective extension of ASP is known as
ASP _modulo difference logic_ (ASP(DL)) [19] and it is natively implemented by
the _Clingo-dl_ system222https://github.com/potassco/clingo-dl. Difference
constraints have a normal form $x-y\leq k$ where $x$ and $y$ are variables and
$k$ is a constant. For the purposes of this work, the domain of integers is
sufficient, since the relative orders attribute values are significant rather
than their magnitudes. It is also worth noting that satisfiability can be
checked very efficiently for sets of difference constraints, e.g., by using
the Bellman-Ford algorithm [4].

Without going into syntactic details of rules used in logic programs, we
recall that the _answer sets_ of a logic program $P$ with a vocabulary
$\tau_{P}$ are sets of atoms $S\subseteq\tau_{P}$ satisfying a fix-point
equation $S=\mathrm{LM}({P}^{S})$ where (i) $\mathrm{LM}(\cdot)$ denotes the
_least model_ of the positive program given as argument and (ii) the reduct
${P}^{S}$ is obtained by partially evaluating the negative conditions of $P$
with respect to $S$ [14]. In view of definitions in Section 2, the answer sets
${S_{1}}{,}\ldots{,}{S_{n}}$ of $P$ induce a $\tau_{P}$-model
$\langle{\\{{1}{,}\ldots{,}{n}\\}},{V_{P}}\rangle$ such that
$V_{P}(a)=\\{{i}\mid{a\in S_{i}}\\}$ for $a\in\tau_{P}$. If difference
constraints $(x-y\leq k)$ are incorporated into a logic program $P$, they are
treated as special _difference atoms_ $p(x,y,k)$ in $\tau_{P}$ forming the
respective subsignature $\tau_{d}\subseteq\tau_{P}$ of the program $P$. If
$\tau_{d}\neq\emptyset$, the set of linear inequalities

$\\{{x-y\leq k}\mid{p(x,y,k)\in\tau_{d}\cap
S}\\}~{}\cup~{}\\{{x-y>k}\mid{p(x,y,k)\in\tau_{d}\setminus S}\\}$ (2)

must be additionally satisfied by $S=\mathrm{LM}({P}^{S})$. Note that the
negations $x-y>k$ in (2) can be rewritten in the normal form as $y-x\leq-k-1$,
i.e., they can be treated efficiently as difference constraints. In contrast
with traditional _satisfiability modulo theories_ (SMT) approaches [3], the
equation (2) gives a _strict_ interpretation to difference atoms [19], since
auxiliary atoms $p(x,y,k)$ falsified by $S$ contribute the negation $(x-y>k)$
to the theory being satisfied. This is in line with the minimality of answer
sets in general and _program completion_ [12]: the satisfaction of (2) amounts
to satisfying the completing formula $p(x,y,k)\leftrightarrow(x-y\leq k)$ for
$p(x,y,k)\in\tau_{d}$. When an answer set $S$ is reported to the user, some
concrete values can be derived for variables $x$ and $y$ involved in atoms
$p(x,y,k)\in S\cap\tau_{d}$.

######  Example 1.

One of the datasets studied in this paper has age as an attribute, the actual
ages ranging from 21 to 81 years. Using an integer variable $a$ for an age
limit, we can define its range in terms of two ASP(DL) facts:

$z-a\leq-21$. $a-z\leq 81$.

where $z=0$ is an auxiliary variable (cf. footnote 3). Assuming a particular
data point where the age happens to be $27$, its Booleanization $b$ can be
captured with a rule

$b\leftarrow a-z\leq 27$.

Using auxiliary atoms introduced above, the facts and rules above give rise to
the following answer sets:

$\begin{array}[]{rcl}S_{1}&=&\\{p(z,a,-21),p(a,z,81),p(a,z,27),b\\}\\\
S_{2}&=&\\{p(z,a,-21),p(a,z,81)\\}.\end{array}$

The former can be backed up by assignments $z=0$ and $a=21$, while the latter
by $z=0$ and $a=28$. These minimal non-negative values are given by _Clingo-
dl_ , too.  $\blacksquare$

Our implementation is based on a modular architecture divided into three
layers as illustrated in Figure 1. The topmost layer is responsible for
generating hypotheses, evaluating them with respect to data points, and
calculating the overall accuracy as part of the objective function. The middle
layer encodes the principles for Booleanizing the data and decides which
attribute is the target of explanation, and which attributes may contribute to
the actual explanations. The bottom layer extracts the values of the
attributes from the real data as a preprocessing step. Although the
architecture in Figure 1 looks _stratified_ at the first glance, it does not
have a unique answer set as typical for stratified programs. In particular,
the two highest layers involve choices that create a search space for
optimization. The layers are further detailed in the respective Sections
4.1–4.3 below. By changing logic programs associated with the two top layers,
it is also easy to take different hypothesis spaces into consideration and to
change the strategy for Booleanizing attributes.

Hypothesis Space … | at(3) | at(2) | and | neg | at(1) | or
---|---|---|---|---|---|---
val/3, expl/1, excl/1Dynamic Booleanization

    &diff{p(I)} >= M :- minval(a(I),M), idx(I).
    &diff{p(I)} <= M :- maxval(a(I),M), idx(I).

    val(D,a(I),0) :- &diff{p(I)} > V, rd(D,I,V), idx(I).
    val(D,a(I),1) :- &diff{p(I)} <= V, rd(D,I,V), idx(I).

rd/3Raw Data Figure 1: Three Layers of the Implementation

###  4.1 Generating the hypothesis space

The topmost layer has access to the underlying data set in Boolean form only.
The data is abstracted using a 3-argument predicate val(D,A,V) where D is the
identifier of the data point, A the index of the attribute, and V the value,
i.e., either 0 or 1. Besides this view to the data, the hypothesis generator
receives a target attribute any denied attributes expressed using predicates
expl/1 and excl/1, respectively. Based on this interface, the generator is
completely unaware whether the _original_ data is Boolean or not and its
implementation can be kept independent of this aspect.

In this research, we use the generator from [18] revised to respect the
interface just described. The generator represents formulas as sequences of
attributes and logical operators following the idea of _reverse Polish
notation_. For instance, the sequence $p_{3},p_{2},\land,\neg,p_{1},\lor$
represents the formula $p_{1}\lor\neg(p_{2}\land p_{3})$ that can be
reconstructed as well as evaluated using a stack, e.g., for any data point.
The encoding deploys a parameter l for the maximum length of the
sequence/explanation. The dots “…” in the first cell (see Fig. 1) illustrate
the fact that the sequence may start by unused symbol positions whenever a
shorter (optimal) explanation is feasible. We refer the reader to [18] for
further details of the encoding, but recall that the objective function
maximizes the accuracy (1) as the primary criterion and minimizes the length
of the formula as the secondary criterion. In this way, the parameter l acts
just as an upper bound and the user need not provide the exact length of
explanations in advance.

It should also be emphasized that the search implemented by a conflict-driven
learning solver (such as _Clingo-dl)_ is far from exhaustive brute-force
search. The search space is limited by several factors such as attribute
values admitted by the underlying data set, the value of the objective
function, as well as conflicts learned during the search. In addition, the
encoding of [18] deploys syntactic symmetry constraints to reduce the
hypothesis space.

###  4.2 Dynamic Booleanization of attributes

The values of Boolean attributes and one-hot encoded categorical attributes
can be forwarded unchanged through the val/3 predicate. The dynamic approach
is applied to other numeric attributes that are made available in terms of the
rd(D,I,V) predicate with arguments D for a data point, I for the index of an
attribute, and V for an integer value, hence assuming some mapping from actual
values to the integer domain (see Section 4.3). Then, given an attribute a(I)
indexed by I, its minimum and maximum values can be determined and recorded as
respective facts using predicates minval/2 and maxval/2.

In the dynamic approach based on a single pivot, the first two rules in Figure
1 choose a value p(I) in this range.333Note that the constant bounds $l\leq
x\leq u$ for a variable $x$ can be expressed by $z-x\leq-l$ and $x-z\leq u$ in
the normal form, using an auxiliary variable $z$ standing for “zero”.  The
latter two rules compare the actual value V at a data point D against p(I) and
Booleanize V as 0, if V is smaller than p(I), and as 1, otherwise. Although
difference atoms are exploited in a corner case (cf. footnote 3), the benefit
is that there is no need to sort the values of the attribute a(I) beforehand:
the Booleanized value can be decided based on the pivot p(I). It should be
emphasized that the choice of pivots is completely dynamic: any values are
fine as long as the required Booleanization of the data can be realized, able
to justify the current hypothesis in terms of accuracy and length—in harmony
with the idea of SMT.

The generalization for the interval-based Booleanization is obtained by
introducing two pivots, namely l(I) and u(I), instead of p(I) for an attribute
a(I). The rules in Figure 1 must be refined accordingly. While the minimum
value of an attribute a(I) is the least value for l(I), its maximum value
gives the greatest value for u(I). Moreover, the difference constraint
&diff{u(I) - l(I)} >= 0 is needed to make l(I) smaller than or equal to u(I).
Based on their selected values, the value V of the attribute a(I) is
Booleanized as 1 if V falls in the range from l(I) to u(I) including endpoints
and 0, otherwise.

###  4.3 Treating non-integer numeric values

Given a _raw_ data file, e.g., as a comma-separated values (CSV) file, the
predicate rd/3 can be computed procedurally. We leave the mapping to integers
at the user’s discretion. For instance, non-integer values can be multiplied
by suitable constants and rounded off or truncated at the decimal point. For
the dynamic approach based on pivots, it is more important to preserve the
order than the magnitudes when considering the values of a particular
attribute.

##  5 Experiments

$20$$40$$60$$80$$100$BankMarketingBreastCancerDiabetesGermanCreditHeartDiseaseHepatitisStudentDropout$96.6$$78.0$$81.6$$87.2$$89.1$$97.1$$76.5$$73.4$$79.5$$83.0$$87.5$$89.4$$97.2$$77.3$$73.7$$79.2$$82.1$$86.2$$88.8$$95.4$$74.2$$69.7$$75.9$$83.8$$78.6$$88.8$$95.5$$73.4$$67.3$$75.5$$82.2$$85.8$$88.9$$95.9$$73.5$$66.9$$71.3$$81.5$$84.9$%Interval
FSM Random forestPivot FSM XGBoostMedian FSM UCI repo

Figure 2: The average test accuracies and standard deviations obtained for
each dataset with the pivoted FSM, random forests, XGBoost and the median FSM.
Where available, we include also accuracies reported as Baseline Model
Performance in the UCI repository, although these are not directly comparable
due to the unknown technical particularities behind the UCI numbers.

We compare the interval, pivot and median variants of the formula-size method
empirically to random forests and XGBoost. The two latter methods are very
commonly used state-of-the-art methods for tabular data [24, 15]. We first
present our experimental setup in detail, then discuss the results both in
terms of accuracy and the formulas obtained, and finally report the runtimes
of the methods.

###  5.1 Experimental setup

We test all five methods on seven binary classification tasks with numeric
attributes originating from the UCI machine learning repository:
BankMarketing, BreastCancer, Diabetes, GermanCredit, HeartDisease, Hepatitis
and StudentDropout. The dataset StudentDropout was originally a ternary
classification task, but we converted it into a binary one. The dataset
Diabetes seems to no longer be available from the UCI repository; it can be
downloaded from https://github.com/susanli2016/Machine-Learning-with-
Python/blob/master/diabetes.csv.

The datasets we use contain both categorical and numeric attributes. For
categorical attributes we perform a simple one-hot Booleanization. For each
categorical attribute $X$ and category $a$ of $X$, we create a Boolean
attribute $p^{X}_{a}$, which is true in the points of category $a$ w.r.t. the
attribute $X$.

Note that at least the implementations of random forests and XGBoost we use
(see the end of this section for details) only support numeric data as input.
This means that no feature can be truly categorical as the user is forced to
fix an order for the categories by giving them numeric names. We feel that the
one-hot encoding described above is closer to a truly categorical attribute
than an arbitrary order.

Furthermore, we have experimentally observed that on the datasets that we use
random forests and XGBoost seem to obtain better results on one-hot encoded
data than on numerically encoded categorical data. Thus we use one-hot
encoding as a global preprocessing step for all methods before the
experiments.

We move on to numeric attributes. As dynamic discretization is an important
feature of our method and the reference methods work directly with numeric
attributes, we do not Booleanize numeric attributes beforehand. We do,
however, conduct a small preprocessing step to transform floating point
attributes into integers for easier compatibility with ASP. We simply multiply
the numbers by a suitable power of ten to make them integers. This
preprocessing is not done in the case of random forests and XGBoost, which
receive the numeric attributes as such.

All missing values have been removed from the datasets. In the case of the
Hepatitis dataset, we removed two columns with many missing values. The
diabetes dataset does not have missing values, but it has features which take
values that are inconsistent with background knowledge (e.g., the dataset
contains rows where the value of blood pressure is zero). We removed all rows
with such values.

We use ten-fold cross-validation for all methods and all datasets. That is, we
randomly split the original data $W$ into ten equal (when possible) parts. For
each such 10 percent part, we feed the 90 percent complement into the method
as the full training data $W_{0}$ and set aside the 10 percent part itself as
holdout data $W_{\mathrm{hold}}$. After the method has obtained a final
classifier using $W_{0}$, we record the accuracy of that classifier on
$W_{\mathrm{hold}}$. We report the average and standard deviation of the
accuracy over the ten 90/10 splits given in the beginning.

As the formula size based methods are computationally quite intensive, we
employ a timeout for our runs. This raises the possibility that our early
stopping trigger is never met before a run times out. In this circumstance, we
record the best validation accuracy obtained before timeout and look back at
the smallest formula length that achieves the best accuracy. In a large
majority of our experiments, the early stopping was triggered before the
computation timed out. The exceptions are the dataset StudentDropout and some
individual splits scattered among the other datasets.

For all three formula-size methods, the runs were conducted on a Linux cluster
featuring Intel Xeon 2.40 GHz CPUs with 8 CPU cores per run, employing a
timeout of 72 hours per instance and a memory limit of 64 GB.

For random forests we use the implementation of Scikit-learn. We use nested
cross-validation for tuning the hyperparameters for both random forests and
XGBoost. As for the formula-size methods, we use a $70/30$-split. For
hyperparameter optimization, we use Optuna [1] with 100 trials. The
hyperparameter spaces used for random forests and XGBoost are the same as in
[15]. For the readers convenience, we also report the used hyperparameter
spaces in the Appendix A. We ran these experiments simply on a laptop.

Dataset |  Example formula |  Average length
---|---|---
BM | $p\lor q_{[794,2769]}$ | 2.8
BC | $p_{[1,5]}\land q_{[1,6]}\land r_{[1,4]}$ | 5.5
D | $p_{[25,60]}\land q_{[128,196]}$ | 4.3
GC | $p\lor q_{[4,30]}$ | 3.1
HD | $p_{[0.0,2.6]}\land(q\lor(r\land s_{[176,240]}))$ | 2.9
H | $p_{[2.9,6.4]}$ | 2.3
SD | $p_{[0,2]}\lor\neg q$ | 4.1
Table 1: An example formula obtained by the interval method from the first
90/10 split of each dataset as well as the average length of the ten formulas
obtained for all ten 90/10 splits.

###  5.2 Results

Dataset |  Example formula |  Average length
---|---|---
BM | $(p\lor q_{\geq 767})\land\neg r$ | 4.0
BC | $\neg(p_{\geq 6}\lor q_{\geq 7}\lor r_{\geq 5})$ | 6.5
D | $q_{\geq 158}\lor(p_{\geq 29}\land r_{\geq 27.6}\land s_{\geq 109})$ | 2.9
GC | $\neg r_{\geq 10876}$ | 3.0
HD | $\neg t\lor(q\land r)$ | 4.4
H | $q$ | 2.8
SD | $\neg q\land(p_{\geq 2}\lor\neg r_{\geq 1})$ | 6.9
Table 2: An example formula obtained by the pivot method from the first 90/10 split of each dataset as well as the average length of the ten formulas obtained for all ten 90/10 splits. Dataset |  Example formula |  Average length
---|---|---
BM | $q_{m}\land(s\lor p)$ | 6.6
BC | $\neg((q_{m}\land s_{m})\lor(p_{m}\land t_{m}\land u_{m}))$ | 8.8
D | $s_{m}\land p_{m}\land(t_{m}\lor q_{m})$ | 4.2
GC | $s\lor\neg(t\land u)$ | 5.5
HD | $q$ | 4.2
H | $q$ | 2.5
SD | $\neg(p_{m}\lor(q\land(s\lor t_{m})))$ | 7.1
Table 3: An example formula obtained by the median method from the first 90/10
split of each dataset as well as the average length of the ten formulas
obtained for all ten 90/10 splits. Propositions that were obtained from
medians of numeric attributes are marked with the subscript $m$.

The results of our experiments are summarized in Figure 2 and Tables 1, 2 and
3. In Figure 2 we report the average accuracy and standard deviation of each
method over ten-fold cross-validation. We can see that the formula size
methods give accuracies comparable to those of random forests and XGBoost for
every dataset. The best performance was obtained in the BankMarketing dataset,
where the difference is less than 1 percentage point, while the largest gap of
4 percentage points is found in the GermanCredit dataset.

To measure interpretability, we report the average length of the final formula
for each dataset in Table 1 for the interval method, in Table 2 for the pivot
method and in Table 3 for the median method, along with an example formula for
each dataset. The example formulas were taken from the first 90/10 split of
each ten-fold cross-validation. We see that the formulas we obtain are
typically very short with all three methods with the largest average length
being 8.8 for the median method on the BreastCancer dataset. Keeping in mind
that formula length includes all connectives as well as proposition symbols,
this means that all formulas we obtain are easily human interpretable.

Comparing the three formula size methods, we can observe from Tables 1, 2 and
3 an antitone relationship between the expressive power of a method and the
average lengths of formulas produced by it: the more expressive the method is,
the shorter the formulas it produces. This makes sense, as the expressive
power of a method is related to how quickly the early stopping criteria is
met. Furthermore, the more expressive methods are also computationally more
intensive, so the timeouts can occur at smaller formula lengths, again
shortening the formulas.

From Figure 2 we can also observe that on these datasets, the greater
expressive power of a method does not always translate to better accuracy on
the holdout data. Indeed, the median method, which has the lowest expressive
power, obtains on several datasets accuracies that are similar to or slightly
better than the ones obtained via the alternative formula size based methods.
With no early stopping or cross-validation, the more expressive methods would
clearly give a better accuracy on the training data. Therefore the explanation
for the lower accuracies of the more expressive methods must relate to these
testing measures.

When utilizing the same formula length, the interval and pivot methods have a
higher susceptibility to overfitting to the training data in contrast to the
median method. Although the higher expressive power of the interval and pivot
methods may result in shorter formulas and hence improved interpretability, it
also heightens the likelihood of generating formulas excessively tailored for
the training data, thus lowering holdout accuracy.

Another possible phenomenon is related to the fact that the final classifier
is formed using the union $W$ of the training data and the validation data.
When forming this classifier, the interval $[r,r^{\prime}]$ or the pivot value
$r$ is optimized again, possibly leading to overfitting to $W$ at this final
stage. Possible ways to address these issues in the future would be to
consider different early stopping criteria or to migrate more information from
the validated classifier to the final one.

For the StudentDropout dataset, the interval and pivot methods have obtained
significantly better accuracies than the median method. This is because in
this dataset, the ability of the more expressive methods to choose a cutoff
point other than the median is very useful. The attribute $p$, that occurs in
the example formulas of all methods in Tables 1, 2 and 3 for the
StudentDropout dataset, corresponds to the number of approved curricular units
a student has completed in the 2nd semester. For this attribute, the pivot and
interval methods choose a pivot point of 2, whereas the median of this
attribute is 5. Although StudentDropout is the only dataset among those we
studied where this behavior is clearly seen, it is nevertheless clear that for
many real-world datasets the ability to use dynamic predicates is essential
for good accuracy.

###  5.3 Runtimes

The average training time of each formula-size method is reported in Table 4.
The runtimes, while considerably large, are still very feasible to run on
powerful computation setups that are available today. We do not report the
runtimes of random forests and XGBoost as they were all run on a laptop and
took less than five minutes.

We see that in the majority of cases, the more computationally intensive
method has a higher runtime. However, there are also cases, where this
relationship is reversed. At least one example of this can be seen for any two
of the three methods; we use the pivot and median methods as an example here.
For these two methods, the reversed relationship occurs on the BankMarketing
dataset. There are at least two possible reasons for this. Firstly, the more
expressive pivot method could trigger early stopping at a lower formula
length, thus saving runtime compared to the median method. Secondly, due to
its computational intensity, the pivot method can reach timeout sooner, while
the median method can keep performing runs for increasing formula lengths and
accumulate more runtime.

Dataset |  Interval FSM |  Pivot FSM |  Median FSM
---|---|---|---
BM | 281120 | 98403 | 249263
BC | 95868 | 4089 | 2193
D | 190758 | 61903 | 19
GC | 103844 | 198504 | 146400
HD | 120394 | 3817 | 139
H | 26722 | 119 | 47
SD | 375093 | 480024 | 260895
Table 4: The average training time in seconds per split of each formula-size
method for each dataset.

##  6 Conclusion

In this work, we introduce a method that extracts immediately human
interpretable classifiers from tabular data using dynamic discretizations that
are determined on the fly while optimizing the accuracy and the length of the
classifier. Based on comparing three variants of the method against state-of-
the art methods for classifying tabular data, viz. random forests and XGBoost,
we conclude that explanations obtained by our dynamic methods considered are
competitively accurate while much easier to interpret.

As regards future work, the layered architecture described in Section 4 opens
up several avenues for further research. When it comes to the syntactic form
of explanations, it is interesting to consider further Boolean connectives and
specific normal forms for the classifiers. Yet further strategies for the
discretization of attributes can be obtained as refinements of the
pivot/interval-based approaches of this work.

Acknowledgments. Tomi Janhunen, Antti Kuusisto, Masood Feyzbakhsh Rankooh and
Miikka Vilander were supported by the Academy of Finland consortium project
_Explaining AI via Logic_ (XAILOG), grant numbers 345633 (Janhunen) and 345612
(Kuusisto). Antti Kuusisto and Miikka Vilander were also supported by the
Academy of Finland project _Theory of computational logics_ , grant numbers
324435, 328987 (to December 2021) and 352419, 352420 (from January 2022). The
author names of this article have been ordered on the basis of alphabetical
order.

## References

  * [1] Takuya Akiba, Shotaro Sano, Toshihiko Yanase, Takeru Ohta, and Masanori Koyama.  Optuna: A next-generation hyperparameter optimization framework.  In Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, page 2623–2631, 2019.
  * [2] Elaine Angelino, Nicholas Larus-Stone, Daniel Alabi, Margo Seltzer, and Cynthia Rudin.  Learning certifiably optimal rule lists for categorical data.  Journal of Machine Learning Research, 18(234):1–78, 2018.
  * [3] Clark Barrett and Cesare Tinelli.  Satisfiability Modulo Theories, pages 305–343.  Springer, Cham, 2018.
  * [4] R. Bellman.  On a routing problem.  Quarterly of Applied Mathematics, 16:78–90, 1958.
  * [5] Dimitris Bertsimas and Jack Dunn.  Optimal classification trees.  Machine Learning, 106(7):1039–1082, 2017.
  * [6] Anselm Blumer, Andrzej Ehrenfeucht, David Haussler, and Manfred K. Warmuth.  Occam’s razor.  Inf. Process. Lett., 24(6):377–380, 1987.
  * [7] Vadim Borisov, Tobias Leemann, Kathrin Seßler, Johannes Haug, Martin Pawelczyk, and Gjergji Kasneci.  Deep neural networks and tabular data: A survey.  IEEE Transactions on Neural Networks and Learning Systems, 2022.
  * [8] Leo Breiman.  Random forests.  Machine Learning, 45(1):5–32, 2001.
  * [9] Gerhard Brewka, Thomas Eiter, and Miroslaw Truszczynski.  Answer set programming at a glance.  Communications of the ACM, 54(12):92–103, 2011.
  * [10] California-OAG.  CCPA regulations: Final regulation text, 2021.
  * [11] Tianqi Chen and Carlos Guestrin.  Xgboost: A scalable tree boosting system.  In Balaji Krishnapuram, Mohak Shah, Alexander J. Smola, Charu C. Aggarwal, Dou Shen, and Rajeev Rastogi, editors, Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, San Francisco, CA, USA, August 13-17, 2016, pages 785–794. ACM, 2016.
  * [12] K. Clark.  Negation as failure.  In Logic and Data Bases, pages 293–322. Plenum Press, 1978.
  * [13] EU-GDPR.  Regulation (EU) 2016/679 of the European parliament and of the council, 2016.
  * [14] M. Gelfond and V. Lifschitz.  The stable model semantics for logic programming.  In Proceedings of ICLP’88, pages 1070–1080, 1988.
  * [15] Leo Grinsztajn, Edouard Oyallon, and Gael Varoquaux.  Why do tree-based models still outperform deep learning on typical tabular data?  In Thirty-sixth Conference on Neural Information Processing Systems Datasets and Benchmarks Track, 2022.
  * [16] Robert C. Holte.  Very simple classification rules perform well on most commonly used datasets.  Machine Learning, 11(1):63–90, April 1993.
  * [17] Alexey Ignatiev, Nina Narodytska, and João Marques-Silva.  Abduction-based explanations for machine learning models.  In The Thirty-Third AAAI Conference on Artificial Intelligence, AAAI 2019, The Thirty-First Innovative Applications of Artificial Intelligence Conference, IAAI 2019, The Ninth AAAI Symposium on Educational Advances in Artificial Intelligence, EAAI 2019, Honolulu, Hawaii, USA, January 27 - February 1, 2019, pages 1511–1519. AAAI Press, 2019.
  * [18] Reijo Jaakkola, Tomi Janhunen, Antti Kuusisto, Masood Feyzbakhsh Rankooh, and Miikka Vilander.  Short Boolean formulas as explanations in practice.  In Sarah Alice Gaggl, Maria Vanina Martinez, and Magdalena Ortiz, editors, Logics in Artificial Intelligence - 18th European Conference, JELIA 2023, Dresden, Germany, September 20-22, 2023, Proceedings, volume 14281 of Lecture Notes in Computer Science, pages 90–105. Springer, 2023.
  * [19] T. Janhunen, R. Kaminski, M. Ostrowski, S. Schellhorn, P. Wanko, and T. Schaub.  Clingo goes linear constraints over reals and integers.  Theory and Practice of Logic Programming, 17(5–6):872–888, 2017.
  * [20] João Marques-Silva.  Logic-based explainability in machine learning.  In Leopoldo E. Bertossi and Guohui Xiao, editors, Reasoning Web. Causality, Explanations and Declarative Knowledge - 18th International Summer School 2022, Berlin, Germany, September 27-30, 2022, Tutorial Lectures, volume 13759 of Lecture Notes in Computer Science, pages 24–104. Springer, 2022.
  * [21] Tom M. Mitchell.  Generalization as search.  Artif. Intell., 18(2):203–226, 1982.
  * [22] Cynthia Rudin.  Stop explaining black box machine learning models for high stakes decisions and use interpretable models instead.  Nature Machine Intelligence, 1(5):206–215, 2019.
  * [23] Andy Shih, Arthur Choi, and Adnan Darwiche.  A symbolic approach to explaining Bayesian network classifiers.  In Jérôme Lang, editor, IJCAI, pages 5103–5111, 2018.
  * [24] Ravid Shwartz-Ziv and Amitai Armon.  Tabular data: Deep learning is not all you need.  Inf. Fusion, 81:84–90, 2022.

##  Appendix A Appendix

Here we report the hyperparameter spaces for random forest and XGBoost.
Hyperparameters not listed were kept at their default values.

Random forest

\- max_depth: None, 2, 3, 4

\- n_estimators: Integer sampled from $[9.5,3000.5]$ using log domain

\- criterion: gini, entropy

\- max_features: sqrt, log2, None, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9

\- min_samples_split: 2, 3

\- min_samples_leaf: Integer sampled from $[1.5,50.5]$

\- bootstrap: True, False

\- min_impurity_decrease: 0.0, 0.01, 0.02, 0.05

XGBoost

\- max_depth: Integer sampled from [1,11]

\- n_estimators: Integer sampled from [100, 5900] with step-size being 200

\- min_child_weight: Float sampled from [1.0, 100.0] using log domain

\- subsample: Float sampled from [0.5, 1.0]

\- learning_rate: Float sampled from [1e-5, 0.7] using log domain

\- colsample_bylevel: Float sampled from [0.5, 1.0]

\- gamma: Float sampled from [1e-8, 7.0] using log domain

\- reg_lambda: Float sampled from [1.0, 4.0] using log domain

\- reg_alpha: Float sampled from [1e-8, 100.0] using log domain


regions of low noise levels. It also increases the converged area in regions
of noise level higher than $-0.22dB$. This shows that CGN improves the
convergence properties when applied to the VW material model.

###### Mooney-Rivlin

The second test that will be made is on the MR material model. Here, ten
samples were taken from the test range $[0.3,0.9]$ to generate the test data.
The test range was changed to represent the compressive range of the MR model.
The results of the test are shown in Fig. 4.11.

Figure 4.11: Figure showing comparison between CGN and GN when tested on the
Mooney-Rivlin material model. The CGN method produced a larger coverage of the
test range and decreased the number of iterations needed to reach a result
across the test range.

There it can be seen that CGN provides an increased convergence area. However,
the iteration count at low distances has increased. This shows that CGN
improves the convergence area when applied to the MR material model.

###  4.6 Conclusion

This chapter showed how a non-linear FEM model can be derived and how the two
corrected methods can be used to create forward and inverse models. The test
result for forwards models revealed that the ENR method improves the
convergence behaviour of the FEM model when the materials exhibits hyper-
elastic behaviour. This chapter, furthermore, showed the extrapolation of a
strategy for selecting an adequate value of the modification parameter $c$,
which allows ENR to converge for configurations which cause NR to fail. The
CGN was also shown to provide an improvement over the standard GN method in
highly non-linear systems in either total coverage of the test range or in
integration count or both depending on the system.

##  Chapter 5 Discussion & Conclusion

###  5.1 Discussion

The aim of this paper was to propose new methods which can improve the
convergence of forward and inverse finite element models. The paper has shown
that the two proposed methods possess improved convergence properties when
compared to their original counterparts. These improved convergence properties
can then translate to improved convergence of the overarching forward and
inverse models that this paper began with. When constructing a forward model
of an object that is made of a material that is hyper-elastic, the usage of
the extended Newton-Raphson method can greatly improve the convergence
properties of the model. The degree to which this improvement will be realised
is based on the choice of the parameter $c$, for which a strategy has been
proposed. The improvement added by ENR can be twofold. It can, firstly, reduce
the iteration count needed to find an answer and, secondly, it can expand the
convergence basin of the problem allowing for more extreme configurations to
be tested.

When creating an inverse model, the same improvements can be seen when using
the corrected Gauss-Newton method. It was shown in multiple examples that the
method yielded a wider convergence range and a reduced iteration count. The
CGN method also proved to be more resilient to noise in the data. This can
make the CGN a vital tool when creating inverse models using data that has a
high degree of noise.

Future development of the methods proposed by this paper can be concentrated
on finding alternative formulations of the non-linear modification function. A
few different functions were tested in the development of this paper. However,
their performance was inadequate due to their insufficient continuity. The
non-linear modification can also be applied to other root-finding methods such
as BFGS. During the course of this project, an attempt was made to create a
multivariate version of the secant method which does not require a large
number of starting choices. This was achieved by taking the main diagonal of
the approximation of the Jacobian matrix. A similar strategy has been proven
to work in quasi-Newtonian methods [4], so pursuing a simplified multivariate
version of the secant method may yield valuable results. A line search
algorithm can be applied to the CGN step to further improve the convergence
behaviour. It is known that the standard GN method can approach quadratic
convergence with an appropriate line search. Hence, if an appropriate line
search is used with CGN better than quadratic convergence could be achieved
due to CGN’s improved rate.

###  5.2 Conclusion

In conclusion, this paper has covered the improvement of the Newton-Raphson
root-finding algorithm through a non-linear modification aimed to reduce the
non-linearity of a function. It also covered how the Gauss-Newton step can be
corrected to produce improved convergence behaviour. Finally, the paper showed
how to derive a non-linear version of FEM and apply the two improved
approaches to forward and inverse models.

This paper showed that the ENR method can perform better than NR in a variety
of highly non-linear functions. When ENR was applied to a FEM model, the paper
showed that for the same initial starting point as NR, ENR was able to
converge to a solution in fewer steps. For some configurations in which NR
failed to converge, ENR was able to produce an answer from the same starting
point. This paper, additionally, showed a bracketing behaviour exhibited by
ENR when it comes to the choice of $c$. This led to the development of a
strategy for selecting $c$ which produced consistent results.

The paper also demonstrated that the CGN method extends the convergence basin
of GN in a variety of highly non-linear functions, yet performs identically in
the linear case. CGN was also shown to reduce the iteration count needed to
reach an answer in some cases. When applied to the two material models
examined by this paper, CGN was, again, able to produce a larger convergence
basin and a reduced step count. This suggests that it can be used as an
alternative to GN when dealing with highly non-linear systems, such as those
presented in this paper.

Both improved methods require extra computational resources. However, their
convergence properties may make them a suitable choice for the analysis of
plastics, soft tissues and other materials which exhibit similar hyper-elastic
behaviour.

## Bibliography

  * [1] Niels Henrik Abel  “Démonstration de l’impossibilité de la résolution algébrique des équations générales qui passent le quatrième degré.”  In _Œuvres Complètes de Niels Henrik Abel_ 2, 1826
  * [2] Rohan Abeyaratne  “Continuum Mechanics, Volume II of Lecture Notes on the Mechanics of Solids”  Rohan Abeyaratne, 1988, 2012
  * [3] Ankush Aggarwal and Sanjay Pant  “Beyond Newton: A New Root-Finding Fixed-Point Iteration for Nonlinear Equations”  In _Algorithms_ 13.78, 2020
  * [4] Neculai Andrei  “Diagonal Approximation of the Hessian by Finite Differences for Unconstrained Optimization”  In _Journal of Optimization Theory and Applications_ 185, 2020  DOI: 10.1007/s10957-020-01676-z
  * [5] T Belytschko, WK Liu, B Moran and K Elkhodary  “Nonlinear Finite Elements for Continua and Structures”  John Wiley & Sons, Ltd, 2014
  * [6] Ake Bjorck  “9. Nonlinear Least Squares Problems”  In _Numerical Methods for Least Squares Problems_ SIAM, 1996, pp. 339–358  DOI: 10.1137/1.9781611971484.ch9
  * [7] BARBARA BLASCHKE, ANDREAS NEUBAUER and OTMAR SCHERZER  “On convergence rates for the Iteratively regularized Gauss-Newton method”  In _IMA Journal of Numerical Analysis_ 17.3, 1997, pp. 421–436  DOI: 10.1093/imanum/17.3.421
  * [8] Jared L. Deutsch and Clayton V. Deutsch  “Latin hypercube sampling with multidimensional uniformity”  In _Journal of Statistical Planning and Inference_ 142.3, 2012, pp. 763–772  DOI: https://doi.org/10.1016/j.jspi.2011.09.016
  * [9] Jacob Fish and Ted Belytschko  “A First Course in Finite Elements”  In _A First Course in Finite Elements_ John Wiley & Sons, Ltd, 2007  DOI: https://doi.org/10.1002/9780470510858.oth1
  * [10] A. Galántai  “The theory of Newton’s method” Numerical Analysis 2000. Vol. IV: Optimization and Nonlinear Equations  In _Journal of Computational and Applied Mathematics_ 124.1, 2000, pp. 25–44  DOI: https://doi.org/10.1016/S0377-0427(00)00435-0
  * [11] Charles R. Harris et al.  “Array programming with NumPy”  In _Nature_ 585.7825 Springer ScienceBusiness Media LLC, 2020, pp. 357–362  DOI: 10.1038/s41586-020-2649-2
  * [12] Glyn James  “Advanced modern engineering mathematics”  Pearson Prentice Hall, 2011, pp. 723–724802–814
  * [13] R.. Kellogg, T.. Li and J. Yorke  “A Constructive Proof of the Brouwer Fixed-Point Theorem and Computational Results”  In _SIAM Journal on Numerical Analysis_ 13.4 Society for IndustrialApplied Mathematics, 1976, pp. 473–483  URL: http://www.jstor.org/stable/2156239
  * [14] Adam B. Levy  “Attraction in Numerical Minimization”  Springer, Cham, 2018, pp. 33–38
  * [15] P…. Martins, R.. Natal Jorge and A… Ferreira  “A Comparative Study of Several Material Models for Prediction of Hyperelastic Properties: Application to Silicone-Rubber and Soft Tissues”  In _Strain_ 42.3, 2006, pp. 135–147  DOI: https://doi.org/10.1111/j.1475-1305.2006.00257.x
  * [16] Pedro Martins, Renato Natal Jorge and Antonio Ferreira  “A comparative study of several material models for prediction of hyper-elastic properties: Application to Silicone-Rubber and Soft Tissues”  In _Strain_ 42, 2012
  * [17] George Matsaglia and George P.. Styan  “Equalities and Inequalities for Ranks of Matrices”  In _Linear and Multilinear Algebra_ 2.3 Taylor & Francis, 1974, pp. 269–292  DOI: 10.1080/03081087408817070
  * [18] Yue Mei, Daniel E. Hurtado, Sanjay Pant and Ankush Aggarwal  “On improving the numerical convergence of highly nonlinear elasticity problems”  In _Computer Methods in Applied Mechanics and Engineering_ 337, 2018, pp. 110–127  DOI: https://doi.org/10.1016/j.cma.2018.03.033
  * [19] M. Mooney  “A Theory of Large Elastic Deformation”  In _Journal of Applied Physics_ 11.9, 1940, pp. 582–592  DOI: 10.1063/1.1712836
  * [20] Sahil Moza  “sahilm89/lhsmdu: Latin Hypercube Sampling with Multi-Dimensional Uniformity (LHSMDU)”  Zenodo, 2020  DOI: 10.5281/zenodo.3929531
  * [21] K Naidu et al.  “COMPARATIVE STUDY OF MATERIALS FOR REINFORCEMENT IN BRIDGE PIER BEARINGS”  In _International Journal of Civil Engineering and Technology_ 10, 2019, pp. 912–919
  * [22] Jorge Nocedal and Stephen J. Wright  “Numerical Optimization”  Springer, New York, NY, 2006  DOI: https://doi.org/10.1007/978-0-387-40065-5
  * [23] Joanna M. Papakonstantinou and Richard A. Tapia  “Origin and Evolution of the Secant Method in One Dimension”  In _The American Mathematical Monthly_ 120.6 Mathematical Association of America, 2013, pp. 500–518  URL: http://www.jstor.org/stable/10.4169/amer.math.monthly.120.06.500
  * [24] R. Penrose  “A generalized inverse for matrices”  In _Mathematical Proceedings of the Cambridge Philosophical Society_ 51.3 Cambridge University Press, 1955, pp. 406–413  DOI: 10.1017/S0305004100030401
  * [25] Endre Süli and David F. Mayers  “An Introduction to Numerical Analysis”  Cambridge University Press, 2003  DOI: 10.1017/CBO9780511801181
  * [26] Piotr Szurgott and Lukasz Jarzebski  “Selection of a Hyper-Elastic Material Model - a Case Study for a Polyurethane Component”  In _Latin American Journal of Solids and Structures_ 16, 2019  DOI: 10.1590/1679-78255477
  * [27] Guorong Wang, Yimin Wei and Sanzheng Qiao  “Structured Matrices and Their Generalized Inverses”  In _Generalized Inverses: Theory and Computations_ Singapore: Springer Singapore, 2018, pp. 225–231  DOI: 10.1007/978-981-13-0146-9_6