11institutetext: The University of Manchester, Alliance Manchester Business
School, Manchester M15 6PB, UK
11email<EMAIL_ADDRESS><EMAIL_ADDRESS>Schlumberger Cambridge Research, Cambridge CB3 0EL, UK

# Cooperative Multi-Agent Search on Endogenously-Changing Fitness Landscapes

Chin Woei Lim 11 Richard Allmendinger 11 0000-0003-1236-3143 Joshua Knowles 22
0000-0001-8112-6112 Ayesha Alhosani 11 Mercedes Bleda 11 0000-0002-0197-0936

###### Abstract

We use a multi-agent system to model how agents (representing firms) may
collaborate and adapt in a business ‘landscape’ where some, more influential,
firms are given the power to shape the landscape of other firms. The
landscapes we study are based on the well-known NK model of Kauffman, with the
addition of ‘shapers’, firms that can change the landscape’s features for
themselves and all other players. Our work investigates how firms that are
additionally endowed with cognitive and experiential search, and the ability
to form collaborations with other firms, can use these capabilities to adapt
more quickly and adeptly. We find that, in a collaborative group, firms must
still have a mind of their own and resist direct mimicry of stronger partners
to attain better heights collectively. Larger groups and groups with more
influential members generally do better, so targeted intelligent cooperation
is beneficial. These conclusions are tentative, and our results show a
sensitivity to landscape ruggedness and “malleability” (i.e. the capacity of
the landscape to be changed by the shaper firms). Overall, our work
demonstrates the potential of computer science, evolution, and machine
learning to contribute to business strategy in these complex environments.

###### Keywords:

Cooperative learning NK models Endogenously-changing landscape Shaping
Searching Adaptation

##  1 Introduction

Most non-trivial social systems are inherently challenging to gauge due to the
potential complexity arising from interactions at both individual and
collective levels [8]. Especially in the business context, the mechanics of
interaction between competing firms (agents) are often based on rather coarse
simplifications and an incomplete understanding of the business landscape. The
sophistry embedded within the interplays of businesses, difficult to
appreciate from the outside, produces counterintuitive resultant behaviours
[7].

Firms compete by developing new strategies, technologies or business models
all of which involve solving complex problems and making a high number of
interdependent choices. To solve these problems, managers need to search their
firms’ business landscapes and find a combination of these choices that allows
them to outperform their competitors. Bounded rational managers cannot easily
identify the optimal combination, and tend to engage in sequential search
processes [13] and, via trial and error, learn and find what combinations are
possible and perform well. Effective search for well-performing solutions in a
business landscape is thus a source of competitive advantage for companies.

Conceptually, a business landscape dictates the effectiveness of a firm’s
search strategy by assigning them a fitness, which typically represents the
level of return or performance. The active revision of a firm’s choices is
crucial in maintaining its competitive advantage, growth and profitability
when competing against other firms on a business landscape. Such revisions are
normally in the form of research and development of any aspect of a firm in
order to find better choices (strategies, methods, and/or products), leading
it towards a better path, and to higher local peaks on the landscape.
Generalising, firms improve their performance by adapting to the business
landscape within which they operate. However, actual business landscapes are
dynamic, and they tend to change not only exogenously as a result of external
factors (changes in government policies and regulations, in demographic and
social trends, etc) but also due to the behaviour and strategies of the firms
competing within them. Firms simply do not limit themselves to only adapting
and accepting the state of their current environment as it is. Capable firms
might be able to _shape_ the business landscape to their advantage (in
addition to _search_ the landscape) [6, 17]. A quintessential example of this
phenomenon was when Apple introduced the iPhone and swiftly shook the
environment in its favour, demolishing Nokia, which was the incumbent cell-
phone market leader at that time.

Management research has used the NK model [9, 11] introduced in the management
literature by [12] to build simulation models to represent business
landscapes, and study different factors that influence the effectiveness of
companies’ search processes (see [2] for a review). Despite the usefulness of
these models, most of them consider that the business landscape within which
companies compete does not change or it changes in an exogenous manner, i.e.
due to factors external to the companies. They thus do not account for the
influence that endogenous changes rooted in firms’ behaviour have in the
business landscape, and in turn in the performance of firms within it. The
first simulation model that has analysed companies search effectiveness when
business landscapes change endogenously was proposed by Gavetti et al in 2017
[10]. The authors extend the NK model to consider two types of firms (agents):
agents that search a landscape only (referred to as _searchers_) and agents
that can both search and shape the landscape (_shapers_). Consequently,
searching firms need to adapt (search) on a landscape that is being shaped
(changed) by the shaping firms. In other words, shapers have the power to
change a landscape endogenously, while searchers perceive these changes as
exogenous. Since all agents (shapers and searchers) search the same landscape,
a change in the landscape (caused by a shaper) affects all agents. The study
of Gavetti et al. [10] focused on studying how the impact of different levels
of landscape ruggedness and complexity, and the proportion of shapers vs
searchers operating in it, affect the performance of both types of firms. In
real contexts, firms do not only compete within a business landscape but in
many cases competing firms try to improve their performance by cooperating,
i.e. via coopetition. Coopetition is the act of cooperation between competing
companies. Businesses that engage in both competition and cooperation are said
to be in coopetition. This paper extends Gavetti et al [10] by allowing firms
to cooperate and analyses how cooperation influence firms performance in
endogenously changing business landscapes. This is achieved by incorporating
cognitive and experiential search into the adaptation process.

The next section details the traditional NK model and the adapted version of
that model by Gavetti et al. [10] (also referred to NKZE model), and explains
the search rules. Section 3 proposes a cooperative approach with learning, and
Section 4 provides details about the experimental study and then analyzes the
proposed approach for different configurations of the simulated changing (and
competitive) business environment. Finally, Section 5 concludes the paper, and
discusses limitations of the work and areas of future research.

##  2 Preliminaries

###  2.1 Kauffman’s NK(C) Model

The NK model of Kauffman [11] is a mathematical model of a tunably rugged
fitness landscape. The ruggedness is encapsulated by the size of the landscape
and the number of local optima, which are controlled by the parameters, $N$
and $K$, respectively. Formally, in an NK model, the fitness $f(\vec{x})$ of
an agent (firm) at location $\vec{g}=(g_{1},\dots,g_{N})$,
$g_{i}\in\\{0,1\\}$, on the landscape can be defined as

$f(\vec{g})=\frac{1}{N}\sum_{i=1}^{N}f_{i}(g_{i},g_{i_{1}},\dots,g_{i_{K}}),$
(1)

where $g_{i}$ is the $i$th (binary) decision variable, and the fitness
contribution $f_{i}$ of the $i$th variable depends on its own value, $g_{i}$,
and $K$ other variable values, $g_{i_{1}},\dots,g_{i_{K}}$. The parameter $K$
has a range of 0 $\leq$ $K$ $\leq$ $N-1$ that determines how many other $K$
different $g_{i}$’s will be affecting each $g_{i}$ when computing fitness. The
relationships between $g_{i}$’s are determined randomly and recorded in an
_interaction matrix_ that shall be left unchanged. The function
$f_{i}:\\{0,1\\}^{K+1}\rightarrow\mathbf{R}$ assigns a value drawn from the
uniform distribution in the range [0,1] to each of its $2^{K+1}$ inputs. The
values ${i_{1}},\dots,i_{K}$ are chosen randomly (without replacement) from
$\\{1,\dots,N\\}$. Increasing the parameter $K$ results in more variables
interacting with each other, and hence a more rugged (epistatic) landscape.
The two extreme cases, $K=0$ and $K=N-1$, refer to the scenarios where the
fitness contributions $f_{i}$ depend only on $g_{i}$ (i.e. each $f_{i}$ can be
optimized independently) and all variables, $g_{1},\dots,g_{N}$, respectively
(maximum ruggedness).

Taking an arbitrary firm with a search policy string of $\vec{g}=(011101)$, we
can calculate the fitness contribution ($f_{i}$) of $g_{1}$ by forming a
temporary string with the $g_{i}$’s that are related to itself by referring to
the interaction matrix. Let us assume that, in this example, $g_{1}$ was
initially and randomly determined to be related to $g_{2},g_{4}$ and $g_{5}$.
Since $g_{1}=0$, $g_{2}=1$, $g_{4}=1$ and $g_{5}=0$, the string formed shall
be $(0110)$. The fitness contribution can then be extracted from the fitness
matrix by taking the value from $6^{th}$ row (0110 in decimal), and the
$i^{th}$ column ($1^{st}$ column in this case). Understandably, the fitness
contributions of subsequent $g_{i}$’s are calculated similarly.

Kauffmann later extended the NK model to introduce coupled landscapes (a.k.a
NKC model) [11], which allows multiple species to exist on different
landscapes, and interact through a phenomena of niche construction.

###  2.2 Gavetti et al.’s NKZE Model

The conventional NK model allows firms to continually adapt on a fixed
landscape until they reach some local or global optima. Additionally, the
action of any firm has no consequence on other competing firms. However, in a
realistic and dynamic business environment, the introduction of disruptive
technologies and concepts can often drastically restructure the business
landscape, thereby needing competing firms to re-strategize towards a new goal
or face obsoleteness.

Gavetti et al. [10] introduced the concept of shapers, which have the ability
to modify the business context to their own advantage on top of the standard
agents (hereinafter known as searchers) in the baseline NK model. They then
studied the effects of different levels of shaping on the performance of
shapers themselves and on searchers. Unsurprisingly, shaper firms and the
level of shaping have great effects on the performance of their competitors as
landscape restructuring always tend to undermine competitor performance.
However, a high level of shaping coupled with a great number of shapers were
found to be highly non-beneficial for shapers and searchers alike, as constant
landscape restructuring changes the objective too fast and too much, thus
rendering local search obsolete and causing massive performance instabilities.

The key feature of shapers is their ability to influence the business context
(hopefully) to its own advantage, and as a side-product, alter the fitness of
their competitors. To achieve this, Gavetti et al. [10] extends the NK model
with an additional $Z$ (binary) decision variables,
$\vec{e}=(e_{1},\dots,e_{Z}),e_{i}\in\\{0,1\\}$. Here, $\vec{e}$ is referred
to as the shape policy string and is globally shared by all firms, differently
to the search policy string, $\vec{g}$, which is controlled by each firm
(agent) independently. Similarly to the parameter $K$ in the NK model, Gavetti
et al. [10] use a parameter $E$ to interlink the shape and search policy
strings: each $g_{i}$ is related to $E$ randomly sampled $e_{i}$’s or
$e_{i_{1}}\dots,e_{i_{E}}$. Such relationships are also recorded in an
interaction matrix that shall be kept constant throughout a run. Notably, $E$
has a range of 0 $\leq$ $E$ $\leq$ $Z$.

Accordingly (for the added dimensions), Gavetti et al. [10] update the fitness
assignment function to $f_{i}:\\{0,1\\}^{K+1+E}\rightarrow\mathbf{R}$, which
assigns a value drawn from the uniform distribution in the range [0,1] to
$2^{K+1+E}$ inputs. In practice, the fitness contributions are stored in a
matrix of size $2^{K+1+E}\times N$, and the interactions are stored in a
matrix of size $N\times K+1+E$. Now, the evaluation of a firm’s fitness
depends on $E$ too. Taking the search policy string for an arbitrary firm to
be $\vec{g}=(011011)$, global shape policy string to be $\vec{e}=(101000)$,
$K=3$ and $E=3$, an interaction matrix was generated randomly, resulting in
$g_{1}$ being related to $g_{2},g_{4},g_{5},e_{2},e_{3}$ and $e_{6}$. The
fitness of contribution of $g_{1}$ ($f_{1}$) will be determined by firstly
forming a temporary string in the order of
$(g_{1}g_{2}g_{4}g_{5}e_{2}e_{3}e_{6})$, making $(0101010)$. Similarly,
$f_{1}$ will be taken as the $42$nd ($0101010$ in decimal) row and the $1$st
column ($ith$ column) of the fitness contribution matrix.

Gavetti et al. [10] use the following approach to tackle the NKZE Model: At
the beginning of the simulation, a predetermined amount of firms will be
turned into shapers based on a shaper proportion ($\beta$) parameter. Firms
are re-ordered randomly at the beginning of each iteration. More specifically,
all firms are allowed to make an action in accordance with the randomly
determined order in each iteration. Thus, the number of actions within an
iteration would be equal to that of the firm (agent) population. In terms of
action, each (shaping and searching) firm is allowed to make one adaption
move. A searching firm flips one randomly selected search policy bit (keeping
the shaping policy unaltered) and if the resulting policy has a better fitness
than the current policy, then the firm retains the new policy; otherwise, the
firm will stick with the old policy. However, when it is the turn of a shaper,
it has the choice of either adapt as a searcher would without altering the
shape policy string, or randomly mutating a single bit of the shape policy
string and evaluate fitness with its original unmutated search policy string.
A shaper will then pick either choice that is better, or end its turn without
adopting any mutation if both the choices were found to be unfit.

Intuitively, $E$ also corresponds to the level of shaping and the malleability
of the fitness landscape. A higher $E$ means that the globally shared shape
policy string has more influence on fitness contributions, and transitively,
the fitness landscape itself. Under this condition, the extent of fitness
landscape restructuring when a shaper acts on the shape policy string is high.
Thus, the fitness landscape is said to be highly malleable at high $E$.

The NKZE model is a variation of Kaufmann’s NKC model. While both models allow
agents to dynamically change the environment, there are critical differences:
in the NKC model (i) each species operate on a separate landscape, (ii) all
species have the ability to change the landscape, and (iii) each species is
represented by one agent only. Consequently, the two models are designed to
simulate different (simplified) business environments.

##  3 Stealthy and Cooperative Learning

The transient environment of the NKZE model caused traditional myopic “hill-
climbing” adaptation to underperform. Additionally, such myopic practices
hardly capture the rationality of realistic firms. Motivated by this weakness,
we extend the methodology to tackle the NKZE model to provide firms with the
ability to learn from one another (stealthily or cooperatively), potentially
allowing firms to adapt more quickly and adeptly to changing environments.

To implement such ideas, a strategy of exploiting multi-agent search in NKZE
with population-based optimisation techniques, specifically particle swarm
optimisation (PSO) [5] and explicit direct memory genetic algorithm [20], was
implemented. This was done by (i) allowing firms to quickly adapt to the
environment by looking towards excellent firms during the exploration phase
following concepts inspired by PSO (similar to neighborhood search [18]) and
(ii) preserving good solutions in a memory and exploiting them at the end of
exploration. We will refer to this strategy as stealthy global learning
(StealthL). StealthL operates in an idealistic environment where intel
regarding the strategies (and success level) of competing firms is always
freely and readily available without limitation (i.e., globally). However,
such limitation does exist and is inherent to the nature of competition.
Additionally, the NKZE and StealthL model do not share similar dynamics, as
the former had a single-mutation restriction whereas the latter allowed for
very rapid adaptation by mutating multiple elements within policies of firms.
The dynamics of NKZE is more realistic as a change in a firm’s policy takes
time and is limited by resources. A complete or near-complete revamp of
policies continuously is not affordable.

As a result, the StealthL model was modified to allow firms to form
collaboration groups. Swarm intelligence and memory scheme were now restricted
within the boundary of these groups, thus limiting the amount of information a
firm gets. This new model will be referred to as the structured cooperation
(StructC) model. Both StealthL and StructC were compared against the standard
adaptation used in the NKZE (hereinafter known as the standard model) in the
next section. First, we will describe StealthL and StructC in the next two
subsections.

###  3.1 Stealthy Global Learning

Our model of how stealthy learning occurs between firms is based on a simple
information-sharing scheme used in the swarm intelligence method, PSO. This is
augmented with the use of a memory of past policies, a technique reminiscent
of poly-ploid organisms’ storage of defunct (inactive) genetic material
(chromosomes) that can be resurrected quickly under environmental stress.

Swarm Intelligence To implement swarm intelligence, the search policy string
of an arbitrary firm was mutated based on a guiding vector that is unique to
each firm [21]. Descriptively, the guiding vector of an arbitrary firm is
$\vec{P}=(p_{1},\ldots,p_{N})$ and has a length of $N$, matching that of the
search policy string $\vec{g}=(g_{1},\ldots,g_{N})$. Each element $p_{i}$ in
the guiding vector represents the probability of its corresponding element
($g_{i}$) in the search policy string to mutate to 1. Naturally, the
probability of which $g_{i}$ mutates to 0 is given by $1-p_{i}$. All guiding
vectors were randomly initialised with a uniform distribution with [0,1] range
at the beginning. At its turn, the firm first learns towards the search policy
string of the current global best performing firm $({g}_{maxf,t})$ at time $t$
using its guiding vector with a learning rate of $\alpha$ where
$0\leq\alpha\leq 1$:

$\vec{P}_{t+1}=(1-\alpha).\vec{P}_{t}+\alpha.\vec{g}_{maxf,t}\;.$

Subsequently, a string of random variates with length $N$ is generated as
$\vec{R}=(r_{1},\ldots,r_{N})$ using a uniform distribution with [0,1] range.
Finally, the new search policy string is determined as follows:

$g_{i}=\left\\{\begin{matrix}1,r_{i}<p_{i}\\\ 0,r_{i}\geq p_{i}\>.\\\
\end{matrix}\right.$

Note that $p_{i}$ has a range of $0.05\leq p_{i}\leq 0.95$ to allow for 5%
random mutation after convergence. This new adaptation was designed to
facilitate fast landscape adaptation via guided multiple mutations. The single
random mutation of the shape policy string was kept without alteration to
preserve the nature of the landscape-shaping dynamics. Finally, the firm
chooses whether to adopt the mutated policy strings or to remain unchanged as
in the standard NKZE model.

Learning from Experience, a.k.a. Polyploidy In addition to swarm intelligence,
the model was also extended to memorise the aggregated policy (search + shape)
of the best performing firms in each iteration [20]. We limit the size of the
database in which these memories are stored to $\Theta$ agents. To ensure
environmental diversity within the database, newly memorized candidate memory
should have a unique shape policy string. If the shape policy string of the
candidate is already present in a memory, the fitter one will be adopted. At
full capacity, replacement can only happen if an environmentally unique
candidate was better than the worse performing memory. To prevent premature
memory exploitation, a parameter $\varepsilon$ representing the probability of
not exploring the database was initialised to 1 at the beginning of the model.
A decay parameter $\gamma<1$ was then set to reduce $\varepsilon$ at the end
of each iteration using $\varepsilon_{t+1}=\varepsilon_{t}.\gamma$.

At the turn of an arbitrary firm, the firm shall only exploit the memory if a
random number generated is greater than $\varepsilon$ without undergoing any
exploration (searching and/or shaping). The firm then, without hesitation,
adopts the best policies of the best performing memory. Note that a searcher
can only adopt the search policy string of the best performing memory, whilst
a shaper adopts both search and shape policy strings simultaneously.

Relevance to practice Coopetition emphasizes the mixed-motive nature of
relationships in which two or more parties can create value by complementing
each other’s activity [3]. In the stealth model, organizations in the
landscape are all competing to reach the best fitness, but also cooperating by
sharing knowledge with each other. This resembles the scenario of
organizations helping each other to reach a new common goal, such as tech
giants Microsoft and SpaceX working together to explore space technology. In
their collaboration, the organizations work together to provide satellite
connectivity between field-deployed assets and cloud resources across the
globe to both the public and private sector via SpaceX’s Starlink satellite
network [16]. At the same time, they are competing to dominate niche segments.
If this collaboration succeeds, Microsoft and SpaceX will be dominating the
space technology market.

###  3.2 Structured Cooperation

By forming random collaboration groups, firms are now equipped to exchange
landscape knowledge amongst their partners, but have no information outside of
the group. Thus, (i) a firm can only refer to the best performing policy
within its group, and (ii) each group has a separate memory. Mutation-wise,
only the element of the search policy string that has the highest probability
of mutating can mutate. If multiple were present, one will be selected at
random to mutate.

Relevance to practice Collaboration between companies can occur in numerous
ways, one of which is sharing knowledge and expertise. In today’s age,
expertise and information are considered valuable strategic assets for
organizations [15]. The StructC model mimics sharing knowledge among a closed
pre-determined group of companies that falls under the same
management/ownership (a.k.a conglomerate). Examples of this type of
corporations are, Alphabet LLC who owns Google, DeepMind; Amazon who owns
Audible, Amazon Fresh, Ring to list a few. Despite sharing their knowledge and
expertise, collaborating companies work towards one goal while maintaining
their independence and decision making [14]. Also, knowledge sharing becomes
the natural required action for the company to reduce costs and save time, and
improve efficiency [22].

##  4 Experimental study

This section outlines the model and algorithm parameter settings followed by
an analysis of results. Table 1 provides an overview of the default parameter
settings for the the three models to be investigated, Standard Model, StealthL
and StructC. Parameter $N$, $Z$ and $\beta$ were chosen in accordance to [10]
for comparison purposes. We simulate a business environment with a population
of $M=10$ agents (firms) and a maximum group size of $\omega_{max}=4$. Table 2
lists the possible group compositions. The default parameters for the other
parameters were set based on preliminary experimentation such that robust
results are obtained.

$\omega$ = 1 |  $\omega$ = 2 |  $\omega$ = 3 |  $\omega$ = 4
---|---|---|---
| 1 searcher
---
($\beta$ = 0)
| 2 searchers
---
($\beta$ = 0)
| 3 searchers
---
($\beta$ =0)
| 4 searchers
---
($\beta$ =0)
| 1 shaper
---
($\beta$ = 1)
| 1 searcher 1 shaper
---
($\beta$ = 0.5)
| 2 searchers 1 shaper
---
($\beta$ = 0.33)
| 3 searchers 1 shaper
---
($\beta$ = 0.25)
- |  | 2 shapers
---
($\beta$ = 1)
| 1 searcher 2 shapers
---
($\beta$ = 0.67)
| 2 searchers 2 shapers
---
($\beta$ = 0.5)
- | - |  | 3 shapers
---
($\beta$ = 1)
| 1 searcher 3 shapers
---
($\beta$ = 0.75)
- | - | - |  | 4 shapers
---
($\beta$ = 1)

Table 2: StructC group combinations.

###  4.1 Experimental results

We investigate the performance of various models as a function of the problem
complexity. We achieve this by visualizing and analysing the model behaviors
during the search process. All models were validated using 50 runs with each
run using a randomly generated fitness landscapes (same set of landscapes were
used for each model) and lasting for 100 iterations. For StructC, group
compositions were randomly sampled to preserve generality. As a result, it was
ensured that each unique group composition will have appeared 50 times
throughout the whole experiment. Each run takes around $2$ minutes depending
on problem complexity using an Intel i7 (8th gen.) CPU, 8 GB DDR3L RAM.111The
code to replicate these experiment can be downloaded at
https://github.com/BrandonWoei/NK-Landscape-Extensions

Standard vs. Stealthy Global Learning Figures 1\- 3 compare the performance of
searchers and shapers for the Standard Model and StealthL as well as for
different learning rates ($\alpha$). Following observations can be made:

Shapers and searchers in StealthL outperform their counterparts in the
Standard Model significantly for rugged landscape ($K>0$) and regardless of
the learning rate $\alpha$. However, for $K=0$, the Standard Model achieved a
better performance because StealthL is suffering from premature convergence
caused by its weakened random perturbation, a trade-off of guided learning.
Stronger mutation is necessary when peaks on the fitness landscape are rare
(or a single peak is present as is the case for $K=0$) and sufficiently far
from one another since StealthL becomes complacent to the point in which
exploration is inhibited when its corresponding agents all have roughly good
and near solutions.

For StealthL, the performance of searchers and shapers is almost identical
regardless of the level of ruggedness $K$. For the Standard Model, the
performance gap between searchers and shapers depends on $K$ with the most
significant performance gap being observed for an intermediate level of
landscape ruggedness.

The learning rate $\alpha$ has a significant impact on search performance. A
decrease in $\alpha$ leads to a slower convergence but an improved final
performance (if given sufficient optimization time). A high learning rate of
$\alpha>0.8$ leads to premature convergence. Generally, slightly higher
learning rates perform better as the level of landscape ruggedness and
malleability increases. Searchers and shapers are affected in a similar
fashion by a changing learning rate, while the performance gap between
searchers and shapers for a specific learning rate is minimal.

Figure 1: Average best fitness (and 95% confidence interval in shading
calculated using $t$ statistic) obtained by the Standard Model and StealthL at
(left plot) minimal ($K=0$, $E=0$), (middle plot) intermediate ($K=5$, $E=6$),
and (right plot) high ($K=11$, $E=12$) landscape ruggedness and malleability.
StealthL Searcher and Stealth Shaper have almost identical performance.
Figure 2: Varying StealthL learning rate $\alpha$ at extreme malleability
($E=12$) with (left plot) minimal ($K=0$) (middle plot) intermediate ($K=5$)
(right plot) high ($K=12$) landscape ruggedness. The confidence interval of
the average fitness value has been omitted in this and following plots for the
sake of readability. Figure 3: Varying StealthL learning rate $\alpha$ at
extreme ruggedness ($K=11$) with (left plot) minimal ($E=0$) (middle plot)
intermediate ($E=6$) (right plot) high ($E=12$) landscape malleability. Figure
4: Average fitness values achieved by StructC for different group sizes at low
ruggedness and malleability ($K=0$, $E=0$). Figure 5: Average fitness values
achieved by StructC for different group sizes at high ruggedness and
malleability ($K=11$, $E=12$).

###  4.2 StructC Model at Low Landscape Ruggedness and Malleability

Figure 4 visualizes the impact of different group sizes on structured
cooperation for $K=0$, $E=0$. Evidently, being in a bigger group with shapers
helped with a searcher’s performance. Lone searchers suffered from severe
premature convergence with no sign of any improvement. The advantage of having
a shaper was especially prominent in a twin group. However, the relationship
between shaper proportion and searcher performance was highly inconspicuous
and nonlinear. A group of size 3 gave the best searcher performance when it
only had 1 shaper. The further increase in shapers decreased searcher
performance but was still better than the complete searcher group. The results
for size 4 groups were even more astonishing, with searchers only groups
outperforming groups with 1 or 2 shapers in terms of searcher performance.
However, a group with 3 shapers enabled its only searcher to outperform
complete searcher groups by a noticeable degree. Whilst shapers are useful to
searchers, the benefits were found to be highly dependent on group size and
the number of shapers.

Albeit outperforming a lone searcher, a lone shaper suffers a similar
premature convergence problem. For twin groups, a complete shaper group gave
better shaper performance than a mixed group. For size 3 groups, a complete
shaper gave the best shaper performance. In terms of the mixed groups, groups
with 1 searcher were found give better shaper performance than groups with 2
searchers. For size 4 groups, shaper dominated groups gave better shaper
performance than the balanced and searcher dominated groups. Unexpectedly, a
shaper dominated group with 1 searcher gives slightly better shaper
performance than a complete shaper group. Generally, a large group size
dominated by shapers ($>75\%$) was found to be group configuration that
produced excellent searchers and shapers. Predictably, all groups were not
able to surpass the standard and StealthL model performance under all
conditions.

###  4.3 StructC Model at High landscape Ruggedness and Malleability

Figure 5 visualizes the impact of different group sizes on structured
cooperation for $K=11$, $E=12$. At high ruggedness and malleability, lone
searchers had no improvement from its initial fitness. Mixed twin groups were
only able to marginally improve their initial searcher performance. Searcher
performance of searcher dominated twin groups was seen to deteriorate slightly
with time. Such temporal deterioration effects were also seen in all complete
searcher groups regardless of group size. Universally, the increase in group
size and shaper proportion improved searcher performance. At maximum group
size, searcher performance was able to exceed that of the standard model by a
noticeable degree. On the other hand, lone shapers found it difficult to
improve themselves. Similarly, a larger group size with more shapers led to
excellent shaper performance.

##  5 Conclusion

We have used a multi-agent system to model how agents (firms) may collaborate
and adapt in a business ‘landscape’ where some, more influential, firms are
given the power to _shape_ the landscape of other firms. We have found that
shapers outperform searchers under all landscape conditions. However,
excessive landscape reshaping can lead to poor collective performance due to
the instability it introduces. Additionally, both searchers and shapers
perform best, even under dynamic business landscapes, when they can keep their
organisational complexities at a minimum by reducing relationships between
elements of their policies (encoded here by the two parameters $K$ and $E$).
Complex organisations can find it hard to cope with landscape changes
especially when the changes are frequent and substantial. However, we found
that this can be overcome to some extent by allowing for collaboration via
experience-sharing between firms. Whilst mutual learning is beneficial, direct
mimicry of best practices can lead to a reduction in collective knowledge as
firm’ shared inertia hinders exploration, thereby weakening the synergistic
effects of collaboration. Lastly, the positive effects of collaboration were
also found to be at their best when collaboration groups are dominated by
shapers.

Despite having extended NK models, a limitation of our work is that reality is
more complex than the abstraction considered here, particularly in terms of
the kinds of strategy different firms might employ. Moreover, the model could
be made more realistic by, for example, accounting for additional objectives
and constraints, observing delays between evaluating an objective/constraint
function and having its value available [1, 4, 19], and simulating more than
two types of agents as firms may vary widely e.g. in terms of capabilities,
goals, and how they engage with other firms.

## References

  * [1] Allmendinger, R., Handl, J., Knowles, J.: Multiobjective optimization: When objectives exhibit non-uniform latencies. European Journal of Operational Research 243(2), 497–513 (2015)
  * [2] Baumann, O., Schmidt, J., Stieglitz, N.: Effective search in rugged performance landscapes: A review and outlook. Journal of Management 45(1), 285–318 (2019)
  * [3] Bonel, E., Rocco, E.: Coopeting to survive; surviving coopetition. International Studies of Management & Organization 37(2), 70–96 (1997)
  * [4] Chugh, T., Allmendinger, R., Ojalehto, V., Miettinen, K.: Surrogate-assisted evolutionary biobjective optimization for objectives with non-uniform latencies. In: Proceedings of the 2018 Annual Conference on Genetic and Evolutionary Computation. pp. 609–616 (2018)
  * [5] Engelbrecht, A.: Particle swarm optimization. In: Proceedings of the Companion Publication of the 2015 Annual Conference on Genetic and Evolutionary Computation. p. 65–91. GECCO Companion ’15, Association for Computing Machinery, New York, NY, USA (2015). https://doi.org/10.1145/2739482.2756564, https://doi.org/10.1145/2739482.2756564
  * [6] Felin, T., Kauffman, S., Koppl, R., Longo, G.: Economic opportunity and evolution: Beyond landscapes and bounded rationality. Strategic Entrepreneurship Journal 8(4), 269–282 (2014)
  * [7] Forrester, J.: Urban dynamics. IMR; Industrial Management Review (pre-1986) 11(3), 67 (1970)
  * [8] Forrester, J.: Counterintuitive behavior of social systems. Theory and decision 2(2), 109–140 (1971)
  * [9] Ganco, M., Agarwal, R.: Performance differentials between diversifying entrants and entrepreneurial start-ups: A complexity approach. Academy of Management Review 34(2), 228–252 (2009)
  * [10] Gavetti, G., Helfat, C.E., Marengo, L.: Searching, shaping, and the quest for superior performance. Strategy Science 2(3), 194–209 (2017)
  * [11] Kauffman, S.A.: The origins of order: Self-organization and selection in evolution. Oxford University Press, USA (1993)
  * [12] Levinthal, D.: Adaptation on rugged landscapes. Management Science 43(7), 934–950 (1997)
  * [13] March, J., Simon, H.: Organizations. New York: Wiley (1958)
  * [14] Rakowski, N., Patz, M.: An Overview and Analysis of Strategic Alliances on the Example of the Car Manufacturer Renault. GRIN Verlag (2009)
  * [15] Sen, Y.: Knowledge as a valuable asset of organizations: Taxonomy, management and implications. Management Science: Foundations and Innovations p. 29–48 (2019)
  * [16] Sokolowsky, J.: Azure space partners bring deep expertise to new venture (2020), https://news.microsoft.com/transform/azure-space-partners-bring-deep-expertise-to-new-venture/#:~:text=Microsoft%20is%20working%20with%20SpaceX,via%20SpaceX’s%20Starlink%20satellite%20network.
  * [17] Uzunca, B., Rigtering, J.C., Ozcan, P.: Sharing and shaping: A cross-country comparison of how sharing economy firms shape their institutional environment to gain legitimacy. Academy of management discoveries 4(3), 248–272 (2018)
  * [18] Wang, H., Sun, H., Li, C., Rahnamayan, S., Pan, J.s.: Diversity enhanced particle swarm optimization with neighborhood search. Information Sciences 223, 119–135 (2013)
  * [19] Wang, X., Jin, Y., Schmitt, S., Olhofer, M., Allmendinger, R.: Transfer learning based surrogate assisted evolutionary bi-objective optimization for objectives with different evaluation times. Knowledge-Based Systems 227, 107190 (2021)
  * [20] Yang, S.: Explicit memory schemes for evolutionary algorithms in dynamic environments. In: Yang, S., Ong, Y.S., Jin, Y. (eds.) Evolutionary Computation in Dynamic and Uncertain Environments, pp. 3–28. Springer Berlin Heidelberg, Berlin, Heidelberg (2007), https://doi.org/10.1007/978-3-540-49774-5_1
  * [21] Yang, S.: Learning the dominance in diploid genetic algorithms for changing optimization problems. In: Proceedings of the 2nd International Symposium on Intelligence Computation and Applications (2007)
  * [22] Ye, J.: The relations between knowledge sharing among technology centers within a conglomerate, its r&d synergy and performance: A case study of dongfeng motor corporate. In: International Conference on Global Economy and Business Management (GEBM 2019). pp. 43–53 (2019)


# Classifying a frequently repeating fast radio burst, FRB 20201124A, with
unsupervised machine learning

Bo Han Chen1, Tetsuya Hashimoto2, Tomotsugu Goto1, Bjorn Jasper R. Raquel3
<EMAIL_ADDRESS>Yuri Uno2, Seong Jin Kim1, Tiger Y.-Y.
Hsiao4, and Simon C.-C. Ho5
1Institute of Astronomy, National Tsing Hua University, No. 101, Section 2,
Kuang-Fu Road, Hsinchu City 30013, Taiwan
2Department of Physics, National Chung Hsing University, No. 145, Xingda Rd.,
South Dist., Taichung, 40227, Taiwan
3Department of Earth and Space Sciences, Rizal Technological University, Boni
Avenue, Mandaluyong, 1550 Metro Manila, Philippines
4Department of Physics and Astronomy, Johns Hopkins University, Baltimore, MD
21218, USA
5Research School of Astronomy and Astrophysics, The Australian National
University, Canberra, ACT 2611, Australia

(Accepted 2023 March 20. Received 2023 March 17; in original form 2022 October
28)

###### Abstract

Fast radio bursts (FRBs) are astronomical transients with millisecond
timescales. Although most of the FRBs are not observed to repeat, a few of
them are detected to repeat more than hundreds of times. There exist a large
variety of physical properties among these bursts, suggesting heterogeneous
mechanisms of FRBs. In this paper, we conduct a categorisation on the
extremely frequently repeating FRB 20201124A with the assistance of machine
learning, as such techniques have the potential to use subtle differences and
correlations that humans are unaware of to better classify bursts. The
research is carried out by applying the unsupervised Uniform Manifold
Approximation and Projection (UMAP) model on the FRB 20201124A data provided
by Five-hundred-meter Aperture Spherical radio Telescope (FAST). The algorithm
eventually categorises the bursts into three clusters. In addition to the two
categories in previous work based on waiting time, a new way for
categorisation has been found. The three clusters are either high energy, high
frequency, or low frequency, reflecting the distribution of FRB energy and
frequency. Importantly, a similar machine learning result is found in another
frequently repeating FRB20121102A, implying a common mechanism among this kind
of FRB. This work is one of the first steps towards the systematical
categorisation of the extremely frequently repeating FRBs.

###### keywords:

(transients:) fast radio bursts – radio continuum: transients – methods: data
analysis

††pubyear: 2022††pagerange: Classifying a frequently repeating fast radio
burst, FRB 20201124A, with unsupervised machine learning–A

##  1 Introduction

Fast radio bursts (FRBs) are astronomical transients that show sudden
brightening at radio wavelengths (e.g. Lorimer et al. 2007). Most of the FRBs
happen in milliseconds timescale (e.g. Petroff et al., 2016) and occur at
cosmological distances. FRBs are usually classified as either non-repeating or
repeating, where they are defined as one-off bursts and repeated bursts
detected from single FRB sources, respectively. Numerous FRB theories have
been proposed (e.g. Platts et al., 2019) to account for the physical origin of
either non-repeating or repeating FRBs. The proposed hypotheses include
neutron star mergers (Yamasaki et al., 2017), pulsar–black hole interactions
(Bhattacharyya, 2017), active galactic nuclei activities (Vieyro et al.,
2017), etc. In spite of the intensive observations and modelling, so far, the
origin of FRBs remains a mystery.

Most of the repeating FRB sources only have a handful of observed burst events
(Masui & CHIME/FRB Collaboration 2021). However, a highly active repeater
could have more than 1000 bursts detected. The well-known FRB20121102 (Li et
al., 2021) and FRB 20201124A (Xu et al., 2022) are both examples of a highly
active repeater also known as frequently repeating FRBs. This kind of FRBs are
rarely seen, but the multitude of their bursts allows us to perform
categorisation with the aid of machine learning.

FRB 20201124A is the most active FRB known so far (Xu et al., 2022). A radio
observational campaign of Five-hundred-meter Aperture Spherical radio
Telescope (FAST) is dedicated to monitoring FRB 20201124A which provided the
necessary data for this research. The campaign monitored FRB 20201124A from
April 1, 2021 up until June 11, 2021, which cumulatively has 96.9 hours of
observation time. The released catalogue contains 1863 bursts, and each burst
is characterised by at most 27 physical properties (Xu et al., 2022).

These 1863 bursts are detected with a signal-to-noise ratio (S/N) > 7, with
913 bright bursts reaching S/N > 50\. The burst rate of FRB 20201124A evolves
slowly between $5.6\ $\mathrm{h}\mathrm{r}^{-1}$$ to $45.8\
$\mathrm{h}\mathrm{r}^{-1}$$. The burst flux density ranges from 0.005 Jy to
11.5 Jy, and the inferred isotropic luminosity spans from $5\times 10^{37}\
$\mathrm{e}\mathrm{r}\mathrm{g}\ \mathrm{s}^{-1}$$ to $3\times 10^{40}\
$\mathrm{e}\mathrm{r}\mathrm{g}\ \mathrm{s}^{-1}$$ (Xu et al., 2022). The data
provided by the observation is homogeneous, therefore analysing it with
unsupervised machine learning is a feasible direction.

There are already several works that applied machine learning to study FRBs.
Most recently, Hewitt et al. 2022 used a supervised deep-learning algorithm to
classify FRB candidates from the FRB121102 data which is detected by the
Arecibo Telescope. Before that Agarwal (2020) also presented a classification
for FRB candidates from Australian Square Kilometre Array Pathfinder (ASKAP)
detection using deep neural networks. In an earlier work, Wagstaff et al.
(2016) trained and deployed a machine learning classifier that marks each
detection as either a known pulsar, artefact due to interference, or potential
new FRB discovery. Suggesting that unsupervised machine learning is a
promising way to identify misclassified FRB repeaters. For example, Chen et
al. (2021) successfully uncloaked the hidden repeating FRB candidates from
non-repeating FRBs in the CHIME/FRB observation with UMAP, an unsupervised
machine learning model. Thus, motivating subsequent FRB studies such as Kim et
al. 2022 and Hashimoto et al. 2022.

In this paper, an unsupervised dimensionality reduction algorithm is used and
followed by an unsupervised clustering algorithm. The dimensionality reduction
algorithm we utilised is known as Uniform Manifold Approximation and
Projection (UMAP), a technique for dimension reduction based on manifold
learning techniques and ideas from topological data analysis (McInnes et al.,
2018). The UMAP algorithm maps the observational parameters of each FRB to a
2D embedding plane after training on the features of the training samples.
This projection aggregates the FRB samples having similar physical properties
and repels those that do not. As a result, the algorithm provides us with a
way to systematically classify the bursts of FRB 20201124A. Next, an
unsupervised clustering algorithm known as K-means is used to label the
clusters resulting from the projection.

In conclusion, we find that the bursts in FRB 20201124A could be roughly
classified into 3 clusters, rather than the previously known two (Xu et al.,
2022). The classification is mainly dominated by the energy and frequency
parameter of the bursts, while both the arrival time and the waiting time do
not have significant effect on the clustering result. These findings are
different from the bimodal distribution of waiting time proposed in Xu et al.
(2022). This work is a first step towards the systematic categorisation of the
extremely frequently repeating FRBs.

This work is organised as follows. We describe our data composition and model
configuration in Section 2. Our UMAP Model classification result is described
in Section 3. We present the Discussion in Section 4, followed by Conclusions
in Section 5.

##  2 Data composition and model configuration

###  2.1 Sample and Data Selection

The data for FRB 20201124A in this work is provided by the observational
campaign of FAST (Xu et al., 2022). The catalogue records 1863 FRBs at a
frequency range between 1.0 GHz to 1.5 GHz from April 1, 2021 up until June
11, 2021. There are 27 observational parameters available in the catalogue. We
included 8 parameters in the main unsupervised machine learning experiment,
where all bursts must have the measurements of these 8 features. 118 samples
are excluded due to the missing parameter values, thus finally 1745 FRB are
utilised in this paper. We further discuss the details of each observational
parameter in Section 2.1.1.

####  2.1.1 The observational parameters

A total of 8 observational parameters are included in our unsupervised
training. All of these parameters come from the original FRB 20201124A data
provided by FAST. Brief descriptions of the parameters are as follows (see Xu
et al., 2022, for details).

  * •

Barycentrical arrival time (BAT)  : The Barycentrical arrival time of the FRBs
in units of Modified Julian Date (MJD), measured from the centroid of the
best-matched boxcar filter. The value of this parameter is ranged from
59307.33458 to 59360.17901.

  * •

S/N : The signal-to-noise ratio of the FRBs. The value of this parameter is
ranged from 7.14 to 2885.83.

  * •

Peak flux density $\mathrm{(Jy)}$ : The peak flux of the FRBs. The value of
this parameter is ranged from 0.005 to 11.505.

  * •

Fluence $\mathrm{(Jy}\cdot\mathrm{ms)}$ : The integration of the pulse flux
density with respect to the time of the FRBs. The value of this parameter is
ranged from 0.0173 to 67.3181.

  * •

Equivalent width $\mathrm{(ms)}$ : A measure of the time width of FRBs,
defined by the fluence divided by the pulse peak flux density. The value of
this parameter is ranged from 1.556 to 28.36.

  * •

Bandwidth $\mathrm{(MHz)}$ : The bandwidth of the FRBs. The value of this
parameter is ranged from 30.151 to 485.84.

  * •

Peak frequency $\mathrm{(MHz)}$ : The peak frequency of the FRBs. The value of
this parameter is ranged from 1000 to 1500.

  * •

Waiting time $\mathrm{(s)}$ : The time difference between each FRB and its
next one. We impose the log scale on the value for training since it spans
several orders of magnitude. We removed 118 FRB samples at the end of each
observation session because their waiting time cannot be measured. The value
of this parameter ranges from 0.004 to 1325.108 seconds.

It is worth noting that, as the FRBs in this paper originated from a single
source, the cosmological effect is not opposing discrepant effect on them. As
a result, we don’t consider the cosmological effect correction in this paper.

####  2.1.2 The Statistical Information Regarding the Parameters

Our data is comprised of the 1745 FRBs originating from FRB 20201124A, and
this sample is the foothold of our unsupervised machine learning model. In
order to further understand the basic composition of our research, we plot the
distribution of the parameters mentioned in Section 2.1.1, shown in Fig. 1.

Figure 1: This figure shows the observational parameter distributions, which
are the input data we apply for our unsupervised learning. See Section 2.1.1
for the details.

###  2.2 Unsupervised Machine Learning Model and Configuration

We utilised 1745 FRBs provided by the FAST observation campaign, and there are
8 physical parameters for each FRB. Serving as our input data, we managed to
conduct two unsupervised machine learning analyses on this data. The applied
machine learning techniques are UMAP, an algorithm based on topological data
analysis and manifold learning, and K-means, an algorithm for unsupervised
clustering.

The training of UMAP involves two stages. In the first stage, the algorithm
constructs the local connectivity of the data manifold by computing the
distances of each data point’s $k$ nearest neighbours under a relative scale
of the distance to nearest neighbours. As a result, $k$ is one of the UMAP
hyperparameters, referred to as n_neighbors. n_neighbors manipulates the
model’s balance between the local and the global manifold structure of the
data. A large n_neighbors demands UMAP to consider a large cluster of data
points when computing the local connectivity. On the other hand, a small
n_neighbors allows the algorithm to concentrate on the local structure of the
data. We discuss our choice of n_neighbors in Sec. 3.

The second stage of UMAP is to map the FRBs (represented by 8 parameters) to
low-dimensional representations. Thus, we need another hyperparameter
n_components to decide the resulting dimensionality of the reduced dimension.
In our work, we adopt n_components = $2$, which means we project the physical
properties of FRBs onto a 2D plane.

To find a 2-dimensional representation that matches the topological structure
of the 8-dimensional data, UMAP performs a stochastic gradient descent by a
specific cross-entropy function. The gradient descent provides an attractive
force between the points where the scaled distance mentioned above is short in
high dimensional space and provides a repulsive force between the points
whenever the distance is large. In order to prevent the resulting low
dimensional projection from clumping together, another hyperparameter min_dist
is induced to constrain the minimum Euclidean distance between the projected
points. However, this hyperparameter is not affecting our result
significantly. In this paper, min_dist is set to have the default value of
$0.1$.

The result provided by UMAP is followed by the unsupervised clustering method
K-means. It is used to group data into clusters based on similarity. It is one
of the most widely used clustering algorithms and is relatively simple to
implement. The algorithm works by first initializing a set of cluster centres,
then iteratively assigning data points to the closest cluster centre, and
updating the cluster centres to be the mean of the points assigned to that
cluster. This process is repeated until the cluster centres stop changing or a
maximum number of iterations is reached. The main advantage of the K-means
algorithm is its speed and simplicity, but it requires a preset cluster number
n_cluster to determine the total cluster numbers. We discuss our choice of
n_cluster in Sec. 3 and Appendix A.

This section concludes with a simple introduction to the UMAP algorithm and
the hyperparameter setting we use. Readers can refer to McInnes et al. (2018)
for the mathematical background of UMAP.

##  3 Result

###  3.1 UMAP Model Classification Result

Our unsupervised UMAP model makes the projection and then K-means identify the
clusters. We expect the projection of the FRBs to have distinct clusters on
the embedding plane. In order to do that, we need to properly set up the UMAP
hyperparameter n_neighbors. Thus, we evaluate the silhouette score (Rousseeuw,
1987) result from each n_neighbors, presented in Fig. 2. The silhouette score
is a metric for evaluating the performance of a clustering algorithm. The
silhouette score ranges from -1 to 1, with a score of 1 indicating that all
the samples are perfectly matched to their own clusters, and a score of -1
indicating that the samples are poorly matched to their own clusters. The
highest silhouette score is obtained when n_neighbors = 7. Therefore, we use
n_neighbors = 7 in the paper.

K-means clustering algorithm requires a predetermined value for the
hyperparameter n_cluster to decide the number of the resulting clusters. In
order to select the best value, we evaluate the silhouette score of each
n_cluster setting, presented in Fig. 3. In the figure, we see the clustering
has the best silhouette score in n_cluster = 2. However, the embedding plane
obviously has more than two clusters. Therefore, we exclude n_cluster = 2 from
further analysis. Except for n_cluster = 2, n_cluster = 3 provides the best
silhouette score. Therefore, we apply n_cluster = 3, separating the embedding
plane into 3 areas.

We present the result of our unsupervised machine learning in Fig. 4.
Aggregation of the FRBs on the embedding plane implies a similar physical
property. The data points on the plane have been classified by the K-means
clustering algorithm. We denote the three clusters as Cluster 1, Cluster 2 and
Cluster 3. Cluster 2 is large and sparse, including more than half of the
samples. In contrast, the structure of cluster 1 and cluster 3 are denser.
This fact implies the clusters on the right-hand side have similar physical
parameters. The visual inspection suggests that Cluster 1 could be composed of
two sub-clusters. However, due to the restrictions of the clustering
algorithm, we are not able to separate these two sub-clusters.

The clusters presented on the embedding plane suggest that there is indeed
heterogeneity in the behaviour of FRB 20201124A. As a result, we look into the
physical property distribution on the embedding plane in the following section
and figure out which factors contribute to the formation of these clusters.

Figure 2: The silhouette score result from each n_neighbor setting for the
UMAP algorithm. Figure 3: The silhouette score result from each n_cluster
setting for the K-means clustering algorithm. Figure 4: The unsupervised
machine learning result of the FRB 20201124A. The distances among the points
represent the similarities of their physical parameters, but the $x$ and $y$
axes do not have explicit physical implications. 1745 FRBs are projected to
the embedding plane according to their physical properties. Classification is
carried out with the K-means clustering algorithm, and the resulting clusters
are marked in different colours.

###  3.2 Physical Property Distribution of Each cluster

We colour the embedding plane according to the values of the observational
parameters as shown in Fig. 5. In the figure, we see S/N, peak flux, fluence,
bandwidth and peak frequency present clear correlation with each other,
implying these are the factors that significantly influenced the machine
learning results. The distributions of BAT, waiting time and equivalent width
are almost random, indicating these parameters have very little to no effect
on the result.

According to Fig. 5, the parameters can be classified into 3 groups:

  * •

Parameter group A: S/N, peak flux and fluence

  * •

Parameter group B: bandwidth and peak frequency

  * •

Parameter group C: BAT, waiting time and equivalent width

Parameter group A are the factors highlighting cluster 3. Due to these
parameters, bursts in cluster 3 generally have higher S/N, peak flux and
fluence. These factors are energy-related. Thus, cluster 3 is the high-energy
cluster. Parameter group B accentuates the wavelength property of the FRBs.
cluster 1 has high bandwidth and peak frequency, while cluster 2 is exactly
the opposite. Intriguingly, cluster 3 includes FRBs with both low and high
values of bandwidth and peak frequency. Lastly, parameter group C consists of
the parameters not influential for the training. BAT, waiting time and
equivalent width are the time factors, suggesting our dimensionality reduction
result almost is time-independent.

Figure 5: The value distribution of each FRB on the embedding plane. The
numerical scale is the same as Fig. 4. The range of the color bar in each
panel is limited to $\pm 3$ standard deviations from the average for clear
visualization.

##  4 Discussion

###  4.1 Comparing with the UMAP result on multi-FRB source (CHIME/FRB)
catalogue.

The analysis in this paper mainly focuses on the 1745 FRBs of FRB 20201124A
which are bursts that originates from a single repeating FRB source. In this
section, we compare it to the analysis result based on the FRB catalogue with
mostly different sources. The catalogue we referred to here was released by
Canadian Hydrogen Intensity Mapping Experiment Fast Radio Burst (CHIME/FRB)
Project. It includes 536 FRB sources at a frequency range between 400 and 800
MHz from July 25, 2018 to July 1, 2019 (Masui & CHIME/FRB Collaboration 2021).

The analysis done for the CHIME catalogue is also based on UMAP, which is
originally presented in the work of Chen et al. (2021). Their work identified
the hidden repeating FRBs from the currently one-off FRBs. We transplant the
analysis and present it in Fig. 6. In order to provide a straightforward
comparison, the x and y-axes of the two figures (Fig. 4 and Fig. 6) are
presented on the same scale. In comparison with our result in Fig. 4, the
clusters in Fig. 6 are more distinct and dense.

The FRBs in Fig. 4 all originated from one single source, while Fig. 6 are the
FRBs aggregated from 535 sources. As a result, we speculate the different
clumping level comes from the different levels of heterogeneity of the
underlying data. When the data includes only one single source, the data
points are more similar thus the clusters are less distinct. In contrast, when
the FRB sources are diverse, the UMAP model can efficiently divide the
significantly different FRBs and group FRBs that has some degree of
similarity.

In summary, from this comparison, we can infer that the heterogeneity among
the FRBs of different sources is higher than that of a single source. The FRBs
within FRB 20201124A have different patterns, but the prominence of the
pattern difference is definitely lower than the one from the multi-source
catalogue.

Figure 6: The UMAP analysis based on the CHIME/FRB catalogue, in which FRBs
are originated from mostly different sources.

###  4.2 Comparing with the UMAP result on another frequently repeating FRB

FRB 20201124A is not the only observed frequently repeating FRB. FRB20121102
(Li et al., 2021) is a source detected to have more than 1000 bursts as well.
Raquel et al. (2022 submitted) carried out the UMAP analysis for FRB20121102,
which is quite similar to this work. In their work, they included
Barycentrical arrival time, dispersion measure, time width, bandwidth, peak
flux density, fluence, and energy of the FRBs as the input parameters.
Comparing our work, their analysis is presented in Fig. 7. Their result is
produced by setting the hyperparameter values for n_neighbors $=6$ and
min_dist $=0$ and clustered by the HDBSCAN algorithm. In this section, we are
comparing the best UMAP result of our source and theirs. Since the source data
are different, their hyperparameter setting is different from ours. Instead of
doing a comparison on identical hyperparameter settings and clustering
algorithms, it would be more common in the machine learning community to
compare two results in their respective best scenarios (e.g. Badia et al.
2020).

We see that despite the source is different, the UMAP analysis produced
results which are similar in nature. Both of them have three main clusters,
appearing in 1 large, sparse cluster and 2 small, dense clusters. Moreover, in
the lower right corner of both panels, there is a smaller cluster not
specified by the clustering algorithm. In summary, the two results are highly
matched. The embedding planes reflect the similarity of the FRB pattern. As a
result, this fact could be evidence that the frequently repeating FRBs are
powered by a similar mechanism.

Figure 7: The UMAP result comparison between FRB 20201124A and FRB20121102.
The x and y-axes of the two figures are presented on a similar scale.

##  5 Conclusions

In this work, we applied unsupervised machine learning analysis to the 1745
bursts of FRB 20201124A. For each burst event, we included Barycentrical
arrival time, signal-to-noise ratio, peak flux, fluence, equivalent width, and
bandwidth of signal as the input data for the algorithm which then projects
the data points onto a 2D plane based on their similarity. Our work resulted
in three main conclusions:

(i) Our UMAP model divides the FRBs into 3 clusters. The classification is
based on either energy or frequency factors. The 3 clusters are either high
energy, high frequency, or low frequency. Interestingly, time-related factors
did not have a significant effect on our results.

(ii) We compare our result with the work of (Chen et al., 2021) which also has
UMAP result based on multi-FRB source data. The multi-FRB source data result
has a clumpier, distinct embedding plane. We infer that the heterogeneity
among the FRBs from different sources is higher than that of a single source.

(iii) We compare our result with the UMAP result based on another frequently
repeating FRB, FRB20121102 (Raquel et al. 2022 submitted). Impressively, their
result is very similar to our result, including the scale of the clusters and
the relative position of the clusters. This fact suggests that these two
frequently repeating FRBs are powered by similar mechanisms.

## Acknowledgements

We sincerely thank the anonymous referee for many insightful comments, which
improved the paper significantly. We sincerely thank Dr H. Xu (Kavli Institute
for Astronomy and Astrophysics, Peking University) for assisting with data
collection. TG acknowledges the support of the National Science and Technology
Council of Taiwan through grants 108-2628-M-007-004-MY3 and
111-2123-M-001-008-. TH acknowledges the support of the National Science and
Technology Council of Taiwan through grants 110-2112-M-005-013-MY3,
110-2112-M-007-034-, and 111-2123-M-001-008-. This work used high-performance
computing facilities operated by the Center for Informatics and Computation in
Astronomy (CICA) at National Tsing Hua University. This equipment was funded
by the Ministry of Education of Taiwan, the National Science and Technology
Council of Taiwan, and National Tsing Hua University.

## Data Availability

The utilised data and the UMAP analysis result are available in the article’s
online supplementary material.

## References

  * Agarwal (2020) Agarwal D., 2020, PhD thesis, West Virginia University, doi:https://doi.org/10.33915/etd.7827
  * Badia et al. (2020) Badia A., Piot B., Kapturowski S., Sprechmann P., Vitvitskyi A., Guo D., Blundell C., 2020, Agent57: Outperforming the Atari Human Benchmark
  * Bhattacharyya (2017) Bhattacharyya S., 2017, arXiv e-prints, p. arXiv:1711.09083
  * Chen et al. (2021) Chen B. H., Hashimoto T., Goto T., Kim S. J., Santos D. J. D., On A. Y. L., Lu T.-Y., Hsiao T. Y.-Y., 2021, Monthly Notices of the Royal Astronomical Society, 509, 1227
  * Hashimoto et al. (2022) Hashimoto T., et al., 2022, MNRAS, 511, 1961
  * Hewitt et al. (2022) Hewitt D. M., et al., 2022, MNRAS, 515, 3577
  * Kim et al. (2022) Kim S. J., Hashimoto T., Chen B. H., Goto T., Ho S. C. C., Hsiao T. Y.-Y., Wong Y. H. V., Yamasaki S., 2022, MNRAS, 514, 5987
  * Li et al. (2021) Li D., et al., 2021, Nature, 598, 267
  * Lorimer et al. (2007) Lorimer D. R., Bailes M., McLaughlin M. A., Narkevic D. J., Crawford F., 2007, Science, 318, 777
  * Masui & CHIME/FRB Collaboration (2021) Masui K., CHIME/FRB Collaboration 2021, in American Astronomical Society Meeting Abstracts. p. 325.01
  * McInnes et al. (2018) McInnes L., Healy J., Melville J., 2018, arXiv preprint arXiv:1802.03426
  * Petroff et al. (2016) Petroff E., et al., 2016, Publ. Astron. Soc. Australia, 33, e045
  * Platts et al. (2019) Platts E., Weltman A., Walters A., Tendulkar S., Gordin J., Kandhai S., 2019, Physics Reports, 821, 1
  * Rousseeuw (1987) Rousseeuw P. J., 1987, Journal of Computational and Applied Mathematics, 20, 53
  * Vieyro et al. (2017) Vieyro F. L., Romero G. E., Bosch-Ramon V., Marcote B., del Valle M. V., 2017, Astronomy and Astrophysics, 602
  * Wagstaff et al. (2016) Wagstaff K. L., et al., 2016, Publications of the Astronomical Society of the Pacific, 128, 084503
  * Xu et al. (2022) Xu H., et al., 2022, Nature, 609, 685
  * Yamasaki et al. (2017) Yamasaki S., Totani T., Kiuchi K., 2017, Publications of the Astronomical Society of Japan, 70

##  Appendix A Additional clustering result

In the main text, we apply the K-means algorithm with n_cluster = 3 to cluster
the result provided by the unsupervised machine learning result. Nevertheless,
it is not the only way to present the cluster. In this section, we present
another result using K-means but set n_cluster = 4, and the other result with
HDBSCAN, a clustering algorithm that does not require a predetermined total
cluster count. The hyperparameter setting for HDBSCAN is min_cluster_size =
120 and min_samples = 25. The two results are presented in Fig. 8 and Fig. 9,
respectively.

Figure 8: K-means clustering result with n_cluster = 4. Figure 9: HDBSCAN
clustering result with min_cluster_size = 120 and min_samples = 25.

In Fig. 8, we find that the n_cluster = 4 case does not identify the possible
sub-cluster at the lower right. Instead, it divides the biggest apparently
single cluster (cluster 2 in Fig. 4) into two. Furthermore, n_cluster = 4 has
a worse silhouette score than n_cluster=3 in Fig. 3. As a result, we decided
to use n_cluster = 3 for further analysis in this research.

In Fig. 9, we see HDBSCAN is offering a result with significantly more
clusters. The cluster count in HDBSCAN is automatically determined, thus we
have set min_cluster_size small to identify the sub-clusters. As a result, the
possible sub-clusters are mostly identified, at the cost of splitting the
largest cluster into several pieces. This cluster result could work as an
alternative when an analysis with a more detailed classification of
FRB20121102A is needed.