# Explanations for Monotonic Classifiers

Joao Marques-Silva  Thomas Gerspacher  Martin Cooper  Alexey Ignatiev  Nina
Narodytska

###### Abstract

In many classification tasks there is a requirement of monotonicity.
Concretely, if all else remains constant, increasing (resp. decreasing) the
value of one or more features must not decrease (resp. increase) the value of
the prediction. Despite comprehensive efforts on learning monotonic
classifiers, dedicated approaches for explaining monotonic classifiers are
scarce and classifier-specific. This paper describes novel algorithms for the
computation of one formal explanation of a (black-box) monotonic classifier.
These novel algorithms are polynomial in the run time complexity of the
classifier and the number of features. Furthermore, the paper presents a
practically efficient model-agnostic algorithm for enumerating formal
explanations.

Machine Learning, ICML

##  1 Introduction

Monotonicity is an often required constraint in practical applications of
machine learning. Broadly, a monotonicity constraint requires that increasing
(resp. decreasing) the value of one or more features, while keep the other
features constant, will not cause the prediction to decrease (resp. increase).
Monotonicity has been investigated in the context of classification (Cano et
al., 2019), including neural networks (Sill, 1997; Magdon-Ismail & Sill, 2008;
Bonakdarpour et al., 2018; Sivaraman et al., 2020; Liu et al., 2020), random
forests (Bartley et al., 2019) and rule ensembles (Bartley et al., 2018),
decision trees (Ben-David et al., 1989; Ben-David, 1995), decision lists
(Potharst & Bioch, 2000) and decision rules (Verbeke et al., 2017), support
vector machines (Bartley et al., 2016), nearest-neighbor classifiers
(Duivesteijn & Feelders, 2008), among others (Fard et al., 2016; Gupta et al.,
2016; You et al., 2017; Bonakdarpour et al., 2018). Monotonicity has been
studied in bayesian networks (van der Gaag et al., 2004; Shih et al., 2018),
active learning (Barile & Feelders, 2012) and, more recently, in fairness
(Wang & Gupta, 2020).

To a much lesser extent, monotonicity has also been studied from the
perspective of explainability, with one recent example being the study of the
explainability of monotonic bayesian networks (Shih et al., 2018). This work
proposes to compile different families of bayesian networks, including naive
bayes and monotonic networks, into a decision diagram, which can then be used
for computing PI-explanations222Given some feature space point $\mathbf{v}$, a
PI-explanation is a subset-minimal subset of features which, the assignment of
the corresponding coordinate value in $\mathbf{v}$, is sufficient for the
prediction.. Approaches based on an intermediate (knowledge) compilation step
are characterized by two main drawbacks, namely their worst-case complexity,
which is exponential both in time and in the size of the representation, but
also the fact that these approaches are not model-agnostic, i.e. some formal
logic representation of the model must be known and reasoned about. Clearly,
model-agnostic heuristic approaches, which include LIME (Ribeiro et al.,
2016), SHAP (Lundberg & Lee, 2017), or Anchor (Lundberg & Lee, 2017), can also
be applied to explaining monotonic classifiers. However, these approaches do
not readily exploit monotonicity, and both the theoretical and practical
performance may be discouraging333In fact, there are recent negative results
on the tractability of exact SHAP learning (Van den Broeck et al., 2020)..
Furthermore, heuristic approaches offer no formal guarantees of rigor, e.g. an
Anchor explanation may be consistent with points in feature space for which
the model’s prediction differ from the target prediction (Ignatiev, 2020).

On a more positive note, recent work proposed polynomial-time exact algorithms
for computing PI-explanations explanations of different classes of classifiers
(Marques-Silva et al., 2020), namely linear and naive bayes classifiers. These
results were complemented by the observation that, for ML models related with
some classes of knowledge representation languages, PI-explanations can also
be computed in polynomial time (Audemard et al., 2020).

This paper extends these initial results to the case of monotonic classifiers,
in a number of ways. First, the paper proposes model-agnostic algorithms for
computing PI-explanations and contrastive explanations (Miller, 2019) for
_any_ monotonic ML model. Second, the complexity of the proposed algorithms is
shown to be polynomial on the time required to run the (black-box) monotonic
classifier and the number of features. Third, the paper proposes an algorithm
for the iterative enumeration of formal explanations444The term formal
explanation is used in contrast with heuristic explanation (Ribeiro et al.,
2016; Lundberg & Lee, 2017; Ribeiro et al., 2018) and it will be defined
precisely in Section 2.. (This algorithm is worst-case exponential, but it is
shown to be remarkably efficient in practice.)

The paper is organized as follows. Section 2 introduces the notation and
definitions used in the rest of the paper. Section 3 details algorithms for
computing one or more formal explanations of monotonic classifiers. Section 4
summarizes initial experiments, which confirm the scalability of the proposed
algorithms. The paper concludes in Section 5.

##  2 Preliminaries

#### Classification problems.

A classification problem is defined on a set of features (or attributes)
${\mathcal{F}}=\\{1,\ldots,N\\}$ and a set of classes
${\mathcal{K}}=\\{c_{1},c_{2},\ldots,c_{M}\\}$. Each feature
$i\in{\mathcal{F}}$ takes values from a domain $\mathbb{D}_{i}$. Domains are
ordinal and bounded, and each domain can be defined on boolean, integer or
real values. If $x_{i}\in\mathbb{D}_{i}$, then $\lambda(i)$ and $\mu(i)$
denote respectively the smallest and largest values that $x_{i}$ can take,
i.e. $\lambda(i)\leq{x_{i}}\leq\mu(i)$. Feature space is defined as
$\mathbb{F}=\mathbb{D}_{1}\times{\mathbb{D}_{2}}\times\ldots\times{\mathbb{D}_{N}}$.
The notation $\mathbf{x}=(x_{1},\ldots,x_{N})$ denotes an arbitrary point in
feature space, where each $x_{i}$ is a variable taking values from
$\mathbb{D}_{i}$. Moreover, the notation $\mathbf{v}=(v_{1},\ldots,v_{N})$
represents a specific point in feature space, where each $v_{i}$ is a constant
representing one concrete value from $\mathbb{D}_{i}$. An _instance_ (or
example) denotes a pair $(\mathbf{v},c)$, where $\mathbf{v}\in\mathbb{F}$ and
$c\in{\mathcal{K}}$. (We also use the term _instance_ to refer to
$\mathbf{v}$, leaving $c$ implicit.) An ML classifier $\mathbb{C}$ is
characterized by a _classification function_ $\kappa$ that maps feature space
$\mathbb{F}$ into the set of classes ${\mathcal{K}}$, i.e.
$\kappa:\mathbb{F}\to{\mathcal{K}}$.

#### Monotonic classification.

Given two points in feature space $\mathbf{a}$ and $\mathbf{b}$,
$\mathbf{a}\leq\mathbf{b}$ if $a_{i}\leq{b_{i}}$, for all
$i\in\\{1,\ldots,N\\}$. A set of classes
${\mathcal{K}}=\\{c_{1},\ldots,c_{M}\\}$ is _ordered_ if it respects a total
order $\preccurlyeq$, with
$c_{1}\preccurlyeq{c_{2}}\preccurlyeq\ldots\preccurlyeq{c_{M}}$. An ML
classifier $\mathbb{C}$ is fully monotonic if the associated classification
function is monotonic, i.e.
$\mathbf{a}\leq\mathbf{b}\Rightarrow\kappa(\mathbf{a})\preccurlyeq\kappa(\mathbf{b})$555The
paper adopts the classification of monotonic classifiers proposed in earlier
work (Daniels & Velikova, 2010).. Throughout the paper, when referring to a
monotonic classifier, this signifies a fully monotonic classifier. In
addition, the interaction with a classifier is restricted to computing the
value of $\kappa(\mathbf{v})$, for some point $\mathbf{v}\in\mathbb{F}$, i.e.
the classifier will be viewed as a black-box.

######  Example 1 (Running example).

Let us consider a classifier for predicting student grades. We assume that the
classifier has learned the following formula (after being trained with grades
of students from different cohorts):

$\begin{array}[]{rcl}S&=&\max\left[0.3\times{Q}+0.6\times{X}+0.1\times{H},R\right]\\\\[2.0pt]
M&=&\textnormal{ite}(S\geq 9,A,\textnormal{ite}(S\geq
7,B,\textnormal{ite}(S\geq 5,C,\\\\[1.5pt] &&\textnormal{ite}(S\geq
4,D,\textnormal{ite}(S\geq 2,E,F)))))\\\ \end{array}$

$S$, $Q$, $X$, $H$ and $R$ denote, respectively, the final score, the marks on
the quiz, the exam, the homework, and the mark of an optional research
project. Each mark ranges from 0 to 10\. (For the optional mark $R$, the final
mark is 0 if the student opts out.) The final score is the largest of the two
marks, as shown above. Moreover, the final grade $M$ is defined using an ite
(if-then-else) operator, and ranges from $A$ to $F$. As a result, $Q$, $X$,
$H$ and $R$ represent the features of the classification problem, respectively
numbered 1, 2, 3 and 4, and so ${\mathcal{F}}=\\{1,2,3,4\\}$. Each feature
takes values from $[0,10]$, i.e. $\lambda(i)=0$ and $\mu(i)=10$. The set of
classes is ${\mathcal{K}}=\\{A,B,C,D,E,F\\}$, with
${F}\preccurlyeq{E}\preccurlyeq{D}\preccurlyeq{C}\preccurlyeq{B}\preccurlyeq{A}$.
Clearly, the complete classifier (that given the different marks computes a
final grade) is monotonic. Moreover, we will we consider a specific point of
feature space representing student $s_{1}$, $(Q,X,H,R)=(10,10,5,0)$, with a
predicted grade of $A$, i.e. $\kappa(10,10,5,0)=A$.

#### Abductive and contrastive explanations.

We now define formal explanations. Prime implicant (PI) explanations (Shih et
al., 2018) denote a minimal set of literals (relating a feature value $x_{i}$
and a constant $v_{i}$ from its domain $\mathbb{D}_{i}$) that are sufficient
for the prediction666PI-explanations are related with abduction, and so are
also referred to as abductive explanations (AXp) (Ignatiev et al., 2019). More
recently, PI-explanations have been studied from a knowledge compilation
perspective (Audemard et al., 2020).. Formally, given
$\mathbf{v}=(v_{1},\ldots,v_{N})\in\mathbb{F}$ with $\kappa(\mathbf{v})=c$, a
PI-explanation (AXp) is any minimal subset
${\mathcal{X}}\subseteq{\mathcal{F}}$ such that,

$\forall(\mathbf{x}\in\mathbb{F}).\left[\bigwedge\nolimits_{i\in{{\mathcal{X}}}}(x_{i}=v_{i})\right]\operatorname*{\rightarrow}(\kappa(\mathbf{x})=c)$
(1)

AXp’s can be viewed as answering a ’Why?’ question, i.e. why is some
prediction made given some point in feature space. A different view of
explanations is a contrastive explanation (Miller, 2019), which answers a ’Why
Not?’ question, i.e. which features can be changed to change the prediction. A
formal definition of contrastive explanation is proposed in recent work
(Ignatiev et al., 2020). Given $\mathbf{v}=(v_{1},\ldots,v_{N})\in\mathbb{F}$
with $\kappa(\mathbf{v})=c$, a CXp is any minimal subset
${\mathcal{Y}}\subseteq{\mathcal{F}}$ such that,

$\exists(\mathbf{x}\in\mathbb{F}).\bigwedge\nolimits_{j\in{\mathcal{F}}\setminus{\mathcal{Y}}}(x_{j}=v_{j})\land(\kappa(\mathbf{x})\not=c)$
(2)

Building on the results of R. Reiter in model-based diagnosis (Reiter, 1987),
(Ignatiev et al., 2020) proves a minimal hitting set (MHS) duality relation
between AXp’s and CXp’s, i.e. AXp’s are MHSes of CXp’s and vice-versa.

######  Example 2 (AXp’s & CXp’s).

As can be readily observed (from the expression for $M$ in Example 1), as long
as $Q$ and $X$ take value 10, the prediction will be $A$, _independently_ of
the values given to $H$ and $R$. Hence, given $(Q,X,H,R)=(10,10,5,0)$, one AXp
is $\\{1,2\\}$. Moreover, to obtain a different prediction, it suffices to
allow a suitable change of value in $Q$ (or alternatively in $X$). Hence,
given $(Q,X,H,R)=(10,10,5,0)$, one CXp is $\\{1\\}$ (and another is
$\\{2\\}$). As can be observed, $\\{1,2\\}$ is the only MHS of
$\\{\\{1\\},\\{2\\}\\}$ and vice-versa. These are the only AXp’s and CXp’s for
the example instance.

#### Boolean satisfiability (SAT).

SAT is the decision problem for propositional logic. The paper uses standard
notation and definitions e.g. (Biere et al., 2009). A propositional formula is
defined on a set $U$ of boolean variables, where the domain of each variable
$u_{i}\in{U}$ is $\\{0,1\\}$. We consider conjunctive normal form (CNF)
formulas, where a formula is a conjunction of clauses, each clause is a
disjunction of literals, and a literal is a variable $u_{i}$ or its negation
$\neg{u_{i}}$. CNF formulas and SAT reasoners are used in Section 3.2.

##  3 Explanations for Monotonic Classifiers

This section describes three algorithms. The first algorithm serves to compute
one AXp (and is referred to as $\mathsf{findAXp}$). Its complexity is
polynomial in the run time complexity of the classifier. The second algorithm
serves to compute one CXp (and is referred to as $\mathsf{findCXp}$). It has
the same polynomial complexity as $\mathsf{findAXp}$. The third algorithm
shows how to use SAT reasoners for iteratively enumerating AXp’s or CXp’s.
This algorithm is inspired by earlier work (Liffiton et al., 2016), but with
key observations that minimize the number of times a SAT reasoner is called.
This algorithm is based on the other two algorithms, and is described in
Section 3.2.

One key property of the three algorithms is that, besides knowing that the
classifier is monotonic, _no_ additional information about the classifier is
required. Indeed, the algorithms described in this section only require
running the classifier for specific points in feature space. Thus, and
similarly to LIME (Ribeiro et al., 2016), SHAP (Lundberg & Lee, 2017) or
Anchor (Ribeiro et al., 2018), the algorithms proposed in this section are
model-agnostic. However, and in contrast also with LIME, SHAP or Anchor, the
proposed algorithms compute rigorously defined AXp’s, CXp’s, and also serve
for the enumeration of explanations.

###  3.1 Finding One AXp and One CXp

The two algorithms $\mathsf{findAXp}$ and $\mathsf{findCXp}$ (shown as
Algorithm 1 and Algorithm 2) share a number of common concepts, while solving
different problems. These concepts are summarized next.

Algorithm 1 Finding one AXp –
$\mathsf{findAXp}({\mathcal{F}},{\mathcal{S}},\mathbf{v})$

1: $\mathbf{v}_{\textnormal{L}}\leftarrow(v_{1},\ldots,v_{N})$

2: $\mathbf{v}_{\textnormal{U}}\leftarrow(v_{1},\ldots,v_{N})$ // Ensures:
$\kappa(\mathbf{v}_{\textnormal{L}})=\kappa(\mathbf{v}_{\textnormal{U}})$

3:
$({\mathcal{C}},{\mathcal{D}},{\mathcal{P}})\leftarrow({\mathcal{F}},\emptyset,\emptyset)$

4: for all $i\in{\mathcal{S}}$ do

5:
$(\mathbf{v}_{\textnormal{L}},\mathbf{v}_{\textnormal{U}},{\mathcal{C}},{\mathcal{D}})\leftarrow\mathsf{FreeAttr}(i,\mathbf{v},\mathbf{v}_{\textnormal{L}},\mathbf{v}_{\textnormal{U}},{\mathcal{C}},{\mathcal{D}})$

6: end for// Require:
$\kappa(\mathbf{v}_{\textnormal{L}})=\kappa(\mathbf{v}_{\textnormal{U}})$,
given ${\mathcal{S}}$

7: for all $i\in{\mathcal{F}}\setminus{\mathcal{S}}$ do // Loop inv.:
$\kappa(\mathbf{v}_{\textnormal{L}})=\kappa(\mathbf{v}_{\textnormal{U}})$

8:
$(\mathbf{v}_{\textnormal{L}},\mathbf{v}_{\textnormal{U}},{\mathcal{C}},{\mathcal{D}})\leftarrow\mathsf{FreeAttr}(i,\mathbf{v},\mathbf{v}_{\textnormal{L}},\mathbf{v}_{\textnormal{U}},{\mathcal{C}},{\mathcal{D}})$

9: if
$\kappa(\mathbf{v}_{\textnormal{L}})\not=\kappa(\mathbf{v}_{\textnormal{U}})$
then // If invariant broken, fix it

10:
$(\mathbf{v}_{\textnormal{L}},\mathbf{v}_{\textnormal{U}},{\mathcal{D}},{\mathcal{P}})\leftarrow\mathsf{FixAttr}(i,\mathbf{v},\mathbf{v}_{\textnormal{L}},\mathbf{v}_{\textnormal{U}},{\mathcal{D}},{\mathcal{P}})$

11: end if

12: end for

13: return ${\mathcal{P}}$

Algorithm 2 Finding one CXp –
$\mathsf{findCXp}({\mathcal{F}},{\mathcal{S}},\mathbf{v})$

1: $\mathbf{v}_{\textnormal{L}}\leftarrow(\lambda(1),\ldots,\lambda(N))$

2: $\mathbf{v}_{\textnormal{U}}\leftarrow(\mu(1),\ldots,\mu(N))$// Ensures:
$\kappa(\mathbf{v}_{\textnormal{L}})\not=\kappa(\mathbf{v}_{\textnormal{U}})$

3:
$({\mathcal{C}},{\mathcal{D}},{\mathcal{P}})\leftarrow({\mathcal{F}},\emptyset,\emptyset)$

4: for all $i\in{\mathcal{S}}$ do

5:
$(\mathbf{v}_{\textnormal{L}},\mathbf{v}_{\textnormal{U}},{\mathcal{C}},{\mathcal{D}})\leftarrow\mathsf{FixAttr}(i,\mathbf{v},\mathbf{v}_{\textnormal{L}},\mathbf{v}_{\textnormal{U}},{\mathcal{C}},{\mathcal{D}})$

6: end for// Require:
$\kappa(\mathbf{v}_{\textnormal{L}})\not=\kappa(\mathbf{v}_{\textnormal{U}})$,
given ${\mathcal{S}}$

7: for all $i\in{\mathcal{F}}\setminus{\mathcal{S}}$ do // Loop inv.:
$\kappa(\mathbf{v}_{\textnormal{L}})\not=\kappa(\mathbf{v}_{\textnormal{U}})$

8:
$(\mathbf{v}_{\textnormal{L}},\mathbf{v}_{\textnormal{U}},{\mathcal{C}},{\mathcal{D}})\leftarrow\mathsf{FixAttr}(i,\mathbf{v},\mathbf{v}_{\textnormal{L}},\mathbf{v}_{\textnormal{U}},{\mathcal{C}},{\mathcal{D}})$

9: if
$\kappa(\mathbf{v}_{\textnormal{L}})=\kappa(\mathbf{v}_{\textnormal{U}})$ then
// If invariant broken, fix it

10:
$(\mathbf{v}_{\textnormal{L}},\mathbf{v}_{\textnormal{U}},{\mathcal{D}},{\mathcal{P}})\leftarrow\mathsf{FreeAttr}(i,\mathbf{v},\mathbf{v}_{\textnormal{L}},\mathbf{v}_{\textnormal{U}},{\mathcal{D}},{\mathcal{P}})$

11: end if

12: end for

13: return ${\mathcal{P}}$

The two algorithms iteratively update three sets of features (${\mathcal{C}}$,
${\mathcal{D}}$ and ${\mathcal{P}}$) and two points in feature space
($\mathbf{v}_{\textnormal{L}}$ and $\mathbf{v}_{\textnormal{U}}$). Using these
variables, the two algorithms maintain two invariants. The first invariant is
that ${\mathcal{C}}$, ${\mathcal{D}}$ and ${\mathcal{P}}$ form a partition of
${\mathcal{F}}$, and represent respectively the candidate, dropped and picked
sets of features (with the picked features denoting those that are included
either in an AXp or an CXp). The second invariant serves to ensure that the
selected set of features satisfies (1) (for $\mathsf{findAXp}$) or (2) (for
$\mathsf{findCXp}$). Maintaining this invariant, requires iteratively updating
two points
$\mathbf{v}_{\textnormal{L}}=(v_{\textnormal{L}_{1}},\ldots,v_{\textnormal{L}_{N}})$
and
$\mathbf{v}_{\textnormal{U}}=(v_{\textnormal{U}_{1}},\ldots,v_{\textnormal{U}_{N}})$,
denoting respectively lower and upper bounds on the class values that can be
obtained given the features that are allowed to take any value in their
domain.

Feat. | Initial values | Changed values | Predictions | Dec. | Resulting values
---|---|---|---|---|---
$\mathbf{v}_{\textnormal{L}}$ | $\mathbf{v}_{\textnormal{U}}$ | $\mathbf{v}_{\textnormal{L}}$ | $\mathbf{v}_{\textnormal{U}}$ | $\kappa(\mathbf{v}_{\textnormal{L}})$ | $\kappa(\mathbf{v}_{\textnormal{U}})$ | $\mathbf{v}_{\textnormal{L}}$ | $\mathbf{v}_{\textnormal{U}}$
1 | $(10{,}10{,}5{,}0)$ | $(10{,}10{,}5{,}0)$ | $(0{,}10{,}5{,}0)$ | $(10{,}10{,}5{,}0)$ | $C$ | $A$ | ✓ | $(10{,}10{,}5{,}0)$ | $(10{,}10{,}5{,}0)$
2 | $(10{,}10{,}5{,}0)$ | $(10{,}10{,}5{,}0)$ | $(10{,}0{,}5{,}0)$ | $(10{,}10{,}5{,}0)$ | $E$ | $A$ | ✓ | $(10{,}10{,}5{,}0)$ | $(10{,}10{,}5{,}0)$
3 | $(10{,}10{,}5{,}0)$ | $(10{,}10{,}5{,}0)$ | $(10{,}10{,}0{,}0)$ | $(10{,}10{,}5{,}0)$ | $A$ | $A$ | ✗ | $(10{,}10{,}0{,}0)$ | $(10{,}10{,}10{,}0)$
4 | $(10{,}10{,}0{,}0)$ | $(10{,}10{,}10{,}0)$ | $(10{,}10{,}0{,}0)$ | $(10{,}10{,}10{,}10)$ | $A$ | $A$ | ✗ | $(10{,}10{,}0{,}0)$ | $(10{,}10{,}10{,}10)$

Table 1: Execution of algorithm for finding one AXp

#### Finding one AXp.

We detail below the main steps of algorithm $\mathsf{findAXp}$ (see Algorithm
1). (Lines 4 to 5 are used for enumerating explanations, and so we assume
${\mathcal{S}}=\emptyset$ for now.) The main goal of $\mathsf{findAXp}$ is to
find a _maximal_ set of features ${\mathcal{D}}$ which are allowed to take
_any_ value, i.e. that are _free_. For such a set ${\mathcal{D}}$, the set of
features that remain fixed to the value specified in $\mathbf{v}$, i.e.
${\mathcal{P}}={\mathcal{F}}\setminus{\mathcal{D}}$, is a minimal set of
(picked) features that is sufficient for the prediction, as intended. The
different sets used by the algorithm are initialized in 3. (As noted earlier,
the sets ${\mathcal{C}}$, ${\mathcal{D}}$ and ${\mathcal{P}}$ form a partition
of ${\mathcal{F}}$, and ${\mathcal{C}}=\emptyset$ upon termination.)

For $\mathsf{findAXp}$, the second invariant of the algorithm is that
$\kappa(\mathbf{v}_{\textnormal{L}})=\kappa(\mathbf{v}_{\textnormal{U}})$,
i.e. by allowing the features in ${\mathcal{P}}\cup{\mathcal{C}}$ to take the
corresponding value in $\mathbf{v}$, the value of the prediction is guaranteed
not to change.

The use of the second invariant
$\kappa(\mathbf{v}_{\textnormal{L}})=\kappa(\mathbf{v}_{\textnormal{U}})$ is
justified by the following result.

######  Proposition 1.

If $\kappa(\mathbf{v}_{\textnormal{L}})=\kappa(\mathbf{v}_{\textnormal{U}})$,
then it holds that,
$\forall(\mathbf{x}\in\mathbb{F}).[\mathbf{v}_{\textnormal{L}}\leq\mathbf{x}\leq\mathbf{v}_{\textnormal{U}}]\rightarrow[\kappa(\mathbf{x})=\kappa(\mathbf{v})]$.

The algorithm starts by enforcing the second invariant as the result of
executing lines 1 and 2.

Moreover, $\mathsf{findAXp}$ analyzes one feature at a time. Starting from the
set ${\mathcal{C}}$ of candidate features (in 7), the algorithm iteratively
picks a feature $i$ from ${\mathcal{C}}$ and makes a decision about whether to
drop the feature from the explanation. The first step is to assume that the
feature $i$ can indeed be allowed to take any value. This is done in 8, by
calling the following function $\mathsf{FreeAttr}$:

$\mathbf{v}_{\textnormal{L}}\leftarrow(v_{\textnormal{L}_{1}},\ldots,\lambda(i),\ldots,v_{\textnormal{L}_{N}})$

$\mathbf{v}_{\textnormal{U}}\leftarrow(v_{\textnormal{U}_{1}},\ldots,\mu(i),\ldots,v_{\textnormal{U}_{N}})$

$({\mathcal{A}},{\mathcal{B}})\leftarrow({\mathcal{A}}\setminus\\{i\\},{\mathcal{B}}\cup\\{i\\})$

return
$(\mathbf{v}_{\textnormal{L}},\mathbf{v}_{\textnormal{U}},{\mathcal{A}},{\mathcal{B}})$

where ${\mathcal{A}}$ is replaced by ${\mathcal{C}}$ and $B$ is replaced by
${\mathcal{D}}$, and so feature $i$ is moved from ${\mathcal{C}}$ to
${\mathcal{D}}$. In addition, the value of $i$ is now allowed to range from
$\lambda(i)$ (in $\mathbf{v}_{\textnormal{L}}$) to $\mu(i)$ (in
$\mathbf{v}_{\textnormal{U}}$),

The next step of the algorithm (in 9) is to decide whether allowing $i$ to
take any value breaks the invariant
$\kappa(\mathbf{v}_{\textnormal{L}})=\kappa(\mathbf{v}_{\textnormal{U}})$. If
the invariant is not broken, then the algorithm moves to analyze the next
feature (in 7). However, if the invariant is broken, then the the feature
cannot take any value, and so it must be fixed to the corresponding value in
$\mathbf{v}$. This is done by calling (in 10) the following function
$\mathsf{FixAttr}$:

$\mathbf{v}_{\textnormal{L}}\leftarrow(v_{\textnormal{L}_{1}},\ldots,v_{i},\ldots,v_{\textnormal{L}_{N}})$

$\mathbf{v}_{\textnormal{U}}\leftarrow(v_{\textnormal{U}_{1}},\ldots,v_{i},\ldots,v_{\textnormal{U}_{N}})$

$({\mathcal{A}},{\mathcal{B}})\leftarrow({\mathcal{A}}\setminus\\{i\\},{\mathcal{B}}\cup\\{i\\})$

return
$(\mathbf{v}_{\textnormal{L}},\mathbf{v}_{\textnormal{U}},{\mathcal{A}},{\mathcal{B}})$

where ${\mathcal{A}}$ is replaced by ${\mathcal{D}}$ and ${\mathcal{B}}$ is
replaced by ${\mathcal{P}}$, and so feature $i$ is moved from ${\mathcal{D}}$
to ${\mathcal{P}}$. In addition, the value of $i$ is once again fixed to the
corresponding value in $\mathbf{v}$. After analyzing all features, the
algorithm $\mathsf{findAXp}$ terminates (in 13) by return the (minimal) set of
features ${\mathcal{P}}$ that are fixed to their value in $\mathbf{v}$. It is
immediate to conclude that each feature is analyzed once, and that for each
feature, the classifier is invoked twice. Given the discussion above, we
conclude that,

######  Theorem 1.

Given a monotonic classifier, an instance $\mathbf{v}$ with prediction
$c=\kappa(\mathbf{v})$, Algorithm 1 computes one AXp in linear time in the
running time complexity of the classifier.

We illustrate the operation of $\mathsf{findAXp}$, with an example.

######  Example 3.

Given the monotonic classifier from Example 1, and the concrete case of
student $s_{1}$, with $(Q,X,H,R)=(10,10,5,0)$ and predicted mark $A$, we show
how one PI-explanation can computed. (In settings with more than one AXp,
changing the order of how features are analyzed, may results in a different
explanation being obtained.) For each feature $i$, $1\leq{i}\leq 4$,
$\lambda(i)=0$ and $\mu(i)=10$. Moreover, features are analyzed in order:
$\langle 1,2,3,4\rangle$; the order is arbitrary. The algorithm’s execution is
summarized in Table 1. As can be observed, features 1 and 2 are kept as part
of the PI-explanation (decision is ✓in 9, i.e. invariant is broken and
features are kept), whereas features 3 and 4 are dropped from the PI-
explanation (decision is ✗, i.e. invariant holds). As a result, the PI-
explanation for the grade of student $s_{1}$ is $\\{1,2\\}$, which denotes
that as long as $(Q=10)\land(X=10)$, the prediction will be $A$.

Feat. | Initial values | Changed values | Predictions | Dec. | Resulting values
---|---|---|---|---|---
$\mathbf{v}_{\textnormal{L}}$ | $\mathbf{v}_{\textnormal{U}}$ | $\mathbf{v}_{\textnormal{L}}$ | $\mathbf{v}_{\textnormal{U}}$ | $\kappa(\mathbf{v}_{\textnormal{L}})$ | $\kappa(\mathbf{v}_{\textnormal{U}})$ | $\mathbf{v}_{\textnormal{L}}$ | $\mathbf{v}_{\textnormal{U}}$
1 | $(0{,}0{,}0{,}0)$ | $(10{,}10{,}10{,}10)$ | $(10{,}0{,}0{,}0)$ | $(10{,}10{,}10{,}10)$ | $E$ | $A$ | ✗ | $(10{,}0{,}0{,}0)$ | $(10{,}10{,}10{,}10)$
2 | $(10{,}0{,}0{,}0)$ | $(10{,}10{,}10{,}10)$ | $(10{,}10{,}0{,}0)$ | $(10{,}10{,}10{,}10)$ | $A$ | $A$ | ✓ | $(10{,}0{,}10{,}0)$ | $(10{,}10{,}10{,}10)$
3 | $(10{,}0{,}0{,}0)$ | $(10{,}10{,}10{,}10)$ | $(10{,}0{,}5{,}0)$ | $(10{,}10{,}5{,}10)$ | $E$ | $A$ | ✗ | $(10{,}0{,}5{,}0)$ | $(10{,}0{,}5{,}10)$
4 | $(10{,}0{,}5{,}0)$ | $(10{,}10{,}5{,}10)$ | $(10{,}0{,}5{,}0)$ | $(10{,}10{,}5{,}0)$ | $E$ | $A$ | ✗ | $(10{,}0{,}5{,}0)$ | $(10{,}10{,}5{,}0)$

Table 2: Execution of algorithm for finding one CXp

#### Finding one CXp.

The two algorithms $\mathsf{findAXp}$ and $\mathsf{findCXp}$ are organized in
a similar way. (This in part results from the fact that AXps are minimal
hitting sets of CXps and vice-versa (Ignatiev et al., 2020).) We briefly
explain the differences when computing a CXp (see Algorithm 2). (Lines 4 to 5
are used for enumerating explanations, and so we assume
${\mathcal{S}}=\emptyset$ for now.)

The main goal of $\mathsf{findCXp}$ is to find a _maximal_ set of features
${\mathcal{D}}$ that are only allowed to take the value specified in
$\mathbf{v}$, i.e. that are _fixed_. For such a set ${\mathcal{D}}$, the set
of features that are allowed to take any value, i.e.
${\mathcal{P}}={\mathcal{F}}\setminus{\mathcal{D}}$, is a minimal set that, by
being allowed to take any value in their domain, suffices for allowing the
prediction to change, as intended. The different sets used by the algorithm
are initialized in 3.

For $\mathsf{findCXp}$, the second invariant of the algorithm is that
$\kappa(\mathbf{v}_{\textnormal{L}})\not=\kappa(\mathbf{v}_{\textnormal{U}})$,
i.e. by allowing the features in ${\mathcal{P}}\cup{\mathcal{C}}$ to take any
value, the value of the prediction does not change. The algorithm starts by
enforcing the second invariant as the result of executing lines 1 and 2.

The use of the second invariant
$\kappa(\mathbf{v}_{\textnormal{L}})\not=\kappa(\mathbf{v}_{\textnormal{U}})$
is justified by the following result.

######  Proposition 2.

If
$\kappa(\mathbf{v}_{\textnormal{L}})\not=\kappa(\mathbf{v}_{\textnormal{U}})$,
then it holds that,
$\exists(\mathbf{x}\in\mathbb{F}).[\mathbf{v}_{\textnormal{L}}\leq\mathbf{x}\leq\mathbf{v}_{\textnormal{U}}]\land[\kappa(\mathbf{x})\not=\kappa(\mathbf{v})]$.

Similarly to $\mathsf{findAXp}$, $\mathsf{findCXp}$ analyzes one feature at a
time. Starting from the set ${\mathcal{C}}$ of candidate features (in 7), the
algorithm iteratively picks a feature $i$ from ${\mathcal{C}}$ and makes a
decision about whether to drop the feature from the explanation. The first
step is to assume that the feature $i$ can indeed be fixed to the
corresponding value in $\mathbf{v}$. This is done in 8, by calling the
following function $\mathsf{FixAttr}$, where ${\mathcal{A}}$ is replaced by
${\mathcal{C}}$, and $B$ is replaced by ${\mathcal{D}}$, and so feature $i$ is
moved from ${\mathcal{C}}$ to ${\mathcal{D}}$. In addition, the value of $i$
is now fixed to its value in $\mathbf{v}$.

The next step of the algorithm (in 9) is to decide whether fixing the value of
$i$ breaks the invariant
$\kappa(\mathbf{v}_{\textnormal{L}})\not=\kappa(\mathbf{v}_{\textnormal{U}})$.
If the invariant is not broken, then the algorithm moves to analyze the next
feature (in 7). However, if the invariant is broken, then the feature cannot
be fixed, and so it must be allowed to take any value from its domain. This is
done by calling (in 10) the following function $\mathsf{FreeAttr}$, with
${\mathcal{A}}$ replaced by ${\mathcal{D}}$ and ${\mathcal{B}}$ replaced by
${\mathcal{P}}$, and so feature $i$ is moved from ${\mathcal{D}}$ to
${\mathcal{P}}$. In addition, the value of $i$ is once again allowed to take
any value from its domain. After analyzing all features, the algorithm
$\mathsf{findCXp}$ terminates (in 13) by returning the (minimal) set of
features ${\mathcal{P}}$ that are allowed to take any value from their domain.
It is immediate to conclude that each feature is analyzed once, and that for
each feature, the classifier is invoked twice. Given the discussion above, we
conclude that,

######  Theorem 2.

Given a monotonic classifier, an instance $\mathbf{v}$ with prediction
$c=\kappa(\mathbf{v})$, Algorithm 2 computes one CXp in linear time in the
running time complexity of the classifier.

We illustrate the operation of $\mathsf{findCXp}$, with an example.

######  Example 4.

For the running example (see Examples 1, 2 and 3), for instance
$\mathbf{v}_{0}=(10,10,5,0)$ with prediction $A$, we illustrate the
computation of one CXp. The algorithm’s execution is summarized in Table 2.
(When computing one CXp, a feature is _kept_ (decision is ✓) if it is declared
free, and it is _dropped_ (decision is ✗) if it must be fixed.) As can be
observed, a contrastive explanation is: $\\{2\\}$, i.e. there is an assignment
to feature 2 (i.e. to $X$), which guarantees a change of prediction when the
other features are kept to their values. For example, by setting $X=0$ (and
keeping the remaining values fixed), the value of the prediction changes.

#### Complexity.

As can be readily concluded from Algorithm 1 and Algorithm 2, the algorithms
execute in linear time in the number of features. However, in each iteration
of the algorithm, the classifier is invoked twice, for finding the predicted
classes for $\mathbf{v}_{\textnormal{L}}$ and for
$\mathbf{v}_{\textnormal{U}}$. We will represent the time required by the
classifier as ${\mathcal{T}}_{\mathbb{C}}$, and so the overall run time of
each algorithm is
${\mathcal{O}}(|{\mathcal{F}}|\times{\mathcal{T}}_{\mathbb{C}})$.

###  3.2 Enumerating Explanations

We first show that for monotonic classifiers, the enumeration of explanations
with polynomial-time delay is computationally hard.

######  Theorem 3.

Determining the existence of $\lfloor N/2\rfloor{+}1$ AXp’s (or CXp’s) of a
monotonic $N$-feature classifier is NP-complete.

(The proof is included in the appendix.) Since the enumeration of AXp’s and
CXp’s with polynomial delay is unlikely, we describe in this section how to
use SAT reasoners for the enumeration of AXp’s and CXp’s of a monotonic
classifier. (Although we prove the algorithm to be sound and complete, the
algorithm necessarily has leeway in selecting the order in which AXp’s and
CXp’s are listed.)

The algorithm uses the following propositional representation:

  1. 1.

The algorithm will iteratively add clauses to a CNF formula ${\mathcal{H}}$.
The clauses in ${\mathcal{H}}$ account for the AXp’s and CXp’s already
computed, and serve to prevent their repetition.

  2. 2.

Formula ${\mathcal{H}}$ is defined on a set of variables $u_{i}$,
$1\leq{i}\leq{n}$, where each $u_{i}$ denotes whether feature $i$ is declared
free ($u_{i}=1$) or is alternatively declared fixed ($u_{i}=0$).

The algorithm proposed in this section requires exactly one call to a SAT
reasoner before computing one explanation (either AXp/CXp), and one additional
call to decide that all explanations have been computed. As a result, the
number of calls to a SAT reasoner is
$|\textnormal{AXp}|+|\textnormal{CXp}|+1$. Furthermore, the size of the
formula grows by one clause after each AXp or CXp is computed. In practice,
for a wide range of ML settings, both the number of variables and the number
of clauses are well within the reach of modern SAT reasoners.

######  Proposition 3.

Let $\mathbf{v}$ be a point in feature space, let
$\kappa(\mathbf{v})=c\in{\mathcal{K}}$, and let
${\mathcal{Z}}\subseteq{\mathcal{F}}$. Then, either (1) (on page 1) holds,
with ${\mathcal{X}}={\mathcal{Z}}$, or (2) (also on page 1) holds, with
${\mathcal{Y}}={\mathcal{F}}\setminus{\mathcal{Z}}$, but not both.

Proposition 3 essentially states that, given a set ${\mathcal{Z}}$ of
features, if these are fixed, and the others are allowed to take any value
from their domains, then either the prediction never changes, or there exists
an assignment to the non-fixed features, which causes the prediction to
change. The approach for enumerating AXp’s and CXp’s is shown in Algorithm 3.

Algorithm 3 Enumeration of AXp’s and CXp’s

1: ${\mathcal{H}}\leftarrow\emptyset$// ${\mathcal{H}}$ defined on set $U$

2: repeat

3: $(\mathsf{outc},\mathbf{u})\leftarrow\mathsf{SAT}({\mathcal{H}})$

4: if $\mathsf{outc}=\textbf{true}{}$ then

5:
$\mathbf{v}_{\textnormal{L}}\leftarrow(v_{\textnormal{L}_{1}},\ldots,v_{\textnormal{L}_{N}}),\>\textnormal{s.t.}\>v_{\textnormal{L}_{i}}\leftarrow\textnormal{ite}(u_{i},\lambda(i),v_{i})$

6:
$\mathbf{v}_{\textnormal{U}}\leftarrow(v_{\textnormal{U}_{1}},\ldots,v_{\textnormal{U}_{N}}),\>\\!\textnormal{s.t.}\;\\!v_{\textnormal{U}_{i}}\leftarrow\textnormal{ite}(u_{i},\mu(i),v_{i})$

7: if
$\kappa(\mathbf{v}_{\textnormal{L}})=\kappa(\mathbf{v}_{\textnormal{U}})$ then

8: ${\mathcal{S}}\leftarrow\\{i\in{\mathcal{F}}\,|\,u_{i}=1\\}$//
${\mathcal{F}}\setminus{\mathcal{S}}\supseteq$ some AXp

9:
${\mathcal{P}}\leftarrow\mathsf{findAXp}({\mathcal{F}},{\mathcal{S}},\mathbf{u})$

10: $\mathsf{reportAXp}({\mathcal{P}})$

11:
${\mathcal{H}}\leftarrow{\mathcal{H}}\cup\\{(\lor_{i\in{\mathcal{P}}}u_{i})\\}$

12: else

13: ${\mathcal{S}}\leftarrow\\{i\in{\mathcal{F}}\,|\,u_{i}=0\\}$//
${\mathcal{F}}\setminus{\mathcal{S}}\supseteq$ some CXp

14:
${\mathcal{P}}\leftarrow\mathsf{findCXp}({\mathcal{F}},{\mathcal{S}},\mathbf{u})$

15: $\mathsf{reportCXp}({\mathcal{P}})$

16:
${\mathcal{H}}\leftarrow{\mathcal{H}}\cup\\{(\lor_{i\in{\mathcal{P}}}\neg{u_{i}})\\}$

17: end if

18: end if

19: until $\mathsf{outc}=\textbf{false}{}$

The algorithm starts in 1 by initializing the CNF formula ${\mathcal{H}}$
without clauses (these will be added as the algorithm executes). The main loop
(from 2 to 19) is executed while the formula ${\mathcal{H}}$ is satisfiable.
This is decided with a call to a SAT reasoner (in 3). Any satisfying
assignment to the formula ${\mathcal{H}}$ partitions the features into two
sets: one denoting the features that can take any value (with $u_{i}=1$) and
another denoting the features that take the corresponding value in
$\mathbf{v}$ (with $u_{i}=0$). (The assignment effectively identifies a set
${\mathcal{Z}}\subseteq{\mathcal{F}}$, of fixed features, and thus we can
invoke Proposition 3.) In 5 and 6, the algorithm creates
$\mathbf{v}_{\textnormal{L}}$ and $\mathbf{v}_{\textnormal{U}}$. For a fixed
feature $i$, both $\mathbf{v}_{\textnormal{L}}$ and
$\mathbf{v}_{\textnormal{U}}$ are assigned value $v_{i}$. For a free feature
$i$, $\mathbf{v}_{\textnormal{L}}$ and $\mathbf{v}_{\textnormal{U}}$ are
respectively assigned to $\lambda(i)$ and $\mu(i)$. Let ${\mathcal{Z}}$ denote
the set of fixed features. In 7, we check in which case of Proposition 3
applies.

If $\kappa(\mathbf{v}_{\textnormal{L}})=\kappa(\mathbf{v}_{\textnormal{U}})$,
then we know that the invariant of Algorithm 1 holds. Moreover,
${\mathcal{F}}\setminus{\mathcal{Z}}$ is a subset of an AXp. Hence, we set
${\mathcal{S}}={\mathcal{F}}\setminus{\mathcal{Z}}$ as the _seed_ for
$\mathsf{findAXp}$. This is shown in lines 8 and 9. After reporting the
computed AXp, represented by the set of features ${\mathcal{P}}$, we prevent
the same AXp from being computed again by requiring that at least one of the
fixed features must be free in future satisfying assignments of
${\mathcal{H}}$. This is represented as a positive clause.

######  Proposition 4.

In the case
$\kappa(\mathbf{v}_{\textnormal{L}})=\kappa(\mathbf{v}_{\textnormal{U}})$, set
${\mathcal{S}}$ is such that, for any previously computed AXp, at least one
feature will be included in ${\mathcal{S}}$ (as a free literal). Since
$\mathsf{findAXp}$ only grows ${\mathcal{S}}$, then the Algorithm 3 does not
repeat AXp’s.

Moreover, if
$\kappa(\mathbf{v}_{\textnormal{L}})\not=\kappa(\mathbf{v}_{\textnormal{U}})$,
then we know that the invariant of Algorithm 2 holds. Moreover,
${\mathcal{Z}}$ is a subset of a CXp. Hence, we set
${\mathcal{S}}={\mathcal{Z}}$ as the _seed_ for $\mathsf{findCXp}$. This is
shown in lines 13 and 14. After reporting the computed CXp, represented by the
set of features ${\mathcal{P}}$, we prevent the same CXp from being computed
by requiring that at least one of the free features must be free in future
satisfying assignments of ${\mathcal{H}}$. This is represented as negative
clause.

######  Proposition 5.

In the case
$\kappa(\mathbf{v}_{\textnormal{L}})\not=\kappa(\mathbf{v}_{\textnormal{U}})$,
set ${\mathcal{S}}$ is such that, for any previously computed CXp, at least
one feature will be included in ${\mathcal{S}}$ (as a fixed literal). Since
$\mathsf{findCXp}$ only grows ${\mathcal{S}}$, then the Algorithm 3 does not
repeat CXp’s.

Given the above, and since the number of AXp’s and CXp’s (being subsets of
${\mathcal{F}}$) is finite, then we have,

######  Theorem 4.

Algorithm 3 is sound and complete for the enumeration of AXp’s and CXp’s.

Step | ${\mathcal{H}}$ | $\mathbf{u}~{}/~{}\mathsf{out}$ | $\mathbf{v}_{\textnormal{L}}$ | $\mathbf{v}_{\textnormal{U}}$ | $\kappa(\mathbf{v}_{\textnormal{L}})$ | $\kappa(\mathbf{v}_{\textnormal{U}})$ | AXp | CXp | Clause added
---|---|---|---|---|---|---|---|---|---
1 | $\emptyset$ | $(0,0,1,0)$ | $(10,10,0,0)$ | $(10,10,10,0)$ | $A$ | $A$ | $\\{1,2\\}$ | – | $\begin{array}[]{l}(u_{1}\lor{u_{2}})\end{array}$
2 | $\begin{array}[]{l}(u_{1}\lor{u_{2}})\end{array}$ | $(1,0,0,1)$ | $(0,10,5,0)$ | $(10,10,5,10)$ | $C$ | $A$ | – | $\\{1\\}$ | $\begin{array}[]{l}(\neg{u_{1}})\end{array}$
3 | $\begin{array}[]{l}(u_{1}\lor{u_{2}})\\\ (\neg{u_{1}})\end{array}$ | $(0,1,1,0)$ | $(10,0,0,0)$ | $(10,10,10,0)$ | $E$ | $A$ | – | $\\{2\\}$ | $\begin{array}[]{l}(\neg{u_{2}})\end{array}$
4 | $\begin{array}[]{l}(u_{1}\lor{u_{2}})\\\ (\neg{u_{1}}),(\neg{u_{2}})\end{array}$ | UNSAT | – | – | – | – | – | – | –

Table 3: Execution of enumeration algorithm

######  Example 5.

Building on earlier examples, we summarize the main steps of the SAT oracle-
based algorithm for enumerating AXp and CXp explanations. Table 3 illustrates
one execution of the proposed algorithm. There are 1 AXp’s and 2 CXp’s.
(Regarding the call to the SAT oracle, the satisfying assignments shown are
intended to be as arbitrary as possible, given the existing constraints; other
satisfying assignments could have been picked.) For each computed AXp, we add
to ${\mathcal{H}}$ one positive clause. In this example, we add the clause
$(u_{1}\lor{u_{2}})$, since the AXp is $\\{1,2\\}$. By adding this clause, we
guarantee that features 1 and 2 will not both be deemed fixed by subsequent
satisfying assignments of ${\mathcal{H}}$. Similarly, for each computed CXp,
we add to ${\mathcal{H}}$ one negative clause. For the example, the clauses
added are $(\neg{u_{1}})$ for CXp $\\{1\\}$, and $(\neg{u_{2}})$ for CXp
$\\{2\\}$. In both cases, the added clause guarantees that feature 1 (resp. 2)
will not be deemed free by subsequent satisfying assignments of
${\mathcal{H}}$. One additional observation is that the number of SAT oracle
calls matches the number of AXp’s plus the number of CXp’s and plus one final
call to terminate the algorithm’s execution. For step 4 of the algorithm, it
is easy to conclude that ${\mathcal{H}}$ is unsatisfiable, as intended.

###  3.3 Related Work

The algorithms for computing one AXp or one CXp for a monotonic classifier are
novel. However, the insight of analyzing elements of a set (i.e. features in
our case) to find a minimal set respecting some property has been studied in a
vast number of settings (e.g. (Chinneck, 2008) for an overview). The proposed
solution for reasoning about ordinal features that can take boolean, integer
or real values, represents another aspect of novelty. In the case of monotonic
classifiers, we obtain a running time that is linear in the running time
complexity of the classifier. This result applies in the case of _any_
monotonic classifier, and so it improves significantly over the worst-case
exponential time and space approach proposed in earlier work (Shih et al.,
2018), for the concrete case of monotonic bayesian networks. The algorithm for
enumerating AXp’s and CXp’s for a monotonic classifier is also novel. However,
it is inspired by the MARCO algorithm for the analysis of inconsistent logic
theories (Liffiton et al., 2016). Although MARCO can be optimized in different
ways, Algorithm 3 can be related with its most basic formulation. Since
computing one AXp or one CXp can be achieved in polynomial time (conditioned
by the classifier run time complexity), then our approach guarantees that
exactly one SAT reasoner call is required for each computed minimal set (i.e.
AXp or CXp in our case).

##  4 Experiments

The objective of this section is to illustrate the scalability of both the
algorithms for finding one explanation, but also the algorithm for enumerating
explanations. The tool MonoXP implements the algorithms 1, 2 and 3. As
observed in recent work, most monotonic classifiers are not publicly available
(Cano et al., 2019). We analyze two publicly available classifiers, and
describe two experiments. The first experiment evaluates MonoXP for explaining
two recently proposed tools, COMET777Available from
https://github.com/AishwaryaSivaraman/COMET. (Sivaraman et al., 2020) and
monoboost888Available from
https://github.com/chriswbartley/monoboost. (Bartley et al., 2018). COMET is
run on the Auto-MPG999https://archive.ics.uci.edu/ml/datasets/auto+mpg.
dataset studied in earlier work (Sivaraman et al., 2020), with the choice
justified by the time the classifier takes to run. monoboost is run on a
monotonic dataset with two classes (as required by the tool) (Bartley et al.,
2018). We use a monotonic subset (PimaMono) of the Pima
dataset101010https://sci2s.ugr.es/keel/dataset.php?cod=21.. A second
experiment compares MonoXP with Anchor (Ribeiro et al., 2018), both in terms
of the number of calls to the classifier and running time, but also in terms
of the quality of the computed explanations111111It should be underlined that
neither Anchor (Ribeiro et al., 2018), LIME (Ribeiro et al., 2016) nor SHAP
(Lundberg & Lee, 2017) can enumerate explanations, neither can these tools
compute heuristic contrastive explanations., namely accuracy and size. This
second experiment also considers two datasets. The first dataset is
BankruptcyRisk (Greco et al., 1998) (which is monotonic if one instance is
dropped). For this dataset, the monotonic decision tree classifier proposed in
earlier work is used (Potharst & Bioch, 2000). The second dataset is PimaMono,
and the classifier used is the one obtained with monoboost (as in the first
experiment). All experiments were run on a MacBook Pro, with a 2.4GHz quad-
core i5 processor, and 16 GByte of RAM, running MacOS Big Sur. For each
dataset, we either pick 100 instances, randomly selected, or the total number
of instances in the dataset, in case this number does not exceed 100.

Dataset/Tool | #Inst. | Avg. # expl. | Avg. AXp sz | Avg. CXp sz | Avg. classif. time | Avg. run time | % classif. time
---|---|---|---|---|---|---|---
AutoMPG/CMT | 100 | 2.35 | 1.49 | 1.02 | 105.90s | 105.92s | 99.99%
PimaMono/MBT | 69 | 9.09 | 1.27 | 3.36 | 16.285s | 16.360s | 99.54%
(a) Assessing MonoXP on the Auto-MPG and PimaMono datasets, using resp. COMET
or monoboost as the classifier

Dataset | #Inst. | Anchor |  MonoXP (AXp) | % diff
---|---|---|---|---
Avg. Xp sz | Avg. time | # Cls calls | Avg. # Xp | Avg. Xp sz | Avg. Xp time | # Cls calls
B. Risk | 39 | 2.18 | 0.11s | 1217 | 1.03 | 2.0 | 0.009s | 24 | 64.1
PimaMono | 69 | 1.26 | 11.2s | 2967 | 5.64 | 1.27 | 1.8s | 16 | 18.8
(b) Assessing MonoXP and Anchor on the Bankruptcy Risk and Pima Mono datasets

Table 4: Results of running MonoXP

###  4.1 Cost of Computing Explanations

We run MonoXP on a neural network classifier envelope implemented with COMET
for the Auto-MPG dataset, and on a tree ensemble obtained with monoboost for
the PimaMono dataset. (Since the running times of COMET can be significant,
this experiment does not consider a comparison with the heuristic explainer
Anchor (Ribeiro et al., 2018). As shown below, Anchor calls the classifier a
large number of times, and that would imply unwieldy running times.)

4(a) shows the results of running MonoXP using COMET as a monotonic envelope
on the Auto-MPG dataset, and monoboost on the PimaMono dataset. As can be
observed, the explanation sizes are in general small, which confirms the
interpretability of computed AXp’s and CXp’s. As a general trend, CXp’s are on
average smaller than AXp’s for Auto-MPG, but larger than AXp’s for PimaMono.
Moreover, the classification time completely dominates the total running time
(i.e. resp. 99.99% and 99.54% of the time is spent running the classifier,
independently of the classifier used). These results offer evidence that the
time spent on computing explanations is in general negligible for monotonic
classifiers. For both datasets, and for the instances considered, it was
possible to enumerate all AXp and CXp explanations, with negligible
computational overhead.

###  4.2 Comparison with Anchor

This section compares MonoXP with Anchor, using two pairs of classifiers and
datasets, i.e. a monotonic decision tree for BankruptcyRisk and monoboost for
PimaMono.

4(b) shows the results of running Anchor and MonoXP on the BankruptcyRisk and
the PimaMono datasets. MonoXP is significantly faster than Anchor (more than 1
order magnitude in the first case, and more than a factor of 5 in the second
case). The justification is the much smaller number of calls to the classifier
required by MonoXP than by Anchor. (While for Anchor the number of calls to
the classifier can be significant, for MonoXP, each AXp is computed with at
most a linear number of calls to the classifier. Thus, unless the number of
features is very substantial, MonoXP has a clear performance edge over
Anchor.) Somewhat surprisingly, over all instances, the average size of AXp’s
computed by MonoXP is smaller than that of Anchor for the BankruptcyRisk
dataset. For the PimaMono dataset, the average size is almost the same. These
results suggest that formally defined explanations need not be significantly
larger than the ones computed with heuristic approaches. Furthermore, for
64.1% (resp. 18.8%) of the instances, Anchor identifies an explanation that
does not hold across all points of feature space, i.e. there are points in
feature space for which the explanation of Anchor holds, but the ML model
makes a different prediction121212Similar observations have been reported
elsewhere (Ignatiev, 2020).. Observe that since MonoXP computes all AXp’s, we
can be certain about whether the explanation of Anchor is a correct
explanation.

##  5 Conclusions & Discussion

This paper proposes novel algorithms for computing a single PI or contrastive
explanation for a monotonic classifier. In contrast with earlier work (Shih et
al., 2018), the complexity of the proposed algorithms is polynomial on the
number of features and the time it takes the monotonic classifier to compute
its predictions. As the experiments demonstrate, for simple ML models, the
algorithm achieves one order of magnitude speed up when compared with a well-
known heuristic explainer (Ribeiro et al., 2018), achieving better quality
explanations of similar size. In contrast, for complex ML models, the
experiments confirm that the running time is almost entirely spent on the
classifier. Furthermore, the paper proposes a SAT-based algorithm for
enumerating PI and contrastive explanations. As the experimental results show,
the use of a SAT solver for enumerating PI and contrastive explanations incurs
a negligible overhead.

One possible criticism of the work is that SAT solvers are used for guiding
the enumeration of explanations. This involves solving an NP-complete decision
problem for each computed explanation, and so it might pose a scalability
concern. (One alternative would be to consider explicit enumeration of
candidate explanations, as proposed in the earlier works on model based
diagnosis (Reiter, 1987; Greiner et al., 1989; Wotawa, 2001).) However, for
classification problems with tens to hundreds of features and targeting
thousands to tens of thousands explanations (and this far exceeds currently
foreseen scenarios), the use of modern SAT reasoners (capable of solving
problems with hundreds of thousands of variable and millions of clauses) can
hardly be considered a limitation. Another possible criticism of this work is
that full monotonicity is required. We conjecture that _full_ monotonicity is
necessary for tractable explanations (conditioned by the classifier run time
complexity). Addressing partial monotonicity (Daniels & Velikova, 2010) is a
subject of future research.

#### Acknowledgments.

This work was supported by the AI Interdisciplinary Institute ANITI, funded by
the French program “Investing for the Future – PIA3” under Grant agreement no.
ANR-19-PI3A-0004, and by the H2020-ICT38 project COALA “Cognitive Assisted
agile manufacturing for a Labor force supported by trustworthy Artificial
intelligence”.

## References

  * Audemard et al. (2020) Audemard, G., Koriche, F., and Marquis, P.  On tractable XAI queries based on compiled representations.  In _KR_ , pp. 838–849, 2020.
  * Babin & Kuznetsov (2011) Babin, M. A. and Kuznetsov, S. O.  Enumerating minimal hypotheses and dualizing monotone boolean functions on lattices.  In _FCA_ , pp. 42–48, 2011.
  * Barile & Feelders (2012) Barile, N. and Feelders, A.  Active learning with monotonicity constraints.  In _SIAM ICDM_ , pp. 756–767, 2012.
  * Bartley et al. (2016) Bartley, C., Liu, W., and Reynolds, M.  Effective monotone knowledge integration in kernel support vector machines.  In _ADMA_ , pp. 3–18, 2016.
  * Bartley et al. (2018) Bartley, C., Liu, W., and Reynolds, M.  A novel framework for constructing partially monotone rule ensembles.  In _ICDE_ , pp. 1320–1323, 2018.
  * Bartley et al. (2019) Bartley, C., Liu, W., and Reynolds, M.  Enhanced random forest algorithms for partially monotone ordinal classification.  In _AAAI_ , pp. 3224–3231, 2019.
  * Ben-David (1995) Ben-David, A.  Monotonicity maintenance in information-theoretic machine learning algorithms.  _Mach. Learn._ , 19(1):29–43, 1995.
  * Ben-David et al. (1989) Ben-David, A., Sterling, L., and Pao, Y.  Learning, classification of monotonic ordinal concepts.  _Comput. Intell._ , 5:45–49, 1989.
  * Biere et al. (2009) Biere, A., Heule, M., van Maaren, H., and Walsh, T. (eds.).  _Handbook of Satisfiability_ , 2009. IOS Press.
  * Bonakdarpour et al. (2018) Bonakdarpour, M., Chatterjee, S., Barber, R. F., and Lafferty, J.  Prediction rule reshaping.  In _ICML_ , pp. 629–637, 2018.
  * Cano et al. (2019) Cano, J. R., Gutiérrez, P. A., Krawczyk, B., Wozniak, M., and García, S.  Monotonic classification: An overview on algorithms, performance measures and data sets.  _Neurocomputing_ , 341:168–182, 2019.
  * Chinneck (2008) Chinneck, J. W.  _Feasibility and Infeasibility in Optimization:: Algorithms and Computational Methods_.  Springer Science & Business Media, 2008.
  * Daniels & Velikova (2010) Daniels, H. and Velikova, M.  Monotone and partially monotone neural networks.  _IEEE Trans. Neural Networks_ , 21(6):906–917, 2010.
  * Duivesteijn & Feelders (2008) Duivesteijn, W. and Feelders, A.  Nearest neighbour classification with monotonicity constraints.  In _ECML/PKDD_ , pp. 301–316, 2008.
  * Fard et al. (2016) Fard, M. M., Canini, K. R., Cotter, A., Pfeifer, J., and Gupta, M. R.  Fast and flexible monotonic functions with ensembles of lattices.  In _NeurIPS_ , pp. 2919–2927, 2016.
  * Greco et al. (1998) Greco, S., Matarazzo, B., and Slowinski, R.  A new rough set approach to evaluation of bankruptcy risk.  In _Operational tools in the management of financial risks_ , pp. 121–136. Springer, 1998.
  * Greiner et al. (1989) Greiner, R., Smith, B. A., and Wilkerson, R. W.  A correction to the algorithm in Reiter’s theory of diagnosis.  _Artif. Intell._ , 41(1):79–88, 1989.
  * Gupta et al. (2016) Gupta, M. R., Cotter, A., Pfeifer, J., Voevodski, K., Canini, K. R., Mangylov, A., Moczydlowski, W., and Esbroeck, A. V.  Monotonic calibrated interpolated look-up tables.  _J. Mach. Learn. Res._ , 17:109:1–109:47, 2016.
  * Ignatiev (2020) Ignatiev, A.  Towards trustable explainable AI.  In _IJCAI_ , pp. 5154–5158, 2020.
  * Ignatiev et al. (2019) Ignatiev, A., Narodytska, N., and Marques-Silva, J.  Abduction-based explanations for machine learning models.  In _AAAI_ , pp. 1511–1519, 2019.
  * Ignatiev et al. (2020) Ignatiev, A., Narodytska, N., Asher, N., and Marques-Silva, J.  On relating ’why?’ and ’why not?’ explanations.  _CoRR_ , abs/2012.11067, 2020.  URL https://arxiv.org/abs/2012.11067.
  * Liffiton et al. (2016) Liffiton, M. H., Previti, A., Malik, A., and Marques-Silva, J.  Fast, flexible MUS enumeration.  _Constraints An Int. J._ , 21(2):223–250, 2016\.
  * Liu et al. (2020) Liu, X., Han, X., Zhang, N., and Liu, Q.  Certified monotonic neural networks.  In _NeurIPS_ , 2020.
  * Lundberg & Lee (2017) Lundberg, S. M. and Lee, S.  A unified approach to interpreting model predictions.  In _NeurIPS_ , pp. 4765–4774, 2017.
  * Magdon-Ismail & Sill (2008) Magdon-Ismail, M. and Sill, J.  A linear fit gets the correct monotonicity directions.  _Mach. Learn._ , 70(1):21–43, 2008.
  * Marques-Silva et al. (2020) Marques-Silva, J., Gerspacher, T., Cooper, M. C., Ignatiev, A., and Narodytska, N.  Explaining naive bayes and other linear classifiers with polynomial time and delay.  In _NeurIPS_ , 2020.
  * Miller (2019) Miller, T.  Explanation in artificial intelligence: Insights from the social sciences.  _Artif. Intell._ , 267:1–38, 2019.
  * Potharst & Bioch (2000) Potharst, R. and Bioch, J. C.  Decision trees for ordinal classification.  _Intell. Data Anal._ , 4(2):97–111, 2000.
  * Reiter (1987) Reiter, R.  A theory of diagnosis from first principles.  _Artif. Intell._ , 32(1):57–95, 1987.
  * Ribeiro et al. (2016) Ribeiro, M. T., Singh, S., and Guestrin, C.  ”Why should I trust you?”: Explaining the predictions of any classifier.  In _KDD_ , pp. 1135–1144, 2016.
  * Ribeiro et al. (2018) Ribeiro, M. T., Singh, S., and Guestrin, C.  Anchors: High-precision model-agnostic explanations.  In _AAAI_ , pp. 1527–1535, 2018.
  * Shih et al. (2018) Shih, A., Choi, A., and Darwiche, A.  A symbolic approach to explaining bayesian network classifiers.  In _IJCAI_ , pp. 5103–5111, 2018.
  * Sill (1997) Sill, J.  Monotonic networks.  In _NIPS_ , pp. 661–667, 1997.
  * Sivaraman et al. (2020) Sivaraman, A., Farnadi, G., Millstein, T. D., and den Broeck, G. V.  Counterexample-guided learning of monotonic neural networks.  In _NeurIPS_ , 2020.
  * Van den Broeck et al. (2020) Van den Broeck, G., Lykov, A., Schleich, M., and Suciu, D.  On the tractability of SHAP explanations.  _CoRR_ , abs/2009.08634, 2020.  URL https://arxiv.org/abs/2009.08634.
  * van der Gaag et al. (2004) van der Gaag, L. C., Bodlaender, H. L., and Feelders, A. J.  Monotonicity in bayesian networks.  In _UAI_ , pp. 569–576, 2004.
  * Verbeke et al. (2017) Verbeke, W., Martens, D., and Baesens, B.  RULEM: A novel heuristic rule learning approach for ordinal classification with monotonicity constraints.  _Appl. Soft Comput._ , 60:858–873, 2017.
  * Wang & Gupta (2020) Wang, S. and Gupta, M. R.  Deontological ethics by monotonicity shape constraints.  In _AISTATS_ , pp. 2043–2054, 2020.
  * Wotawa (2001) Wotawa, F.  A variant of Reiter’s hitting-set algorithm.  _Inf. Process. Lett._ , 79(1):45–51, 2001.
  * You et al. (2017) You, S., Ding, D., Canini, K. R., Pfeifer, J., and Gupta, M. R.  Deep lattice networks and partial monotonic functions.  In _NeurIPS_ , pp. 2981–2989, 2017.

##  Appendix A Appendix

######  Theorem 3.

Determining the existence of $\lfloor N/2\rfloor{+}1$ AXp’s (or CXp’s) of a
monotonic $N$-feature classifier is NP-complete.

###### Proof.

We say that a CNF is trivially satisfiable if some literal occurs in all
clauses. Clearly, SAT restricted to non-trivial CNFs is still NP-complete. Let
$\Phi$ be a not trivially-satisfiable CNF on variables $x_{1},\ldots,x_{k}$.
Let $N=2k$. Let $\tilde{\Phi}$ be identical to $\Phi$ except that each
occurrence of a negative literal $\overline{x_{i}}$ ($1\leq i\leq k$) is
replaced by $x_{i+k}$. Thus $\tilde{\Phi}$ is a CNF on $N$ variables each of
which occur only positively. Define the boolean classifier $\kappa$ by
$\kappa(x_{1},\ldots,x_{N})=1$ if $x_{i}=x_{i+k}=1$ for some
$i\in\\{1,\ldots,k\\}$ or $\tilde{\Phi}(x_{1},\ldots,x_{N})=1$ (and 0
otherwise). To show that $\kappa$ is monotonic we need to show that
$\mathbf{a}\leq\mathbf{b}\Rightarrow\kappa(\mathbf{a})\leq\kappa(\mathbf{b})$.
This follows by examining the two cases in which $\kappa(\mathbf{a})=1$: if
$a_{i}{=}a_{i+k}\land\mathbf{a}{\leq}\mathbf{b}$, then $b_{i}{=}b_{i+k}$,
whereas, if $\tilde{\Phi}(\mathbf{a}){=}1\land\mathbf{a}{\leq}\mathbf{b}$,
then $\tilde{\Phi}(\mathbf{b})=1$ (by positivity of $\tilde{\Phi}$), so in
both cases $\kappa(\mathbf{b})=1{\geq}\kappa(\mathbf{a})$.

We first consider AXp’s. Clearly $\kappa(\mathbf{1})=1$. There are $N/2$
obvious AXp’s of this prediction, namely $(i,i{+}k)$ ($1{\leq}i{\leq}k$).
These are minimal by the assumption that $\Phi$ is not trivially satisfiable.
Suppose that $\Phi(\mathbf{u}){=}1$. Let ${\mathcal{X}}_{\mathbf{u}}$ be
$\\{i\,|\,1{\leq}i{\leq}k\land
u_{i}{=}1\\}\cup\\{i{+}k\,|\,1{\leq}i{\leq}k\land u_{i}{=}0\\}$. Then (some
subset of) ${\mathcal{X}}_{\mathbf{u}}$ is an AXp of the prediction
$\kappa(\mathbf{1}){=}1$. The converse also holds. Thus, determining whether
$\kappa(\mathbf{1}){=}1$ has more than $N/2$ AXp’s is equivalent to testing
the satisfiability of $\Phi$. NP-completeness follows from the fact that
$\lfloor N/2\rfloor{+}1$ AXp’s are a polytime verifiable certificate.

The proof for CXp’s is similar. Clearly $\kappa(\mathbf{0})=0$. Again, there
are $N/2$ obvious CXp’s of this prediction, namely $(i,i{+}k)$
($1{\leq}i{\leq}k$) and (some subset of) ${\mathcal{X}}_{\mathbf{u}}$ is a CXp
iff $\Phi(\mathbf{u}){=}1$. Thus, determining whether $\kappa(\mathbf{0}){=}0$
has more than $N/2$ CXp’s is equivalent to testing the satisfiability of
$\Phi$, from which NP-completeness again follows. ∎

In the case of AXp’s, the theorem follows from a result on boolean monotone
functions (Babin & Kuznetsov, 2011), but for clarity of exposition we
preferred to give a direct proof.


††thanks: Electronic address<EMAIL_ADDRESS>
# Higher Order Corrections to Antisymmetric Spherical Solutions in $f(R)$
gravity

Garrett Gould  Department of Physics, University of California, Santa Barbara

###### Abstract

The static, spherical solutions that exhibit an antisymmetry between the
temporal and radial coordinates in $f(R)$ gravity theories are presented. I
present the constraint for this antisymmetry and show that pure $R^{2}$ models
produce these solutions. An asymptotic expansion approach is taken to find the
solutions, presenting first and second order corrections to the Schwarzschild
metric. The null geodesics and photon sphere are determined for these
corrections. Furthermore, I show that electromagnetic fields produce this
antisymmetry as well. The metric solutions are analyzed in the context of
black holes, bringing forward potentially a new class of black holes in
modified gravity theories.

##  1 Introduction

Despite the high success of general relativity since its inception, it is
accepted to be a low energy effective theory Donoghue (1994); Burgess (2004),
and either lacks detailed description or fails to describe select
gravitational phenomena, such as dark energy Peebles and Ratra (2003), cosmic
inflation Tsujikawa (2003); Debono and Smoot (2016); Vachaspati and Trodden
(1999), and dark matter Srivastava et al. (2023); Capozziello et al. (2006);
Arun et al. (2017); Srednicki (2000). This issue suggests work leading towards
describing the higher energy regime of gravitational interactions, reaching a
theory of strong gravity. Previous work in modified gravity theories, such as
$f(R)$ gravity theory, show that these said phenomena can be described Carroll
et al. (2004). In this paper, I analyze the static, spherical solutions in
$f(R)$ gravity theories, and pay special attention to those models that have
an antisymmetry in the metric between the temporal and radial coordinates, an
antisymmetry that is present in the Schwarzschild solution.

###  1.1 $f(R)$ gravity

The Einstein-Hilbert action, which gives rise to the Einstein field equations,
is

$S=\frac{1}{16\pi}\int\differential^{4}x\;\sqrt{-g}R$ (1.1)

The lagrangian for this action can be generalized to take on any arbitrary
function of $R$, $f(R)$. In the metric formalism of $f(R)$ gravity, the action
is Sotiriou and Faraoni (2010); Faraoni (2008),

$S=\frac{1}{16\pi}\int\differential^{4}x\;\sqrt{-g}f(R)$ (1.2)

Minimizing the action gives the $f(R)$ field equations,

$f^{\prime}(R)R_{\mu\nu}-\frac{1}{2}f(R)g_{\mu\nu}+(g_{\mu\nu}\Box-\nabla_{\mu}\nabla_{\nu})f^{\prime}(R)=8\pi
T_{\mu\nu}$ (1.3)

It is useful when determining the dynamics of the equations and solving for
the metric tensor by examining the trace of the field equations. For $f(R)$
models, this is,

$Rf^{\prime}(R)-2f(R)+3\Box f^{\prime}(R)=8\pi T$ (1.4)

##  2 Vacuum Spherical Solutions in $f(R)$ gravity

The general form of a static, spherical metric is Carroll (2022)

$\differential s^{2}=-e^{2a(r)}\differential t^{2}+e^{2b(r)}\differential
r^{2}+r^{2}\differential\Omega^{2}$ (2.1)

Upon which, for the general relativistic case, the unique vacuum solution is
the Schwarzschild metric, solved by both $R=0$ and $R_{\mu\nu}=0$ Wald (2009).
There exists a symmetry between $R_{tt}$ and $R_{rr}$, upon which for $R=0$,
then $a(r)=-b(r)$, giving rise to anti-symmetry between the temporal and
radial coordinates Carroll (2022).
In the simplest nontrivial $f(R)$ model, where $f(R)=R^{2}$, this anti-
symmetry is present. In general, for any choice of $f(R)$, the Ricci tensor in
vacuum is,

$R_{\mu\nu}=\frac{1}{f^{\prime}(R)}\left(\frac{1}{2}f(R)g_{\mu\nu}+\nabla_{\mu}\nabla_{\nu}-g_{\mu\nu}\Box
f^{\prime}(R)\right)$ (2.2)

As with the general relativistic case, the setup for determining that the
temporal and radial components of the Ricci tensor vanish independently
involves setting $e^{2(a-b)}R_{rr}+R_{tt}=0$ Carroll (2022), which is true in
vacuum for $f(R)=R$. In general, this expression in vacuum is,

$e^{2(a-b)}R_{rr}+R_{tt}=\frac{1}{f^{\prime}(R)}\nabla_{r}\nabla_{r}f^{\prime}(R)$
(2.3)

Where all non-radial derivatives of $f^{\prime}(R)$ go to zero due to the
scalar curvature being a single variable function of the radial coordinate.
For any static, spherical metric in $f(R)$ gravity, the scalar curvature must
be time-independent, a result guaranteed by the extension of Birkhoff’s
theorem to $f(R)$ gravity Ravindranath et al. (2018); Capozziello and Sáez-
Gómez (2012).
In order for the symmetry between the temporal and radial metric coefficients
to present itself in $f(R)$ gravity models, the right hand side of equation
(2.3) must go to zero.

$\displaystyle\frac{1}{f^{\prime}(R)}\nabla_{r}\nabla_{r}f^{\prime}(R)$
$\displaystyle=\frac{1}{f^{\prime}(R)}\left(\partial_{r}^{2}f^{\prime}(R)+\Gamma_{rr}^{r}\partial_{r}f^{\prime}(R)\right)$
(2.4)
$\displaystyle=\frac{1}{f^{\prime}(R)}\left(\partial_{r}^{2}f^{\prime}(R)+\partial_{r}b\partial_{r}f^{\prime}(R)\right)$
(2.5)

So, for any case of $f(R)$ such that
$\partial_{r}^{2}f^{\prime}(R)+\partial_{r}b\partial_{r}f^{\prime}(R)=0$, the
metric will have $a(r)=-b(r)$.

###  2.1 Pure $R^{2}$ solutions

Consider the simplest nontrivial choice of $f(R)$, which is $f(R)=\alpha
R^{2}$, where $\alpha$ is a normalization constant. I will work in units such
that $\alpha$ = 1 to simplify all equations. The trace equation for this model
simplifies down to,

$\Box R=0$ (2.6)

Expanding the d’Alembert operator for this metric,

$\Box R=g^{\mu\nu}\nabla_{\mu}\nabla_{\nu}=g^{rr}\nabla_{r}\nabla_{r}$ (2.7)

Where $\Box R$ is reduced to derivatives of the radial coordinate due to the
scalar curvature being time-independent and spherically symmetric. Since this
equation is zero, then the $R^{2}$ model in vacuum must satisfy the condition
in equation (2.5), and therefore the metric will be anti-symmetric between the
temporal and radial coordinate.
In this model, even in vacuum, there should exist non-zero values of $R$ due
to the second-order homogeneous equation in (2.7). This can also be seen by
solving for the Ricci tensor in (1.3).

$R_{\mu\nu}=\frac{1}{f^{\prime}(R)}\left(\frac{1}{2}f(R)g_{\mu\nu}+\nabla_{\mu}\nabla_{\nu}f^{\prime}(R)-g_{\mu\nu}\Box
f^{\prime}(R)\right)$ (2.8)

So in the case of $f^{\prime}(R)\neq 1$, the Ricci tensor will have non-zero
values in vacuum. Because of the anti-symmetry between $R_{tt}$ and $R_{rr}$,
it will be best to solve for $R$ by solving both the trace equation (2.6) and
solving for the angular Ricci tensor,

$R_{\theta\theta}=\frac{1}{4}r^{2}R$ (2.9)

Where $\Box R=0$ was substituted to simplify the equation. The primary aim of
this section is to find the metric for this specific model. Rewriting the
metric in (2.1) such that $e^{2a(r)}=A(r)$,

$\differential s^{2}=-A(r)\differential t^{2}+\frac{1}{A(r)}\differential
r^{2}+r^{2}\differential\Omega^{2}$ (2.10)

In terms of the metric coefficient $A(r)$, the scalar curvature $R$ is Carroll
(2022),

$R(r)=-\frac{2}{r^{2}}\left(A(r)+2rA^{\prime}(r)-1\right)-A^{\prime\prime}(r)$
(2.11)

Upon which the homogeneous solution to this equation, as well as the angular
Ricci tensor equation is the Schwarzschild metric. So, for non-zero scalar
curvature values, there should exist a metric such that,

$A(r)=1-\frac{2M}{r}+\epsilon(r)$ (2.12)

Where $\epsilon(r)$ is the inhomogeneous coupled solution between the angular
Ricci tensor and the scalar curvature. This extra term can be considered to be
the "correction" to the Schwarzschild metric arising from considering a
modified gravity theory such as $f(R)=R^{2}$. In terms of the metric in
(2.10), the $R_{\theta\theta}$ component is Carroll (2022),

$\displaystyle R_{\theta\theta}$ $\displaystyle=1-A(r)-rA^{\prime}(R)$ (2.13)
$\displaystyle=-\epsilon(r)-r\epsilon^{\prime}(r)$ (2.14)

Setting (2.9) and (2.14) equal to each other, I get the second differential
equation to solve for this case,

$-\epsilon(r)-r\epsilon^{\prime}(r)=\frac{1}{4}r^{2}R(r)$ (2.15)

Since the highest order term of $R(r)$ in terms of the metric will come from
$A^{\prime\prime}(r)$, and $R(r)$ vanishes for $A(r)\sim 1/r^{2}$, then I can
establish a constraint on the possible choices of $R(r)$,

$\lim_{r\rightarrow\infty}R(r)r^{2}=0$ (2.16)

A constraint that is required for all asymptotically flat spacetimes, and also
shown in Pretel and Duarte (2022). Since $R(r)$ must now be asymptotically
flat and approach flatness faster than the growth of $r^{2}$, a possible
solution of $R(r)$ is an asymptotic expansion about $r=0$,

$R(r)=\sum_{n=2}^{\infty}\frac{b_{n}(2M)^{n-2}}{r^{n}}$ (2.17)

Where powers of the Schwarzschild radius, $2M$, are included to preserve
units. From (2.3) and (2.11), a choice for the perturbation to the metric is
another asymptotic expansion about $r=0$,

$\epsilon(r)=\sum_{n=0}^{\infty}\frac{a_{n}(2M)^{n}}{r^{n}}$ (2.18)

Now, to solve for the coefficients of the $R$ expansion, I solve (2.7) for a
static, spherical metric. This equation becomes,

$A\partial_{r}^{2}R+\frac{2}{r}\partial_{r}R-\partial_{r}A\partial_{r}R=0$
(2.19)

Because of the highly non-linear nature of equation (2.15), a series solution
approach will provide a great approximation to the solutions, especially an
asymptotic series approach. Direct substitution of equation (2.13), (2.14)
into the differential equation (2.15) gives,

$\displaystyle\sum_{n=2}^{\infty}\frac{n(n-1)b_{n}(2M)^{n-2}}{r^{n+2}}-\sum_{n=2}^{\infty}\frac{n^{2}b_{n}(2M)^{n-1}}{r^{n+3}}+\sum_{n=2}^{\infty}\frac{(2M)^{n-4}}{r^{n}}\left(S_{n}^{(1)}-S_{n}^{(2)}\right)=0$
(2.20)

Where,

$\displaystyle S_{n}^{(1)}$
$\displaystyle=\sum_{l=0}^{n}(n-l)(n-l+1)a_{l-2}b_{n-l}$ (2.21) $\displaystyle
S_{n}^{(2)}$ $\displaystyle=\sum_{l=0}^{n}(l-2)(n-l)a_{l-2}b_{n-l}$ (2.22)

To account for the nonlinear terms of the differential equation, I employed
the Cauchy product (cite here), which gives rise to the finite sums in
equations (2.17), (2.18). Since the powers and the starting index of the sums
are not aligned, the first few terms of the first and third infinite series in
(2.16) need to be extracted and evaluated. This gives the three equations,

$\displaystyle S_{2}^{(1)}-S_{2}^{(2)}$ $\displaystyle=0$ (2.23)
$\displaystyle S_{3}^{(1)}-S_{3}^{(2)}$ $\displaystyle=0$ (2.24)
$\displaystyle S_{4}^{(1)}-S_{4}^{(2)}+2b_{2}$ $\displaystyle=0$ (2.25)

Expanding out these equations,

$\displaystyle 3a_{-1}b_{1}+10a_{-2}b_{2}=0$ (2.26) $\displaystyle
2a_{0}b_{1}+8a_{-1}b_{2}+18a_{-2}b_{3}=0$ (2.27) $\displaystyle
a_{1}b_{1}+2b_{2}+6a_{0}b_{2}+15a_{-1}b_{3}+28a_{-2}b_{4}=0$ (2.28)

Examining the first equation, since $a_{-1}=a_{-2}=0$, this is statement that
the expression is true, and $b_{2}$ is undefined from this equation. Moving on
to the second equation, $b_{2}$ and $b_{3}$ are also undefined, as well as
$a_{0}$. So, the third equation will be needed to determine these
coefficients. Solving for $b_{2}$ gives,

$b_{2}=\frac{-a_{1}b_{1}-15a_{-1}b_{3}-28a_{-2}b_{4}}{2+6a_{0}}$ (2.29)

Upon which $b_{2}=0$. Now, since the $a_{0}$ term in the expansion of
$\epsilon(r)$ will be a constant, and will not contain any powers of $r$, then
it is assumed that $a_{0}=0$, otherwise the constant term of the metric will
be not equal to one, which is a requirement for the metric to be
asymptotically flat. Since this is a constant, a transformation of the
coordinates could also take care of this problem, but I can set $a_{0}=0$
without loss of generality. It will be shown later that this is concluded as
well when relating $b_{n}$ to $a_{n}$ via the $R_{\theta\theta}$ equation.
Since $b_{2}=2$, the recursion relation for higher terms starts for $n>2$, and
is

$b_{n+1}=\frac{1}{n(n+1)}(n^{2}b_{n}+S_{n+3}^{(2)}-S_{n+3}^{(1)})$ (2.30)

All coefficients in the series can be expressed by the preceding coefficient.
Since the non-zero series coefficients start at $b_{3}$, and this coefficient
can take on any numerical value without loss of generality, the series
expression of $R$ can be expanded in terms of $b_{3}$

$\displaystyle
R(r)=\frac{(2M)}{r^{3}}b_{3}-\frac{3}{4}\frac{(2M)^{2}}{r^{4}}b_{3}(a_{1}-1)+\frac{(2M)^{3}}{r^{5}}\left(\frac{3}{5}b_{3}(a_{1}-1)^{2}-\frac{3}{10}a_{2}b_{3}\right)+O\left(\frac{1}{r^{6}}\right)$
(2.31)

Where $a_{1},a_{2}$ are determined by solving for $\epsilon(r)$ using (2.11).
Taking the first order of $R$ and solving (2.11) gives $a_{1}$,

$a_{1}=-\frac{b_{3}}{4}\ln(r/r_{0})$ (2.32)

Where $r_{0}$ is a length scale parameter introduced when solving the
differential equation in order to keep the argument of the logarithm
dimensionless. Higher terms of $a_{n}$ are solved to be,

$a_{n}=\frac{1}{4(n-1)}b_{n+2},\;\;n>1$ (2.33)

Upon which, these coefficients give the corrections to the metric. So, as
previously stated, the $a_{0}$ term in the expansion of $\epsilon(r)$ must be
zero. Relabeling $b_{3}=\xi$, the higher order corrected metric is now,

$\displaystyle
A(r)=1-\frac{2M}{r}-\frac{M}{2}\frac{\xi\ln(r/r_{0})}{r}+\frac{3}{16}\frac{M^{2}}{r^{2}}\left(\xi^{2}\left(1+\ln(\frac{r}{r_{0}})\right)+4\xi\right)+O\left(\frac{1}{r^{3}}\right)$
(2.34)

Since higher-order corrections to the metric than stated in equation (2.34)
will contain higher powers of $\xi$, in order for these additional terms to be
considered corrections and not diverge, the constraint $|\xi|<1$ is needed.
Since the logarithm term will grow at a slower rate than the $r^{-n}$ terms,
no further constraints are needed on the logarithm terms in order to ensure
convergence.

##  3 Vacuum spherical geodesics

The first application of the results derived in Section 2.1 is the effects of
the corrections on the geodesics, specifically the corrections to the stable
orbits. In general, the geodesics will be determined by Carroll (2022),

$\derivative[2]{x^{\mu}}{\tau}+\Gamma_{\lambda\nu}^{\mu}\derivative{x^{\lambda}}{\tau}\derivative{x^{\nu}}{\tau}$
(3.1)

For any metric coefficient $A(r)$, the orbit potential will be,

$V(r)=\frac{1}{2}A(r)\left(\frac{L^{2}}{r^{2}}+\kappa\right)$ (3.2)

Where $\kappa=-g_{\mu\nu}U_{\mu}U_{\nu}$ is the normalization constant of the
four-velocity. In the case of null geodesics, this will always be zero.
To first order, the potential in the $R^{2}$ model is,

$V(r)=\frac{1}{2}\left(1-\frac{2M}{r}-\frac{M}{2}\frac{\xi\ln(r/r_{0})}{r}\right)\left(\frac{L^{2}}{r^{2}}+\kappa\right)$
(3.3)

Figure 1: The potentials for massless particles are plotted, varying in $\xi$.
The maxima decrease with $\xi$, but the location of maxima, where the stable
orbital radii are, increase.

The stable orbital radii can be determined by minimizing the potential. In the
case of massless particles, where $\kappa=0$,

$r_{c}=3M\left(1-\frac{\xi}{12}+\frac{\xi\ln(r_{c}/r_{0})}{4}\right)$ (3.4)

Because of the logarithm term, the solution will involve the Lambert W
function Corless et al. (1996). The photon sphere radius in the first order
corrected metric is,

$r_{c}=3M-\frac{3}{4}M\xi W_{0}(z)$ (3.5)

Where $W_{0}$ is the principal branch of the Lambert W function, and

$z=-\frac{4r_{0}}{3M\xi}\exp(\frac{1}{3}-\frac{4}{\xi})$ (3.6)

Figure 2: The general relativistic and the first order corrected photon sphere
radii are plotted. The photon sphere due to the first order corrections is
increased slightly as a function of $\xi$. This corresponds with 1, where the
maxima of the potentials shifted towards larger $r$ at larger values of $\xi$.

Because of the logarithm term, all higher order corrections require numerical
evaluations of the photon sphere radius, as there will not be a closed form
solution.

##  4 $f(R)$ spherical solutions with matter sources

When introducing a matter field, the Ricci tensor in (2.2) becomes,

$R_{\mu\nu}=\frac{1}{f^{\prime}(R)}\left(8\pi
T_{\mu\nu}+\frac{1}{2}f(R)g_{\mu\nu}+\nabla_{\mu}\nabla_{\nu}f^{\prime}(R)-g_{\mu\nu}\Box
f^{\prime}(R)\right)$ (4.1)

So, for $R^{2}$ models, the trace equation is,

$\Box R=\frac{8\pi}{3}T$ (4.2)

And the angular Ricci tensor component is,

$R_{\theta\theta}=\frac{1}{2R}\left(8\pi
T_{\theta\theta}+\frac{1}{2}r^{2}R^{2}-\frac{8\pi}{3}r^{2}T\right)$ (4.3)

For the anti-symmetry between the temporal and radial components of the metric
to exist when introducing a source, the equation (2.3) is generalized to
include source terms,

$e^{2(a-b)}R_{rr}+R_{tt}=\frac{1}{f^{\prime}(R)}\nabla_{r}\nabla_{r}f^{\prime}(R)+\frac{8\pi}{f^{\prime}(R)}(g_{tt}g^{rr}T_{rr}+T_{tt})$
(4.4)

So, for this anti-symmetry to continue over when introducing a matter field,
there needs to be the same anti-symmetry in the energy-momentum tensor, that
is,

$g^{rr}T_{rr}=-g^{tt}T_{tt}$ (4.5)

Matter fields that automatically satisfy this relation are those with a
traceless energy-momentum tensor, such as for electromagnetic fields.

###  4.1 $R^{2}$ corrections with electromagnetic fields

Since the trace of the electromagnetic energy-momentum tensor, then the trace
equation for $R^{2}$ remains the same as vacuum, and therefore $R(r)$ is
unchanged. But, the metric will gain terms corresponding to the electric and
magnetic fields. This is a similar process to that of the Reissner-Nordström
metric, where the electric field is introduced as a second term to the
Schwarzschild metric Wald (2009); Carroll (2022).
To solve for the metric, I employ equation (2.15), but with a source term
introduced,

$-\epsilon(r)-r\epsilon^{\prime}(r)=\frac{1}{4}r^{2}R+\frac{4\pi}{R}T_{\theta\theta}$
(4.6)

The angular component of the electromagnetic energy-momentum tensor is
Matyjasek (1997); Anghinoni et al. (2022),

$\displaystyle T_{\theta\theta}$ $\displaystyle=-r^{4}\sigma_{\theta\theta}$
(4.7) $\displaystyle=\frac{1}{2}\frac{Q_{0}^{2}}{16\pi^{2}}$ (4.8)

Where $\sigma_{ij}$ is the Maxwell stress tensor, and
$Q_{0}^{2}=(Q_{e}^{2}+Q_{m}^{2})$ represents the total electric and magnetic
charge of the system. Although $R(r)$ exhibits the same solutions in the case
of an electromagnetic field source, I transform the powers of $2M$ in the
expansion of the metric and the scalar curvature to a general variable
$\mathcal{R}$, upon which the coordinate singularity radius for each specific
case would be substituted. In the static case, $\mathcal{R}=2M$, and for the
non-rotating charged case,
$\mathcal{R}=\frac{1}{2}\left(2M\pm\sqrt{4M^{2}-Q_{0}^{2}/\pi}\right)$. So,
the expansion of $R(r)$ becomes,

$R(r)=\frac{\mathcal{R}}{r^{3}}\xi-\frac{3}{4}\frac{\mathcal{R}^{2}}{r^{4}}\xi(a_{1}-1)+\frac{\mathcal{R}^{3}}{r^{5}}\left(\frac{3}{5}\xi(a_{1}-1)^{2}-\frac{3}{10}a_{2}\xi\right)+O\left(\frac{1}{r^{6}}\right)$
(4.9)

For the case of the electromagnetic field source, the $a_{1}$ coefficient of
the metric expansion is,

$a_{1}=-\frac{1}{4}\mathcal{R}\xi\ln(r/r_{0})-\frac{Q^{2}(r)}{32\pi\xi\mathcal{R}^{2}}$
(4.10)

Where $Q(r)=Q_{0}r^{2}$. The metric is then,

$A(r)=1-\frac{2M}{r}-\frac{M}{2}\mathcal{R}\xi\frac{\ln(r/r_{0})}{r}-\frac{MQ^{2}(r)}{16\pi\xi\mathcal{R}^{2}r}+O\left(\frac{1}{r^{2}}\right)$
(4.11)

Because $Q(r)$ will overpower the denominator at large $r$, the metric will
not be asymptotically flat. This is due to the $1/R$ term in the angular Ricci
tensor equation. A power series solution for $R(r)$ would overcome this, but
then would be significantly overpowered by the first term on the right-hand-
side of the angular Ricci tensor equation. So, for an antisymmetric metric
with an electromagnetic source, the spacetime would be highly unstable,
rendering the $R^{2}$ model unphysical in the case of an antisymmetric metric.

##  5 Modified Black Holes

Even though the electromagnetic case was rendered to be unphysical for the
specific metric considered in this paper, a point mass $M$ can still be
considered. Because of the corrections to the metric in equation (2.34), the
dynamics of the event horizons and the coordinate singularities will change,
in comparison to the Schwarzschild metric. Consider the first order
corrections to the metric,

$A(r)=1-\frac{2M}{r}-\frac{M}{2}\frac{\xi\ln(r/r_{0})}{r}$ (5.1)

There are two things to note here. First, when $r=r_{0}$, the classical
Schwarzschild metric is produced, regardless of the choice of $\xi$. Second,
the slope of the light cone behaves differently than in the general
relativistic case. The slope is Carroll (2022),

$\frac{dt}{dr}=\pm\left(1-\frac{2M}{r}-\frac{M}{2}\frac{\xi\ln(r/r_{0})}{r}\right)^{-1}$
(5.2)

Which diverges not only at a new coordinate singularity location larger than
the Schwarzschild radius, depending on the value of $\xi$. But, the slope of
the light cone flips direction again once reaching an inflection point at the
maximum point of the slope. Transforming to a modified version of Eddington-
Finkelstein coordinates for infalling geodesics to represent this behavior
Finkelstein (1958), the metric is,

$\differential
s^{2}=-\left(1-\frac{2M}{r}-\frac{M}{2}\frac{\xi\ln(r/r_{0})}{r}\right)\differential
v^{2}+\differential v\differential r+\differential r\differential
v+r^{2}\differential\Omega^{2}$ (5.3)

The maximum of the slope of the light cone occurs at,

$r_{\text{max}}=r_{0}\exp(1-\frac{4}{\xi})$ (5.4)

At this point, the light cone slope rapidly increases, until crossing a second
coordinate singularity that arises from the first order correction. Upon
crossing this point, the light cone inverts its behavior back to when it was
outside the first coordinate singularity, with a rapidly decreasing slope. The
second coordinate singularity occurs at,

$r=-\frac{1}{2}M\xi W_{0}\left(z_{\text{EH}}\right)$ (5.5)

Where the argument of the Lambert W function is,

$z_{\text{EH}}=-\frac{2r_{0}}{M\xi}e^{-4/\xi}$ (5.6)

Figure 3: Plot of the slope of the light cone. To show the inflection point.
As with previous figures, the units are such that $M=r_{0}=1$. To emphasize
the location of the inflection point, $\xi=3$, even though it is constrained
to $|\xi|<1$ from Section 1. The two coordinate singularities are highlighted
as red dashed lines. Figure 4: Diagram of the light cones when crossing the
two coordinate singularities $r_{1}$ and $r_{2}$, as well as the inflection
point. The diagram is on a v-r plot, indicating infalling geodesics.

##  6 Conclusion

The static, spherical metric solutions in $f(R)=R^{2}$ that exhibit an
antisymmetry between the temporal and radial coordinates were presented by an
asymptotic series expansion due to the singularity at $r=0$. To ensure each
higher order term does not rapidly diverge, a constraint on the expansion
coefficient $\xi$ is set such that $|\xi|<1$. Because of this constraint, all
deviations away from general relativity due to the corrections from $R^{2}$
gravity are small. This presents challenges when testing for these
corrections. For example, at the limit of the $\xi$ constraint, the photon
sphere radius increases by $M/40$, a very small deviation away from general
relativity.
When introducing a matter source, the only viable energy-momentum tensor that
would conserve the antisymmetric metric was an electromagnetic field, due to
the traceless nature of its energy-momentum tensor. But, after further
analysis, these results and the corrections that the electromagnetic fields
produced onto the metric in $R^{2}$ gravity were shown to either be unphysical
or highly unstable, indicating that a different approach to finding the scalar
curvature and metric in the presence of an electromagnetic field is needed.
The condition of the antisymmetry between the radial and temporal coordinate
may need to be relaxed in order to determine these solutions. In the case of
black holes, in the first order correction, a second horizon presented itself.
With the constraint on $\xi$, this second horizon would be very close to the
curvature singularity at the center. There exists a point of inflection, where
the tilting of the light cones due to the crossing of the first event horizon
reverts, and the slopes of the light cone rapidly increases as the test
particle falls towards the second horizon. Upon crossing the second horizon,
the light cone points outwards, behaving inversely to the outer event horizon,
preventing any infalling matter from reaching the center.
It is possible that this second horizon provides a "shell" around the
curvature singularity at $r=0$. Test particles that fall inwards inside the
black hole would be trapped, since they would have a strong gravitational
field acting inwards towards the center, but upon crossing the second horizon,
the test particles would experience an inverse force acting outwards,
launching the particles outside of the second horizon, trapping them around
this region. This would create a shell of matter that surrounds the
singularity, never reaching it, and only space-like geodesics crossing the
second horizon could reach the singularity, inversely to the first event
horizon.
Taking this into consideration, a new class of black holes, specifically black
holes in modified gravity theories like $f(R)$ gravity, may be imminent.
Further development of the solutions in $f(R)$ gravity models, such as other
choices of models, and other modified theories of gravity are needed. Such an
example is $f(R,T)$ gravity, where the energy-momentum trace is directly
included into the Einstein-Hilbert action Harko et al. (2011). The further
study of spherical solutions, and potentially other rotating models such as
Kerr-Newman black holes, may be a route in finding further deviations away
from general relativity that present themselves in modified gravity theories.
These deviations are likely a great route in testing and falsifying these
models of gravity to determine what is needed to seek for a effective theory
of gravity that provides better insight than general relativity.

## References

  * Donoghue (1994) J. F. Donoghue, Phys. Rev. D 50, 3874 (1994), URL https://link.aps.org/doi/10.1103/PhysRevD.50.3874.
  * Burgess (2004) C. P. Burgess, Living Reviews in Relativity 7 (2004).
  * Peebles and Ratra (2003) P. J. E. Peebles and B. Ratra, Rev. Mod. Phys. 75, 559 (2003), URL https://link.aps.org/doi/10.1103/RevModPhys.75.559.
  * Tsujikawa (2003) S. Tsujikawa, arXiv:hep-ph/0304257 (2003).
  * Debono and Smoot (2016) I. Debono and G. F. Smoot, Universe 2 (2016), ISSN 2218-1997, URL https://www.mdpi.com/2218-1997/2/4/23.
  * Vachaspati and Trodden (1999) T. Vachaspati and M. Trodden, Phys. Rev. D 61, 023502 (1999), URL https://link.aps.org/doi/10.1103/PhysRevD.61.023502.
  * Srivastava et al. (2023) Y. Srivastava, G. Immirzi, J. Swain, O. Panella, and S. Pacetti, The European Physical Journal C 83 (2023).
  * Capozziello et al. (2006) S. Capozziello, V. F. Cardone, and A. Troisi, Journal of Cosmology and Astroparticle Physics 2006, 001–001 (2006).
  * Arun et al. (2017) K. Arun, S. Gudennavar, and C. Sivaram, Advances in Space Research 60, 166–186 (2017).
  * Srednicki (2000) M. Srednicki, The European Physical Journal C 15, 143–144 (2000).
  * Carroll et al. (2004) S. M. Carroll, V. Duvvuri, M. Trodden, and M. S. Turner, Phys. Rev. D 70, 043528 (2004), URL https://link.aps.org/doi/10.1103/PhysRevD.70.043528.
  * Sotiriou and Faraoni (2010) T. P. Sotiriou and V. Faraoni, Rev. Mod. Phys. 82, 451 (2010), URL https://link.aps.org/doi/10.1103/RevModPhys.82.451.
  * Faraoni (2008) V. Faraoni, arXiv preprint arXiv:0810.2602 (2008).
  * Carroll (2022) S. M. Carroll, _Spacetime and geometry: An introduction to general relativity_ (Cambridge University Press, 2022).
  * Wald (2009) R. M. Wald, _General relativity_ (Univ. of Chicago Press, 2009).
  * Ravindranath et al. (2018) P. J. Ravindranath, Y. Aditya, D. R. Reddy, and M. V. Subba Rao, The European Physical Journal Plus 133 (2018).
  * Capozziello and Sáez-Gómez (2012) S. Capozziello and D. Sáez-Gómez, Annalen der Physik 524, 279 (2012), eprint https://onlinelibrary.wiley.com/doi/pdf/10.1002/andp.201100244, URL https://onlinelibrary.wiley.com/doi/abs/10.1002/andp.201100244.
  * Pretel and Duarte (2022) J. M. Pretel and S. B. Duarte, Classical and Quantum Gravity 39, 155003 (2022).
  * Corless et al. (1996) R. M. Corless, G. H. Gonnet, D. E. Hare, D. J. Jeffrey, and D. E. Knuth, Advances in Computational Mathematics 5, 329–359 (1996).
  * Matyjasek (1997) J. Matyjasek, Phys. Rev. D 55, 809 (1997), URL https://link.aps.org/doi/10.1103/PhysRevD.55.809.
  * Anghinoni et al. (2022) B. Anghinoni, G. Flizikowski, L. Malacarne, M. Partanen, S. Bialkowski, and N. Astrath, Annals of Physics 443, 169004 (2022), ISSN 0003-4916, URL https://www.sciencedirect.com/science/article/pii/S0003491622001506.
  * Finkelstein (1958) D. Finkelstein, Phys. Rev. 110, 965 (1958), URL https://link.aps.org/doi/10.1103/PhysRev.110.965.
  * Harko et al. (2011) T. Harko, F. S. N. Lobo, S. Nojiri, and S. D. Odintsov, Phys. Rev. D 84, 024020 (2011), URL https://link.aps.org/doi/10.1103/PhysRevD.84.024020.