# Mining Potentially Explanatory Patterns via Partial Solutions

A proposed methodology and experimental results

GianCarlo A.P.I. Catalano<EMAIL_ADDRESS>,  Alexander E.I. Brownlee
<EMAIL_ADDRESS>,  David Cairns<EMAIL_ADDRESS>University of
StirlingComputing Science and MathematicsStirlingUK ,  John McCall
<EMAIL_ADDRESS>Robert Gordon UniversityNational Subsea CentreAberdeenUK
and  Russell Ainslie<EMAIL_ADDRESS>BT TechnologyApplied
ResearchIpswichUK

###### Abstract.

Genetic Algorithms have established their capability for solving many complex
optimization problems. Even as good solutions are produced, the user’s
understanding of a problem is not necessarily improved, which can lead to a
lack of confidence in the results. To mitigate this issue, explainability aims
to give insight to the user by presenting them with the knowledge obtained by
the algorithm.

In this paper we introduce Partial Solutions in order to improve the
explainability of solutions to combinatorial optimization problems. Partial
Solutions represent beneficial traits found by analyzing a population, and are
presented to the user for explainability, but also provide an explicit model
from which new solutions can be generated. We present an algorithm that
assembles a collection of Partial Solutions chosen to strike a balance between
high fitness, simplicity and atomicity. Experiments with standard benchmarks
show that the proposed algorithm is able to find Partial Solutions which
improve explainability at reasonable computational cost without affecting
search performance.

Genetic Algorithms, Explainable AI (XAI), combinatorial optimization problems

††ccs: Computing methodologies Genetic algorithms††ccs: Theory of computation
Models of computation

##  1\. Introduction

Genetic Algorithms (GA) are able to solve many kinds of optimization problems
by producing high-quality results, which promotes their adoption even in
critical applications (Xu et al., 2019; Grace et al., 2018). A problem
presented by GAs, but also by Machine Learning, is that they act as black
boxes which humans struggle to trust, despite the apparent quality of the
outputs (Barredo Arrieta et al., 2020). This mistrust is not unjustified: a
GA’s results can be sub-optimal for many reasons, especially due to stochastic
behavior, incorrect fitness function implementations, unintended biases and
design issues (Xu et al., 2019).

This is not a small problem, since user trust becomes increasingly important
as GAs are introduced into high-stake applications (Xu et al., 2019). The
solution to this issue is explainability, where the user is presented with the
knowledge synthesized by the algorithm to solve the task. This can be used to
debug the implementation, provide documentation for auditing purposes,
manually adjusting and improve understanding of computer-generated solutions
(Xu et al., 2019).

###  1.1. Optimization problems

The present paper will concentrate on combinatorial optimization problems,
which can be particularly complex to solve despite their relatively simple
construction. Many have been shown to be computationally hard, such as Graph
Coloring (Hartmanis, 2006), Knapsack (Vance, 2006) and resource allocation
problems (Rabiner, 1984). As such, these often have to be approached using
metaheuristics such as GAs, but their internal mechanisms are often difficult
for humans to follow, suggesting the need for techniques similar to XAI for
Machine Learning (Xu et al., 2019).

There are many open questions which apply to optimization problems, that are
also particularly suited for explainability:

  * •

Why is a solution good, or better than another?

  * •

What characterizes good solutions?

  * •

How can the global optima be detected?

The approach we propose offers one way to address these questions. We
introduce the concept of a Partial Solution (PS), an explicit decomposition of
a positive trait identified in high-fitness solutions. A collection of PSs is
assembled and acts as a model that is then sampled to generate new solutions.
In doing so, the algorithm is inherently explainable, offering the ability for
the end user to identify key components of the recommended solutions.

These Partial Solutions offer a novel way of providing both global and local
explanations, but at the same time improve the search efficiency when solving
the problem. We test whether it is possible to consistently find these Partial
Solutions for some benchmark functions which have known underlying structures,
in Section 5.2.1 and Section 5.2.2. Afterwards, the ability to solve problems
efficiently is compared against a traditional GA and UMDA in Section 5.2.3.

##  2\. Related work

We now summarize relevant research in the EC and optimization literature
related to explainability.

###  2.1. Models within GAs and EDAs

Approaches which use a model to represent the problem, including surrogate
fitness functions such as neural networks and decision trees, or the
probabilistic models of Estimation of Distribution Algorithms (EDA), have been
a popular target for developing explainability functionality. This is because
these models are designed to represent abstract information in order to
accurately generalize the solution space. This knowledge can also be presented
to the user for explainability purposes (Bacardit et al., 2022), drawing on
the techniques developed for explainability in Machine Learning.

This paper makes use of related concepts from EDAs, as they present a fertile
avenue for explainability applications that can be adapted to metaheuristics.
Although traditional GAs lack explicit models, such structures can be added.
For example, Brownlee (2016) approximates the fitness function using a Markov
Network, which can be then presented to the user, whereas Oh and Ahn (2023)
use a marginal product model for linkage learning, and Yu and Goldberg (2006)
use a hierarchical linkage model. Some methods can be argued to be using the
population itself as a model, such as (Fyvie et al., 2021a, b), where the GA
is analyzed in terms of the trajectory within the search space that the
population takes during the evolution process. Similar to how an EDA generates
its model from a population, the present work will analyze a population and
produce a set of PSs, which then act as a model from which new solutions are
generated. On the other hand, since this “model” is only built once and is not
iterative, our approach is not a fully fledged EDA. Section 3.2 will give a
formal definition of the proposed concepts, and Section 3.4 will review the
novelty of this approach in more detail.

###  2.2. Backbones and building blocks

As discussed in (Bacardit et al., 2022), backbones are related to
explainability because they give insight into patterns found in solutions.
Backbones mostly relate to SAT problems, and consist of sub-configurations of
values which are shared among the global optima. Mathematically, backbones and
PSs share some similarities because they both represent commonalities within
sets of solutions, but PSs are made with explainability in mind, and thus they
are generated using different quality metrics and used very differently.

Both concepts are reminiscent of schemata, introduced by Holland in some early
theoretical work to formalize the convergence properties of GAs (Holland,
1992; White, 2014), but these theories were later shown to not match observed
results (Thornton, 1997). It should be noted that these results are limited to
Holland’s conjectures, and do not relate to the explicit construction of
Partial Solutions or their use in explainability, which will be the central
topics of this paper.

##  3\. Partial Solutions

###  3.1. Informal Introduction

To understand the behavior of the fitness function of a combinatorial
optimization problem, one can look at patterns associated with higher fitness
solutions. These may be higher-order interactions between groups of variables,
e.g. ‘all variables in the group must be the same value’, but in the present
paper the focus is on patterns where certain variables are fixed and others
can vary (indicated using *). For example 0**1 applies to 0001, 0011, 0101 and
0111.

A collection of such patterns, found by analyzing a population, provides a
succinct and expressive representation of what makes a solution “good” (as
seen in Figure 1).

Figure 1. The positive traits in a collection of full solutions
$\text{P}_{\text{Ref}}$ can be described by a collection of Partial Solutions
$X^{\psi}$

The items in $X^{\psi}$ represent beneficial traits that were observed in the
population $X$. While $X_{3}$ is the global optimum, the rest of the
population contains small imperfections which allow the algorithm to determine
what makes $X_{3}$ good.

In a real problem, PSs can offer useful insight into what variables are
linked, and what values they should take. In an employee-allocation problem,
for instance, these sub-configurations might indicate that it is beneficial to
assign certain groups of employees to specific tasks. A solution such as
1111111 might have its high fitness explained by the presence of 111*** and
*****1, although 111001 is even better because of the presence of ***00*.

Additionally, PSs might show the presence of positive traits that cannot
coexist (e.g., **00 and *11*), which can explain the presence of many near-
global optima. The ability to find such disagreeing PSs can help in explaining
situations where generally positive traits are not necessarily found in the
global optima.

Partial solutions are an explicit decomposition of the problem, which has two
benefits:

  * •

Partial Solutions are interpretable:

    * –

Global explanations: describing the fitness landscape by finding simple and
atomic Partial Solutions which are associated with high-fitness.

    * –

Local explanations: describing a solution by pointing out the Partial
Solutions it contains.

  * •

Partial Solutions assist the search process, because they act as a model that
describes high-fitness regions of the solution space (similar to Probabilistic
Graphical Models) from which solutions may be constructed (see Section 4.2).

###  3.2. Formal definition

Table 1. Notation Symbol | Meaning
---|---
$F$ | Fitness function to be maximized
$F^{\psi}$ | Fitness function for partial solutions
$n$ | The number of parameters in the solutions
$a,b,\dots,x$ | Variables used to denote a _full_ solution
$x_{1},x_{2},\dots,x_{n}$ | The parameters of $x$, a full solution
$a^{\psi},b^{\psi},\dots,x^{\psi}$ | Variable used to denote a _Partial_ Solution
$x^{\psi}_{1},x^{\psi}_{2},\dots,x^{\psi}_{n}$ | The parameters of $x^{\psi}$, a Partial Solution
$P_{\text{Ref}}$ | A collection of evaluated full solutions

Table 1 summarizes our notation. A Partial Solution $x^{\psi}$ can be
described as “a sub-configuration of parameter values which is associated with
high-fitness solutions within a reference population”, where that reference
population will be denoted as $P_{\text{Ref}}$.

In the same way that a full solution is a tuple of values
$x=(x_{1},x_{2},...,x_{2})$, a Partial Solution is
$x^{\psi}=(x^{\psi}_{1},x^{\psi}_{2},...x^{\psi}_{n})$ where each
$x^{\psi}_{i}$ can take any value that $x_{i}$ can, as well as “*”, the “any”
value. A solution $x$ is said to “contain” the Partial Solution $x^{\psi}$
when they agree on the fixed values found in $x^{\psi}$.

(1)
$\text{contains}(x,x^{\psi})\Leftrightarrow(x^{\psi}_{i}=*)\vee(x^{\psi}_{i}=x_{i})\
\forall\ i\in[1,n]$

Conversely, it is possible to find all the solutions in $P_{\text{Ref}}$ which
contain the Partial Solution $x^{\psi}$:

(2) $\text{observations}_{P_{\text{Ref}}}(x^{\psi})=\\{x\in P_{\text{Ref}}\ |\
\text{contains}(x,x^{\psi})\\}$

For an optimization problem, within all of the possible PSs there is a small
subset which is preferable for explainability, the catalog . These are the PSs
which are optimal with respect to three metrics described in the next section:
simplicity, mean fitness and atomicity.

###  3.3. Partial Solution metrics

It is preferable for a PS to have many of its parameters take the * value,
since fewer fixed parameters improve interpretability. We define this metric
to be:

(3) $\text{simplicity}(x^{\psi})=|\\{i\ \in\ [1,n]|\ x^{\psi}_{i}=*\\}|$

The simplest PS possible has no fixed variables, and consists entirely out of
*’s. Other than being simple, a PS also needs to be associated with high-
fitness solutions, which can be quantified by finding the mean fitness of its
observations (cf. Equation 2).

(4)
$\text{meanFitness}(x^{\psi})=\text{average}(F(\text{observations}_{P_{\text{Ref}}}(x^{\psi}))\
$

In order to have PSs assist with explainability, the previous two metrics are
not sufficient. Assume $x^{\psi}$ and $y^{\psi}$ are optimal in different
metrics, and let $z^{\psi}$ be the union of their fixed values. $z^{\psi}$ is
likely to also be high-scoring, since it will have similar observations to its
parents, and its simplicity might not be significantly worse than $x^{\psi}$
and $y^{\psi}$, even though it is redundant compared to them.

This effect was described in Holland’s Building Block Hypothesis as beneficial
to the GA process, but for our purposes this is actually undesirable, since it
gets in the way of interpretability due to all the redundancy present. To
solve this problem, it will be necessary to define a third metric, which will
encapsulate the idea that a PS should be “irreducible”: atomicity.

For a Partial Solution $x^{\psi}$, its fixed-value elements $x^{\psi}_{j}$
each contribute in some way to its average fitness. We define a normalized
fitness function, remapped in the range [0, 1]:

(5) $\displaystyle\text{F}_{\text{Norm}}(x)=\frac{\text{F}(x)-m}{S}$
$\displaystyle(\text{where}\ m=\min_{x\ \in\ P_{\text{Ref}}}F(x),S=$
$\displaystyle\sum_{x\ \in\ P_{\text{Ref}}}F(x)-m)$

In order to measure the contribution of the variable $x^{\psi}_{k}$ within
$x^{\psi}$, define the following modifications of $x^{\psi}$:

  * •

isolate: where only $x^{\psi}_{k}$ is fixed

  * •

exclude: where $x^{\psi}_{k}$ is *

And we define their “benefit” to be the sum of the normalized fitnesses of
their observations.

(6) $\displaystyle\text{isolate}(x^{\psi},k)=y^{\psi}\ s.t.\
y^{\psi}_{i}=\left\\{\begin{matrix}x^{\psi}_{i}&\text{if}\ i=k\\\
*&\text{otherwise}\end{matrix}\right.$ (7)
$\displaystyle\text{exclude}(x^{\psi},k)=z^{\psi}\ s.t.\
z^{\psi}_{i}=\left\\{\begin{matrix}x^{\psi}_{i}&\text{if}\ i\neq k\\\
*&\text{otherwise}\end{matrix}\right.$ (8)
$\displaystyle\text{benefit}(a^{\psi})=\sum_{x\ \in\
\text{observations}_{P_{\text{Ref}}}(a^{\psi})}F_{\text{Norm}(x)}$

Finally, we define the contribution of the variable $x^{\psi}_{k}$ as:

(9) $\displaystyle\text{contribution}(x^{\psi},k)$ $\displaystyle=p_{AB}\
\cdot\ \log\left(\frac{p_{AB}}{p_{A*}\ \cdot\ p_{*B}}\right)$
$\displaystyle\text{where}\ p_{AB}$ $\displaystyle=\text{benefit}(x^{\psi}),$
$\displaystyle p_{A*}$
$\displaystyle=\text{benefit}(\text{isolate}(x^{\psi},k)),$ $\displaystyle
p_{*B}$ $\displaystyle=\text{benefit}(\text{exclude}(x^{\psi},k))$

This is derived from mutual information, as used in DSMGA-II (Hsu and Yu,
2015). Continuing the example from Section 3.1, the Partial Solution 111**1
will have its contributions calculated at $k\in\\{1,2,3,6\\}$, and since
independent groups are present the contribution scores will be low. For
111***, the contributions are calculated at $k\in\\{1,2,3\\}$, which will all
be greater than those in the previous example.

In the same way that a chain is as strong as its weakest link, the atomicity
of $x^{\psi}$ is defined as the minimum of its contributions, which ensures
that all of the fixed parameters in the PS depend on each other and are
essential.

(10)
$text{atomicity}(x^{\psi})=\text{min}\left\\{\left.\text{contribution}(x^{\psi},k))\
\right|\begin{subarray}{c}\ k\ \in\ [1,\ n]\\\ \ x^{\psi}_{k}\neq\
\texttt{*}\end{subarray}\ \right\\}$

Each of these metrics is necessary in order to characterize the PSs which are
best suited for explainability:

  * •

Simplicity prioritizes more understandable PSs

  * •

Mean-fitness promotes beneficial PSs

  * •

Atomicity prevents redundant PSs

###  3.4. Novelty

Having defined Partial Solutions, we can now review related existing concepts
to demonstrate where our novelty lies.

Schemata and backbones are similarly structured to Partial Solutions, but they
are used very differently. Schemata are abstractions that were used for
theoretical work on GAs, and are implicitly sampled, i.e. not meant to be
explicitly instantiated. In contrast, our work does instantiate them, similar
to other studies which use them to improve search efficiency (Fletcher and
Wennekers, 2017; Goldberg et al., 2003), whereas in this work they are
designed mainly for explainability.

Backbones relate to SAT problems where multiple satisfactory solutions share a
common sub-configuration. For most problems there is only one backbone, and it
is preferred for it to have many fixed variables in order to shrink the search
space (Slaney and Walsh, 2001), whereas in this work the aims are reversed to
improve interpretability.

Since in this work the PS catalog is treated as a model, EDAs are also an
adjacent topic. In this sense, there are some related works: Linkage learning
(Baluja and Davies, 1997; Hsu and Yu, 2015) and Family-of-Subsets (Thierens
and Bosman, 2011a) aim to find which variables are interacting in the fitness
function, and are then assembled into structures such as sets or cliques
(Pelikan et al., 2002). These structures only determine which variables are
interacting, but their optimal value assignments are not specified, whereas
PSs contain both. In other words, PSs are instances of the cliques / linked
sets found by the aforementioned approaches.

##  4\. Overall system

The system that we propose consists of:

  * •

PS Miner: obtains the PS catalog from a reference population, described in
Section 4.1 and Algorithm 1

  * •

Pick & Merge: an algorithm which forms solutions by combining elements from
the catalog , described in Section 4.2 and Algorithm 2.

The following (one pass) sequence of operations combines them:

  1. (1)

Generate a reference population $P_{\text{Ref}}$, and evaluate it

  2. (2)

Apply the PS Miner on $P_{\text{Ref}}$ to generate the PS catalog

  3. (3)

Apply pick_and_merge on the PS catalog to obtain full solutions

The process of producing solutions offers the explainability benefits of PSs
and maintains the ability to solve the original optimization problem (see
Figure 2), all of which will be tested in Section 5.

Figure 2. The proposed system, which uses the PS catalog for both
explainability and to generate full solutions.

By using Partial Solutions as an intermediate, the system is able to generate
full solutions and explainability, specifically local explanations (when
comparing against full solutions) and global explanations.

###  4.1. Finding Partial Solution Catalogs

The algorithm to find the PS catalog will be implemented using the metrics in
Section 3.2, but the design of the search method involves some non-trivial
design choices, which we will now describe.

####  4.1.1. Reference Population

$P_{\text{Ref}}$ is used to evaluate the meanFitness and atomicity of PSs, and
remains fixed throughout the PS Mining process. It’s important for
$P_{\text{Ref}}$ to cover enough of the search space, and in the present work
this was done using random uniform sampling. Section 5.2.2 investigates the
necessary size of $P_{\text{Ref}}$ and whether it might be beneficial to
evolve it.

Table 2. PS local search variants Variant |  get_local |  get_init
---|---|---
Simplification-only |  Simplifications of $x^{\psi}$ |  returns the top individuals in $P_{\text{Ref}}$ as Partial Solutions
Specialization-only |  Specializations of $x^{\psi}$ |  returns the universal PS: **…**
Full-local |  simplifications $\cup$ specializations | The union of the above two methods

####  4.1.2. Multiple objectives and $F^{\psi}$

To simplify the system, the three objectives are aggregated into a single
objective $F^{\psi}$. This is obtained by remapping each metric into the range
$[0,1]$ and taking the mean, as seen in Algorithm 1.

Note that an evaluation of a full solution using $F$ is a different operation
from evaluating a Partial Solution with $F^{\psi}$, and the latter uses pre-
calculated fitness values from $P_{\text{Ref}}$ for efficiency.

####  4.1.3. Exclusion archive

A problem to be tackled when searching for the catalog is that the PS search
space is large and contains many local optima clusters, which hinder the
search process. Their presence is due to the fact that local optima for $F$
also act as fully-set PS local optima for meanFitness, and their neighborhoods
contain many similar individuals but with ‘*’ values. In order to assist the
algorithm, an exclusion archive based method is proposed, where the items
selected in every iteration are put into an archive and are not allowed to
reappear in the population again. In another sense, it marks which regions
have been visited and prevents the population from revisiting them. The
archive stores all of the best solutions seen throughout the iterations of the
algorithm, and upon termination these are evaluated and the top scoring PSs
are returned. This approach is implemented as shown in Algorithm 1, and it was
tested as part of Section 5.2.1.

1

2

3

4

5Def  _remap(_values_)_:

6 return $\frac{values-\min(values)}{\max(values)-\min(values)}$

7

8Def  _$F^{\psi}$(_$X^{\psi}$_)_:

9 $M\leftarrow\textnormal{{{remap}}}(\texttt{meanFitness}(X^{\psi}))$

10 $S\leftarrow\textnormal{{{remap}}}(\texttt{simplicity}(X^{\psi}))$

11 $A\leftarrow\textnormal{{{remap}}}(\texttt{atomicity}(X^{\psi}))$

12 return $\\{\texttt{avg}(m,s,a)\ \texttt{for}\
m,s,a\in\texttt{zip}(M,S,A)\\}$

13

14Def  _top(_$X^{\psi}$ , quantity_)_:

15
$\text{sorted}\leftarrow\text{sort}(X^{\psi},\text{key}=\textnormal{{$F^{\psi}$}})$

16 return sorted[:quantity]

17

19

18 get_init and get_local are discussed in Section 4.1

20

21

22Def  _mine_ps(_get_init , get_local, pop_size, qty_ret_)_:

23
$X^{\psi}\leftarrow\textnormal{{{get\\_init}(}}\textnormal{\emph{}}\textnormal{{)}}$

24 $\text{archive}\leftarrow\\{\\}$

25

26 while  _termination_criteria_not_met_ do

27 $\text{selected}\leftarrow\texttt{TournamentSelection}(X^{\psi})$

28

29 $\text{localities}\leftarrow\bigcup_{x^{\psi}\in\
\text{seleted}}\text{get\\_local}(x^{\psi})$

30

31 $\text{archive}\leftarrow\text{archive}\cup\text{selected}$

32

33
$X^{\psi}\leftarrow(X^{\psi}\cup\text{localities})\smallsetminus\text{archive}$

34

35 $X^{\psi}\leftarrow\textnormal{{{top}}}(X^{\psi},\texttt{pop\\_size})$

36

37

38 return $\textnormal{{{top}}}(\texttt{archive},\texttt{qty\\_ret})$

39

Algorithm 1 Archive-based PS Miner

1

2

3

4Def  _weighted_random_choice(_$X^{\psi}$_)_:

/* samples a random element of $X^{\psi}$ using $F^{\psi}$ as weights */

5

6

7Def  _merge_from(_$X^{\psi}$ , limit_)_:

8 $\texttt{available}\leftarrow\texttt{copy}(X^{\psi})$

9 $x^{\psi}\leftarrow\texttt{**..**}$

10 $\texttt{added}\leftarrow 0$

11 while  _$\texttt{available}\neq\varnothing$ &
$\texttt{added}<\texttt{limit}$_ do

12
$y^{\psi}\leftarrow\textnormal{{{weighted\\_random\\_choice}(}}\textnormal{\emph{{available}}}\textnormal{{)}}$

13
$\texttt{available}\leftarrow\texttt{available}\smallsetminus\\{y^{\psi}\\}$

14 if  _$\texttt{mergeable}(x^{\psi},y^{\psi})$_ then

15 $x^{\psi}\leftarrow\texttt{merge}(x^{\psi},y^{\psi})$

16 $\texttt{added}\leftarrow\texttt{added}+1$

17

18

19 return $x^{\psi}$

20

21

22Def  _fill_gaps(_$x^{\psi}$_)_:

23 for  _$x^{\psi}_{i}\in x^{\psi}$_ do

24 if  _$x^{\psi}=*$_ then

25
$x^{\psi}_{i}\leftarrow\texttt{random}.\texttt{randrange}(\texttt{cardinalities}[i])$

26

27

28 return $x^{\psi}$

29

30Def  _generate_via_pick_and_merge(_$X^{\psi}$ , merge_limit_)_:

31 return fill_gaps(_merge_from(_$X^{\psi}$ , merge_limit_) _)

32

Algorithm 2 Pick-and-Merge algorithm

####  4.1.4. Local search

Knowing that the PS search space might contain many clustered optima, it is
beneficial to include a form of local search. Given a PS $x^{\psi}$, its local
neighborhood consists of

  * •

Simplifications: one of the fixed variables is replaced by *

  * •

Specializations: one of the * is replaced by a fixed value

Local search may use either or both of these, which determine the
implementation of get_local and get_init in Algorithm 1. The three possible
implementations will be tested as part of Section 5.2.1.

###  4.2. Pick & Merge algorithm

The PS catalog obtained via Algorithm 1 can be used as a model to construct
high-quality full solutions. Algorithm 2 simply picks items from the PS
catalog and merges them when possible. At the end, any remaining unfixed
parameters are filled with random values, so that a full solution can be
returned. Let $X^{\psi}$ be the catalog found for a problem, and define the
following operations, which are used in Algorithm 2:

(11) $\displaystyle\text{mergeable}(x^{\psi},y^{\psi})$
$\displaystyle=\neg\exists\ i\ s.t.\ (x^{\psi}_{i}\neq*\neq y^{\psi}_{i}\wedge
x^{\psi}\neq y^{\psi})$ (12) $\displaystyle\text{merge}(x^{\psi},y^{\psi})$
$\displaystyle=z^{\psi}\ s.t.\
z^{\psi}_{i}=\left\\{\begin{matrix}x^{\psi}_{i}&\text{if}\
x^{\psi}_{i}\neq*\\\ y^{\psi}_{i}&\text{if}\ y^{\psi}_{i}\neq*\\\
*&\text{otherwise}\end{matrix}\right.$ (13)
$\displaystyle\text{to\\_full}(x^{\psi})$ $\displaystyle=x\ s.t.\
x_{i}=x^{\psi}_{i}\ \forall\ i\in[1,n]$ $\displaystyle(\text{defined only
when}\ \neg\exists\ x^{\psi}_{i}=\ *)$

The parameter merge_limit used in pick_and_merge determines how many PSs
should be contained in a generated solution, which for this paper was set to
$\lceil\sqrt{n}\rceil$. Note how this algorithm does not use the original
fitness function, as it is assumed that $X^{\psi}$ is representative of
beneficial traits.

##  5\. Experiments

This section will present three research questions, RQ1, RQ2 and RQ3, answered
by the testing rounds T1, T2 and T3 respectively.

  * •

RQ1: Which algorithm parameters are best suited for finding the PS catalog ?

  * •

RQ2: How is the ability to extract Partial Solutions from a reference
population affected by its size and the amount of generations it has been
evolved for?

  * •

RQ3: Can Partial Solutions be used to obtain high-fitness full solutions? Our
goal is increasing explainability, but the algorithm must still be able to
solve the problem.

T1 and T2 will check which algorithms can produce the PS catalog effectively,
which we assume will help explainability. For testing, we’ll use benchmark
problems with known catalogs, defined in Section 5.1. T3 will test the
system’s ability to solve the original optimization problem by checking how
often it is able to find the global optima. The algorithms and the parameters
being tested will be contained in Section 5.2, along with the description of
the testing rounds. Finally, the results are presented and discussed in
Section 6.

###  5.1. Benchmark Problems

Here we draw upon benchmarks previously used by the EDA community, as they
have clearly defined structures that map well to graphical models, and present
well-defined “target” PSs which we aim for our algorithms to find. The
benchmark problems will need to be challenging enough to benefit from
explainability, and hence “simpler” ones such as OneMax and BinVal will not be
considered. Instead, problems which have dependencies between variables will
be chosen, and to prevent any biases due to proximity in linked variables, the
input variables will have their positions shuffled.

The Royal Road problem assigns a fitness to the input bit strings by dividing
them into disjoint adjacent groups of $k$ variables, and counting how many of
them are entirely consisting of 1’s (Mitchell et al., 1992). The target PSs
here are the individual groups, with their values all being set to 1. In this
study $k=4$ (as can be seen in Figure 3), with 20 parameters in total, meaning
that there are 5 target PSs.

Figure 3. Comparison of the RR and RRO problems, where RRO has more complex
global optima

The main difference between the two problems is that the Partial Solutions in
RRO can ”overlap”, meaning that some of them might not be able to coexist in
the same full solution. This has the interesting effect of often producing
many distinguished global optima, whereas in RR there is only one single
global optima, composed of the union of all the items in the PS catalog

A property of the Royal Road problem is that the target PSs are disjoint,
which might be too trivial for the algorithm to solve. In order to increase
the problem complexity, a variant called Royal Road with Overlaps (RRO) will
be introduced, in which the groups are allowed to have overlaps and use 0’s
instead, as seen in Figure 3.

More formally, RRO is described by $k$ (the size of each group), $q$ (the
number of groups) and $l$ (the amount of bits in the total input). When
$l<q\cdot k$ the groups will be forced to overlap, which increases problem
complexity. Let $G^{\psi}$ be the targets, in the form of PSs, then the
fitness function will be:

(14) $F_{RRO}(x)=|\\{g^{\psi}\in G^{\psi}\ \text{s.t.}\ \text{contains}\
(x,g^{\psi})\\}|$

In this paper, $k=4$, $q=5$ and $l=15$ (guaranteeing overlap), with the groups
being randomized for every run. A full solution containing all of the groups
will often not be possible, and thus the task of finding PSs does not align
with finding just the building blocks present in the global optima alone.

Another popular benchmark is Trap-k, in which each group of k bits is assigned
a fitness which increases monotonically with the amount of 0’s, but achieves
its maximum when it contains only 1’s, as seen in Equation 16 (Goldberg, 1989;
Thierens and Bosman, 2011b).

(15) $\displaystyle
F_{\text{Trap-k}}(x)=\sum_{i=0}^{(l/\kappa)-1}F^{\text{sub}}_{\text{Trap-k}}\left(\sum_{j=\kappa\cdot
i}^{(\kappa+1)\cdot i-1}x_{j}\right)$ (16) $\displaystyle
F^{\text{sub}}_{\text{Trap-k}}(u)=\begin{cases}k&\text{if }u=k,\\\
k-u-1&\text{if }u<k.\end{cases}$

The problem was devised to be challenging for GAs, due to its deceptive
fitness function, and it’s been considered here to check whether deceptive PSs
are a problem for our algorithms. Here, $k=5$ and $l=25$, so that there will
be 5 deceptive groups. The target PSs correspond to these groups, with all
their parameters being set to 1.

###  5.2. Setup

Table 3. PS Mining configurations Algorithm | Parameters
---|---
Own method | $\begin{cases}\text{Population\ Size}\ \in\\{50,100,150\\}\\\ \text{Local\ Search}\ \in\begin{Bmatrix}\begin{matrix}\text{Simplification-only},\\\ \text{Specialization-only},\\\ \text{Full-local}\end{matrix}\\\ \end{Bmatrix}\\\ \text{Uses\ Archive}\ \in\\{\text{True},\text{False}\\}\end{cases}$
GA | $\begin{cases}\text{One point mutation, rate of 0.075}\\\ \text{Two point crossover, rate of 0.7}\end{cases}$
Hill Climber | $\begin{cases}\text{One point mutation, rate of 0.075}\\\ \text{Repeated trials, with each trial ending}\\\ \text{when a local optima is reached}\end{cases}$
Table 4. Reference Population GA Parameters Option | Value
---|---
Selection | Tournament (size = 2), Elite = 2
Mutation | One point mutation, rate of 0.075
Crossover | Two point crossover, rate of 0.7

When devising a PS mining algorithm, there are many potential sources of
variability, including:

  * •

The algorithm itself, due to selection operators

  * •

The reference population, being randomly generated

  * •

The problem definition, in the case of RRO

Due to this, whenever data is gathered about an algorithm configuration being
executed on a certain problem, the run is executed 100 times, with a freshly
generated reference population (and problem definition for RRO). Additionally,
the algorithms will be restricted to returning at most 50 PSs for every run
(qty_ret in Algorithm 1).

####  5.2.1. T1

Many configurations of algorithm parameters will be compared in their ability
to find PSs efficiently. This will be done by measuring the number of
$F^{\psi}$ evaluations needed in order to find all the target PSs for the
benchmark problems. The algorithms and the parameters that will be compared
are shown in Table 3.

In each run, $P_{\text{Ref}}$ will consist of $10^{4}$ uniformly randomly
generated full solutions. The algorithms will then be executed with the
termination criterion being that they have found all of the target PSs, or
that they have gone over $10^{5}$ $F^{\psi}$ evaluations. Once the run
terminates, it is considered successful only if the output catalog contains
all of the target PSs.

####  5.2.2. T2

This round will contain a short analysis of how the ability to find the PS
catalog is affected by the size of $P_{\text{Ref}}$ and the number of
generations it has been evolved for. Populations of sizes
$\in\\{100,200,500,1K,2K,5K,10K\\}$ are evolved for generations
$\in\\{0,10,20,50,100,200\\}$. The evolution process is done using a
traditional GA, as described in Table 4 (derived from (Brownlee et al., 2010;
Shakya et al., 2005)), and the success criteria is the same as T1.

####  5.2.3. T3

The ability for the system to produce high-quality full solutions will be
tested by applying pick_and_merge as defined in Section 4.2. For evaluation
budgets $\in$ $\\{$1K, 5K, 10K, 15K, 20K, 25K, 30K$\\}$, three methods to
generate full solutions will be compared: a traditional GA (as defined for
T2), UMDA (both with population size of 150) and the method suggested in
Section 4. The latter will use the settings selected from T1 to produce the
catalog $X^{\psi}$ via the PS Miner algorithm, by consuming a percentage
($\in$ 10%, 20%, $\dots$, 90%) of the evaluation budget with $F^{\psi}$
evaluations. The remaining evaluations are spent on $F$ to evaluate solutions
and construct $P_{\text{Ref}}$. The catalog obtained is passed to
pick_and_merge to produce full solutions. For each of these three methods, a
run is considered successful if it produces a global optima for the problem.

##  6\. Results

###  6.1. T1

Table 5 shows the success rate for each PS Mining algorithm configuration,
with the clear winner having the following properties:

• method: Archive-based • population size: 150 • variant: Specialization-only

This configuration will be used in the tests that follow. In general,
configurations with the archive enabled were more successful than their
counterparts, with the specialization-only variants being the best across all
problems. Additionally, the GA and Hill Climber were less successful,
particulary with regard to the Trap-k problem.

Table 5. Results for T1: success rates for different algorithm configurations (in bold are figures ¿ 90%) Search method | Problem
---|---
Local search |  Pop. size |  Uses archive | RR | RRO | Trap-k
Simplification-only | 50 |  $\left\\{\begin{matrix}\\\ \\\ \\\ \\\ True\\\ \\\ \\\ \\\ \\\ \end{matrix}\right.$ | 0% | 1% | 0%
Simplification-only | 100 | 0% | 1% | 0%
Simplification-only | 150 | 0% | 5% | 0%
Specialization-only | 50 | 64% | 61% | 47%
Specialization-only | 100 | 96% | 93% | 86%
Specialization-only | 150 | 100% | 94% | 93%
Full-local | 50 | 0% | 36% | 0%
Full-local | 100 | 0% | 27% | 0%
Full-local | 150 | 0% | 31% | 0%
Simplification-only | 50 |  $\left\\{\begin{matrix}\\\ \\\ \\\ \\\ False\\\ \\\ \\\ \\\ \\\ \end{matrix}\right.$ | 0% | 0% | 0%
Simplification-only | 100 | 0% | 0% | 0%
Simplification-only | 150 | 0% | 0% | 0%
Specialization-only | 50 | 61% | 12% | 1%
Specialization-only | 100 | 95% | 85% | 10%
Specialization-only | 150 | 100% | 93% | 30%
Full-local | 50 | 0% | 0% | 0%
Full-local | 100 | 0% | 6% | 0%
Full-local | 150 | 0% | 18% | 0%
GA | 50 | - | 0% | 0% | 0%
GA | 100 | - | 0% | 0% | 0%
GA | 150 | - | 0% | 0% | 0%
Hill-climber | - | - | 28% | 31% | 0%

###  6.2. T2

The data in Table 6 suggests that the best $P_{\text{Ref}}$ to find the
catalog is unevolved, with its minimum size = $10^{4}$ for the tested
problems. A possible explanation for this effect is the phenomenon of
“hitchhiking”, where sub-optimal traits spread throughout a population due to
their co-occurrence with optimal traits. This effect was hypothesized by
Mitchell et al. (1993) upon observing similar results for Royal functions,
where methods which rely on schemata (such as (Brown et al., 2002)) perform
worse when applied on pre-evolved populations.

Table 6. PS catalog success rate for T2, in bold are figures ¿ 90% (success is finding all of the target Partial Solutions) RR | Gen’s | 0 | 10 | 20 | 50 | 100 | 200 |
---|---|---|---|---|---|---|---
$P_{\text{Ref}}$ Size |  |  |  |  |  |  |
100 | 0% | 0% | 0% | 0% | 0% | 0% |
200 | 0% | 0% | 0% | 0% | 0% | 0% |
500 | 52% | 23% | 22% | 12% | 23% | 23% |
1000 | 90% | 77% | 80% | 75% | 80% | 82% |
2000 | 100% | 99% | 98% | 99% | 98% | 99% |
5000 | 100% | 100% | 100% | 100% | 100% | 100% |
10000 | 100% | 100% | 100% | 100% | 100% | 100% |

RRO Gen’s | 0 | 10 | 20 | 50 | 100 | 200
---|---|---|---|---|---|---
$P_{\text{Ref}}$ Size |  |  |  |  |  |
100 | 4% | 0% | 0% | 0% | 0% | 0%
200 | 22% | 8% | 9% | 4% | 5% | 4%
500 | 66% | 55% | 47% | 50% | 45% | 46%
1000 | 83% | 82% | 77% | 75% | 81% | 78%
2000 | 89% | 88% | 89% | 87% | 87% | 84%
5000 | 89% | 86% | 88% | 91% | 90% | 89%
10000 | 94% | 89% | 89% | 89% | 89% | 89%

Trap-k Gen’s | 0 | 10 | 20 | 50 | 100 | 200
---|---|---|---|---|---|---
$P_{\text{Ref}}$ Size |  |  |  |  |  |
100 | 0% | 0% | 0% | 0% | 0% | 0%
200 | 0% | 0% | 0% | 0% | 0% | 0%
500 | 0% | 0% | 0% | 0% | 0% | 0%
1000 | 0% | 0% | 0% | 0% | 0% | 0%
2000 | 2% | 0% | 0% | 0% | 0% | 0%
5000 | 38% | 11% | 4% | 3% | 3% | 6%
10000 | 93% | 51% | 60% | 53% | 49% | 60%

Figure 4. (Relating to T3) Mean fitness obtained by dedicating some of the
evaluation budget to PS catalog mining.

For all of the problems, the highest mean fitness is obtained by using 20% to
25% of the budget to PS evaluations.

###  6.3. T3

Table 7 and Figure 4 show three important facts:

  1. (1)

The PS catalog mediated approach is able to consistently find the global
optima for all of the problems, when between 20% and 60% being the recommended
proportion of the budget to dedicate to $F^{\psi}$ evaluations.

  2. (2)

for the RR and RRO problems, the GA and the PS catalog mediated approach
perform similarly when around 50% of the evaluations are used for $F^{\psi}$.

  3. (3)

For Trap-k, only the proposed method was able to find the global optima
consistently, whereas the GA and UMDA almost always got stuck in the local
optima.

Table 7. Global optima success rate for T3, in bold are figures ¿ 90% (success is finding a global optima) RR | Eval’s | 1K | 5K | 10K | 15K | 20K | 25K | 30K
---|---|---|---|---|---|---|---
P&M |  |  |  |  |  |  |
% of eval’s for PS catalog | 10% | 9% | 71% | 100% | 100% | 100% | 100% | 100%
20% | 8% | 100% | 100% | 100% | 100% | 100% | 100%
30% | 6% | 100% | 100% | 100% | 100% | 100% | 100%
40% | 3% | 100% | 100% | 100% | 100% | 100% | 100%
50% | 3% | 100% | 100% | 100% | 100% | 100% | 100%
60% | 2% | 100% | 100% | 100% | 100% | 100% | 100%
70% | 1% | 98% | 100% | 100% | 100% | 100% | 100%
80% | 0% | 88% | 98% | 100% | 100% | 100% | 100%
90% | 0% | 17% | 79% | 95% | 100% | 100% | 100%
GA | 13% | 100% | 100% | 100% | 100% | 100% | 100%
UMDA | 0% | 3% | 1% | 0% | 0% | 0% | 0%

RRO Eval’s | 1K | 5K | 10K | 15K | 20K | 25K | 30K
---|---|---|---|---|---|---|---
P&M |  |  |  |  |  |  |
% of eval’s for PS catalog | 10% | 86% | 96% | 96% | 97% | 98% | 97% | 98%
20% | 84% | 95% | 97% | 99% | 99% | 98% | 97%
30% | 85% | 97% | 97% | 98% | 98% | 98% | 96%
40% | 81% | 97% | 97% | 99% | 96% | 98% | 96%
50% | 84% | 98% | 98% | 98% | 96% | 97% | 97%
60% | 76% | 96% | 98% | 99% | 98% | 98% | 97%
70% | 77% | 94% | 92% | 97% | 96% | 97% | 96%
80% | 68% | 92% | 93% | 99% | 95% | 96% | 96%
90% | 56% | 81% | 84% | 96% | 96% | 96% | 96%
GA | 98% | 100% | 99% | 100% | 99% | 98% | 98%
UMDA | 73% | 77% | 69% | 68% | 63% | 73% | 71%

Trap-k Eval’s | 1K | 5K | 10K | 15K | 20K | 25K | 30K
---|---|---|---|---|---|---|---
P&M |  |  |  |  |  |  |
% of eval’s for PS catalog | 10% | 0% | 0% | 0% | 53% | 71% | 82% | 81%
20% | 0% | 0% | 13% | 39% | 97% | 99% | 100%
30% | 0% | 0% | 5% | 79% | 99% | 99% | 99%
40% | 0% | 0% | 21% | 87% | 95% | 99% | 100%
50% | 0% | 0% | 14% | 74% | 92% | 99% | 100%
60% | 0% | 0% | 6% | 52% | 89% | 96% | 99%
70% | 0% | 0% | 4% | 38% | 72% | 89% | 93%
80% | 0% | 0% | 0% | 11% | 24% | 58% | 71%
90% | 0% | 0% | 0% | 0% | 4% | 3% | 8%
GA | 0% | 0% | 0% | 0% | 0% | 0% | 0%
UMDA | 0% | 0% | 0% | 0% | 0% | 0% | 0%

These results give evidence that the proposed methodology can solve
optimization problems with similar performance to a GA, improving on Trap-k,
while adding explainability.

##  7\. Conclusion

The algorithms presented in this paper describe a novel technique for
explainability, by using Partial Solutions. The PS catalog provides a global
explanation of the fitness landscape by presenting beneficial sub-
configurations, while local explanations are obtained by checking which PSs
are found in a solution. Moreover, the catalog acts similar to a model in an
EDA, from which new solutions can be sampled and thus solving the original
optimization problem. The testing rounds showed the following results:

  * •

RQ1: The algorithm which is most apt for finding the PS catalog in the
benchmarks was the Specialization-only local search, using the exclusion
archive.

  * •

RQ2: To obtain the PS catalog , it is better for $P_{\text{Ref}}$ to be a
large unevolved population

  * •

RQ3: The system is consistently able to get the global optima for the
benchmark problems by first generating a PS catalog and then applying
pick_and_merge. This suggests that the PS catalog can capture important
information about the fitness landscape. The recommended proportion of
$F^{\psi}$ evaluations is between 20% and 60%.

It should be noted that the PS Miner algorithms implemented in this work are
not as efficient as they could be, since $F^{\psi}$ involves complex
calculations, and experimentation indicates this to be between 50 and 100
times slower than $F$ for $|P_{\text{Ref}}|=10^{4}$. Therefore the
explainability of Partial Solutions still comes at a computational cost,
although this is might be partially countered by an improved chance at finding
global optima (as seen for Trap-k). The source code for the sytem can be
viewed at (Catalano, 2024).

##  8\. Future work

The PS catalog represents beneficial traits for the fitness function, and thus
they can act as a model that could be used to implement a surrogate fitness
function (similar to (Brownlee, 2009)) or improved GA operators (similar to
Optimal Mixing in (Thierens and Bosman, 2011b)). Overall, the system could be
readapted to be a full EDA, for example by feeding the outputs of
pick_and_merge into $P_{\text{Ref}}$.

The problems discussed in this document are admittedly simpler than real-world
problems, and therefore it would be useful to see how well the PSs perform
with some inherently inelegant problems, such as rostering and scheduling
tasks. In particular, there will be cases where interpretable PSs are not
always present or difficult to find. This will be investigated in problems
without clean-cut cliques, where the explainability of PSs will be explored
more thoroughly.

As mentioned in Section 3.2, a PS may be generalized to higher-order
relationships. This will both greatly improve the descriptive power of PSs,
but also greatly increase the complexity of the search methods. $F^{\psi}$
could also be extended by adding more metrics, such as novelty, robustness (to
changes in the parameter values with respect to fitness), or constraint
satisfaction.

###### Acknowledgements.

Work funded under a PhD project supported by The Data Lab and BT Group plc.
Many thanks go to the patience of my supervisors.

## References

  * (1)
  * Bacardit et al. (2022) Jaume Bacardit, Alexander E. I. Brownlee, Stefano Cagnoni, Giovanni Iacca, John McCall, and David Walker. 2022.  The Intersection of Evolutionary Computation and Explainable AI. In _Proceedings of the Genetic and Evolutionary Computation Conference Companion_ (Boston, Massachusetts) _(GECCO ’22)_. Association for Computing Machinery, New York, NY, USA, 1757–1762.  https://doi.org/10.1145/3520304.3533974
  * Baluja and Davies (1997) Shumeet Baluja and Scott Davies. 1997.  Using Optimal Dependency-Trees for Combinational Optimization. In _Proceedings of the Fourteenth International Conference on Machine Learning_ _(ICML ’97)_. Morgan Kaufmann Publishers Inc., San Francisco, CA, USA, 30–38.
  * Barredo Arrieta et al. (2020) Alejandro Barredo Arrieta, Natalia Díaz-Rodríguez, Javier Del Ser, Adrien Bennetot, Siham Tabik, Alberto Barbado, Salvador Garcia, Sergio Gil-Lopez, Daniel Molina, Richard Benjamins, Raja Chatila, and Francisco Herrera. 2020.  Explainable Artificial Intelligence (XAI): Concepts, taxonomies, opportunities and challenges toward responsible AI.  _Information Fusion_ 58 (2020), 82–115.  https://doi.org/10.1016/j.inffus.2019.12.012
  * Brown et al. (2002) D.F. Brown, A.B. Garmendia-Doval, and J.A.W. McCall. 2002.  Markov Random Field Modelling of Royal Road Genetic Algorithms. In _Artificial Evolution_ , Pierre Collet, Cyril Fonlupt, Jin-Kao Hao, Evelyne Lutton, and Marc Schoenauer (Eds.). Springer Berlin Heidelberg, Berlin, Heidelberg, 65–76.
  * Brownlee (2016) Alexander E.I. Brownlee. 2016.  Mining Markov Network Surrogates for Value-Added Optimisation. In _Proceedings of the 2016 on Genetic and Evolutionary Computation Conference Companion_ (Denver, Colorado, USA) _(GECCO ’16 Companion)_. Association for Computing Machinery, New York, NY, USA, 1267–1274.  https://doi.org/10.1145/2908961.2931711
  * Brownlee (2009) Alexander Edward Ian Brownlee. 2009.  _Multivariate Markov networks for fitness modelling in an estimation of distribution algorithm._ Ph. D. Dissertation. ”Robert Gordon University”.
  * Brownlee et al. (2010) Alexander E. I. Brownlee, Olivier Regnier-Coudert, John A. W. McCall, and Stewart Massie. 2010.  Using a Markov network as a surrogate fitness function in a genetic algorithm. In _IEEE Congress on Evolutionary Computation_. IEEE, New York City, USA, 1–8.  https://doi.org/10.1109/CEC.2010.5586548
  * Catalano (2024) Giancarlo Catalano. 2024.  _PS Assisted Explainability_.  https://github.com/Giancarlo-Catalano/PS_Minimal_Showcase
  * Fletcher and Wennekers (2017) Jack McKay Fletcher and Thomas Wennekers. 2017.  A natural approach to studying schema processing.  arXiv:1705.04536 [cs.NE]
  * Fyvie et al. (2021a) Martin Fyvie, John A.W. McCall, and Lee A. Christie. 2021a.  Non-deterministic solvers and explainable AI through trajectory mining.  _SICSA XAI Workshop 2021_ 2894 (6 2021), 75–78.  https://rgu-repository.worktribe.com/output/1395881https://rgu-repository.worktribe.com/output/1395881.abstract
  * Fyvie et al. (2021b) Martin Fyvie, John A. W. McCall, Lee A. Christie, Alexander E. I. Brownlee, and Manjinder Singh. 2021b.  Towards explainable metaheuristics: Feature extraction from trajectory mining.  _Expert Systems_ n/a, n/a (2021), e13494.  https://doi.org/10.1111/exsy.13494 arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1111/exsy.13494
  * Goldberg (1989) David E Goldberg. 1989.  Genetic algorithms and Walsh functions: Part 2, Deception and its analysis.  _Complex systems_ 3 (1989), 153–171.
  * Goldberg et al. (2003) David E. Goldberg, Kumara Sastry, and Yukio Ohsawa. 2003.  _Discovering Deep Building Blocks for Competent Genetic Algorithms Using Chance Discovery via KeyGraphs_.  Springer, Berlin, Heidelberg, Tokyo, Japan, 276–301.  https://doi.org/10.1007/978-3-662-06230-2_19
  * Grace et al. (2018) Katja Grace, John Salvatier, Allan Dafoe, Baobao Zhang, and Owain Evans. 2018.  Viewpoint: When Will AI Exceed Human Performance? Evidence from AI Experts.  _Journal of Artificial Intelligence Research_ 62 (7 2018), 729–754.  https://doi.org/10.1613/JAIR.1.11222
  * Hartmanis (2006) Juris Hartmanis. 2006.  Computers and Intractability: A Guide to the Theory of NP-Completeness (Michael R. Garey and David S. Johnson).  _https://doi.org/10.1137/1024022_ 24 (7 2006), 90–91.  Issue 1. https://doi.org/10.1137/1024022
  * Holland (1992) John H. Holland. 1992.  _Adaptation in Natural and Artificial Systems: An Introductory Analysis with Applications to Biology, Control, and Artificial Intelligence_.  The MIT Press, Massachussets, USA.  https://doi.org/10.7551/mitpress/1090.001.0001
  * Hsu and Yu (2015) Shih-Huan Hsu and Tian-Li Yu. 2015.  Optimization by pairwise linkage detection, incremental linkage set, and restricted/back mixing: DSMGA-II. In _Proceedings of the 2015 Annual Conference on Genetic and Evolutionary Computation_. Association for Computing Machinery, New York, NY, USA, 519–526.
  * Mitchell et al. (1992) Melanie Mitchell, Stephanie Forrest, and John Holland. 1992.  The Royal Road for Genetic Algorithms: Fitness Landscapes and GA Performance.  _European Conference on Artificial Life_ 1 (11 1992), 23–33.
  * Mitchell et al. (1993) Melanie Mitchell, John H. Holland, and Stephanie Forrest. 1993.  When will a genetic algorithm outperform hill climbing?. In _Proceedings of the 6th International Conference on Neural Information Processing Systems_ (Denver, Colorado) _(NIPS’93)_. Morgan Kaufmann Publishers Inc., San Francisco, CA, USA, 51–58.
  * Oh and Ahn (2023) Sanghoun Oh and Chang Wook Ahn. 2023.  Evolutionary Approach for Interpretable Feature Selection Algorithm in Manufacturing Industry.  _IEEE Access_ 11 (2023), 46604–46614.  https://doi.org/10.1109/ACCESS.2023.3274490
  * Pelikan et al. (2002) Martin Pelikan, David E Goldberg, and Fernando G Lobo. 2002.  A survey of optimization by building and using probabilistic models.  _Computational optimization and applications_ 21 (2002), 5–20.
  * Rabiner (1984) L. Rabiner. 1984.  Combinatorial optimization:Algorithms and complexity.  _IEEE Transactions on Acoustics, Speech, and Signal Processing_ 32 (12 1984), 1258–1259.  Issue 6. https://doi.org/10.1109/TASSP.1984.1164450
  * Shakya et al. (2005) Siddhartha Shakya, John McCall, and Deryck Brown. 2005.  Using a Markov Network Model in a Univariate EDA: An Empirical Cost-Benefit Analysis. In _Proceedings of the 7th Annual Conference on Genetic and Evolutionary Computation_ (Washington DC, USA) _(GECCO ’05)_. Association for Computing Machinery, New York, NY, USA, 727–734.  https://doi.org/10.1145/1068009.1068130
  * Slaney and Walsh (2001) John Slaney and Toby Walsh. 2001.  Backbones in Optimization and Approximation. In _Proceedings of the 17th International Joint Conference on Artificial Intelligence - Volume 1_ (Seattle, WA, USA) _(IJCAI’01)_. Morgan Kaufmann Publishers Inc., San Francisco, CA, USA, 254–259.  https://doi.org/10.5555/1642090.1642125
  * Thierens and Bosman (2011a) Dirk Thierens and Peter AN Bosman. 2011a.  Optimal mixing evolutionary algorithms. In _Proceedings of the 13th annual conference on Genetic and evolutionary computation_. Association for Computing Machinery, New York, NY, USA, 617–624.
  * Thierens and Bosman (2011b) Dirk Thierens and Peter A.N. Bosman. 2011b.  Optimal Mixing Evolutionary Algorithms. In _Proceedings of the 13th Annual Conference on Genetic and Evolutionary Computation_ (Dublin, Ireland) _(GECCO ’11)_. Association for Computing Machinery, New York, NY, USA, 617–624.  https://doi.org/10.1145/2001576.2001661
  * Thornton (1997) Chris Thornton. 1997.  The building block fallacy.  _Complexity International_ 4 (1997), 1–1.
  * Vance (2006) Pamela H. Vance. 2006.  Knapsack Problems: Algorithms and Computer Implementations (S. Martello and P. Toth).  _https://doi.org/10.1137/1035174_ 35 (7 2006), 684–685.  Issue 4. https://doi.org/10.1137/1035174
  * White (2014) David White. 2014.  An Overview of Schema Theory.  arXiv:1401.2651 [cs.NE]
  * Xu et al. (2019) Feiyu Xu, Hans Uszkoreit, Yangzhou Du, Wei Fan, Dongyan Zhao, and Jun Zhu. 2019.  Explainable AI: A Brief Survey on History, Research Areas, Approaches and Challenges.  _Lecture Notes in Computer Science (including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics)_ 11839 LNAI (2019), 563–574.  https://doi.org/10.1007/978-3-030-32236-6_51/FIGURES/12
  * Yu and Goldberg (2006) Tian-Li Yu and David E. Goldberg. 2006.  Conquering Hierarchical Difficulty by Explicit Chunking: Substructural Chromosome Compression. In _Proceedings of the 8th Annual Conference on Genetic and Evolutionary Computation_ (Seattle, Washington, USA) _(GECCO ’06)_. Association for Computing Machinery, New York, NY, USA, 1385–1392.  https://doi.org/10.1145/1143997.1144210


DESY-23-135

# Scalar dark matter production from the inflaton

M. Pierre  Deutsches Elektronen-Synchrotron DESY, Notkestr. 85, 22607 Hamburg,
Germany

###### Abstract

A curved space-time is known to act as a source for particle production in
environments where gravity plays a significant role. We explore this effect in
a minimal setup of cosmic inflation on the production of a scalar dark matter
candidate during and after the inflationary stage of the universe. We consider
the production of dark matter via direct coupling to the inflaton field or
from pure gravitational interactions. Cosmological constraints from structure
formation and dark matter isocurvature perturbations are discussed. A new
analytical expression for the isocurvature power spectrum is provided.

##  1 Dark matter production during and after inflation

We consider a minimal setup where a scalar dark matter (DM) candidate $\chi$
couples to a scalar inflaton field $\phi$ via a quartic coupling $\sigma$
whose general action reads

$\mathcal{S}\;=\;\int\mathop{}\\!\mathrm{d}^{4}x\sqrt{-g}\left[-\frac{1}{2}M_{P}^{2}R+\frac{1}{2}(\partial_{\mu}\phi)^{2}-V(\phi)+\frac{1}{2}(\partial_{\mu}\chi)^{2}-\frac{1}{2}m_{\chi}^{2}\chi^{2}-\dfrac{\sigma}{2}\phi^{2}\chi^{2}\right]\,.$
(1)

For concreteness, we consider the T-model [1] inflationary potential
$V(\phi)=6\lambda M_{P}^{4}\tanh^{2}(\phi/\sqrt{6}M_{P})$ that can be
normalized for a nominal choice of $N_{*}=55$ $e$-folds for the fiducial CMB
scale crossing $k_{*}=0.05\leavevmode\nobreak\ \text{Mpc}^{-1}$. The choice
$\lambda\simeq 2\times 10^{-11}$ corresponds to a spectral tilt $n_{s}\simeq
0.963$, tensor-to-scalar ratio $r\simeq 0.004$ and amplitude of scalar-
perturbations power spectra $A_{S^{*}}$ compatible with current CMB
measurements. Close to the minimum $\phi\ll M_{P}$, the potential can be
approximated by $V(\phi)\simeq(m_{\phi}^{2}/2)\phi^{2}$ where the inflaton
mass can be expressed as $m_{\phi}=\sqrt{2\lambda}M_{P}\simeq 1.6\times
10^{13}\leavevmode\nobreak\ \text{GeV}$. Introducing the rescaled DM field
$X\equiv a\chi$, varying the action (1) with respect to $X$ and moving to
Fourier space, we obtain the equation of motion expressed in terms of the
conformal time $\mathop{}\\!\mathrm{d}\eta=\mathop{}\\!\mathrm{d}t/a$ with $a$
being the scale factor and $t$ the cosmic time [2]

$(\partial_{\eta}^{2}+\omega_{k}^{2})X_{k}=0,\quad\text{with}\quad\omega_{k}^{2}=k^{2}+a^{2}m_{\rm{eff}}^{2}\quad\text{and}\quad
m_{\rm{eff}}^{2}=m_{\chi}^{2}+\sigma\phi^{2}+\frac{1}{6}R\,,$ (2)

where the effective DM mass $m_{\rm{eff}}$ is time dependent. In this work we
consider only scalars minimally coupled to gravity. A non-minimal DM coupling
to gravity can manifest itself as an extra term in the effective mass squared.
Extension of this discussion to non-minimal couplings can be found in [3].

During inflation, the quasi-homogeneous inflaton field $\phi(t)$ strongly
affects the background and space-time curvature. The corresponding Ricci
scalar $R$ causes the squared-mass to become negative and induces a tachyonic
excitation of light scalar fields for wavelengths stretched on super-horizon
scales, resulting in copious particle production. After the end of inflation
$a>a_{\text{end}}$, the fast inflaton oscillations about the quadratic minimum
are also responsible for particle production, induced via direct coupling
$\sigma$ or by pure gravitational interactions, i.e. via $\lambda$, in the
form of the Ricci scalar in Eq. (2).

Considering initially the absence of any DM particle corresponds to the Bunch-
Davies vacuum initial condition, given by
$\lim\limits_{\eta\rightarrow-\infty}X_{k}(\eta)=\frac{1}{\sqrt{2k}}e^{-ik\eta}$.
The comoving number density of produced scalar DM, $n_{\chi}$, can be computed
using the following expression [4, 5]:

$n_{\chi}\left(\frac{a}{a_{\rm{end}}}\right)^{3}\;=\;\int_{k_{0}}^{\infty}\mathop{}\\!\mathrm{d}k\frac{k^{3}}{2\pi^{2}}f_{\chi}(k,t)\leavevmode\nobreak\
\leavevmode\nobreak\ {\rm{with}}\leavevmode\nobreak\ \leavevmode\nobreak\
f_{\chi}(k,t)\,\equiv\,\frac{1}{2\omega_{k}}\left|\omega_{k}X_{k}-iX^{\prime}_{k}\right|^{2}\,,$
(3)

with $f_{\chi}$ being the Phase Space Distribution (PSD). By defining the
rescaled dimensionless comoving momentum $q=(k/m_{\phi})(a/a_{\text{end}})$,
one can distinguish modes with $q<1$ and $q>1$ that are respectively excited
during or after inflation.

Modes excited during inflation. For $q<1$, as discussed further on the phase
space distribution scales as $f_{\chi}\sim
q^{2m_{\text{eff}}^{2}/(3H^{2})-3}$. If $\sigma=0$, the DM production is
purely gravitational and the PSD scales as $f_{\chi}\propto q^{-3}$ for light
scalar DM $m_{\rm{\chi}}\ll H$. As $\sigma$ increases, the tachyonic
excitation of the DM scalar field becomes suppressed, resulting in a blue-tilt
of the spectrum at small $q$. A numerical evaluation of the PSD is represented
on the left panel of Fig. 1.

Modes excited after the end of inflation. Such modes with $q>1$ never
experience a tachyonic excitation as they remain sub-horizon during inflation.
For mild values of $0<\sigma/\lambda<10^{2}$, both production via direct or
gravitational couplings, i.e. induced via $\sigma$ or $\lambda$, result in a
PSD scaling as $f_{\chi}\propto q^{-9/2}$ and increasing with
$(\sigma/\lambda)^{2}$. For larger couplings $\sigma/\lambda>10^{2}$, strong
parametric resonances induced by the oscillating mass term in Eq. (2) result
in the formation of quasi-stochastic peaks in the PSD which has to be
evaluated numerically.

Figure 1: Dark matter phase space distribution (left) and isocurvature power
spectrum (right) for selected values of coupling ratio $\sigma/\lambda$. The
colour of each line codes the value of the coupling ratio $\sigma/\lambda$ and
is identical on both left and right panels.

##  2 Structure formation constraints

Light DM candidates produced from a state out of thermal equilibrium may
possess a sizable pressure component. Such pressure can suppress overdensities
on galactic scales during the matter domination era, resulting in a cutoff in
the matter power spectrum $\mathcal{P}(k)$ for scales $k$ larger than the
free-streaming horizon wavenumber
$k_{\text{H}}(a)\,=\,\left(\int_{0}^{a}k_{\text{FS}}^{-1}(\tilde{a})\mathop{}\\!\mathrm{d}\log\tilde{a}\right)^{-1}$,
function of the free-streaming wave number $k_{\text{FS}}$ and the expansion
history of the universe. Measurements from the Lyman-$\alpha$ forest
absorption-lines constrain this cutoff scale to
$k_{\text{H}}(a=1)>15\,h\,\text{Mpc}^{-1}$. It can be translated into a bound
on a generic warm dark matter (WDM) candidate initially thermally coupled to
the SM bath
$m_{\text{WDM}}>m_{\text{WDM}}^{\text{Ly}\mbox{-}\alpha}\;\sim\;3\leavevmode\nobreak\
\text{keV}$. The free-streaming wave number
$k_{\text{FS}}\,=\,\sqrt{9\mathcal{H}^{2}/(10w_{\chi})}$ is uniquely a
function of the conformal Hubble rate $\mathcal{H}\equiv aH$ and the DM
equation of state $w_{\chi}$ which can be approximated as

$w_{\chi}\,\simeq\,\dfrac{T_{\star}^{2}}{3m_{\chi}^{2}}\dfrac{\langle
q^{2}\rangle}{a^{2}}\,,\qquad\text{with}\qquad\langle
q^{2}\rangle\,\equiv\,\dfrac{\int\mathop{}\\!\mathrm{d}q\,q^{4}f_{\chi}(q)}{\int\mathop{}\\!\mathrm{d}q\,q^{2}f_{\chi}(q)}\,,\qquad\text{and}\qquad
T_{\star}\,\equiv\,m_{\phi}\left(\dfrac{a_{\text{end}}}{a_{0}}\right)\,,$ (4)

Since only the equation of state controls the cutoff scale in the matter power
spectrum, one can translate the WDM bound to our scenario, by using the
procedure introduced in [6], via

$m_{\chi}\,>\,m_{\chi}^{\text{Ly}\mbox{-}\alpha}\;=\;m_{\rm
WDM}^{\text{Ly}\mbox{-}\alpha}\left(\frac{T_{\star}}{T_{\rm
WDM,0}}\right)\sqrt{\frac{\langle q^{2}\rangle}{\langle q^{2}\rangle_{\rm
WDM}}}\,,$ (5)

where $\langle q^{2}\rangle_{\rm WDM}\simeq 12.93$ $T_{\rm WDM,0}$ is the WDM
temperature saturating the dark matter abundance.

Constraints for $\sigma/\lambda>1$. In this case, the second moment $\langle
q^{2}\rangle$ is UV dominated and depends on the reheating temperature
$T_{\text{reh}}$ through $\langle q^{2}\rangle\simeq
2.43\sqrt{a_{\text{reh}}/a_{\text{end}}}$. This scaling allows to cancel out
the $T_{\text{reh}}$ dependence on the dark matter energy-density at the
present epoch, which asymptotes to $m_{\chi}\,>\,32.4\leavevmode\nobreak\
\text{eV}$, almost independently of $\sigma/\lambda$ for $\sigma/\lambda>1$.
However, parametric resonances induce minor deviations from this relation as
represented numerically in Fig. 2.

Constraints for $\sigma/\lambda<1$. In this case, the PSD and corresponding
second moment become sensitive to the IR part. In that case, the constraint
has to be evaluated numerically. As the coupling $\sigma/\lambda$ decreases,
irrespectively of $\sigma/\lambda\ll 1$ the Lyman-$\alpha$ bound asymptotes to
$m_{\chi}\,>2\times 10^{-4}\leavevmode\nobreak\ \text{eV}$. The full numerical
evaluation of the constraints are depicted in Fig. 2, showing a good agreement
with the analytical approximations provided in this section.

##  3 Isocurvature perturbations

The presence of a light scalar spectator field during inflation can lead to
large isocurvature perturbations [7], in conflict with CMB measurements. The
rapid increase in DM energy density is essentially driven by the quadratic
fluctuations that substantially contribute to the variance
$\langle\chi^{2}\rangle$ [7, 5]. Initially assumed to vanish, the DM
misalignment remains $\langle\chi\rangle=0$ throughout inflation. The DM
inhomogeneities do not directly affect the curvature perturbation, and they
can be treated as pure isocurvature fluctuations in the comoving gauge [7].
The second-order contribution to the isocurvature power spectrum is given by
[7, 5]
$\mathcal{P}_{\mathcal{S}}(k)\;=\;k^{3}/(2\pi^{2}\rho_{\chi}^{2})\int\mathop{}\\!\mathrm{d}^{3}\bm{x}\
\langle\delta\rho_{\chi}(\bm{x})\delta\rho_{\chi}(0)\rangle
e^{-i\bm{k}\cdot\bm{x}}\,,$ where $\rho_{\chi}$ and $\delta\rho_{\chi}$ denote
the DM energy density and its fluctuation, respectively. The current
constraints on the isocurvature power spectrum from ${\textit{P}lanck}$ impose
$\mathcal{P}_{\mathcal{S}}(k_{*})\lesssim 8.3\times 10^{-11}$ for the pivot
scale $k_{*}=0.05\,{\rm Mpc}^{-1}$.

To derive an analytical expression for the isocurvature power spectrum, we
consider for simplicity a finite duration of inflation $N_{\text{tot}}=76.5$
from $N_{i}=N_{\text{end}}-N_{\text{tot}}$ until $N_{\text{end}}$ which
corresponds to IR and UV cutoffs $k_{\text{IR}}=a_{i}H\simeq 10^{-10}k_{*}$
and $k_{\text{UV}}=a_{\rm{end}}H\simeq 10^{24}k_{*}$. We further assume a
constant Hubble rate $H$ and dark matter mass. By neglecting time derivatives
of the mode function, the isocurvature power spectrum can be expressed as [8]

$\mathcal{P}_{\mathcal{S}}(k)\;=\;\frac{k^{2}}{(2\pi)^{4}\langle\chi^{2}\rangle^{2}}\int_{0}^{\infty}\mathop{}\\!\mathrm{d}p\
p\int_{|k-p|}^{k+p}\mathop{}\\!\mathrm{d}q\ q\
|\chi_{p}|^{2}|\chi_{q}|^{2}\,,\qquad\text{with}\qquad\langle\chi^{2}\rangle\,\equiv\,\int_{k_{\text{IR}}}^{k_{\text{UV}}}\,\mathcal{P}_{\chi}\,\dfrac{\mathop{}\\!\mathrm{d}k}{k}\,,$
(6)

with $\mathcal{P_{\chi}}\,\equiv\,k^{3}/(2\pi^{2})|\chi_{k}|^{2}$. Assuming
$m_{\chi}^{2}/(3H^{2})<1$ , the mode function on super-horizon scales can be
approximated by $|\chi_{k}(k\ll
aH)|\,\simeq\,(2^{1/3}aH^{1/3})^{-3/2}(k/(aH))^{-\nu}\,,$ with $\nu\simeq
3/2-\beta$ and $\beta\equiv m_{\chi}^{2}/(3H^{2})$. The corresponding power
spectrum $\mathcal{P_{\chi}}\,=H^{2}/(4\pi^{2})(k/(aH))^{2\beta}$ is scale
invariant for $\beta=0$ and blue-tilted when $\beta>0$. Defining
$\Delta\hat{N}\equiv N_{\text{end}}-\hat{N}$ with $\hat{N}$ being the number
of $e$-folds at horizon crossing for a given scale $k=a(\hat{N})H(\hat{N})$,
we can generalize the formula derived in [8] to arbitrary scales $k$ and time-
dependent DM mass

$\mathcal{P}_{\mathcal{S}}(k)\;\simeq\;\dfrac{4\sigma^{2}\hat{\phi}^{4}(N_{\text{tot}}-\Delta\hat{N})}{9\hat{H}^{4}}\,e^{-4\Delta\hat{N}\sigma\hat{\phi}^{2}/(3\hat{H}^{2})}\,\left(1-e^{-2N_{\textrm{tot}}\sigma\hat{\phi}^{2}/(3\hat{H}^{2})}\right)^{-2}\,,$
(7)

where the $\,\hat{}\,$ notation means that background quantities have to be
evaluated at horizon crossing. Generalization to any inflaton-induced DM mass
$m_{\text{eff}}(\phi)$ can be done by substituting
$\sigma\hat{\phi}^{2}\rightarrow m_{\text{eff}}^{2}(\hat{\phi})$. From Eq.
(7), the isocurvature constraint
$\mathcal{P}_{\mathcal{S}}(k_{*})<\beta_{\text{iso}}\mathcal{P}_{\mathcal{R}}(k_{*})\simeq
10^{-11}$ corresponds to $\sigma/\lambda\,>\,0.02$ which is in good agreement
with the fully numerical result shown in Fig. 2 left-panel of Ref. [8] and
depicted in Fig. 2. The available parameter space is represented as a white
region in Fig. 2 while coloured region are excluded. See Ref. [2] for details.

Figure 2: The white space represents the allowed parameter space for the DM
mass and its coupling to the inflaton, in which the relic abundance constraint
can be saturated. Forbidden regions correspond to the overproduction of
isocurvature (green) and the suppression of small structure (red). Requiring
the reheating temperature to be above that required for successful Big Bang
Nucleosynthesis excludes the blue region. Ensuring that perturbative reheating
can be achieved excludes the orange region.

## Acknowledgments

I acknowledge support by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) under Germany’s Excellence Strategy – EXC 2121 “Quantum
Universe” – 390833306.

## References

## References

  * [1] R. Kallosh and A. Linde, JCAP 10 (2013), 033 doi:10.1088/1475-7516/2013/10/033 [arXiv:1307.7938 [hep-th]].
  * [2] M. A. G. Garcia, M. Pierre and S. Verner, Phys. Rev. D 107 (2023) no.4, 043530 doi:10.1103/PhysRevD.107.043530 [arXiv:2206.08940 [hep-ph]].
  * [3] M. A. G. Garcia, M. Pierre and S. Verner, [arXiv:2305.14446 [hep-ph]].
  * [4] L. Kofman, A. D. Linde and A. A. Starobinsky, Phys. Rev. D 56 (1997), 3258-3295 doi:10.1103/PhysRevD.56.3258 [arXiv:hep-ph/9704452 [hep-ph]].
  * [5] S. Ling and A. J. Long, Phys. Rev. D 103 (2021) no.10, 103532 doi:10.1103/PhysRevD.103.103532 [arXiv:2101.11621 [astro-ph.CO]].
  * [6] G. Ballesteros, M. A. G. Garcia and M. Pierre, JCAP 03 (2021), 101 doi:10.1088/1475-7516/2021/03/101 [arXiv:2011.13458 [hep-ph]].
  * [7] D. J. H. Chung, E. W. Kolb, A. Riotto and L. Senatore, Phys. Rev. D 72 (2005), 023511 doi:10.1103/PhysRevD.72.023511 [arXiv:astro-ph/0411468 [astro-ph]].
  * [8] M. A. G. Garcia, M. Pierre and S. Verner, Phys. Rev. D 107 (2023) no.12, 123508 doi:10.1103/PhysRevD.107.123508 [arXiv:2303.07359 [hep-ph]].