11institutetext: Dipartimento di Fisica, Università Sapienza, P.le Aldo Moro
2, 00185, Roma, Italy 22institutetext: Wits Centre for Astrophysics, School of
Physics, University of the Witwatersrand, 1 Jan Smuts Avenue, Johannesburg
2050, South Africa

# Sunyaev Zel’dovich effect in galaxy clusters cavities: thermal or non-
thermal origin?

Paolo Marchegiani 1122<EMAIL_ADDRESS>
###### Abstract

Several galaxy clusters host X-ray cavities, often filled with relativistic
electrons emitting in the radio band. In the cluster MS 0735.6+7421 the
cavities have been detected through the Sunyaev Zel’dovich (SZ) effect, but it
has not been possible to determine if this effect is thermal (produced by a
very high temperature gas filling the cavity) or non-thermal (produced by the
relativistic electrons that produce the diffuse radio emission detected in the
cavity). In this paper we discuss the role of the density of the high
temperature gas inside the cavities in determining whether the dominant SZ
effect is the thermal or the non-thermal one, and how it can be possible to
distinguish between the two possibilities, discussing the role of observations
at higher energy bands.

##  1 Introduction

Several relaxed galaxy clusters show the presence inside their Intra Cluster
Medium (ICM) of couples of cavities in the X-ray emission, usually located in
opposite directions compared to the central Brightest Cluster Galaxy
Fabian2012 ; Gitti2012 . The origin of these cavities is probably due to the
mechanical power produced by AGN jets that is converted to heat during the
expansion of the lobe in the ICM McNamara2007 . Many X-ray cavities host a
diffuse radio emission Birzan2020 , showing that relativistic electrons and
magnetic fields are present in them. However, it has been found that the
energy stored in these relativistic electrons, determined under the minimum
energy assumption, is well below the energy necessary to inflate the cavity
Ito2008 , indicating that another component inside the cavity should be
dominant in energy. This dominant component can be a high temperature thermal
gas with a temperature of the order of $kT\sim 100$ keV or higher according to
hydrodynamic simulation Sternberg2009 ; Prokhorov2012 , or a population of
relativistic protons.

While radio observations are suitable to constrain the high energy part of the
relativistic electrons spectrum and the magnetic field, they are not able to
provide information on the low energy part of the electrons spectrum, and on
the thermal gas inside the cavity. The Sunyaev-Zel’dovich effect (SZE)
produced inside the cavity, i.e. the distortion of the Cosmic Microwave
Background (CMB) spectrum by inverse Compton scattering (ICS) Sunyaev1972 ,
has instead been proposed as a probe to constrain the low energy electrons or
the thermal gas properties Colafrancesco2005 ; Pfrommer2005 ; Prokhorov2010 .
With this kind of information it is possible to constrain the energetic
content of the different components inside the cavity, because the low energy
part of the electrons spectrum is the dominant one in determining the total
energy stored in them Colafrancesco2011 .

The galaxy cluster MS 0735.6+7421 hosts a couple of X-ray cavities
McNamara2005 ; Gitti2007 ; Vantyghem2014 , where a diffuse radio emission
filling the cavities is also observed Birzan2008 ; Birzan2020 . These cavities
have also been observed through the SZE at 30 GHz using the Combined Array for
Research in Millimeter-wave Astronomy (CARMA) interferometer Abdulla2019 ,
which detected a clear deficit in the SZ signal in the direction of the
cavities. From these observations, it has not been possible to determine if
the dominant component inside the cavities is thermal or non-thermal, because
the data were compatible with both a non-thermal electrons population having a
low value of the minimum moment ($p_{1}\sim 1-10$), and with a thermal
population with very high temperature of $kT\sim 1000$ keV, or even higher
(see fig.8 in Abdulla2019 ).

In a recent paper Marchegiani2021 we discussed this problem, pointing up that
these two different populations, i.e. the non-thermal electrons and the high
temperature gas, are linked each other because of the effect of the Coulomb
interactions. In particular, a higher density of the thermal gas implies a
higher rate of the Coulomb losses, which have the effect of reducing the
density of the low energy non-thermal electrons and therefore of flattening
their spectrum, and as a consequence of reducing the intensity of the non-
thermal SZE. Therefore, this effect should be taken into account when
interpreting the results of observations of the SZE inside the cluster
cavities in terms of thermal or non-thermal populations.

In this paper, we summarize the results obtained in Marchegiani2021 , and
expand further the discussion about how it can be possible to discriminate
between the thermal and the non-thermal origin of the SZE inside the cavities,
also considering the expected SZE in spectral bands with higher energies than
the microwave one.

##  2 Modeling the spectrum of non-thermal electrons

The time evolution of the non-thermal electrons spectrum $N_{e}(p,t)$, written
as a function of the time $t$ and the normalized momentum $p=\beta\gamma$, is
given the following equation Schlickeiser2002 :

$\frac{\partial N_{e}(p,t)}{\partial t}=\frac{\partial}{\partial
p}\left[\left(-\frac{2}{p}D_{pp}+\sum_{i}b_{i}(p)\right)N_{e}(p,t)+D_{pp}\frac{\partial
N_{e}(p,t)}{\partial p}\right],$ (1)

where $D_{pp}$ is the diffusion term in the momentum space, associated to
Fermi-II acceleration processes, which we model as $D_{pp}=\chi p^{2}/4$ so
that the characteristic acceleration time is
$\tau_{acc}=p^{2}/4D_{pp}=\chi^{-1}$ Brunetti2007 . The energy loss term
$b_{i}(p)$ is given by the sum of the Coulomb losses, $b_{Coul}\propto n_{th}$
Gould1972 , where $n_{th}$ is the density of the thermal gas, of the radiative
losses (synchrotron and inverse Compton scattering),
$b_{rad}\propto(B^{2}+B^{2}_{ICS})p^{2}$ Winner2019 , and of the adiabatic
expansion losses, $b_{ad}\propto(dV/dt)p$, where we assume that the volume is
changing with a Sedov-like expansion $V(t)\propto(t/t_{0})^{6/5}$ Ensslin2001
.

We assume an initial electrons spectrum with a power law shape,
$N_{e}(p,t_{0})=k_{0}p^{-s}$ with $s=2.7$, and initial values of the magnetic
field and the thermal density so that at the present time they are equal
respectively to the equipartion value $B=4.7$ $\mu$G Birzan2008 and to a
reference value chosen in the range $10^{-6}-10^{-3}$ cm-3. We therefore
calculate the electrons spectrum after a time of the order of the age of the
bubble, estimated to be $\sim 160$ Myr from X-ray observations Vantyghem2014
and MHD simulations Ehlert2019 , and choose the electrons normalization
$k_{0}$ and the Fermi-II acceleration parameter $\chi$ in order to match the
observed spectrum of the radio emission, with the best parameters values
resulting to be $k_{0}=1.5\times 10^{-3}$ cm-3 and $\chi=3.5\times 10^{-16}$
s-1 Marchegiani2021 .

In Fig.1, left panel, the resulting radio spectrum, compared with the observed
data Birzan2008 ; Birzan2020 , is shown for four values of the thermal density
inside the cavity, while the electrons spectra for the same values of the
thermal density are shown in the right panel. It is interesting to see that,
while the value of the thermal density does not impact too much on the shape
of the radio spectrum (in fact the relative differences between different
models are smaller than 10% in the whole frequency range between 100 MHz and
10 GHz), it heavily affects the spectrum of the electrons at low energies. As
shown in the next Section, this fact has strong consequences on the SZE
produced by these electrons.

Figure 1: Radio spectrum after 160 Myr from the injection (left panel) and the
corresponding non-thermal electrons spectra (right panel) for four values of
the thermal density: $10^{-6}$ (solid line), $10^{-5}$ (dashed line),
$10^{-4}$ (dot-dashed line), and $10^{-3}$ (three dots-dashed line) cm-3.
Figures are from Marchegiani2021 .

##  3 The non-thermal Sunyaev-Zel’dovich effect

Using the electrons spectra calculated in the previous Section, it is possible
to calculate the non-thermal SZE, using the full relativistic formalism
Wright1979 ; Colafrancesco2003 , where the intensity of the SZE is given by
$\Delta I(x)=\tau[J_{1}(x)-I_{0}(x)]$. Here $x=h\nu/(k_{B}T_{0})$ is the
frequency normalized at the CMB temperature $T_{0}$, $I_{0}(x)$ is the
unperturbed CMB spectrum, $\tau$ is the electrons optical depth, proportional
to the electrons density integrated along the line of sight, and $J_{1}(x)$ is
calculated by convolving the CMB spectrum with the single scattering
redistribution function, which is given by the convolution of the electrons
spectrum with the Compton scattering redistribution function for an electron
with a given energy Ensslin2000 .

The resulting SZE spectra for the four values of the thermal density that were
used in the previous Section are shown in the first panel of Fig.2. As it is
possible to see, smaller values of the thermal density imply reduced Coulomb
losses, and therefore a higher density of the non-thermal electrons at low
energies and as a consequence a more intense SZE.

Figure 2: First panel: non-thermal SZE produced by the electrons with the
spectra shown in Fig.1, with the same meaning of the line styles. Other
panels: thermal SZE inside the cavities (solid lines) compared with the non-
thermal SZE calculated for the corresponding values of the thermal density
(dashed lines) for values of the temperature of 1000, 1500, and 2000 keV. The
dotted line in last three panels indicates the zero level. Figures are from
Marchegiani2021 .

The SZE produced by the non-thermal electrons inside the cavity can be
compared with the thermal SZE produced by the high temperature gas inside the
same cavity. Since we don’t have observational information about the
properties of this gas, we can assume several values for its temperature, and
assume for its density the values for which its pressure is equal to the one
of the external ICM at the location of the cavity, $P_{cav}=6\times 10^{-11}$
erg cm-3 Gitti2007 , i.e. $n_{th}=P_{cav}/(kT)$. Using the full relativistic
formalism, we calculate the thermal SZE for these temperatures and the optical
depths derived from these values of the density, and the non-thermal SZE for
the Coulomb losses rates provided by the same values of the thermal density.
In the last three panels of Fig.2 the results for three values of the
temperature are shown.

As it is possible to see, values of the temperature until $\sim 1500$ keV
allow to have quite high values of the thermal density, and this fact
increases the intensity of the thermal SZE and at the same time decreases the
one of the non-thermal SZE because of the Coulomb losses. For higher
temperatures, we can expect that the non-thermal SZE should be dominant, while
for lower temperatures the thermal effect can be dominant, even if it is
possible that the thermal density is lower than the value assumed here, so it
is not possible to state whether in this case the dominant effect is
definitely the thermal one.

##  4 Thermal or non-thermal SZE?

In this Section we discuss how it can be possible to determine if the dominant
SZE inside the cavity is the thermal or the non-thermal one, adding some other
considerations as compared to the previous paper Marchegiani2021 , to which we
refer instead for other details on the models, on the properties of the
electrons populations used in the calculations, on the comparison between the
SZE inside the cavity and the one of the surrounding ICM, and on the
observational challenges to be faced in order to derive this kind of
information.

As discussed in the previous paper, the angular size of the cavities in MS
0735.6+7421 is of the order of 1 arcmin; therefore instruments with angular
resolution of the order of 10 arcsec like NIKA2 are well suited to resolve the
SZ signal inside the cavities, and to distinguish it from the one coming from
the rest of cluster. While NIKA2 can provide information at the frequencies of
150 and 260 GHz, a wider spectral coverage could be useful to better
distinguish between a thermal and a non-thermal origin of the SZE; this can be
done with instruments operating at 90 GHz, like Mustang-2 at the Green Bank
Telescope (GBT) and Mistral at the Sardinia Radio Telescope (SRT), and at
higher frequencies, like the planned space telescope Millimetron Kardashev2014
, with its planned angular resolution of a few arcsec and a spectral coverage
until a few THz.

However, also a full spectral frequency coverage of the SZE inside the cavity
might not be sufficient to determine if its origin is thermal or non-thermal.
In fact, in Fig.2 we can see that for temperatures in the range 1000 – 1500
keV the thermal and the non-thermal SZE are expected to have very similar
spectral shapes (note for example that the crossover frequency in the top
right panel of Fig.2 is basically the same for the two effects).

A possible solution might be to look at the SZE emission at frequencies higher
than the THz, where it should reach its maximum, depending on the properties
of the electrons population. For example, it has been suggested that the Hard
X-ray band can be good to distinguish between the ICS emission of the non-
thermal electrons and the bremsstrahlung of the high temperature thermal gas,
because the efficiency of the ICS process is much higher than the
bremsstrahlung one Prokhorov2012 . In order to check this possibility, in
Fig.3 we compare the thermal and non-thermal SZE spectra calculated at high
frequencies, from the THz until $10^{19}$ Hz (corresponding to a photon energy
of $\sim 40$ keV), for values of the gas temperature inside the cavity of 1000
and 1500 keV.

Figure 3: SZE spectra at high frequencies for the thermal gas inside the
cavity (solid lines) with temperatures of 1000 (left panel) and 1500 keV
(right panel), compared with the SZE spectra of the non-thermal electrons
(dashed lines), and with the bremsstrahlung of the high temperature gas inside
the cavity (dotted line) and of the external ICM along the line of sight of
the cavity center (long dashed lines).

As it is possible to see, the spectra of the thermal ad the non-thermal SZE
are quite similar until frequencies of $\nu\sim 10^{14}$ Hz, while they start
to diverge at higher frequencies, with the thermal SZE falling down quite
abruptly, and the non-thermal SZE continuing instead with a power-law shape
until the Hard X-ray band, with a progressive steepening due to the
corresponding high energy steepening of the electrons spectrum, visible in the
right panel of Fig.1. Therefore, a hypothetical detection of the SZE inside
the cavity at frequencies $>10^{14}$ Hz, i.e. from the near infrared band to
higher frequencies, should be attributed to the non-thermal electrons rather
than to the high temperature gas. Such an effect should appear as a positive
diffuse emission spatially coincident with the cavity.

In Fig.3 we also show an estimate of the surface brightness in the direction
of the cavity center due to the thermal bremsstrahlung of the high temperature
gas inside the cavity and of the external ICM present along the line of sight.
This last quantity has been calculated assuming the presence of a gas with
constant temperature of 5.5 keV and with the double beta radial profile
derived from X-ray observations (see Marchegiani2021 for details), extended
until the virial radius of 2.23 Mpc Gitti2007 , along the line of sight at the
distance of 150 kpc from the cluster center, and excluding from the
integration along the line of sight the cavity region, approximated as a
sphere with radius of 100 kpc. The bremsstrahlung produced by the high
temperature gas inside the cavity results to be much lower than the one from
the ICM, as expected because of its lower density, higher temperature, and
smaller size. The bremsstrahlung produced by the ICM is instead dominant on
the SZE produced by the non-thermal electrons at frequencies higher than $\sim
2\times 10^{15}$ Hz, i.e. in the far ultraviolet band.

Therefore there is a spectral region, going from the near infrared to the far
ultraviolet, where the non-thermal SZE should be the dominant diffuse
component inside the cavity. Should a diffuse emission spatially coincident
with the cavity in these spectral bands be detected, it should therefore be
attributed to the non-thermal SZE. In particular, the most favorable spectral
region results to be at a wavelength around 400 nm, i.e. between the near
ultraviolet and the visible spectral regions. Instruments that can be suitable
to detect such an emission might be UVOT aboard the Swift satellite, or the
WFC3 camera aboard the Hubble Space Telescope; in order to establish the
effective possibility to detect this signal, a dedicated study would be
required.

## Acknowledgements

The Author thanks the organizers and the participants of the “mm Universe
@NIKA2” Conference for very interesting and stimulating talks and discussions.
The Author also thanks the Referee for useful comments that helped to improve
the quality of the paper, and A. Maselli for his help and suggestions.

## References

  * (1) Fabian A.C., 2012, ARA&A, 50, 455
  * (2) Gitti M., Brighenti F., McNamara B.R. 2012, AdAst, 2012, 950641
  * (3) McNamara B. R., Nulsen, P. E. J., 2007, ARA&A, 45, 117
  * (4) Birzan L. et al., 2020, MNRAS, 496, 2613
  * (5) Ito H., Kino M., Kawakatu N., et al., 2008, ApJ, 685, 828
  * (6) Sternberg A., Soker N., 2009, MNRAS, 398, 422
  * (7) Prokhorov D.A., Moraghan A., Antonuccio-Delogu V., Silk J., 2012, MNRAS, 425, 1753
  * (8) Sunyaev R.A., Zel’dovich Ya.B., 1972, CoASP, 4, 173
  * (9) Colafrancesco S., 2005, A&A, 435, L9
  * (10) Pfrommer C., Enßlin, T.A., Sarazin C.L., 2005, A&A, 430, 799
  * (11) Prokhorov D.A., Antonuccio-Delogu V., Silk J., 2010, A&A, 520, A106
  * (12) Colafrancesco S., Marchegiani P., 2011, A&A, 535, A108
  * (13) McNamara B.R., Nulsen P.E.J., Wise M.W., et al., 2005, Nature, 433, 45
  * (14) Gitti M., McNamara B.R., Nulsen P.E.J., Wise M.W., 2007, ApJ, 660, 1118
  * (15) Vantyghem A.N., McNamara B.R., Russell H.R., et al., 2014, MNRAS, 442, 3192
  * (16) Birzan L., McNamara B.R., Nulsen P.E.J., et al., 2008, ApJ, 686, 859
  * (17) Abdulla Z. et al., 2019, ApJ, 871, 195
  * (18) Marchegiani P., 2021, MNRAS, 503, 4183
  * (19) Schlickeiser R., 2002, Cosmic Ray Astrophysics, Springer-Verlag, Berlin
  * (20) Brunetti G., Lazarian A., 2007, MNRAS, 378, 245
  * (21) Gould R.J., 1972, Physica, 60, 145
  * (22) Winner G., Pfrommer C., Girichidis P., Pakmor R., 2019, MNRAS, 488, 2235
  * (23) Enßlin T.A., Gopal-Krishna, 2001, A&A, 366, 26
  * (24) Ehlert K., Pfrommer C., Weinberger R., Pakmor R., Springel V., 2019, ApJ, 872, L8
  * (25) Wright E.L., 1979, ApJ, 232, 348
  * (26) Colafrancesco S., Marchegiani P., Palladino E., 2003, A&A, 397, 27
  * (27) Enßlin T.A., Kaiser C.R., 2000, A&A, 360, 417
  * (28) Kardashev N.S. et al., 2014, Physics-Uspekhi, 57, 12


[acronym]long-short

#  toki: A Build- and Test-Platform for Prototyping and Evaluating Operating
System Concepts in Real-Time Environments ††thanks:  This work was supported
by the European Union (EU) under the Horizon 2020 program, projects TAPPS
(Trusted Apps for open CPSs) and 5GCroCo (5G Cross Border Control), and the
German Federal Ministry of Economics and Technology (BMWi) under the Smart
Service World program, project PASS (Platform for Automotive Apps Guaranteeing
Security and Safety).

Oliver Horst  fortiss GmbH – Research Institute of the Free State of Bavaria
Guerickestr. 25, 80805 Munich, Germany
Email<EMAIL_ADDRESS>Uwe Baumgarten  Technical University of Munich
Department of Informatics, Germany
Email<EMAIL_ADDRESS>
###### Abstract

Typically, even low-level os concepts, such as resource sharing strategies and
predictability measures, are evaluated with Linux on PC hardware. This leaves
a large gap to real industrial applications. Hence, the direct transfer of the
results might be difficult. As a solution, we present _toki_ , a prototyping
and evaluation platform based on FreeRTOS and several open-source libraries.
toki comes with a unified build- and test-environment based on Yocto and Qemu,
which makes it well suited for rapid prototyping. With its architecture chosen
similar to production industrial systems, toki provides the ground work to
implement early prototypes of real-time systems research results, up to
technology readiness level 7, with little effort.

##  I Introduction

Currently, most applied real-time systems research prototypes are developed
and evaluated on top of Linux on PC hardware. This leaves a large gap between
real industrial applications in that field and the prototype. In case of low-
level os concepts concerning, context switch times, resource sharing, intra-
node communication, and predictability, the drawn conclusions could even be
void due to the completely different nature of the industrial platform.
Furthermore, we see a lack in practical examinations of which latency and how
much temporal predictability, in the sense of [1], is achievable with certain
configurations (software architectures and predictability measures). Hence, we
see the need to ease constructing early prototypical implementations of
research results on relevant hardware in relevant environments.

On the one hand, Linux seems to be a good choice here. It is available for a
wide variety of hardware platforms, provides excellent third-party library
support, and can fulfill real-time requirements to some extent. However, it is
a complex task to configure Linux in a way that its influence on low-level
performance benchmarks is negligible or at least predictable. Moreover,
integrating own concepts into the Linux kernel requires detailed knowledge
about the kernel sources and its concepts. Hence, even though specialized
distributions, such as LITMUSRT [2], can drastically reduce the configuration
effort, the integration complexity remains as main issue.

Industrial-grade rtos, on the other hand, provide considerably less influence
on performance measurements, but at the costs of usability. They are either
delivered as minimal systems, like FreeRTOS [3], which lack tooling support
and provide not much more than a scheduler, or as sophisticated platforms
tailored to specific industrial fields, such as AUTOSAR [4], which come along
with royalty fees, tooling, and complex, configurable software stacks.

As an exception, the GENODE OS framework [5] comes with toolchains for several
hardware platforms and provides a convenient, production ready rtos.
Unfortunately, GENODE’s overall orientation to safety and security severely
hinders rapid prototyping. Its micro-kernel based design requires that all
adaption and changes to the os obey the strict isolation mantra, which is
conceptually challenging and time consuming.

Therefore, we see the need for a minimal, yet flexible real-time system
framework that provides comfort similar to commercial platforms, but comes
without a complex software architecture and strict inherent design concepts.
Ideally, the framework should be fully based on open-source software, provide
an integrated standard C-library, a target toolchain, configuration tools, and
an emulation environment for testing. A plus would be the possibility to
safety certify the system.

Accordingly, we present _toki_ , a flexible and configurable os framework,
close to production industrial systems, but without the hassle of complex
software architectures.

##  II The Build- and Test-Platform Toki

toki, Japanese for “when an action occurs”, is a build- and test-environment
constructed around FreeRTOS or SafeRTOS [6], respectively. toki’s main goal is
to compose an easy to use rapid prototyping platform to develop and evaluate
new operating-system- and intra-node-communication-concepts for cyber-physical
systems on commodity microcontrollers. Therefore, toki combines the following
open-source projects into the unified architecture and build-system
illustrated in Fig. 1: FreeRTOS [3] (v10.0.0), newlib [7] (v3.0.0),
libXil/libXilPm [8] (v2019.1), and lwIP [9] (v2.0.3). These base components
are supplemented by the ESFree Scheduling Library [10], HST User Mode
Scheduler [11], and TACLe benchmarks [12].

These components were chosen under the following aspects: _(i)_ the simplicity
to modify and extend their code base, _(ii)_ their ability to closely mimic
the software stacks of comparable industrial systems, _(iii)_ the flexibility
of their design regarding rapid prototyping, and _(iv)_ compatible software
licenses.

Figure 1:  Sketch of toki’s architecture and build-system, showing the
included components (3,4), differentiated among integrated (3) and newly
created (4) components, and their relations. Given a deployment configuration
(1), Yocto (2) builds all included components and collates them into a SD card
image (5), which can be deployed and tested on the hardware or in Qemu (6).

At platform-level, the toki build-system (Fig. 1) utilizes Yocto [13] to
provide a self-contained build-environment that builds all components of toki
together with their dependencies (a GCC cross-compiler). At project-level, we
ensured proper include prefixes for all components, by restructuring their
sources and adding a CMake [14] based build-system, where needed. This
reorganization is conducted by a dedicated Python script, to allow nearly
automatic upstream pulls.

To allow software developers to benefit from the code insight and debugging
capabilities of their favorite IDE, single-core deployments of toki can be
built by CMake. Production-ready images of multi-core deployments, on the
other hand, can only be built by Yocto. toki deployment configurations, in
general, are specified either manually or with assistance of the toki-config
tool, and guide Yocto in the build process.

Besides the changes required for the build-environment, we contribute the
following new features to FreeRTOS on ARMv8: amp boot support, a newlib
syscall interface, a memguard [15] implementation, and software-tracing
support via ARM’s system trace macrocell [16].

The toki build-environment is seamlessly integrated with a virtual test-
environment. Hence, all images built by Yocto can either be tested on the real
hardware or its emulated counterpart. The emulation is handled by a tailored
Qemu [17], built by Yocto, which enables a fully virtual development cycle.

Currently, toki is solely tested and ready-to-run on the Xilinx Zynq
UltraScale+ MPSoC platform [18], a contemporary ARMv8 multi-core SoC with an
integrated FPGA. This platform was selected, because of _(i)_ Xilinx’s
extensive software support for it, including a bare-metal driver kit, and
_(ii)_ the included FPGA, which allows us to increase the accuracy of our
performance evaluations. Nevertheless, toki was designed with portability in
mind; hence, it can easily be ported to other target platforms, the STM32Cube
by STmicro [19].

In the future, we plan to conduct measurements to compare the interrupt
handling latency of toki with other software stacks and extend toki with
FreeRTOS+POSIX [20], a communication middleware, and ptp support.

Feel free to try out the latest version of toki, by downloading it from:
https://git.fortiss.org/toki

##  III Demonstration

We will demonstrate toki through two showcases:

  1. 1.

A video of the TAPPS project’s [21] final demonstrator, showing the on-the-fly
installation of an application into the real-time critical control path of the
throttle control of a production electric motorcycle, realized with toki.

  2. 2.

A live-demo of toki’s configuration, build, simulation, and deployment cycle
on the example of memory benchmarks deployed to distinct cores and regulated
by memguard.

The first showcase focuses on the real-time capabilities and applicability to
industrial use-cases, and the second on the configuration flexibility and
measurement capabilities of toki.

## Acknowledgment

We would like to thank all contributors to toki for their work, namely: Martin
Jobst, Johannes Wiesböck, Dorel Coman, Ulrich Huber, Mahmoud Rushdi, Tuan Tu
Tran, Firas Trimech, Raphael Wild, Andreas Ruhland, and Dhiraj Gulati.

## References

  * [1] B. Sun _et al._ , “Definitions of predictability for cyber physical systems,” _J. of Syst. Architecture - Embedded Syst. Des._ , vol. 63, pp. 48–60, 2016.
  * [2] J. M. Calandrino _et al._ , “LITMUS^RT : A Testbed for Empirically Comparing Real-Time Multiprocessor Schedulers,” in _27th IEEE Int. Real-Time Syst. Symp. (RTSS)_ , Dec. 2006, pp. 111–126.
  * [3] https://www.freertos.org.
  * [4] AUTOSAR, “Classic platform release 4.3.1,” AUTOSAR Std., Dec. 2017.
  * [5] N. Feske, _GENODE Foundations – Operating System Framework 19.05_. GENODE Labs, 2019.
  * [6] https://www.highintegritysystems.com/safertos/.
  * [7] https://sourceware.org/newlib/.
  * [8] https://github.com/Xilinx/embeddedsw.
  * [9] https://savannah.nongnu.org/projects/lwip/.
  * [10] R. Kase, “Efficient Scheduling Library for FreeRTOS,” Master’s thesis, KTH Information and Communication Technology, 2016.
  * [11] F. E. Páez _et al._ , “FreeRTOS user mode scheduler for mixed critical systems,” in _6th Argentine Conf. on Embedded Syst. (CASE)_ , Aug. 2015.
  * [12] H. Falk _et al._ , “TACLeBench: A Benchmark Collection to Support Worst-Case Execution Time Research,” in _16th Int. Workshop on Worst-Case Execution Time Analysis (WCET)_ , 2016, pp. 2:1–2:10.
  * [13] https://www.yoctoproject.org.
  * [14] https://cmake.org.
  * [15] H. Yun _et al._ , “MemGuard: Memory bandwidth reservation system for efficient performance isolation in multi-core platforms,” in _19th IEEE Real-Time and Embedded Technol. and Appl. Symp. (RTAS)_ , Apr. 2013.
  * [16] _CoreSight System Trace Macrocell – Technical Reference Manual_ , r0p1 ed., ARM Limited, Dec. 2010.
  * [17] https://www.qemu.org.
  * [18] https://www.xilinx.com/products/silicon-devices/soc/zynq-ultrascale-mpsoc.html.
  * [19] https://www.st.com/content/st_com/en/stm32cube-ecosystem.html.
  * [20] https://www.freertos.org/FreeRTOS-Plus/FreeRTOS_Plus_POSIX/.
  * [21] http://www.tapps-project.eu.


[1]Nathan Haut

[1]CMSE, Michigan State University, 426 Auditorium Rd., East Lansing, 48824,
Michigan, USA

2]Computer Science, Michigan State University, 426 Auditorium Rd., East
Lansing, 48824, Michigan, USA

# Active Learning in Genetic Programming: Guiding Efficient Data Collection
for Symbolic Regression

https://orcid.org/0000-0002-3989-6532<EMAIL_ADDRESS>Wolfgang Banzhaf
https://orcid.org/0000-0002-6382-3245<EMAIL_ADDRESS>Bill Punch
https://orcid.org/0000-0002-4493-6877<EMAIL_ADDRESS>*  [

###### Abstract

This paper examines various methods of computing uncertainty and diversity for
active learning in genetic programming. We found that the model population in
genetic programming can be exploited to select informative training data
points by using a model ensemble combined with an uncertainty metric. We
explored several uncertainty metrics and found that differential entropy
performed the best. We also compared two data diversity metrics and found that
correlation as a diversity metric performs better than minimum Euclidean
distance, although there are some drawbacks that prevent correlation from
being used on all problems. Finally, we combined uncertainty and diversity
using a Pareto optimization approach to allow both to be considered in a
balanced way to guide the selection of informative and unique data points for
training.

###### keywords:

Active learning, Genetic programming, Symbolic regression

##  1 Introduction

In applications of data science, the task of collecting and labelling data is
often time-consuming and expensive. In some cases where data doesn’t yet
exist, it may be very expensive to run experiments to gather data, or possibly
it could take long periods of time for experiments to complete. In these
cases, it would be ideal to target specific experiments where maximal
information will be gained, so fewer experiments have to be run to gain the
desired insight into the system of study. In other cases, large masses of data
may already exists, but the process of labelling the data is time-consuming.
Here, it would be ideal to target a subset of samples, that when labelled,
will provide the most information. To achieve these time-savings and cost
reductions we can use machine learning (ML) not only to build models to
describe these systems, but also to predict the information gained by each
training sample. The process of using machine learning to iteratively select
data to best inform machine learning model development is called active
learning.

More specifically, active learning (AL) is a method used in conjunction with
machine learning to actively select new training data with the goal of
selecting data points that will maximally inform the machine learning model
[1]. Various forms of active learning exist, with three types dominating:
pool-based AL, stream-based AL, and membership query synthesis [19]. Figure 1
shows a simple visual representation to compare the three methods of active
learning. Pool-based and stream-based methods both have a set of training
samples to choose from, with the goal of selecting and training on only a
small subset of maximally informative cases. The key difference between pool-
based and stream-based methods is that pool-based methods search over a set of
data points for the ones that are most informative. Steam-based methods differ
by checking each potential training case in order one-by-one and only admit
them to the training set if data points are ”informative”. Membership query
synthesis approaches do not have a set of already existing training samples to
choose from, instead, they search a training space to find and synthesize new
training data points that are expected to maximally inform the machine
learning model. Once synthesized, a new data point is then labelled by the
researcher via experimentation or expert knowledge.

Figure 1: The three main types of active learning: Stream-based, pool-based,
and membership query synthesis are visually demonstrated. Stream-based
approaches, shown on the left, search through the samples one at a time and
either mark them for labelling or skip them. Green indicates a sample is found
to be informative and is marked for labelling, red indicates a sample is
skipped. Pool-based approaches, shown in the middle, assigns an information
score to each potential training sample and the most informative sample is
chosen to be labelled and added to the training set. Membership query
synthesis, shown on the right, searches a space of potential points not yet
collected while maximizing an information measure and selects a point to be
synthesized and labelled that maximizes the information score. The selected
point is indicated by the green circle, while the y-axis of the curve
represents the informativeness measure and the x-axis is representative of the
sample space.

Active learning is a versatile method with uses ranging from effectively sub-
sampling of data from a huge set for training, sampling of data with specific
goals such as to maximize diversity, to guiding experimentation by suggesting
experiments that will be most informative to the researcher in the model
building process. It can be used to focus on interesting samples from large
sets or to expand small data sets while minimizing data collection efforts.
For example, AL has recently been used to explore a space of 16 million
potential catalysts to maximize the conversion rate of methane to methanol,
which without active learning would not have been possible to search
effectively within a reasonable time [15]. Active learning has also been shown
to effectively sub-sample training data for identifying malware-infected PDF
documents [14]. The authors found that when using active learning they could
reduce the training set size to 1/30-th of the original size, while
maintaining the same performance as models trained on the whole set.

For a wide range of machine learning methods active learning approaches have
been developed, e.g. for support vector machines or neural networks. In
support vector machines, for instance, AL has been realized by computing the
distance of all points to the separating hyperplane and selecting the point
nearest the hyperplane to be labelled [11]. For neural networks, one AL
variant has been to select points with minimum difference between the two most
probable predicted labels [16]. This distribution was defined as
$M=P(l1|x)-P(l2|x)$, where $M$ is the margin between the two most probable
labels, $l1$ is the most probably label for input $x$, and $l2$ is the second
most probable label for input $x$.

In this contribution, we apply active learning strategies for genetic
programming used in symbolic regression tasks. The goal is to exploit some of
the features of GP, in particular its reliance on a population of models. More
specifically, we want to utilize uncertainty and diversity measures in a model
population context to accelerate the discovery of models (physics equations in
our study). The idea is to look for disagreement among high-quality
individuals in the population as a guide to locate informative data points to
add to the training set.

##  2 Related Works

Active learning methods for machine learning have shown to be very successful
in applied settings to improve the method of labelling and collecting data
with various machine learning types. AL has recently been demonstrated to
significantly reduce the labelling efforts required for labelling data
associated with identifying heart disease [5]. The authors demonstrated that
they could find more accurate models using fewer data points when compared to
a random point selection strategy.

AL has been applied to genetic programming classification tasks as well. Using
an ensemble of GP models, the models ”vote” on the class of data pairs, and
points are only labelled when the committee of developing models encounters
pairs that can’t be classified [3]. This was found to reduce the total effort
needed to label training points, since only a subset had to be labelled before
finding accurate models. Where GP training sets are large, AL has been
successfully applied by selecting sub-samples to be used for training [13, 2].
In [2] AL is performed by segmenting the data into smaller blocks and training
the models using one randomly selected block at a time using uniform
probability. As training continues, bias is introduced into the probability by
increasing the tendency to select blocks that haven’t been seen in a while, as
well as blocks where the models performed poorly during training. AL for sub-
sampling with genetic programming was found to decrease training times to find
better binary classification models by an order of magnitude [2]. In [13]
subsets were selected by dynamically developing a fitness case topology that
could be used to create minimally related subsets of data. In this context,
the strength of a relationship between two training cases was indicated by the
number of individuals that were able to solve both training cases.

In the discovery of biological networks AL methods have also been employed
successfully [21]. Several different approaches were explored by the authors
for determining which new data points would be maximally informative for a
wide range of machine learning models, including Boolean networks, causal
Bayesian networks, differential equation models, etc. One approach the authors
explored was the maximum difference method in which two best-fit models are
chosen and a new data point is selected where those two best-fit models have
the largest difference in predictions. They also examined entropy score
maximization. In that method a new data point is selected that maximizes an
entropy score, where entropy can be thought of as the amount of information to
be gained by gathering that data point. The entropy score $H_{e}$ is computed
as follows:

$H_{e}=-\sum_{x=1}^{x_{e}}\frac{e_{x}}{|M|}\log_{2}\frac{e_{x}}{|M|}$

where $M$ is the set of Boolean networks, $x_{e}$ is the number of network
states for a given data point, and $e$ is the set of all potential data
points.

In chemical engineering AL has been applied to expedite a reaction screening
process by only selecting a subset of maximally informative experiments to
complete rather than by exhaustively performing all possible experiments [6].
This was done by training neural networks and using them to select a subset of
experiments that maximized the information gain. Maximal information gain was
determined by looking at the standard deviation of an ensemble of neural
networks.

Kotanchek et al. (2009) used genetic programming for active design of
experiments, where models developed by a GP system are used to find optimal
conditions in a system of study. Active design of experiments is an
application of active learning, where it has the goal of designing experiments
that have specific properties or yield maximal information. The authors
proposed to employ ensembles of models from symbolic regression to find
regions of uncertainty in order to gather new data with high information
content. While this method has been proposed for how an active learning method
using model ensembles could be applied to GP for symbolic regression, there
has yet to be any research showing how active learning methods affect the
performance of GP symbolic regression tasks or how the method to quantify
uncertainty affects the quality of points selected for inclusion in the
training data. As well, it is yet to be shown that this idea of selecting an
ensemble from a model population and searching for points of high uncertainty
or disagreement among models is generalizable to any machine learning method
where a population of models is available.

##  3 Methods

We compare two classes of active learning: uncertainty and diversity based.
The implementations are described in detail below. We use two random sampling
methods as a baseline to compare the performance of the active learning
methods. The key features of the GP system we used, StackGP, are also
discussed.

Figure 2: An overview of the iterative active learning approach. It begins
with an initially randomly selected dataset. It then iteratively evolves
models and selects new training points that maximize uncertainty of an
ensemble of models. By maximizing ensemble uncertainty to select new training
samples, points with relatively high information content are added to the
training set each iteration.

###  3.1 Active Learning

Two general types of active learning were implemented to work with StackGP for
the purpose of accelerating the development of models to fit physics data from
the Feynman Symbolic Regression Dataset (Tegmark). The first type of active
learning explored was uncertainty-based, a model-driven approach to active
learning, where an ensemble of diverse, high-quality models from a population
was used to search for regions in the search space where there was high
uncertainty or disagreement between the models. The second type of active
learning explored was diversity-based active learning, where new points are
selected that differ maximally from the points already in the training sample.
This second type of active learning is a data-driven approach rather than a
model-driven approach. The first type of active learning is summarized in
Figure 2.

Both types of active learning methods were implemented to determine how they
each impact the success of evolution in genetic programming symbolic
regression tasks. Several different uncertainty and diversity metrics are
implemented to determine their respective impact on the success of the task.
Success of active learning by maximizing uncertainty would indicate that the
diversity of the population can be utilized to guide the collection of
informative data. Success of diversity sampling would indicate that GP
symbolic regression model development benefits from improved data sampling.

####  3.1.1 Maximizing Uncertainty

Several different uncertainty metrics were explored to determine how different
measures impact the success of active learning, and which approach would
generally work best. As an overview, each approach begins by selecting an
ensemble of models using the same method, then a function that uses the
specific uncertainty metric along with the ensemble and current training set
is created. This function is then fed to an optimizer to search for regions of
relatively high uncertainty. The most uncertain point found is then returned
and selected to be added to the training set. In total, there were 6 different
uncertainty maximization approaches tested which varied in how they quantified
disagreement, whether outlier predictions were considered, and which optimizer
was used. The steps and methods will be described in greater detail below and
the entire process is depicted in Algorithm 2.

Generating the ensemble is the first step in uncertainty-based active
learning. The goals for generating the ensemble were to capture diverse, high-
quality individuals from the population while keeping the size of the ensemble
relatively small so that the computational cost of optimizing uncertainty is
reasonable. The diversity goal is essential to the success of active learning
since disagreement between models is a necessary requirement. The method
chosen to capture both diversity and quality from the model population works
by clustering the training data using the input space and selecting a model
that best fits each cluster, ensuring no model is selected more than once. If
a model is already selected by another cluster, the next best unselected model
is chosen. The minimum number of clusters is set to 3 and the maximum is set
to 10. Thus, 3-10 models are chosen for inclusion in an ensemble. Data
clustering was chosen with the intent to capture diversity by focusing on
models that have biases for different regions of the training space. Quality
in the population would be captured since only models with the best fitness
were selected for each cluster. The algorithm to generate the ensemble is
described in detail in Algorithm 1.

Algorithm 1 Ensemble generation process to select diverse high-quality models.

procedure EnsembleSelect($models$,$trainingData$,$responseData$)

$selectedModels\leftarrow[]$ $\triangleright$ Initialize ensemble

$nClusters\leftarrow min(len(trainingData),10)$ $\triangleright$ Determine
number of clusters

$clusters\leftarrow KMeans(nClusters).fit\\_predict(trainingData)$

for $i=0$;$i++$; $i<nClusters$ do $\triangleright$ Loop over data clusters

$modelErrors\leftarrow computeError(models,clusters[i])$

$sortedModels\leftarrow sortBy(models,modelErrors)$

$j=0$

while $sortedModels[j]$ in $selectedModels$ do $\triangleright$ Find best
unselected model

$j++$

end while

$selectedModels=join(selectedModels,sortedModels[j]$ $\triangleright$ Add to
ensemble

end for

return $selectedModels$ $\triangleright$ Return ensemble

end procedure

The second step of this method is to utilize the specified uncertainty
function with both the current training data and the selected ensemble. The
function is then given to the optimizer with the search space boundaries to
find a point of relatively high uncertainty. In the case that an already
selected point is re-selected, a new search is initiated within a random sub-
region until a unique point is added. This ensures that new information is
added in each iteration to the training set.

The two methods used for optimization were Scipy Optimize’s minimize and
differential evolution (SciPy, SciPy).

In total 5 different uncertainty metrics were used, shown by Equations 1 to 5,
where Equation 5 is used twice, once with Scipy’s minimize function for
optimization, and a second time with Scipy’s differential evolution function
for optimization.

$\Delta=\frac{\text{Std}(\text{EnsembleResponses})}{\text{Mean(Abs(}\text{EnsembleResponses}))}$
(1)

$\Delta=\frac{\text{TrimmedStd}(\text{EnsembleResponses},0.3)}{\text{TrimmedMean(Abs(}\text{EnsembleResponses}),0.3)}$
(2)

$\Delta=\frac{\text{Std}(\text{EnsembleResponses})}{\text{TrimmedMean(Abs(}\text{EnsembleResponses}),0.3)}$
(3)

$\Delta=\text{Std}(\text{EnsembleResponses})$ (4)

$\Delta=\text{DifferentialEntropy}(\text{EnsembleResponses})$ (5)

Algorithm 2 Active Learning Process Using Uncertainty

$TrainingData\leftarrow 3StartingPoints$ $\triangleright$ Generate initial
random training data

$Models\leftarrow RandomModels$ $\triangleright$ Generate initial random
models

$Models\leftarrow Evolve(TrainingData,Models)$ $\triangleright$ Train models
on starting data

while $BestModelError\neq 0$ do $\triangleright$ While perfect model not found

$Ensemble\leftarrow EnsembleSelect(Models)$. $\triangleright$ Select ensemble
of models

$NewPoint\leftarrow MaxUncertainty(Ensemble)$ $\triangleright$ Find point of
max uncertainty

if $NewPoint\subset TrainingData$ then $\triangleright$ If point already
selected

$NewPoint\leftarrow MaxUncertainty(SubSpace(Ensemble))$$\triangleright$ Search
a subspace

end if

$TrainingData\leftarrow Append(TrainingData,NewPoint)$ $\triangleright$ Add
new point

$Models\leftarrow Evolve(TrainingData,Models)$ $\triangleright$ Evolve new
models with new data using best models to seed evolution

end while

####  3.1.2 Point Diversity

A data-driven active learning approach was also explored, aiming to maximize
data diversity rather than maximize ensemble uncertainty. The goal was to
determine if GP evolution for symbolic regression tasks would benefit
significantly from improved sampling of the data for training. Two different
metrics were used to quantify diversity: point distance and point correlation.
Point distance was implemented by measuring both the minimum and average
Euclidean distance to all points in the training set. Point correlation was
defined as the average correlation to all points in the training set. When
selecting a new point, the goal was to either maximize the distance or
minimize the correlation to the current training set.

To minimize the correlation when selecting a new point, Pearson’s $R^{2}$ was
computed between each point and the potential new point. The equation for
computing Pearson’s $R$ is shown in Equation 6. Here $y$ represents the new
training point, $\hat{y}$ represents a point already in the set, and each
instance $i$ represents the value in the ith dimension of the point. The
overall method for computing the joint correlation of a new point to the
training set is summarized in Algorithm 3.

$R=\frac{\sum_{i=1}^{N}(y_{i}-\bar{y})(\hat{y}_{i}-\bar{\hat{y}})}{\sqrt{\sum_{i=1}^{N}(y_{i}-\bar{y})^{2}\times\sum_{i=1}^{N}(\hat{y}_{i}-\bar{\hat{y}})^{2}}}$
(6)

Algorithm 3

1:procedure JointCorrelation($trainingSet$,$newPoint$)

2: $r2Values\leftarrow[PearsonR(trainPt,newPoint)^{2}$ for $trainPt$ in
$trainingSet$] $\triangleright$ $R^{2}$ vals

3: $avgCorr\leftarrow mean(r2Values)$ $\triangleright$ Compute average
correlation

4: Return $avgCorr$

5:end procedure

Algorithm 4 Active Learning Process Using Diversity

$TrainingData\leftarrow 3StartingPoints$ $\triangleright$ Generate initial
random training data

$Models\leftarrow RandomModels$ $\triangleright$ Generate initial random
models

$Models\leftarrow Evolve(TrainingData,Models)$ $\triangleright$ Train models
on starting data

while $BestModelError\neq 0$ do $\triangleright$ While perfect model not found

$NewPoint\leftarrow MaxDiversity(TrainingData)$ $\triangleright$ Find point of
max uncertainty

if $NewPoint\subset TrainingData$ then $\triangleright$ If point already
selected

$NewPoint\leftarrow MaxUncertainty(SubSpace(TrainingData))$$\triangleright$
Search a subspace

end if

$TrainingData\leftarrow Append(TrainingData,NewPoint)$ $\triangleright$ Add
new point

$Models\leftarrow Evolve(TrainingData,Models)$ $\triangleright$ Evolve new
models with new data using best models to seed evolution

end while

####  3.1.3 Benchmark Testing

Each active learning approach was compared on a benchmark set of 35 of the 100
equations from the Feynman Symbolic Regression Dataset [23]. These particular
35 problems were selected since they were thought to be most appropriate for a
study in active learning. In a previous study, 37 other of the 100 equations
were consistently found to need just 3 data points to be solved when using
StackGP [8]. This would render active learning useless in such cases. The
remaining 28 equations generally required all the data points up to 1000 (as
we tested) to reach moderate results, so it did not seem that this type of
active learning, adding one point at a time, would be appropriate for those
problems

###  3.2 StackGP

StackGP is a stack-based genetic programming implementation in Python [8] and
is available here (Haut).

####  3.2.1 Model Structure

Similar to PushGP (Spector), StackGP models use multiple stacks, where the
model evaluation is driven by an operator stack while variables, constants,
and other data types are stored on separate stacks. For symbolic regression
tasks, we have a total of 2 stacks, the operator stack and the
variables/constants stack.

####  3.2.2 Correlation Fitness Function

Unlike many symbolic regression implementations that use (R)MSE as the fitness
function, we employ correlation as the fitness function, together with a
linear scaling post-processing step. This was shown to perform better than
(R)MSE in earlier work [9]. The fitness is optimized during search by first
maximizing $R^{2}$, which is computed using Equation 7, where $N$ is the
number of data points $i$, $y_{i}$ is the target output, and $\hat{y}_{i}$ the
output calculated by the model.

$R=\frac{\sum_{i=1}^{N}(y_{i}-\bar{y})(\hat{y}_{i}-\bar{\hat{y}})}{\sqrt{\sum_{i=1}^{N}(y_{i}-\bar{y})^{2}\times\sum_{i=1}^{N}(\hat{y}_{i}-\bar{\hat{y}})^{2}}}$
(7)

The search is then completed using a post-processing step, which aligns the
resulting models via a simple linear regression step (eq. 8), minimizing

$\operatorname*{argmin}_{a_{0},a_{1}}\sum_{i=1}^{N}(|y_{i}-(a_{1}\hat{y_{i}}+a_{0})|)$
(8)

####  3.2.3 Algorithm

An overview of the algorithm is shown in Algorithm 5. The parameters used to
run the algorithm are shown in Table 1. Note that crossover and mutation calls
in the algorithm are simplified and actually represent applying crossover and
mutation to the correct fractions of models as shown in the parameters.

Crossover is performed using a 2-point crossover operator where two points are
selected in the operator stack of each parent and the operators, along with
the associated variables and constants between the points, are swapped between
the parents. Mutation has several different forms, each occurring with equal
probability: random replacement of a variable, random replacement of an
operator, pushing a random operator to the top of the operator stack and
pushing variables/constants to the second stack when arity is greater than 1,
popping a random number of operators off the operator stack and the correct
number of variables/constants off the second stack, inserting a single
operator at a random position in the stack, 2-point crossover with a random
model, and appending a random operator to the bottom of the operator stack.
There is then a repair mechanism that will push variables and constants to the
top of the second stack if - after mutation - there are not enough items in
the variable/constant stack for the operators.

The tournament selection method used was Pareto tournament selection, where
correlation and complexity were the two objectives. Complexity was measured as
the combined stack lengths.

Parameter | Setting
---|---
Mutation Rate | 79
Crossover Rate | 11
Spawn Rate | 10
Elitism Rate | 10
Crossover Method | 2 Pt.
Tournament Size | 5
Population Size | 300
Selection Rate | 20
Parallel Runs | 4
Generations | 1000
Table 1: StackGP & Active learning Parameter Settings Algorithm 5 StackGP
Search Algorithm

1:procedure Evolve($trainingData$,$models$)

2: for generations 1 to 100 do

3: $models\leftarrow setModelQuality(models,trainingData)$

4: $newPop\leftarrow ElitismSelection(models,20\%)$

5: $models\leftarrow tournamentSelection(models)$

6: $newPop\leftarrow newPop+crossover(models)+mutation(models)$

7: $newPop\leftarrow newPop+randomNewModels$

8: $newPop\leftarrow deleteDuplicates(newPop)$

9: $models\leftarrow newPop$

10: end for

11: $alignedModels\leftarrow alignment(models,trainingData)$

12: Return $alignedModels$

13:end procedure

###  3.3 Random Sampling

As a baseline, we used random sampling of data points from uniform and normal
distributions to determine if an active learning method improves learning
progress over a naive sampling of training data. Uniform random sampling was
chosen since it is a commonly used distribution and would likely be a first
choice for naively sampling data. A normal distribution was selected since
according to the central limit theorem, normal distributions tend to arise in
nature, so a data set sampled from natural processes would likely be a normal
distribution.

To create a fair comparison against the active learning methods, a simple
substitution was made where instead of using active learning to maximize
uncertainty or diversity, a random point was added in each iteration. Beyond
that substitution, the algorithm remains the same.

The normal distribution for each variable was defined using the midpoint
between the sampling bounds as the mean and $1/6$ of the difference between
the upper and lower bounds as the standard deviation. This places 99.8% of the
distribution between the upper and lower bounds of each variable. If a point
is sampled beyond a boundary it is adjusted to be on the boundary instead,
although this is unlikely to occur frequently.

##  4 Results and Discussion

Several different approaches for computing uncertainty and diversity were
compared using the Feynman Symbolic Regression Dataset. We then combine
diversity and uncertainty using a Pareto optimization approach and compare
that multi-objective method to using both uncertainty and diversity alone. The
Pareto approach is then tested on two additional benchmark problems from the
SRBench benchmark set.

###  4.1 Active Learning Uncertainty Sampling

The results of comparing the different uncertainty-based active learning
methods are shown in Figures 3 and 4 and the full table is in the Appendix as
Table 3. Figure 3 uses uniform random sampling as the baseline for comparison,
shown as the blue line in the figure. We also include normally distributed
random sampling for comparison as the red distribution. The results show that
the relative uncertainty measures, where we divide by the mean or trimmed
mean, do not consistently perform better than uniform random sampling. The
non-relative uncertainty measure performed well more consistently with the
methods that use differential entropy performing best. The fact that standard
deviation alone as an uncertainty metric performs consistently well is
appealing since it is very cheap and easy to implement relative to some of the
others. Differential entropy when using differential evolution as the
optimizer performed best. The fact that differential evolution as the
optimizer worked best with differential entropy likely indicates that the
surface is highly non-convex, so differential evolution was better able to
search the uncertainty space.

Figure 3: Comparing Relative Performance of Uncertainty Methods Using Uniform
Random Selection as Baseline. Shown here are the performance differences of AL
uncertainty methods compared to uniform random selection as the baseline (blue
line) and normally distributed random selection (red distribution). We see
that using the relative uncertainty measures where we divided by the mean we
get inconsistent performance, sometimes performing much better than random but
sometimes performing much worse. The non-relative approaches all consistently
perform better than random selection with the methods that use differential
entropy performing best. Using differential entropy with differential
evolution (brown) we observe the best performance. The distributions represent
the median performances of 100 independent runs across all test problems. For
completeness, there is one point not shown for the std/tr. mean approach that
is around -200. Figure 4: Comparing Performance of Uncertainty Methods Against
Uniform Random Selection. Each method is compared to uniform random sampling
and the number of times that the method outperforms and underperforms is
reported. The number of times each method outperforms is shown on the left and
the number of times each method underperforms is shown on the right.
Outperforming means that a method used fewer points than uniform random
sampling. Underperforming means that it required more points. Ties are not
counted but can be easily determined by taking the difference of 35 and the
two values reported. The results show that the methods that use differential
entropy work well most consistently, outperforming more frequently and
underperforming infrequently. We can also see that the relative uncertainty
measures were very inconsistent in their performance.

Figure 4 compares the performance of each method against uniform random
sampling for each problem and displays the number of times each method
outperforms or underperforms random sampling. If a method outperforms random
sampling that means that the method required fewer points to solve a problem.
If a method underperforms random sampling that means that the method required
more points to solve a problem. The results show that the methods using
differential entropy work best, outperforming in the most number of cases and
underperforming in the fewest number of cases. The differential entropy method
that used differential evolution as the optimizer worked better than just
using differential entropy with SciPy Optimize’s minimize function. This
indicates that differential evolution was able to search the uncertainty
surface more effectively. The results also show that the relative uncertainty
methods that divided the mean or trimmed mean were not consistent in their
performance, frequently having a similar number of cases where the methods
outperformed and underperformed.

We see that the relative measures sometimes perform well and sometimes perform
poorly, but on average they are centered around the baseline performance. The
original assumption was that the relative uncertainty measures would be
appealing since it was thought that they would reduce a bias towards selecting
points where the predicted response is larger and thus naturally leads to
wider distributions of the ensemble. This may have been the case occasionally
where those methods did perform much better than uniform random sampling, but
they ware not consistent. Looking at their formulations there is a risk of
selecting points where the mean is near 0 which results in asymptotic behavior
of the uncertainty function.

Considering the results, we also see that of the two random sampling methods,
normally distributed random sampling seems to perform a bit better than
uniform sampling. This indicates that if a researcher does not want to use
active learning to guide their data collection, they would typically be better
off using a normal distribution than a uniform distribution for their samples.

###  4.2 Active Learning Diversity Sampling

The different metrics for determining point diversity were compared to
determine if there are clear differences in what they are measuring and also
to ensure there aren’t any obvious flaws with any of the metrics. When
comparing minimum distance and average distance an initial randomly generated
training set with 3 data points in 3 dimensions was generated. Figure 5 shows
the comparison where new points were selected iteratively to add to the
training set using the minimum distance metric for selection. We can see that
the correlation, $R^{2}$ is actually pretty weak between the two, indicating
they are providing different measures. As well, we recorded the Spearman Rho,
rank-correlation, since that indicates if the methods are ranking points
similarly or not. If methods rank points similarly, then they would likely not
provide unique information if used as a diversity metric. It was found that
the Spearman Rho was 0.44, which means that the two methods are ranking points
differently and could provide unique information.

Figure 5: Comparing minimum Euclidean distance against mean Euclidean distance
as a diversity metric. Here minimum distance is used to select the next point
in the set and both metrics of those points are displayed. We can see that
there is little correlation between the two metrics indicating they provide
different information. The $R^{2}$ between these two metrics on these points
is just 0.37. The Spearman Rho, rank-correlation, is also low at 0.44. Figure
6: Comparing minimum Euclidean distance against mean Euclidean distance as a
diversity metric. Here mean distance is used to select the next point in the
set and both metrics of those points are displayed. We can see that when mean
distance is used to select new points, we get many points with a minimum
distance of 0. This indicates that we are very frequently reselecting points
already in the set. This shows that minimum distance is a better metric than
mean distance.  Figure 7: Comparing minimum Euclidean distance against mean
point correlation as a diversity metric. Here minimizing mean correlation is
used to select the next point in the set and both metrics of those points are
displayed. We can see that there appears to be a weak positive correlation
between the two, indicating that they provide some of the same information but
are not the same, so may have different advantages. It is also promising that
the minimum distance shows that we are not reselecting points already in the
training set. Comparing the metrics for these points we get an $R^{2}$ of 0.35
and a Spearman Rho of 0.33.

To further compare the minimum and mean distance metrics, the analysis was
flipped, such that mean distance was used to select new points and both
metrics were recorded on the selected points. These results are shown in
Figure 6. Here it becomes obvious that mean distance is not a good metric
since the minimum distance metric indicates that we are repeatedly selecting
points already in the set. This is shown by the consistent minimum distance
value of 0 after around 10 iterations. This result led to mean distance being
thrown out as a potential choice of metric.

Minimum distance and correlation were also compared to determine if they
provide unique measure of diversity. The results are shown in Figure 7. For
this analysis, lack of correlation to the training set was used to select new
points and both metrics were recorded. This analysis was slightly different
than the previous ones since for this problem the points were embedded in a 10
dimensional space instead of just 3. The results show that the two metric do
provide unique information since an $R^{2}$ value of 0.35 and a Spearman Rho
value of 0.33 were recorded, which are both low. Since these metrics were
determined to provide unique information without any clear flaws both were
included to be explored, with the one limitation that correlation as a
diversity metric could not be used on problems of less than 3 dimensions.

Figure 8: Comparing Relative Performance of Diversity Methods Using Uniform
Random Selection as Baseline. Shown here are the performance differences of
both the AL diversity methods compared to uniform random selection as the
baseline (blue line) and normally distributed random selection (red
distribution). We see that using minimum distance (green distribution)
performs consistently better than the baseline and correlation (blue
distribution) works best as a diversity metric. The drawback with using
correlation as the diversity metric though is that it requires problems with
more than two dimensions, so the problems with two dimensions are ignored when
using correlation. The distributions represent the median performances of 100
independent runs across all test problems.

The results of comparing the different data diversity-based active learning
methods are summarized in Figures 8 and 9 and the full results are shown in
Table 4 in the Appendix. Figure 8 uses uniform random sampling as the baseline
for comparison, shown as the blue line. We again include normally distributed
random sampling for comparison as the red distribution. We can see that both
diversity metrics have better performance than uniform random sampling, on
average requiring fewer training points to find a solution. We also see that
correlation as a diversity metric performs best, often requiring the least
number of training data points to find a solution. Correlation does have the
disadvantage, though, of not working on the problems with just two dimensions.
Those two problems are not represented in the correlation bar in the chart
since they are not applicable.

Figure 9 shows the number of cases where each method either outperformed or
underperformed when compared to uniform random sampling. We see again that
correlation has the best performance. This indicates that not only does
correlation lead to requiring fewer training points on average, but also
indicates that it most consistently requires fewer points. We see that
distance as a metric requires fewer points than uniform and random sampling,
but is not as consistent as correlation.

Figure 9: Comparing Performance of Diversity Methods Against Uniform Random
Selection. Each method is compared to uniform random sampling and the number
of times that the method outperforms and underperforms is reported. The number
of times each method outperforms is shown on the left and the number of times
each method underperforms is shown on the right. Outperforms means that a
method used fewer points than uniform random sampling. Underperforms means it
required more points. Ties are not counted but can be easily determined by
taking the difference of 35 and the two values reported. The results show that
correlation performed best, underperforming the fewest times and outperforming
the most.

###  4.3 Comparing Diversity, Uncertainty, and Pareto Optimization of Both

Next we explore how the performance compares when using uncertainty and
diversity together to see if there are benefits to considering both for
selecting training data with AL compared to just uncertainty or diversity
alone. For this comparison, we selected one diversity metric and one
uncertainty metric. For the uncertainty metric, we chose differential entropy
since it was shown to be the best performing metric in Figure 3. For the
diversity metric, we chose minimum distance. Although it didn’t perform best,
it is most versatile since it isn’t restricted to problems with more than 2
dimensions. For the combination method, we used a Pareto optimization to find
the points with the best trade-off of both the uncertainty and diversity
metrics from 10,000 randomly generated points each iteration. From the Pareto
front of points that are non-dominated in those two objectives, we ordered
them based on their uncertainty score and selected the median point. Note that
sorting based on uncertainty is just the reverse order of a sort by diversity,
so which objective you choose to sort by shouldn’t have a significant impact.
The only impact would be on cases where an even number of points are on the
front so the point you select isn’t the true median but rather one of the
points near the median. When this occurs, we round down to select the median
point, which would give a slight bias toward uncertainty. By selecting the
median point we are attempting to choose a point that has a relatively good
balance between the two objectives.

The results of this comparison are shown in Figures 10 and 11 with the results
from each problem shown in the Appendix in Table 5. Again in Figure 10, we use
uniform random sampling as the baseline (blue line) and include normally
distributed random sampling for comparison. The results show that all three
methods work better than the baseline and normally distributed random
sampling. Using the uncertainty metric, differential entropy, works slightly
better than using the distance metric, minimum distance. We also see that
there is a benefit to combining both metrics using the Pareto optimization
since we see an improvement in the upper quartile of performance. It is also
interesting to note, as can be seen in Figure 11, that the diversity metric
alone performed worse than uniform random sampling in 8 of the 35 cases,
whereas the uncertainty approach and the Pareto approach only performed worse
in 4 of the cases, demonstrating that the uncertainty and Pareto approaches
offer more consistent improvements. This indicates that it is important to
consider the current models to help guide the AL process. This makes sense
since the goal is to select training points that will best inform the current
model population, using only diversity doesn’t consider the current state of
models, so it is less likely that the training points selected will most
inform those models. Statistical significance tests were also performed and
the number of cases determined to be statistically significant are shown in
the darker regions in the figure. The Mann-Whitney test was used to test for
significance and a threshold of 0.05 was used. The Pareto approach was found
to be statistically significant in 18 of the 20 cases where the Pareto
approach outperformed.

Figure 10: Comparing Relative Performance of Diversity, Uncertainty, and
Pareto Optimization Using Uniform Random Selection as Baseline. Shown here are
the performance differences of AL diversity, uncertainty and Pareto methods
compared to uniform random selection as the baseline (blue line) and normally
distributed random selection (red distribution). We see that using the
diversity metric, minimum distance (green distribution), performs consistently
better than the baseline and the uncertainty metric, DE (blue distribution),
performs a bit better than the diversity method. When using a Pareto
optimization of both diversity and uncertainty we get even better performance.
The distributions represent the median performances of 100 independent runs
across all test problems. For completeness, there is a single point around
-150 for the Pareto approach.  Figure 11: Comparing Performance of Diversity,
Uncertainty, and Pareto Optimization Against Uniform Random Selection. Each
method is compared to uniform random sampling and the number of times that the
method outperforms and underperforms is reported. The number of cases where
the differences are statistically significant is shown in the darker regions.
The number of times each method outperforms is shown on the left and the
number of times each method underperforms is shown on the right. Outperforms
means that a method used fewer points than uniform random sampling.
Underperforms means it required more points. Ties are not counted but can be
easily determined by taking the difference of 35 and the two values reported.
The results show that DE, the uncertainty method works best. The Pareto
approach ties for the least number of underperforming cases, matching DE, and
outperforms between DE and Min. Distance. Statistical significance was
determined using the threshold of 0.05 with the Mann-Whitney test.

Looking at the results, there are two instances where the Pareto approach
performed considerably worse than the uncertainty and diversity approaches.
Those are equations 9 and 71. Table 5 in the Appendix shows that the combined
method performs worse than focusing alone on either diversity or uncertainty
for those two problems. This is likely a result of equations 9 and 71 being
higher dimensional problems with 6 and 5 dimensions, respectively, so the
10,000 randomly generated points don’t sufficiently fill the search space to
find points with high values for both uncertainty and diversity.

Equation 71 was further explored to see if sampling additional points improved
the performance when using the combined diversity uncertainty approach and to
verify that sparse sampling was at least part of the issue as suspected.
Equation 71 was retested using 100,000 randomly sampled points to search for
the best trade-off between diversity and uncertainty. When using 100,000
points the median number of points required to solve the problem decreased to
42 points from 50.5, confirming that better sampling of the space improves the
performance in this higher dimensional problem. The median performance of 42
points is still worse than either of the uncertainty or diversity approaches,
so more points could be used, but increasing the number of points beyond
100,000 begins to make that search rather expensive. Rather than randomly
sampling the points then selecting the Pareto front from those points, an
alternative optimization method, such as NSGA II [4], could be used in future
studies which might be cheaper and likely more effective.

###  4.4 Additional Benchmark Problems

To further test the Pareto AL approach, we selected two problems from a more
recent benchmark set, SRBench [12]. One that is on the easier side for StackGP
and one that is a bit more challenging. The easier problem selected was the
van der Pol oscillator problem, referred to as ”strogatz_vdp1” in SRBench. The
equation for the van der Pol oscillator problem that we are trying to
rediscover is $x^{\prime}=10*(y-(1)/(3)*(x^{3}-x))$. The more challenging
problem was the bar magnet problem, referred to as ”strogatz_barmag1” in
SRBench and the equation for the bar magnet problem that we are trying to
rediscover is $x^{\prime}=0.5*sin(x-y)-sin(x)$. As with the previous problems,
we performed each experiment 100 times and computed the median number of
points to find the solution. The results of those experiments are shown in
Table 2. We can see that the Pareto approach performs significantly better
than randomly sampling from a normal distribution and performs about 27.8%
better than randomly sampling from a uniform distribution on the bar magnet
problem. The performance gains over the normal and uniform distributed
samplings are statistically significant considering a threshold of 0.05 using
the Mann-Whitney test. We computed a p-value of 3.490$*{10}^{-11}$ when
comparing to the normal distribution and 6.481$*{10}^{-6}$ when comparing to
the uniform sampling. We also see better performance on the van der Pol
oscillator, although since it was an easy problem there isn’t as much
opportunity for improvement, so we only see a reduction of a few points. The
performance gains over the normal and uniform distributions are again
statistically significant with a p-value of 2.51$*{10}^{-7}$ when compared
with the results from using normally distributed sampling and a p-value of
4.008$*{10}^{-13}$ when compared with the results from using uniform random
sampling.

SRBench | N. Ran | U. Ran | Pareto AL
---|---|---|---
Problem | Data Pts. | Data Pts. | Data Pts.
Bar Magnet #1 | 51 | 18 | 13
Van der Pol Osc. #1 | 10 | 9 | 7
Table 2: Shown are the median numbers of points needed to solve each equation.
A total of 100 independent trials were performed for each equation. We compare
the active learning method that uses both diversity and uncertainty and
compare the performance against random sampling on two problems from the
SRBench.

##  5 Conclusion

Both uncertainty and diversity metrics for active learning were explored to
see how each metric impacts the success of active learning in genetic
programming. As well, a Pareto approach was defined that allows both diversity
and uncertainty to be considered for active learning. Of the uncertainty
approaches, it was observed that differential entropy performed best. It was
also observed that relative uncertainty functions did not perform well. When
using differential entropy it was found that performance could be boosted by
using differential evolution as the optimizer over Scipy Optimize’s minimize
function. This indicates that the search space is not convex and requires a
good optimizer to find solutions with high uncertainty.

When comparing the data diversity methods, it was found that correlation
performed better than minimum Euclidean distance. Although correlation worked
better, it does not work on cases with 2 dimensions or less. Thus, minimum
Euclidean distance was selected for the Pareto approach. Future
implementations may default to using minimum Euclidean distance for all cases
with 1 or 2 dimensions and using correlation for higher dimensional problems.
Mean distance was considered, but determined to be uninformative due to its
frequency of identifying repeat points.

When comparing the Pareto approach which used both differential entropy and
minimum Euclidean distance to differential entropy, minimum Euclidean
distance, uniform random selection, and normally distributed random selection,
it was found that differential entropy worked best, with the Pareto approach
performing between differential entropy and minimum Euclidean distance.
Looking at individual problems, there were a few cases where the Pareto
approach actually worked better than both differential entropy and minimum
Euclidean distance on their own, indicating potential benefits of combining
the two approaches. For the cases where the Pareto approach did not work as
well, it was identified that the multi-objective optimization strategy may
have been at fault since it relies on randomly generating N points and
selecting the median value in the Pareto front. Better methods such as NSGA-II
could be explored in future studies to see if improved optimization methods
leads to better active learning performance.

Overall, it was found that active learning can be efficiently utilized with
genetic programming to reduce training data requirements. In practice, this
would be useful to apply in scenarios where collecting data or labelling data
is expensive, and model training is relatively cheap. In these scenarios,
active learning could be used to guide data collection and labelling so that
good models can be arrived at using as few data points as possible. This
application has the potential to accelerate data driven research, since it
could lead to finding solutions with fewer resources in less time.

Acknowledgments

Computer support by MSU’s iCER high-performance computing center is gratefully
acknowledged.

Data availability The datasets generated/analysed during the current study are
available from the corresponding author on reasonable request. Code
availability The code for StackGP with active learning can be found here:
https://github.com/hoolagans/StackGP

##  6 Appendix

EQ | U. Ran | N. Ran | Std/Mean | TrStd/TrMean | Std/TrMean | Std | DE | DE (DE)
---|---|---|---|---|---|---|---|---
Num | Data | Data | Data | Data | Data | Data | Data | Data
2 | 54.5 | 97.5 | 50 | 39 | 47 | 53 | 82.5 | 47
3 | $>1000$ | $>1000$ | 876 | 692 | 741 | $>1000$ | $>1000$ | 724.5
4 | 30 | 21.5 | 21.5 | 20 | 23 | 20 | 28 | 19.5
7 | 88.5 | 82 | 23 | 21 | 22.5 | 39 | 52.5 | 35
9 | 120.5 | 155.5 | 150.5 | 73.5 | 359.5 | 100.5 | 153 | 160
10 | 6 | 6 | 6 | 11 | 6 | 7 | 6 | 6
13 | 13 | 14.5 | 15 | 15 | 14 | 14 | 12 | 12
14 | 30.5 | 33 | 28 | 24 | 31 | 23.5 | 24 | 22.5
23 | 8 | 8 | 7 | 8 | 7 | 8 | 7.5 | 7
24 | 49 | 58 | 39 | 29.5 | 31 | 26 | 22 | 28
27 | 30 | 17 | 20 | 13 | 19.5 | 18 | 14 | 15
32 | 17 | 16 | 20 | 18 | 21 | 16 | 18 | 12
35 | 19 | 17 | 17 | 6 | 21 | 18 | 13.5 | 12.5
39 | 10 | 10 | 10 | 12 | 11 | 10 | 9 | 9
41 | 7 | 7 | 7 | 8 | 7 | 8 | 7 | 7
43 | 453 | 82 | 876 | 218 | 202.5 | 144 | 326 | 192.5
47 | 13 | 12 | 14 | 12 | 13 | 13 | 12 | 12
48 | 15.5 | 14 | 18 | 17 | 17.5 | 14 | 13 | 12.5
52 | 9.5 | 9 | 10 | 10 | 9.5 | 10 | 9 | 9
55 | 10 | 10 | 11 | 12 | 10 | 11 | 10 | 9
57 | 30.5 | 31.5 | 25.5 | 27 | 24 | 17 | 23 | 21
60 | 7 | 7 | 7 | 7 | 7 | 7 | 7 | 7
61 | 18.5 | 20 | 20 | 19 | 18 | 18 | 16 | 17
62 | 34.5 | 56.5 | 37.5 | 34.5 | 33 | 34 | 30 | 28.5
63 | 14 | 13 | 15 | 16 | 15.5 | 14 | 13 | 13
66 | 11 | 9 | 15 | 14 | 14 | 15 | 10 | 9
67 | 10.5 | 11 | 11 | 10 | 10 | 10 | 11 | 11
71 | 51.5 | 47 | 34 | 58 | 38.5 | 31 | 30 | 35
83 | 5 | 5 | 5 | 5 | 5 | 5 | 5 | 5
85 | 4 | 4 | 4 | 4 | 4 | 4 | 4 | 4
89 | 5 | 5 | 4 | 5 | 4 | 4.5 | 5 | 5
93 | 8 | 8 | 8 | 8 | 8 | 7 | 8 | 8
95 | 11 | 9 | 12 | 11.5 | 11 | 12 | 10 | 11
98 | 8 | 7 | 9 | 9 | 9 | 9 | 7 | 7.5
99 | 30 | 20 | 31 | 35.5 | 35 | 25 | 24 | 21
Vs. |  |  |  |  |  |  |  |
U. Sampl. | - | - | 19 | 20 | 24 | 27 | 31 | 33
Vs. |  |  |  |  |  |  |  |
N. Sampl. | - | - | 22 | 20 | 20 | 22 | 29 | 30
Table 3: Shown are the median number of points needed to solve each equation. A total of 100 independent trials were performed for each equation. The last row indicates the number of cases where each of the active learning methods matched or performed better than uniform random sampling of training data. Where informative, the minimum number is bolded.  EQ | U. Rand | Pt. Dist | Pt. Corr
---|---|---|---
Num | Data Pts. | Data Pts. | Data Pts.
2 | 54.5 | 44 | -
3 | $>1000$ | $>1000$ | $>1000$
4 | 30 | 19 | 29.5
7 | 88.5 | 35.5 | 60
9 | 120.5 | 210.5 | 102.5
10 | 6 | 6 | 5
13 | 13 | 12 | 14
14 | 30.5 | 24 | 22
23 | 8 | 7 | 8
24 | 49 | 23 | 26.5
27 | 30 | 13.5 | 11.5
32 | 17 | 15.5 | 14
35 | 19 | 13.5 | 15
39 | 10 | 11 | 9
41 | 7 | 8 | 6
43 | 453 | 533.5 | 136.5
47 | 13 | 14.5 | 16
48 | 15.5 | 13 | 13
52 | 9.5 | 10 | 9
55 | 10 | 10 | 10
57 | 30.5 | 28 | 29.5
60 | 7 | 7 | 7
61 | 18.5 | 17.5 | 17.5
62 | 34.5 | 29.5 | 36.5
63 | 14 | 14 | 12
66 | 11 | 12 | 10
67 | 10.5 | 11 | 15
71 | 51.5 | 29 | 30.5
83 | 5 | 5 | 5
85 | 4 | 4 | -
89 | 5 | 5 | 5
93 | 8 | 7.5 | 7
95 | 11 | 10 | 8
98 | 8 | 8 | 7
99 | 30 | 25 | 20.5
Perf. Count | - | 27/35 | 29/33
Table 4: Shown are the median number of points needed to solve each equation. A total of 100 independent trials were performed for each equation. There are 2 equations that have a dash instead of a number and that is because they have only two dimensions, so selecting points with minimal correlation to the rest of the training set is not possible. The approach using uniformly random data points was included in the first column represented as a baseline. The last row indicates the number of cases where each of the point diversity methods matched or performed better than the random approach. EQ | Pt. Dist | Pareto | Pt. Unc.
---|---|---|---
Num | Data Pts. | Data Pts. | Data Pts.
2 | 44 | 36.5 | 82.5
3 | $>1000$ | 501 | $>1000$
4 | 19 | 29 | 28
7 | 35.5 | 48.5 | 52.5
9 | 210.5 | 304 | 153
10 | 6 | 6 | 6
13 | 12 | 12 | 12
14 | 24 | 22 | 24
23 | 7 | 8 | 7.5
24 | 23 | 27 | 22
27 | 13.5 | 19 | 14
32 | 15.5 | 14 | 18
35 | 13.5 | 14 | 13.5
39 | 11 | 10 | 9
41 | 8 | 7 | 7
43 | 533.5 | 314.5 | 326
47 | 14.5 | 12 | 12
48 | 13 | 13 | 13
52 | 10 | 10 | 9
55 | 10 | 9 | 10
57 | 28 | 24 | 23
60 | 7 | 7 | 7
61 | 17.5 | 15 | 16
62 | 29.5 | 21.5 | 30
63 | 14 | 14 | 13
66 | 12 | 10 | 10
67 | 11 | 11 | 11
71 | 29 | 50.5 | 30
83 | 5 | 5 | 5
85 | 4 | 4 | 4
89 | 5 | 5 | 5
93 | 7.5 | 8 | 8
95 | 10 | 13 | 10
98 | 8 | 8 | 7
99 | 25 | 25.5 | 24
Worst Count | 13 | 11 | 8
Best Count | 16 | 19 | 21
Table 5: Shown are the median number of points needed to solve each equation. A total of 100 independent trials were performed for each equation. Here the trade-off between diversity and uncertainty is explored. The second to last row indicates the number of times each approach was the worst of the three approaches. The last row indicates the number of cases where each approach was the best or tied for the best of the three approaches. Minimum point distance was used for the diversity metric and differential entropy was used as the uncertainty metric. EQ | p-value | Significant
---|---|---
Num |  |
2 | 9.18077$*{10}^{-4}$ | Yes
3 | 1.4012505382465074$*{10}^{-11}$ | Yes
4 | 0.612865 | No
7 | 9.606121209101481$*{10}^{-6}$ | Yes
9 | 1.78579$*{10}^{-2}$ | Yes
10 | 3.40099$*{10}^{-5}$ | Yes
13 | 2.38254$*{10}^{-5}$ | Yes
14 | 7.46075$*{10}^{-4}$ | Yes
23 | 3.3175302779045554$*{10}^{-6}$ | Yes
24 | 8.64911$*{10}^{-5}$ | Yes
27 | 0.0354696 | Yes
32 | 1.30021$*{10}^{-3}$ | Yes
35 | 2.90356$*{10}^{-3}$ | Yes
39 | 7.61664$*{10}^{-3}$ | Yes
41 | 5.872949580996268$*{10}^{-7}$ | Yes
43 | 3.38224$*{10}^{-3}$ | Yes
47 | 1.91012$*{10}^{-4}$ | Yes
48 | 1.78267$*{10}^{-5}$ | Yes
52 | 0.433203 | No
55 | 1.98428$*{10}^{-4}$ | Yes
57 | 2.51679$*{10}^{-5}$ | Yes
60 | 6.84394$*{10}^{-3}$ | Yes
61 | 1.39827$*{10}^{-3}$ | Yes
62 | 2.6897830944058467$*{10}^{-20}$ | Yes
63 | 0.0399119 | Yes
66 | 1.55871$*{10}^{-3}$ | Yes
67 | 0.245499 | No
71 | 0.392183 | No
83 | 9.35933923552931$*{10}^{-10}$ | Yes
85 | 1.059722812150408$*{10}^{-12}$ | Yes
89 | 0.140942 | No
93 | 1.54882$*{10}^{-5}$ | Yes
95 | 0.0345741 | Yes
98 | 2.98341$*{10}^{-3}$ | Yes
99 | 5.63906$*{10}^{-3}$ | Yes
Significance Count | 30/35 |
Table 6: Statistical significance of Pareto AL approach vs. uniform random
sampling. We are using a threshold of 0.05 to test for significance. The Mann-
Whitney test was used to test for significance.

## References

  *   * Cohn et al [1996] Cohn DA, Ghahramani Z, Jordan MI (1996) Active learning with statistical models. Journal of Artificial Intelligence Research 4(1):129–145
  * Curry et al [2007] Curry R, Lichodzijewski P, Heywood MI (2007) Scaling genetic programming to large datasets using hierarchical dynamic subset selection. IEEE Transactions on Systems, Man, and Cybernetics, Part B 37(4):1065–1073. 10.1109/TSMCB.2007.896406
  * De Freitas et al [2010] De Freitas J, Pappa GL, da Silva AS, et al (2010) Active learning genetic programming for record deduplication. In: IEEE Congress on Evolutionary Computation, pp 1–8, 10.1109/CEC.2010.5586104
  * Deb et al [2002] Deb K, Pratap A, Agarwal S, et al (2002) A fast and elitist multiobjective genetic algorithm: NSGA-II. IEEE Transactions on Evolutionary Computation 6(2):182–197. 10.1109/4235.996017
  * El-Hasnony et al [2022] El-Hasnony IM, Elzeki OM, Alshehri A, et al (2022) Multi-label active learning-based machine learning model for heart disease prediction. Sensors (Basel) 22(3)
  * Eyke et al [2020] Eyke N, Green W, Jensen K (2020) Iterative experimental design based on active machine learning reduces the experimental burden associated with reaction screening. Reaction Chemistry & Engineering 5:1963–1972
  * Haut [2021] Haut N (2021) StackGP. URL https://github.com/hoolagans/StackGP
  * Haut et al [2022a] Haut N, Banzhaf W, Punch B (2022a) Active learning improves performance on symbolic regression tasks in stackgp. In: Proceedings of the Genetic and Evolutionary Computation Conference Companion. Association for Computing Machinery, New York, NY, USA, GECCO ’22, p 550–553
  * Haut et al [2022b] Haut N, Banzhaf W, Punch B (2022b) Correlation versus RMSE Loss Functions in Symbolic Regression Tasks. In: Banzhaf W, Trujillo L, Winkler S, et al (eds) Genetic Programming Theory and Practice XIX. Springer
  * Kotanchek et al [2009] Kotanchek M, Smits G, Vladislavleva E (2009) Exploiting trustable models via pareto gp for targeted data collection. In: Riolo R, Soule T, Worzel B (eds) Genetic Programming Theory and Practice VI. Springer, p 145–162
  * Kremer et al [2014] Kremer J, Steenstrup Pedersen K, Igel C (2014) Active learning with support vector machines. WIREs Data Mining and Knowledge Discovery 4(4):313–326
  * La Cava et al [2021] La Cava W, Orzechowski P, Burlacu B, et al (2021) Contemporary symbolic regression methods and their relative performance. 10.48550/ARXIV.2107.14351, URL https://arxiv.org/abs/2107.14351
  * Lasarczyk et al [2004] Lasarczyk CW, Dittrich P, Banzhaf W (2004) Dynamic subset selection based on a fitness case topology. Evolutionary Computation 12(2):223 – 242
  * Li et al [2022] Li Y, Wang X, Shi Z, et al (2022) Boosting training for pdf malware classifier via active learning. International Journal of Intelligent Systems 37(4):2803–2821
  * Nandy et al [2022] Nandy A, Duan C, Goffinet C, et al (2022) New strategies for direct methane-to-methanol conversion from active learning exploration of 16 million catalysts. Journal of the American Chemical Society Au 2(5):1200–1213
  * Ren et al [2021] Ren P, Xiao Y, Chang X, et al (2021) A survey of deep active learning. ACM Computing Surveys 54(9)
  * SciPy [Last Access 2023a] SciPy (Last Access 2023a) scipy.optimize.differential_evolution - scipy v1.11.1 manual. URL https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.differential_evolution.html
  * SciPy [Last Access 2023b] SciPy (Last Access 2023b) scipy.optimize.minimize - scipy v1.11.1 manual. URL https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html
  * Settles [2009] Settles B (2009) Active learning literature survey. Computer Sciences Technical Report 1648, University of Wisconsin–Madison
  * Spector [2005] Spector L (2005) PushGP. URL http://faculty.hampshire.edu/lspector/push.html
  * Sverchkov and Craven [2017] Sverchkov Y, Craven M (2017) A review of active learning approaches to experimental design for uncovering biological networks. PLOS Computational Biology 13:1–26
  * Tegmark [2020] Tegmark M (2020) Welcome to the Feynman Symbolic Regression Database! URL https://space.mit.edu/home/tegmark/aifeynman.html#:~:text=As%20opposed%20to%20linear%20regression,any%20combination%20of%20mathematical%20symbols.
  * Udrescu and M. [2020] Udrescu SM, M. T (2020) A physics-inspired method for symbolic regression. Science Advances 6:eaay2631