Astro2020 Science White Paper

Studying Magnetic Fields in Star Formation and the Turbulent Interstellar
Medium

Thematic Areas: $\square$ Planetary Systems $\blacksquare$ Star and Planet
Formation $\square$ Formation and Evolution of Compact Objects $\square$
Cosmology and Fundamental Physics $\square$ Stars and Stellar Evolution
$\square$ Resolved Stellar Populations and their Environments $\square$ Galaxy
Evolution $\square$ Multi-Messenger Astronomy and Astrophysics

Principal Author:

Name: Laura Fissel Institution: NRAO Email<EMAIL_ADDRESS>Phone:
(434)244-6821

Co-authors: Charles L. H. Hull (NAOJ/ALMA), Susan E. Clark (Institute for
Advanced Study), David T. Chuss (Villanova), Philippe André (CEA Saclay),
François Boulanger (Universite Paris Sud/IAS), C. Darren Dowell (JPL/Caltech),
Edith Falgarone (ENS), Brandon Hensley (Princeton), A. Lazarian (U. Wisconsin,
Madison), Giles Novak (Northwestern), Ian Stephens (Harvard CfA), Siyao Xu (U.
Wisconsin, Madison)

Abstract: Understanding the physics of how stars form is a highly-prioritized
goal of modern Astrophysics, in part because star formation is linked to both
galactic dynamics on large scales and to the formation of planets on small
scales. It is well-known that stars form from the gravitational collapse of
molecular clouds, which are in turn formed out of the turbulent interstellar
medium. Star formation is highly inefficient, with one of the likely culprits
being the regulation against gravitational collapse provided by magnetic
fields. Measurement of the polarized emission from interstellar dust grains,
which are partially aligned with the magnetic field, provides a key tool for
understanding the role these fields play in the star formation process.

Over the past decade, much progress has been made by the most recent
generation of polarimeters operating over a range of wavelengths (from the
far-infrared through the millimeter part of the spectrum) and over a range of
angular resolutions (from less than an arcsecond through fractions of a
degree). Future developments in instrument sensitivity for ground-based,
airborne, and space-borne polarimeters operating over range of spatial scales
are critical for enabling revolutionary steps forward in our understanding of
the magnetized turbulence from which stars are formed.

##  1 Introduction

The magnetic field is a key regulator of star formation at all spatial scales
(Elmegreen & Scalo, 2004; Crutcher, 2012). Magnetic fields shape turbulence
over a broad range of length scales in the diffuse interstellar medium (ISM)
and are a critical factor in the formation and evolution of molecular clouds
as well as protostellar cores, envelopes, disks, outflows, and jets. The
efficiency of the conversion of molecular gas to stars during the star-
formation process is notoriously low, due to regulation from a variety of
sources including turbulent gas motions, magnetic fields, and feedback from
young stars (McKee & Ostriker, 2007; Federrath, 2015). In particular, magnetic
fields play an important role in slowing the process of star formation by
inhibiting movement of gas in the direction perpendicular to the field lines.

Over the last two decades a vast amount of work has focused on characterizing
magnetic fields and turbulence in the ISM from galactic to cloud to
protostellar core scales, using interferometers as well as ground-, balloon-,
and space-based single-dish instruments. By measuring polarized thermal
emission from dust grains aligned with respect to the magnetic field, we can
map the magnetic field morphology in star forming regions (Jones et al., 1984;
Lazarian & Hoang, 2007; Andersson et al., 2015). The Planck satellite recently
mapped dust polarization across the entire Milky Way at a best resolution of
5′ (Planck Collaboration Int. XIX et al., 2015). However, Planck lacks the
sensitivity needed to study magnetic fields and turbulence in detail in all
but the very nearest molecular clouds, and cannot resolve dense filamentary
structures and cores in any cloud. In contrast, ground-based polarimeters have
produced high-resolution maps of dense sub-regions within molecular clouds;
however, such instruments are sensitivity-limited as a result of observing
through the atmosphere, and can therefore typically only map very small
($<\,$0.01 deg2), bright regions within molecular clouds (e.g. Dotson et al.
2010; Matthews et al. 2009).

High-sensitivity observations of a large number (of order thousands) of
clouds, in addition to extremely well resolved maps of the well known set of
nearby molecular clouds, will allow us to answer key questions: How do
magnetic fields affect the formation of molecular clouds from the diffuse
interstellar medium? And do magnetic fields influence the formation and
evolution of gravitationally unstable cloud sub-structures, and consequently,
the efficiency with which stars form?

##  2 The diffuse magnetized medium

Molecular clouds condense out of the diffuse ISM, and so the process of star
formation is regulated in its earliest phases by the flow of material into the
molecular phase. Magnetic fields influence structure formation directly,
especially through magnetic pressure, and via their influence on the
acceleration and propagation of cosmic rays (Strong et al., 2007) and the
transfer of heat and matter between the hot, warm, and cold phases of the ISM
(Draine, 2011). In addition to being the fuel for star formation, the diffuse
ISM is also the medium that responds to supernova explosions and other large-
scale injections of energy. Despite its importance, a comprehensive
understanding of the magnetized diffuse ISM is challenging because of its
diverse composition, its sheer expanse, and the multi-scale nature of the
physics that shapes it.

Recent high-resolution observations of the diffuse ISM have revealed a wealth
of complex structure, deeply influenced by magnetized turbulence and the
structure of the ambient magnetic field. The structure of the cold neutral
medium (CNM) in particular is highly anisotropic, largely organized into
filamentary structures that are aligned with the magnetic field (Clark et al.,
2014). A preferred alignment between ISM structures and the magnetic field
traced by dust polarization has recently been measured in both neutral
hydrogen (Clark et al., 2015; Kalberla et al., 2016) and dust emission (Planck
Collaboration Int. XXXII, 2016; Malinen et al., 2016). Sensitive observations
of polarized dust emission toward diffuse regions of the ISM are crucial for
understanding the dynamical importance of the magnetic field for structure
formation at a given scale.

The diffuse ISM is an ideal laboratory for studying magnetohydrodynamic (MHD)
turbulence. Turbulence correlates matter and energy over a range of scales,
traditionally parameterized as a cascade of energy across the wavenumber
spectrum (Elmegreen & Scalo, 2004). Studying turbulence necessitates high
dynamic range observations in order to have robust statistical samples on many
different scales, and dust polarization provides one of the few avenues for
high-dynamic-range magnetic field observations. One statistical measure is the
cross-correlation between dust temperature ($T$) and the dust polarization
field, decomposed into the rotationally invariant $E$\- and $B$-modes. Planck
measured these quantities, with several unexpected results. There is a
positive $TE$ correlation in the diffuse ISM, as well as a non-unity ratio of
$E$-mode to $B$-mode polarization ($EE/BB\sim 2$). The $TE$ and $EE/BB$
correlations are likely due to magnetically induced anisotropy in the dust
distribution (Clark et al., 2015; Planck Collaboration Int. XXXVIII, 2016),
and to MHD turbulence parameters like the Alfvén Mach number (e.g., Caldwell
et al., 2017; Kandel et al., 2017). Perhaps even more surprising, Planck
measured a nonzero $TB$ correlation across a range of scales in the diffuse
ISM (Planck Collaboration 2018 XI, 2018). The nonzero $TB$ measurement could
be a signature of a large-scale helical magnetic field (e.g. Bracco et al.,
2019).

Detailed studies of the diffuse ISM are highly complementary to the search for
inflationary $B$-mode polarization in the cosmic microwave background (CMB),
as both science cases require sensitive polarimetric observations of low-
column-density regions of sky. In addition, a better understanding of the
magnetized ISM allows better removal of the polarized foreground for next-
generation CMB polarization experiments. This is critical for the current
generation of ground and balloon-borne CMB experiments (Henderson et al.,
2016; EBEX Collaboration et al., 2018; Gandilo et al., 2016; Gualtieri et al.,
2018; The BICEP/Keck Collaboration et al., 2018; Galitzki et al., 2018;
Harrington et al., 2016), and will be even more important for future
experiments that will have increased sensitivity (Abitbol et al., 2017; Hanany
et al., 2019; Suzuki et al., 2018).

##  3 The role of magnetic fields in molecular clouds

Current observations suggest that the outer envelopes of molecular clouds can
be supported against gravity by magnetic fields at $\sim$ parsec scales, but
that gravity dominates at the $\lesssim$ 0.1 pc scale of dense cores, and so
these dense structures can collapse to form stars (Crutcher et al., 2010).
Furthermore, Planck observations of the ten nearest (mostly low-mass)
molecular clouds show a change in alignment of cloud structure from parallel
to perpendicular to the magnetic field with increasing column density (Planck
Collaboration Int. XXXII, 2016; Planck Collaboration Int. XXXV, 2016); this
effect can only be reproduced in simulations if the cloud-scale magnetic
energy density is equal to or larger than the energy density of turbulence
(Soler et al., 2013; Chen et al., 2016; Mocz & Burkhart, 2018). More
observations are needed to determine whether this “strong field” case is true
for giant molecular clouds, to measure the physical scale at which molecular
clouds make the transition from being magnetically dominated to turbulence-
and/or gravity-dominated, and to determine whether the cloud magnetic field
strength affects star formation efficiency within molecular clouds.

Figure 1:  500 $\mu$m image of the Vela C molecular cloud from the balloon-
borne BLASTPol instrument; the image has a resolution of 2.5′, or $\sim$ 500
pc at the distance to Vela C (Fissel et al., 2016). Next-generation space-
based single-dish missions will enable even higher resolution and higher
sensitivity images of nearby molecular clouds.

Polarimeters proposed as part of next-generation space observatories, such as
OST, PICO, and SPICA, will have better sensitivity to low surface brightness
dust and a higher spatial dynamic range and could build on the progress made
by Planck, SOFIA, and BLASTPol in the area of detailed studies of nearby
molecular clouds (such as the BLASTPol observations of Vela C: see Figure 1,
from Fissel et al. 2016). Full-sky coverage polarization maps at a resolution
of 1′, would be especially helpful since all maps would have better than 1 pc
resolution out to a distance of 3.4 kpc, yielding detailed maps of magnetic
fields in over 2,000 molecular clouds with 1000s to 100,000s of independent
measurements per cloud. These observations will complement other methods for
probing magnetic fields in molecular clouds such as Zeeman splitting of
molecular lines (Crutcher et al., 2010), analysis of velocity structure
(Lazarian & Yuen, 2018), and optical/near-IR polarization due to selective
extinction of background starlight (Marchwinski et al., 2012).

Observing a large sample of clouds is crucial because measurements of
polarization angle are sensitive to only the magnetic field orientation
projected on the plane of the sky. Polarization maps will look very different
for molecular clouds for different viewing angles; for young clouds versus
evolved clouds; and for low-mass clouds compared with high-mass star forming
regions. With future balloon and space polarimeters providing detailed maps of
magnetic fields in hundreds of molecular clouds it will be possible to compare
cloud magnetization to the efficiency of star formation, measured from near
and far-IR observations of dense cores and protostars with Herschel, Spitzer,
and JWST. With observations of thousands of molecular clouds it should be
possible to determine the role of magnetic fields in star formation as a
function of age and mass.

As a first step toward this goal, balloon-borne sub-millimeter polarimeters
such as the upcoming BLAST-TNG telescope Galitzki et al. (2014) and the PILOT
experiment (Bernard et al., 2016), which operate above $>$99.5$\%$ of the
atmosphere, will be able to map dozens of molecular clouds with $<$ 1′
resolution. However, polarization capabilities on next-generation far-IR
satellites, which can be cooled to $<$10 K and therefore will have orders of
magnitude improvements in sensitivity, will enable much more detailed studies
of the nearest molecular clouds. With 10′′ FWHM resolution (15$\times$ higher
than the BLASTPol map of Vela C shown in Figure 1), such a telescope could
measure the cloud magnetization and MHD turbulence over a wide range of
density scales, from the diffuse cloud envelope, to the scale of dense, $<$0.1
pc width filaments and cores.

##  4 Do magnetic fields affect the collapse of cores and filaments?

Filamentary structures within molecular clouds appear to contain most of the
dense, gravitationally unstable gas and are observed to be the preferred sites
of star formation (André et al., 2014). However, we still do not know if the
magnetic field plays an important role in filament formation and evolution, as
Planck could not spatially resolve filaments, which often have widths
$\sim$0.1 pc (Arzoumanian et al., 2019). Polarization observations of dense
filaments and cores are particularly challenging as dust grains embedded in
filaments are expected to be less well aligned, and thus polarization levels
are expected to be lower (Andersson et al., 2015).

Multi-scale comparisons of the orientation between column density structures
and the magnetic field, indicate that core scale fields may have magnetic
energy density equal to or less than that of turbulence (Hull et al., 2017,
see Fig. 2 (left)). In contrast, in similar studies of cloud scale fields the
magnetic energy appears to be at least as strong as turbulence (Planck
Collaboration Int. XXXV, 2016; Soler et al., 2017). Filaments and cores could
therefore represent the densities at which magnetic energy becomes subdominant
to gas kinetic energy (Chen et al., 2016).

Systematic surveys of magnetic fields in filaments and cores are beginning,
with instruments such as the TolTEC instrument on the LMT ($\lambda_{central}$
= 1.1, 1.4, and 2.1mm; best FWHM resolution of 5′′;Ward-Thompson et al. 2017;
Pattle et al. 2018), POL-2/SCUBA2 on the JCMT ($\lambda_{central}$ = 450 and
850 $\mu m$, with FWHM resolution of 14′′ at 850 $\mu m$), and the HAWC+
polarimeter on SOFIA (50-214 $\mu$m, best FWHM resolution of 6′′; Harper et
al., 2018), all utilizing large focal plane cameras. These instruments provide
the ability to resolve filaments in nearby clouds, but are typically limited
by sensitivity to observe only high column density filaments ($A_{V}\,>10$
mag). A complete survey of fields within all filaments across many clouds,
with resolution to isolate sub-filamentary structures seen in many nearby
clouds (e.g. Hacar et al. 2018) is needed. Such a survey would probe whether
magnetic fields provide significant support against gravitational collapse,
and elucidate the role magnetic fields play in filament and core formation.
With increased sensitivity, statistical quantities such as the histogram of
relative orientations (HRO; Soler et al., 2013) can be utilized to compare
models of clouds, filaments, and cores with data (See Fig. 2 (right)).

Weak magnetic fields can still strongly inhibit protostellar disk formation if
the field is parallel to the angular momentum axis of the disk, a conundrum
known as the magnetic braking catastrophe (Allen et al., 2003; Wurster & Li,
2018). The extent to which magnetic fields limit disk formation, perhaps due
to misalignment between the angular momentum and magnetic field (Joos et al.,
2012), is an important outstanding question. Large surveys of magnetic fields
on filament and core scales with single dish polarimeters, combined
observations of disks properties from ALMA and the EVLA, will be able to
determine whether there is a correlation between disk sizes and masses and
core/cloud magnetization and field geometry. Polarimetry using interferometers
such as ALMA or potentially the ngVLA can trace magnetic fields at even higher
resolution, down to the scales of of protostellar envelopes and protoplanetary
disks (Hull et al., 2018; Hull & Zhang, 2019).

Figure 2: Figures from Hull et al. (2017). (left) A multi-scale polarimetric
view towards the Ser-emb 8 protostar. (right) With increased statistics from
increasingly sensitive polarimeters, statistics such as the histogram of
relative orientations (HRO) can be used to compare models to measurements and
in turn extract physical parameters from polarimetric observations.

##  5 Conclusions

To summarize, progress in understanding the magnetized diffuse ISM and the
role of magnetic fields in formation of stars in molecular clouds requires:

  * •

High sensitivity polarization maps of the entire sky. All-sky surveys of dust
polarization will allow detailed studies of magnetic fields and magnetized
turbulence in the diffuse interstellar medium. Since most of the magnetic
field sampled with polarized dust emission above the Galactic plane is located
within $\sim$200 pc of the sun, high sensitivity surveys with $<$ 1.5′
resolution will resolve scales below 0.1 pc. In addition, all-sky polarization
maps from proposed experiments such as PICO will produce detailed magnetic
field maps of thousands of molecular clouds.

  * •

Detailed polarization maps with $\leq$10′′ resolution over a large spatial
dynamic range. To study the interaction between turbulence, magnetic fields,
feedback, and gravity within molecular clouds, we require detailed
polarization maps that can both resolve magnetic fields on $<<$0.1 pc scales
of cores and filaments, while also having the sensitivity to produce high-
fidelity maps of the emission from low-column-density molecular cloud
envelopes. Balloon-borne polarimeters will soon produce sub-arcmin resolution
maps of many molecular clouds; however, higher sensitivity and resolution maps
are needed from next-generation satellites such as OST (Staguhn et al., 2018)
and SPICA (Rodriguez et al., 2018).

  * •

Next-generation single-dish cameras with high resolution. Polarimeters on
large single-dish ground based telescopes, such as TolTEC on LMT (Austermann
et al., 2018), POL-2 on JCMT (Friberg et al., 2018), NIKA2pol on IRAM (Adam et
al., 2018), and future instruments on proposed facilities like AtLAST
(Mroczkowski & Noroozian, 2018), are needed to map magnetic field in dense
structures, such as in cores and filaments.

  * •

Continued development of polarization capabilities for interferometers.
Improved polarization sensitivity and resolution for ALMA and new facilities
like the ngVLA (Murphy et al., 2018), are need to survey magnetic fields in
protostellar and protoplanetary disks.

## References

  * Abitbol et al. (2017) Abitbol, M.H., Chluba, J., et al., 2017.  Monthly Notices of the Royal Astronomical Society, 471:1126
  * Adam et al. (2018) Adam, R., Adane, A., et al., 2018.  Astronomy & Astrophysics, 609:A115
  * Allen et al. (2003) Allen, A., Li, Z.Y., et al., 2003.  Astrophysical Journal, 599:363
  * Andersson et al. (2015) Andersson, B.G., Lazarian, A., et al., 2015.  Annual Review of Astronomy and Astrophysics, 53:501
  * André et al. (2014) André, P., Di Francesco, J., et al., 2014.  Protostars and Planets VI, 27–51
  * Arzoumanian et al. (2019) Arzoumanian, D., André, P., et al., 2019.  Astronomy & Astrophysics, 621:A42
  * Austermann et al. (2018) Austermann, J.E., Beall, J.A., et al., 2018.  Journal of Low Temperature Physics, 193:120
  * Bernard et al. (2016) Bernard, J.P., Ade, P., et al., 2016.  Experimental Astronomy, 42:199
  * Bracco et al. (2019) Bracco, A., Candelaresi, S., et al., 2019.  Astronomy & Astrophysics, 621:A97
  * Caldwell et al. (2017) Caldwell, R.R., Hirata, C., et al., 2017.  Astrophysical Journal, 839:91
  * Chen et al. (2016) Chen, C.Y., King, P.K., et al., 2016.  Astrophysical Journal, 829:84
  * Clark et al. (2015) Clark, S.E., Hill, J.C., et al., 2015.  Physical Review Letters, 115(24):241302
  * Clark et al. (2014) Clark, S.E., Peek, J.E.G., et al., 2014.  Astrophysical Journal, 789:82
  * Crutcher (2012) Crutcher, R.M., 2012.  Annual Review of Astronomy and Astrophysics, 50:29
  * Crutcher et al. (2010) Crutcher, R.M., Wandelt, B., et al., 2010.  Astrophysical Journal, 725:466
  * Dotson et al. (2010) Dotson, J.L., Vaillancourt, J.E., et al., 2010.  Astrophysical Journal Supplement Series, 186:406
  * Draine (2011) Draine, B.T., 2011.  _Physics of the Interstellar and Intergalactic Medium_
  * EBEX Collaboration et al. (2018) EBEX Collaboration, Abitbol, M., et al., 2018.  The Astrophysical Journal Supplement Series, 239:8
  * Elmegreen & Scalo (2004) Elmegreen, B.G. & Scalo, J., 2004.  Annual Review of Astronomy and Astrophysics, 42:211
  * Federrath (2015) Federrath, C., 2015.  Monthly Notices of the Royal Astronomical Society, 450:4035
  * Fissel et al. (2016) Fissel, L.M., Ade, P.A.R., et al., 2016.  Astrophysical Journal, 824:134
  * Friberg et al. (2018) Friberg, P., Berry, D., et al., 2018.  In _Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series_ , vol. 10708, 107083M
  * Galitzki et al. (2014) Galitzki, N., Ade, P.A.R., et al., 2014.  Journal of Astronomical Instrumentation, 3:1440001
  * Galitzki et al. (2018) Galitzki, N., Ali, A., et al., 2018.  In _Millimeter, Submillimeter, and Far-Infrared Detectors and Instrumentation for Astronomy IX_ , vol. 10708 of _Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series_ , 1070804
  * Gandilo et al. (2016) Gandilo, N.N., Ade, P.A.R., et al., 2016.  In _Millimeter, Submillimeter, and Far-Infrared Detectors and Instrumentation for Astronomy VIII_ , vol. 9914 of _Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series_ , 99141J
  * Gualtieri et al. (2018) Gualtieri, R., Filippini, J.P., et al., 2018.  Journal of Low Temperature Physics, 193:1112
  * Hacar et al. (2018) Hacar, A., Tafalla, M., et al., 2018.  Astronomy & Astrophysics, 610:A77
  * Hanany et al. (2019) Hanany, S., Alvarez, M., et al., 2019.  arXiv e-prints, arXiv:1902.10541
  * Harper et al. (2018) Harper, D.A., Runyan, M.C., et al., 2018.  JAI, 7(2):1840008
  * Harrington et al. (2016) Harrington, K., Marriage, T., et al., 2016.  In _Millimeter, Submillimeter, and Far-Infrared Detectors and Instrumentation for Astronomy VIII_ , vol. 9914 of _Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series_ , 99141K
  * Henderson et al. (2016) Henderson, S.W., Allison, R., et al., 2016.  Journal of Low Temperature Physics, 184:772
  * Hull et al. (2018) Hull, C.L.H., Carrasco-González, C., et al., 2018.  In E. Murphy, ed., _Science with a Next Generation Very Large Array_ , vol. 517 of _Astronomical Society of the Pacific Conference Series_ , 357
  * Hull et al. (2017) Hull, C.L.H., Mocz, P., et al., 2017.  Astrophysical Journal Letters, 842:L9
  * Hull & Zhang (2019) Hull, C.L.H. & Zhang, Q., 2019.  Frontiers in Astronomy and Space Sciences, 6:3
  * Jones et al. (1984) Jones, T.J., Hyland, A.R., et al., 1984.  Astrophysical Journal, 282:675
  * Joos et al. (2012) Joos, M., Hennebelle, P., et al., 2012.  Astronomy & Astrophysics, 543:A128
  * Kalberla et al. (2016) Kalberla, P.M.W., Kerp, J., et al., 2016.  Astrophysical Journal, 821:117
  * Kandel et al. (2017) Kandel, D., Lazarian, A., et al., 2017.  Monthly Notices of the Royal Astronomical Society, 472:L10
  * Lazarian & Hoang (2007) Lazarian, A. & Hoang, T., 2007.  Monthly Notices of the Royal Astronomical Society, 378:910
  * Lazarian & Yuen (2018) Lazarian, A. & Yuen, K.H., 2018.  Astrophysical Journal, 853:96
  * Malinen et al. (2016) Malinen, J., Montier, L., et al., 2016.  Monthly Notices of the Royal Astronomical Society, 460:1934
  * Marchwinski et al. (2012) Marchwinski, R.C., Pavel, M.D., et al., 2012.  Astrophysical Journal, 755:130
  * Matthews et al. (2009) Matthews, B.C., McPhee, C.A., et al., 2009.  Astrophysical Journal Supplement Series, 182:143
  * McKee & Ostriker (2007) McKee, C.F. & Ostriker, E.C., 2007.  Annual Review of Astronomy and Astrophysics, 45:565
  * Mocz & Burkhart (2018) Mocz, P. & Burkhart, B., 2018.  Monthly Notices of the Royal Astronomical Society, 480:3916
  * Mroczkowski & Noroozian (2018) Mroczkowski, T. & Noroozian, O., 2018.  In _Atacama Large-Aperture Submm/mm Telescope (AtLAST)_ , 26
  * Murphy et al. (2018) Murphy, E.J., Bolatto, A., et al., 2018.  In E. Murphy, ed., _Science with a Next Generation Very Large Array_ , vol. 517 of _Astronomical Society of the Pacific Conference Series_ , 3
  * Pattle et al. (2018) Pattle, K., Ward-Thompson, D., et al., 2018.  Astrophysical Journal Letters, 860:L6
  * Planck Collaboration 2018 XI (2018) Planck Collaboration 2018 XI, 2018.  ArXiv e-prints
  * Planck Collaboration Int. XIX et al. (2015) Planck Collaboration Int. XIX, 2015.  Astronomy & Astrophysics, 576:A104
  * Planck Collaboration Int. XXXII (2016) Planck Collaboration Int. XXXII, 2016.  Astronomy & Astrophysics, 586:A135
  * Planck Collaboration Int. XXXV (2016) Planck Collaboration Int. XXXV, 2016.  Astronomy & Astrophysics, 586:A138
  * Planck Collaboration Int. XXXVIII (2016) Planck Collaboration Int. XXXVIII, 2016.  Astronomy & Astrophysics, 586:A141
  * Rodriguez et al. (2018) Rodriguez, L., Poglitsch, A., et al., 2018.  Journal of Low Temperature Physics, 193(3-4):449
  * Soler et al. (2017) Soler, J.D., Ade, P.A.R., et al., 2017.  Astronomy & Astrophysics, 603:A64
  * Soler et al. (2013) Soler, J.D., Hennebelle, P., et al., 2013.  Astrophysical Journal, 774:128
  * Staguhn et al. (2018) Staguhn, J., Amatucci, E., et al., 2018.  In _Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series_ , vol. 10698, 106981A
  * Strong et al. (2007) Strong, A.W., Moskalenko, I.V., et al., 2007.  Annual Review of Nuclear and Particle Science, 57:285
  * Suzuki et al. (2018) Suzuki, A., Ade, P.A.R., et al., 2018.  Journal of Low Temperature Physics, 193:1048
  * The BICEP/Keck Collaboration et al. (2018) The BICEP/Keck Collaboration, :, et al., 2018.  arXiv e-prints, arXiv:1807.02199
  * Ward-Thompson et al. (2017) Ward-Thompson, D., Pattle, K., et al., 2017.  Astrophysical Journal, 842:66
  * Wurster & Li (2018) Wurster, J. & Li, Z.Y., 2018.  Frontiers in Astronomy and Space Sciences, 5:39


# Priors For Symbolic Regression

Deaglan J. Bartlett  CNRS & Sorbonne Université, Institut d’Astrophysique de
Paris (IAP),UMR 7095, 98 bis bd Arago,ParisF-75014France
<EMAIL_ADDRESS>,  Harry Desmond  Institute of Cosmology &
Gravitation, University of Portsmouth,Dennis Sciama BuildingPortsmouthPO1
3FXUK and  Pedro G. Ferreira  Astrophysics, University of Oxford,Denys
Wilkinson Building, Keble RoadOxfordOX1 3RHUK

(2023)

###### Abstract.

When choosing between competing symbolic models for a data set, a human will
naturally prefer the “simpler” expression or the one which more closely
resembles equations previously seen in a similar context. This suggests a non-
uniform prior on functions, which is, however, rarely considered within a
symbolic regression (SR) framework. In this paper we develop methods to
incorporate detailed prior information on both functions and their parameters
into SR. Our prior on the structure of a function is based on a $n$-gram
language model, which is sensitive to the arrangement of operators relative to
one another in addition to the frequency of occurrence of each operator. We
also develop a formalism based on the Fractional Bayes Factor to treat
numerical parameter priors in such a way that models may be fairly compared
though the Bayesian evidence, and explicitly compare Bayesian, Minimum
Description Length and heuristic methods for model selection. We demonstrate
the performance of our priors relative to literature standards on benchmarks
and a real-world dataset from the field of cosmology.

model selection, minimum description length, symbolic regression, equation
learning, data analysis, cosmology, language model

††journalyear: 2023††copyright: acmlicensed††conference: Genetic and
Evolutionary Computation Conference Companion; July 15–19, 2023; Lisbon,
Portugal††booktitle: Genetic and Evolutionary Computation Conference Companion
(GECCO ’23 Companion), July 15–19, 2023, Lisbon, Portugal††price: 15.00††doi:
10.1145/3583133.3596327††isbn: 979-8-4007-0120-7/23/07††ccs: Mathematics of
computing Probabilistic inference problems††ccs: Computing methodologies
Machine learning

##  1\. Introduction

Symbolic regression (SR) is traditionally viewed as a Pareto optimisation
problem, with the aim of finding the symbolic expressions which cannot be made
more accurate without becoming more complex. The task of the SR algorithm –
whether a genetic algorithm (Turing, 1950; David, 1989; Haupt and Haupt, 2004;
Lemos et al., 2022), a deterministic search (Worm and Chiu, 2013; Kammerer et
al., 2021; Rivero and Fernandez-Blanco, 2019; McConaghy, 2011), reinforcement
learning (Sutton and Barto, 2018), a physics-inspired search (Udrescu and
Tegmark, 2020; Udrescu et al., 2020) or an exhaustive scan of functional
parameter space (Bartlett et al., 2022) – is to find candidate functions, fit
them to the data by minimising a loss function, and identify those that
Pareto-dominate the rest. While the loss is typically straightforward to
quantify as the negative of the likelihood, complexity is an inherently
ambiguous notion. A variety of statistics have been constructed to quantify
this (see e.g. (Smits and Kotanchek, 2005; Kommenda et al., 2021; Chen et al.,
2022; Bomarito et al., 2022)), which may be a function of the number and/or
type of operators and parameters (e.g. (Udrescu and Tegmark, 2020; Udrescu et
al., 2020; Cranmer et al., 2020; Cranmer, 2020)) or the degree of nonlinearity
across the domain of interest (Vladislavleva et al., 2009).

Rather than using two separate objectives, recent work (Bartlett et al., 2022;
Guimerà et al., 2020) has focused on combining accuracy and simplicity into a
single goodness-of-fit statistic. There are fundamentally two approaches,
which turn out to be closely related. Utilising the minimum description length
(MDL) principle (Rissanen, 1978; Grünwald and Roos, 2019; Grunwald, 2007;
Lanterman, 2001), Bartlett et al. (2022) derived the codelength (number of
nats) required to send the data with the help of a given function. By putting
both on a common scale of information content, this naturally combines the
preference for accuracy (lower residuals around the function’s expectation)
and simplicity (less information required to convey the structure of the
function and its parameters). Models are penalised which contain many terms, a
large number of operators or a fine-tuning of parameters to fit the data well.
Guimerà et al. (2020) instead considered a Bayesian approach where the
goodness of fit metric is the evidence. In this case the preference for
functional simplicity is encoded in their prior. The choice of Guimerà et al.
is to assign this prior based on a corpus of equations in a given domain such
that functions more similar to those in the corpus are upweighted (said to be
simpler). Each operator is assigned a probability based on its frequency in
the corpus and the operators are then assumed to be independent when
calculating the probability of the full function.

Our goal is to compare the Bayesian and MDL methods and propose two upgrades
to the Bayesian approach, one targeting the parameter prior and one the
operator prior. Our starting point in each case is the method of Guimerà et
al. (2020). They model the accuracy and parameter parts of the evidence using
the Bayesian Information Criterion, an asymptotic Gaussian expansion of the
likelihood that holds in the limit that the number of data points greatly
exceeds the number of free parameters. This is however not a robust
approximation to the evidence without a well-motivated prior: the posterior
and hence evidence scale with the width of the prior (e.g. the range of an
“uninformative” uniform), affecting in general the Bayes factor between two
functions. Their prior on the structure of functions treats each node of the
tree independently, which is in general not desirable behaviour. For example,
consider the two functions $y=\sin\left(\sin\left(x_{0}+x_{1}\right)\right)$
and $y=\sin(x_{0})+\sin(x_{1})$. Both of these functions can be expressed as
trees with the same number of nodes and the same number of each type of
operator, and thus a prior that depends on this would assign equal
probabilities to these functions. However, a physicist would find a sum of
sines more likely a priori than a nested sine, suggesting that a correlation
between a node and its children should be incorporated into the operator
prior. We propose solutions to both problems.

The structure of the paper is as follows. We review the Bayesian and MDL model
selection methods in Sections 2.1 and 2.2 before discussing the choice of
parameter prior in Section 2.3. To incorporate the correlation between
operators in a prior on functions, we introduce a method for using a language
model as a structural prior in Section 2.4 and implement it using a corpus of
well-known scientific equations. We benchmark the models that we consider in
Section 3 before applying them to a real-world example (the expansion history
of the universe) in Section 4. We show that our new methods can lead to
preference for more physically realistic equations than those in Bartlett et
al. (2022). We discuss the context of our work in Section 5, and conclude in
Section 6.

##  2\. Model selection and priors

###  2.1. Bayesian model selection

When performing Bayesian model selection, one begins with a set of candidate
models, $\\{f_{0},f_{1},\ldots\\}$, and their parameters
$\\{\bm{\theta_{0}},\bm{\theta_{1}},\ldots\\}$, which aim to describe a
dataset, $D$. Using Bayes’ theorem, one can compute the probability of one of
these models, $f_{i}$, given $D$ to be

(1) $\begin{split}P\left(f_{i}|D\right)&=\frac{1}{P(D)}\int
P\left(D|f_{i},\bm{\theta}_{i}\right)P\left(\bm{\theta}_{i}|f_{i}\right)P\left(f_{i}\right)\mathrm{d}\bm{\theta}_{i}\\\
&\equiv\frac{P\left(f_{i}\right)}{P(D)}\mathcal{Z}\left(D|f_{i}\right),\end{split}$

where $P(D)$ is the probability of the data (under any model),
$\mathcal{Z}\left(D|f_{i}\right)$ is the evidence,
$P\left(D|f_{i},\bm{\theta}_{i}\right)$ is the likelihood,
$P\left(\bm{\theta}_{i}|f_{i}\right)$ is the prior on the parameters, and
$P\left(f_{i}\right)$ is the prior on the model itself. When comparing models
on the same data one can ignore the constant $P(D)$, and thus one wishes to
find the model that minimises

(2) $-\log P\left(f_{i}|D\right)=-\log
P\left(f_{i}\right)-\log\mathcal{Z}\left(D|f_{i}\right).$

The evidence is often approximated in the case of a uni-modal posterior
distribution by applying Laplace’s approximation about the maximum a
posteriori point,
$\hat{\bm{\theta}}_{i}^{\mathcal{H}}\equiv\operatorname*{arg\,max}_{\bm{\theta}_{i}}\mathcal{H}(D,f_{i},\bm{\theta}_{i})$
for $\mathcal{H}(D,f_{i},\bm{\theta}_{i})\equiv
P\left(D|f_{i},\bm{\theta}_{i}\right)P\left(\bm{\theta}_{i}|f_{i}\right)$, to
obtain

(3)
$\log\mathcal{Z}\left(D|f_{i}\right)\simeq\log\mathcal{\hat{H}}+\frac{p}{2}\log
2\pi-\frac{1}{2}\log\det\hat{\bm{I}}^{\mathcal{H}},$

where
$\mathcal{\hat{H}}=\mathcal{H}(D,f_{i},\hat{\bm{\theta}}_{i}^{\mathcal{H}})$
is the maximum of the posterior,
$\hat{I}^{\mathcal{H}}_{\alpha\beta}=-\partial_{\alpha}\partial_{\beta}\log\mathcal{H}(D,f_{i},\bm{\theta}_{i})|_{\hat{\bm{\theta}_{i}}^{\mathcal{H}}}$
and $p$ is the number of parameters of the model. In the limit of a large
number of data points, $N$, the final term can be approximated by $-p/2\log N$
and the term proportional to $p$ is neglected to obtain (minus one half of)
the Bayesian Information Criterion (BIC) (Schwarz, 1978). Since the full
evidence is typically expensive to compute, we will use the approximation
given by Eq. 3 throughout the paper (but not the simpler BIC approximation).
In all numerical experiments, $\hat{\bm{I}}^{\mathcal{H}}$ is computed via
finite differences with the numdifftools package (Brodtkorb and D’Errico,
2015).

###  2.2. Model selection by minimum description length

An alternative approach to model selection is the minimum description length
(MDL) principle. This was first applied to Genetic Programming by Iba et al.
(1994) and to SR by Udrescu and Tegmark (2020), and refined and formalised by
Bartlett et al. (2022). Under a two-part coding scheme, the amount of
information (measured in nats, since we use natural logarithms) required to
send a dataset $D$ given a hypothesis (trial function) $H$ is

(4) $L(D)=L(H)+L(D|H),$

where $L(H)$ is the number of nats required to send the hypothesis, and
$L(D|H)$ the number needed to send the data given this hypothesis. The optimal
Shannonn–Fano coding scheme gives $L(D|H)=-\log
P\left(D|f_{i},\hat{\bm{\theta}}_{i}^{\mathcal{L}}\right)$ (Cover and Thomas,
1991) for hypothesis $H=f_{i}$, where
$\hat{\bm{\theta}}_{i}^{\mathcal{L}}\equiv\operatorname*{arg\,max}_{\bm{\theta}_{i}}\mathcal{L}(D,f_{i},\bm{\theta}_{i})$
is the maximum likelihood point. One then chooses an encoding for the function
and its parameters, $L(H)$, and selects the model that minimises $L(D)$.

Bartlett et al. choose this encoding such that the description length of a
function comprised of $n$ operators, $p$ parameters $\bm{\theta}$, with
constants $\\{c_{\alpha}\\}$ and expressed as a tree with $k$ nodes is

$\displaystyle L(D)=$ $\displaystyle-\log
P\left(D|f_{i},\hat{\bm{\theta}}_{i}^{\mathcal{L}}\right)+k\log(n)-\frac{p}{2}\log(3)$
$\displaystyle+\sum_{\alpha}\log(c_{\alpha})+\sum_{\alpha}^{p}\left(\frac{1}{2}\log\left(\hat{I}^{\mathcal{L}}_{\alpha\alpha}\right)+\log(|\hat{\theta}_{\alpha}^{\mathcal{L}}|)\right).$

Note that negative constants are encoded via a sign operator that is counted
in $k$ and $n$. The term
$\hat{I}^{\mathcal{L}}_{\alpha\beta}=-\partial_{\alpha}\partial_{\beta}\log
P\left(D|f_{i},\bm{\theta}_{i}\right)|_{\hat{\bm{\theta}_{i}}^{\mathcal{L}}}$
is derived by transmitting a discretised
$\hat{\bm{\theta}}_{\alpha}^{\mathcal{L}}$ in the code, and considering the
optimal discretisation for each of the parameters separately.

An obvious extension to Eq. 2.2 is to discretise the parameters jointly rather
than separately in each dimension, capturing the degeneracy between them in
fitting the data. Using a rectangular lattice with a freely varying
orientation in parameter space replaces
$\sum_{\alpha}^{p}\left(\frac{1}{2}\log\left(\hat{I}^{\mathcal{L}}_{\alpha\alpha}\right)\right)-\frac{p}{2}\log(3)$
by $\frac{1}{2}\log\det\hat{\bm{I}}^{\mathcal{L}}+\frac{p}{2}\log(\nu_{p})$,
where $\log(\nu_{p})=1-\log(3)\,\forall p$ (Takeuchi, 1997). A further
reduction in description length is possible by instead using an optimal
quantising lattice, which gives $\log(\nu_{p})=1+\log(4\kappa_{p})$ (Wallace
and Freeman, 1987, 1992; Lanterman, 2001). Here $\kappa_{p}$ is a constant
relating to the geometry of the optimal lattice in $p$ dimensions, with the
limits $\kappa_{1}=1/12$ and $\kappa_{p}\to(2\pi e)^{-1}$ for $p\gg 1$ (Conway
et al., 2013). For simplicity we will use a rectangular lattice hereafter.

A more fundamental extension is to use prior information on the parameters to
encode them optimally by assigning shorter codelengths to more a priori
probable values. In this case the optimal parameter discretisation is
determined with respect to $\hat{\bm{I}}^{\mathcal{H}}$ rather than
$\hat{\bm{I}}^{\mathcal{L}}$. Combining $\log
P\left(D|f_{i},\bm{\theta}_{i}\right)$ with the parameter encoding $\log
P(\bm{\theta}_{i}|f_{i})$, removing the
$\log(|\hat{\theta}_{\alpha}^{\mathcal{L}}|)$ that derives from integer
encoding and using the rectangular discretisation, the total description
length becomes

(6)
$\tilde{L}(D)=-\log\mathcal{\hat{H}}+\frac{1}{2}\log\det\bm{I}^{\mathcal{H}}+\frac{p}{2}\log(\nu_{p})+k\log(n)+\sum_{\alpha}\log
c_{\alpha}.$

By comparing to Eq. 2 and using the approximation Eq. 3, one can then identify
the prior that is implicitly imposed on the functions by MDL:

(7) $-\log P\left(f_{i}\right)=k\log n+\sum_{\alpha}\log
c_{\alpha}+\frac{p}{2}\log\left(2\pi\nu_{p}\right).$

We see that the first two terms in this expression are also included in
Section 2.2; they penalise equations with many nodes, many different operators
or many or finely tuned constants. These appear as logarithms by construction,
since they are encoded using Risannen’s universal prefix code for the integers
(Rissanen, 1983), but we have neglected the
$\mathcal{O}\left(\log\log(\cdot)\right)$ term that is formally required to
send all nested prefixes. This also explains the appearance of the
$\log|\hat{\theta}_{\alpha}|$ term in Section 2.2: without prior information
on the parameters we have encoded them as integers given the precision
$\Delta_{\alpha}$.

The final term in the functional prior prefers functions with fewer parameters
and is merely the number of these, $p$, multiplied by an $\mathcal{O}(1)$
coefficient ($\sim$0.87 for rectangular quantisation or 0.69-0.87 for optimal
quantisation). These are desirable properties of a prior which prefers
“simple” functions, making it a good choice of SR loss function in the absence
of any prior preference for the type and arrangement of operators in a
function. Each of the model selection metrics discussed here prefer simpler
functions and hence limit or prevent overfitting, although the extent to which
they do this relative to the use of a test–train split deserves further
investigation.

###  2.3. Priors on parameters

To compute Eq. 2, one must choose a prior on the parameters of the function
considered. The typical Bayesian view is that these encode subjective prior
knowledge, but this is problematic in a SR context where one knows nothing a
priori about the thousands or millions of functions one is considering. One
therefore needs an automated method for constructing the prior.

Two common choices are uniform priors on all parameters or a Jeffreys prior
(Jeffreys, 1946) to afford reparameteristaion invariance. However, in the
absence of bounds on the parameter ranges a uniform prior is improper
(integrates to infinity if it is non-zero) and this is often the case for
Jeffreys priors too. This makes Eq. 2 infinite and hence unusable.
Alternatively, one could introduce bounds on the parameters but this simply
scales $\mathcal{Z}(D|f_{i})$ and hence Eq. 2 by an amount related to the
range. Even if one could choose “reasonable” bounds on all parameters, this
rescaling will be different for functions with different numbers of parameters
and hence this method cannot provide a self-consistent model comparison
metric. The typical derivation of the BIC ignores this constant by expanding
the likelihood rather than the posterior, but this appears not to be
acceptable in general. One can alter the loss function of Guimerà et al.
(2020) and hence the function ranking by modifying the uniform prior bounds
without any change to the posteriors. Alternative candidate priors include the
reference prior (Bernardo, 1979) and the maximum-entropy prior. However, the
reference prior, which maximises the expected information gain due to the
data, reduces to the Jeffreys prior for asymptotically normal likelihoods, and
hence is unsuitable as per the above. The maximum-entropy prior may include a
normalisability constraint, but must include other constraints to make this
feasible as the maximum-entropy distribution without any further information
is an (infinite) uniform.

A solution to this problem is afforded by the Partial Bayes Factor (Lempers,
1971). Here, one splits the observed data into two samples: one for training
and another for model comparison. One first runs an inference on the training
sample with a given choice of prior (which does not need to be proper as the
evidence is not required here) to obtain a posterior distribution of the
parameters, which will be normalisable. This posterior distribution is then
used as the prior for the second inference using the remaining data, from
which the evidence is computed and models are compared. Of course, one must
choose a size of the training set and decide which data points go into it.
Berger and Pericchi (1996) suggest, for a chosen training set size $M$, using
all possible training samples and averaging the result to obtain the Intrinsic
Bayes Factor, although which choice of average in unclear and this will
require one to run $N!/(M!(N-M)!)$ inferences if one has $N$ data points.

To alleviate these problems, O’Hagan (1995) introduced the Fractional Bayes
Factor (FBF). Noting that if $M$ and $N$ are large, so that the likelihood of
the test set is approximately the likelihood for the full dataset to the power
$b=M/N$, O’Hagan suggests comparing models based on

(8) $B_{b}=\frac{q_{1}(b)}{q_{2}(b)},\quad q_{i}(b)=\frac{\int
P\left(D|f_{i},\bm{\theta}_{i}\right)P\left(\bm{\theta}_{i}|f_{i}\right)\mathrm{d}\bm{\theta}_{i}}{\int
P\left(D|f_{i},\bm{\theta}_{i}\right)^{b}P\left(\bm{\theta}_{i}|f_{i}\right)\mathrm{d}\bm{\theta}_{i}},$

which is well-defined for an improper prior
$P\left(\bm{\theta}_{i}|f_{i}\right)$. This simply requires a choice of $b$,
or test–train split fraction. The basic requirements are that $b$ must be
between $M_{0}/N$ and 1 (where $M_{0}$ is the minimum permissible training set
size) and $b\to 0$ as $N\to\infty$. For robustness of training, however, one
prefers a larger $b$, where ideally $Nb\to\infty$ as $N\to\infty$. From these
considerations, O’Hagan suggests using $b=N^{-1}\log N$ or $b=N^{-1/2}$, where
the latter is found to be more robust. This is how models are compared in a SR
context by (Bomarito et al., 2023), where $q_{i}$ is computed using a
Sequential Monte Carlo algorithm. One could also vary $b$ to find the
corresponding uncertainty on $q_{i}$, but for simplicity we too choose
$b=N^{-1/2}$.

To apply the Fractional Bayes Factor to the Bayesian or MDL model selection
procedure, we note that Eq. 8 implies that one can take all previous formulae
and replace $P(\bm{\theta}_{i}|f_{i})$ by $P_{b}(\bm{\theta}_{i}|f_{i})$,
where

(9)
$P_{b}\left(\bm{\theta}_{i}|f_{i}\right)=\frac{P\left(\bm{\theta}_{i}|f_{i}\right)}{\int
P\left(D|f_{i},\bm{\theta}_{i}^{\prime}\right)^{b}P\left(\bm{\theta}_{i}^{\prime}|f_{i}\right)\mathrm{d}\bm{\theta}_{i}^{\prime}}.$

In a SR context where this is applied to a great number of equations, one may
want to approximate Eq. 9 for rapid computation. Under the Laplace
approximation, expanding about
$\tilde{\bm{\theta}}_{i}\equiv\operatorname*{arg\,max}_{\bm{\theta}_{i}}\mathcal{H}_{b}(D,f_{i},\bm{\theta}_{i})$,
where we define $\mathcal{H}_{b}(D,f_{i},\bm{\theta}_{i})\equiv
P\left(D|f_{i},\bm{\theta}_{i}\right)^{b}P\left(\bm{\theta}_{i}|f_{i}\right)$,
Eq. 9 becomes

(10) $\log
P_{b}\left(\bm{\theta}_{i}|f_{i}\right)\simeq\log\frac{P\left(\bm{\theta}_{i}|f_{i}\right)}{\mathcal{H}_{b}(D,f_{i},\tilde{\bm{\theta}}_{i})}-\frac{p}{2}\log
2\pi+\frac{1}{2}\log\det\tilde{\bm{I}}^{\mathcal{H}_{b}},$

where
$\tilde{I}^{\mathcal{H}_{b}}_{\alpha\beta}=-\partial_{\alpha}\partial_{\beta}\log\mathcal{H}_{b}(D,f_{i},\bm{\theta}_{i})|_{\tilde{\bm{\theta}_{i}}}$.
To compute the description length (Eq. 6) or the approximation to the evidence
in Eq. 3, one must in general perform two parameter optimisations: first, one
must find $\tilde{\bm{\theta}}_{i}$, then use this value to compute
$\hat{\bm{\theta}}_{i}$ under the new prior, where
$\hat{\bm{\theta}}_{i}=\operatorname*{arg\,max}_{\bm{\theta}_{i}}P\left(D|f_{i},\bm{\theta}_{i}\right)P_{b}\left(\bm{\theta}_{i}|f_{i}\right)$.

In this framework, wide and uniform priors, $P(\bm{\theta}_{i}|f_{i})={\rm
const}$, are acceptable and even desirable, so we adopt them. In this case,
$\hat{\bm{\theta}}_{i}=\tilde{\bm{\theta}}_{i}$ and these are equal to the
maximum likelihood value, so this is a special case where one only needs to
perform a single optimisation (of the likelihood). One also sees that
$b^{-1}\tilde{I}^{\mathcal{H}_{b}}_{\alpha\beta}=\tilde{I}^{\mathcal{H}}_{\alpha\beta}=\tilde{I}^{\mathcal{L}}_{\alpha\beta}$,
so
$b^{-p}\det\tilde{I}^{\mathcal{H}_{b}}=\det\tilde{I}^{\mathcal{H}}=\det\tilde{I}^{\mathcal{L}}$.
Subsituting these results into Eqs. 10 and 3, one finds for a FBF prior with
uniform $P(\bm{\theta}_{i}|f_{i})$,

(11) $\log\mathcal{Z}\left(D|f_{i}\right)\simeq(1-b)\log
P\left(D|f_{i},\hat{\bm{\theta}}_{i}\right)+\frac{p}{2}\log b.$

Note that for $b=N^{-1/2}$ this is almost proportional to the BIC, but with an
extra factor of $1/2$ multiplying the $p\log N$ term and an extra term
$N^{-1/2}\log P\left(D|f_{i},\hat{\bm{\theta}}_{i}\right)$, which becomes
unimportant for large $N$. It is interesting that when using a uniform
training-set prior any dependence on $I^{\mathcal{L}}$ has cancelled upon
application of the FBF, independent of the choice of $b$. This may be
undesirable from a model selection perspective compared to the prior on the
integers used by Bartlett et al. (2022), since it does not penalise finely
tuned parameters.

###  2.4. Priors on functions

We now turn to construction of an optimal prior on the structure of candidate
functions.

####  2.4.1. Formalism

Suppose we have a function represented as a tree, $\bm{T}$, and define
$\bm{\mathcal{D}}_{i}$ to be the set of nodes of depth $i$, i.e.
$\bm{\mathcal{D}}_{0}$ is the root node, $\bm{\mathcal{D}}_{i}$ contains the
nodes of the children of $\bm{\mathcal{D}}_{i-1}$ and $i=0,\ldots,d-1$. By the
chain rule, the probability of this tree is

(12)
$P\left(\bm{T}\right)=\prod_{i=0}^{d-1}P\left(\bm{\mathcal{D}}_{i}|\bm{\mathcal{D}}_{i-1},\ldots,\bm{\mathcal{D}}_{0}\right),$

where
$P\left(\bm{\mathcal{D}}_{i}|\bm{\mathcal{D}}_{i-1},\ldots,\bm{\mathcal{D}}_{0}\right)$
is the conditional probability of obtaining $\bm{\mathcal{D}}_{i}$ given all
the nodes with smaller depths. Our goal is to learn this probability from a
set of equations in order to produce a function prior that encapsulates
domain-specific knowledge concerning the types of equations one expects to
find. We achieve this by training a language model on the equation set to
assign probabilities to operators based on those that precede it in the tree.

We begin with two simplifying assumptions:

  1. (1)

We approximate Eq. 12 as an $n$-gram, such that the probability of a set of
nodes with given depth is only dependent on nodes which have a depth within
some value, $n$.

  2. (2)

We group all nodes at a given depth into siblings of nodes which share a
parent ($\bm{\mathcal{D}}_{i}=\\{s_{i}^{0},\ldots,s_{i}^{n_{i}-1}\\}$ where
$s_{i}^{j}$ is the the $j$th set of siblings for $\bm{\mathcal{D}}_{i}$). We
assume that the probability of each set of siblings only depends on their
ascendants and that all sets of siblings are independent at a given depth.

Under these approximations, one can write

(13)
$P\left(\bm{T}\right)\approx\prod_{i=0}^{d-1}\prod_{j=0}^{n_{i}-1}P\left(s_{i}^{j}|a_{i-1}^{ij},\ldots
a_{i-(n-1)}^{ij}\right),$

where $a_{k}^{ij}$ is the direct ancestor of siblings $s_{i}^{j}$ at depth
$k$.

To evaluate the conditional probabilities of the $n$-grams we use a Katz back-
off model (Katz, 1987), although we note that our approach can easily be
extended to other language models. We choose this model for its simplicity
given its ability to assign a non-zero probability to expressions not in the
training set. We are therefore not simply searching for equations which have
appeared in the corpus, but rather are upweighting equations more similar to
those in the counts and ordering of their operators.

For this model, the conditional probability of the word (operator) $w_{i}$
following the sequence of words (phrase) $w_{i-n+1}\ldots w_{i-1}$ is

(14) $\begin{split}&P\left(w_{i}|w_{i-n+1},\ldots,w_{i-1}\right)\\\
&=\begin{cases}\mbox{$d_{w_{i-n+1},\ldots,w_{i}}\frac{C(w_{i-n+1},\ldots,w_{i})}{C(w_{i-n+1},\ldots,w_{i-1})}$}&C(w_{i-n+1},\ldots,w_{i})>k\\\
\mbox{$\alpha_{w_{i-n+1},\ldots,w_{i-1}}P\left(w_{i}|w_{i-n+2},\ldots,w_{i-1}\right)$}&\text{otherwise},\end{cases}\end{split}$

where $C(\cdot)$ is the number of times a given phrase has been seen in
training. $k$ is a hyperparameter determining how many times one needs to see
a phrase before using the $n$-gram. If a phrase appears fewer than $k$ times
the algorithm “backs off” to the $(n-1)$-gram to evaluate the probability. We
choose $k=0$ for simplicity. We compute the amount of discounting,
$d_{(\cdot)}$, to be $C^{\ast}(\cdot)/C(\cdot)$, where $C^{\ast}(\cdot)$ is
the Good-Turing estimate (GOOD, 1953; Gale and Sampson, 1995) of the number of
counts of a phrase. The left-over probability mass for the $(n-1)$-gram is

(15) $\begin{split}&\beta_{w_{i-n+1},\ldots,w_{i-1}}\\\
&=1-\sum_{w_{i}:C(w_{i-n+1},\ldots,w_{i})>k}d_{w_{i-n+1},\ldots,w_{i}}\frac{C(w_{i-n+1},\ldots,w_{i})}{C(w_{i-n+1},\ldots,w_{i-1})},\end{split}$

and thus the back-off weight is defined to be

(16)
$\alpha_{w_{i-n+1},\ldots,w_{i-1}}=\frac{\beta_{w_{i-n+1},\ldots,w_{i-1}}}{\sum_{w_{i}:C(w_{i-n+1},\ldots,w_{i})\leq
k}P\left(w_{i}|w_{i-n+2},\ldots,w_{i-1}\right)}.$

If no data are available to evaluate these formulae, then one recursively
considers the next-shortest $n$-grams until such data exists.

For any function (whether from the corpus or if we are evaluating the prior
for this function), we first read the mathematical expression using the sympy
package (Meurer et al., 2017) and convert this to a function tree using esr
(Bartlett et al., 2022). We then identify all sets of siblings and the direct
ancestors of these siblings, which are the “phrases” contained in this
function. The set of words used for this phrase are then the set of ancestor
nodes and the terminating sibling. We split each sibling into a left and right
node, and multiply the probability of obtaining the left node by the
conditional probability of obtaining the right node given the left one. Since
the structure of the phrases are different for these two probabilities (one
always ends with a single leaf node, but the other has two leaf nodes), we use
a different back-off model for each and multiply the results.

####  2.4.2. Example - Scientific equations

For the applications of the language model prior considered in this work, we
use a compilation of 161 scientific equations as our training data. This
contains the 100 equations in the Feynman Symbolic Regression Database (FSReD)
(Udrescu and Tegmark, 2020), which were selected from the Feynman Lectures on
Physics (Feynman et al., 2015; Feynman et al., 2006a, b). We also include the
20 “bonus” equations from this database, which were taken from other physics
textbooks (Goldstein et al., 2002; Jackson, 1999; Weinberg, 1972; Schwartz,
2014). Finally, inspired by Guimerà et al. (2020), we augment our training set
with a selection of 41 equations taken from pages linked to Wikipedia’s “List
of scientific equations named after people.” These cover equations found in
topics such as classical mechanics, electromagnetism and special relativity,
contain the operators $\\{+$, $-$, $\times$, $\div$, $\operatorname{pow}$,
$\sqrt{\cdot}$, $\exp$, $\log$, $\sin$, $\cos$, $\arcsin$, $\arccos$,
$\tanh\\}$ and are functions of between one and eleven variables.

The remaining degree of freedom is the length of the $n$-grams used to produce
prior probabilities. In Fig. 1 we investigate the effect of varying $n$ by
constructing a back-off model with the FSReD functions with different values.
We evaluate the prior probabilities of all equations which can be constructed
from the basis $\\{x$, $a$, $\operatorname{inv}$, $+$, $-$, $\times$, $\div$,
$\operatorname{pow}\\}$ up to complexity (number of tree nodes) 10. Here $n=1$
means that we do not consider the effect of any parents, for $n=2$ the
probabilities are conditioned on the parent nodes only, and for $n=3$ we
condition on both the parent and the grand-parent nodes. From the large
scatter when comparing $n=1$ and $n=2$ we learn that there is significant
information in the parent nodes. This is expected (given a node, some children
are more likely than others), illustrating the importance of modelling
correlations in the function prior. There is, however, less scatter in the
probabilities between $n=2$ and $n=3$, hence we choose $n=2$. We find that the
one-dimensional distributions are very similar for all $n$.

Figure 1.  One- and two-dimensional histograms of prior probabilities for all
functions up to complexity 10 using the basis $\\{x$, $a$,
$\operatorname{inv}$, $+$, $-$, $\times$, $\div$, $\operatorname{pow}\\}$. The
priors are based on the compilation of scientific equations described in
Section 2.4.2 and are evaluated using a back-off model with different lengths,
$1\leq n\leq 3$. For the two-dimensional plots, points lying on the red line
have an equal prior for both values of $n$ in the corresponding panel. In the
one-dimensional plots, the black lines are the distributions for the
corresponding $n$ and we plot the histograms for the other values of $n$ in
red for reference.

##  3\. Benchmark tests

We have discussed a range of methods for SR model selection from either a
Bayesian or MDL perspective, involving different choices of prior or using
more heuristic approaches. To study the effect of these choices we apply the
following six model selection procedures to a set of benchmark datasets:

  1. (1)

Likelihood: Choose the model with the maximum likelihood,
$P(D|\bm{\theta}_{i},f_{i})$.

  2. (2)

Score: As a heuristic approach, choose the model which has the best “score”,
defined to be the model on the Pareto front with the largest negative of the
derivative of log-likelihood with respect to complexity (number of nodes in
the tree representation). This is the choice adopted by PySR (Cranmer et al.,
2020; Cranmer, 2020). Note we always report the logarithm of the score due to
the potentially large range of values.

  3. (3)

MDL: Choose the MDL model using the prescription of Bartlett et al. (2022)
(Section 2.2).

  4. (4)

MDL+LM: Choose the MDL model according to Section 2.2, but using the prior on
functions described in Section 2.4 instead of the term $k\log n$.

  5. (5)

MDL+FBF+LM: Choose the MDL model using FBF priors on initially uniform
parameter priors and with a language model function prior.

  6. (6)

Bayes+FBF+LM: Choose the model with the largest Bayesian probability (Eq. 2)
using the approximation to the Bayesian evidence from a FBF procedure with
initially uniform parameter priors (Eq. 11) and using the language model prior
discussed in Section 2.4.

We choose these methods since they encompass a broad range of model selection
approaches, however we leave a more comprehensive survey of different methods
to future work.

To benchmark these methods, we consider the five standard test functions given
in Table 1 which were proposed as SR benchmarks by Uy et al. (2011) and Korns
(2011) and collated by McDermott et al. (2012), where we use domains
comparable to these references. Since our model selection methods are designed
for noisy data, we first generate $10^{5}$ values of $x$ from a uniform
distribution and measure the standard deviation of the resulting $y$ values.
For all our benchmark studies, we add Gaussian noise with a standard
deviation, $\sigma$, equal to half this value. Our likelihood is then

(17) $\log
P\left(D|f_{i},\bm{\theta}_{i}\right)=-\frac{1}{2\sigma^{2}}\sum_{a}\left(y_{a}-f_{i}(x_{a}|\bm{\theta}_{i})\right)^{2}+{\rm
const},$

where we sum over the mock data points $\\{(x_{a},y_{a})\\}$. To ensure our
results are robust to the random seed, for each function and choice of number
of data points, we generate five data sets which we independently analyse.

Since our focus here is on model selection rather than equation generation we
use the esr algorithm so that we are guaranteed to evaluate all functions of a
given complexity, rather than a subset as would be generated using a non-
exhaustive method. We use the operator basis set $\\{x$, $a$, $\sqrt{\cdot}$,
$(\cdot)^{2}$, $\sin$, $\cos$, $+$, $\times$, $-$, $\div$,
$\operatorname{pow}\\}$ and run esr to complexity 8. Although one of the
benchmark functions contains an exponential, we do not include this operator
since the generating function has a lower-complexity variant using the power
operator. We stress, however, that these model selection procedures do not
require an exhaustive search so could equally well be implemented in
conjunction with any other function generation algorithm.

Our results are plotted in Fig. 2, where we show the difference in each metric
between the top-ranked function and the truth (if the truth is not ranked top)
or the truth and the second-best function (if the truth is ranked top). We
consider the “truth” to be any function mathematically equivalent to the
generating function as given in Table 1. This is determined by the function
simplification step of esr and subsequent visual inspection. The model
selection method has successfully recovered the truth for all points above
zero in the plot, and the further above zero the stronger the preference for
the truth compared to any other function. If the truth lies in the top two
functions, we colour the point red, otherwise it is blue.

One sees that assessing functions solely based on their likelihood does not
find the truth to be in the top two equations for any of the cases considered.
This is as one would expect for noisy data: one can make a function
arbitrarily accurate by making it sufficiently complex, and thus the
relatively simple generating equations are not recovered and over-fitting is
preferred. The ‘score’ method performs better, almost always selecting the
truth for the Nguyen-8 function, although in many other cases the truth does
not appear anywhere in the function ranking. This is a consequence of the
score only being defined on the Pareto-front; if the truth does not lie on
this (even if it is infinitesimally far away), then this metric cannot
accommodate it. Perhaps more concerning is that, even if the truth is on the
Pareto-front, this metric does not necessarily prefer the truth as the number
of data points tends to infinity; the likelihood scales approximately linearly
with the number of data points, and thus the overall ranking does not
necessarily change.

Turning our attention to the MDL and Bayesian metrics, we find that these more
often select the truth than the other methods. For the simplest function
(Nguyen-8), the MDL prescription of Bartlett et al. (2022) performs best and
always selects the truth. We find that the methods involving the language
model prior do not perform as well for this equation, with the truth often
appearing as the third best function. However, we note that $y=x^{\theta_{0}}$
is frequently the second-best function in this case, which is a generalisation
of the truth. This is likely a consequence of variables being raised to a
power being more common in our training set than square-roots. For the Korns-1
and Korns-6 functions, the MDL approach outperforms the other methods,
although larger sample sizes are required to recover the truth than for the
Nguyen-8 data. The MDL+LM method is the only other procedure which can select
the truth for any run considered, correctly identifying the generating
equation for Korns-1 when $\sim 10^{4}$ data points are used. All MDL and
Bayesian methods perform similarly well for the Korns-4 dataset, with the
truth almost always being selected, contrary to the other approaches.

The language model methods outperform the other criteria when applied to the
Korns-7 dataset. Since the function prior is based on scientific equations and
we are applying this to a different context (namely functions used as SR
benchmarks), there is no reason to suppose a priori that this would be the
case. However, when one examines the truth, one can imagine a plethora of
scientific applications for such a function. As such, including a prior based
on these equations allows the truth to be selected with fewer data points than
would be required otherwise.

Dataset | Function | Complexity | Domain
---|---|---|---
Nguyen-8 | $\sqrt{x}$ | 2 | $[0,4]$
Korns-1 | $1.57+2.43x$ | 5 | $[-50,50]$
Korns-4 | $-2.13+0.13\sin(x)$ | 6 | $[-50,50]$
Korns-6 | $1.3+0.13\sqrt{x}$ | 6 | $[0,50]$
Korns-7 | $213.81(1-e^{-0.547x})$ | 7 | $[0,50]$
Table 1. Datasets used as benchmarks, taken from (Uy et al., 2011) (Nguyen)
and (Korns, 2011) (Korns), as collated in (McDermott et al., 2012). Figure 2.
Preference for the truth for benchmark tests according to different model-
selection methods as a function of number of data points, when considering
functions up to complexity 8 with the basis $\\{x$, $a$, $\sqrt{\cdot}$,
$(\cdot)^{2}$, $\sin$, $\cos$, $+$, $\times$, $-$, $\div$,
$\operatorname{pow}\\}$. Any point above zero favours the truth. Red points
indicate runs where the truth is in the top two functions, otherwise the point
is blue.

##  4\. Physically reasonable functions

We now apply our language model prior to the context for which it is designed:
to investigate how it affects the ranking of functions we apply it to a real-
world SR problem drawn from the field of cosmology. We reanalyse the results
of Bartlett et al. (2022) where the expansion rate of the Universe, $H$, was
determined as a function of the redshift, $z$ (the fractional change in a
photon’s wavelength due to the expansion of the Universe; a higher redshift is
an earlier time in the history of the Universe) using the esr algorithm. This
is a univariate SR problem, where one wishes to find $y\left(x\equiv
1+z\right)\equiv\left(H\left(z\right)\right)^{2}$, which, according to the
concordance cosmological model (called $\rm\Lambda CDM$), satisfies the
Friedmann equation

(18) $y_{\rm\Lambda CDM}\left(x\right)=\theta_{0}+\theta_{1}x^{3}.$

We note that this equation does not appear in our corpus, so we expect to
‘discover’ a new equation in this case study.

In particular we focus on the analysis of the Pantheon+ sample of 1590
Cepheid-calibrated supernovae (Scolnic et al., 2021; Riess et al., 2022). The
observable here is not $y$, but the “distance modulus” $\mu$ which can be
calculated by computing an integral of $y$ (see (Bartlett et al., 2022)). The
likelihood in this case is a modelled as a multivariate Gaussian

(19) $\begin{split}\log
P\left(D|f_{i},\bm{\theta}_{i}\right)=&-\frac{1}{2}\left(\bm{\mu}_{i}\left(\bm{z}^{\mathrm{o}}\right)-\bm{\mu}^{\mathrm{o}}\right)^{\rm
T}\mathbf{\Sigma}^{-1}\left(\bm{\mu}_{i}\left(\bm{z}^{\mathrm{o}}\right)-\bm{\mu}^{\mathrm{o}}\right)\\\
&-\frac{1}{2}\log\left(\det\left(2\pi\mathbf{\Sigma}\right)\right),\end{split}$

where we use the publicly available covariance matrix, $\bm{\Sigma}$ (Brout et
al., 2022), $\bm{\mu}_{i}$ is the vector of predicted distance moduli at
redshifts $\bm{z}^{\mathrm{o}}$ using equation $f_{i}$ with parameters
$\bm{\theta}_{i}$, and $\bm{\mu}^{\mathrm{o}}$ are the observed values.

In the absence of a prior on functions from the literature, the functions were
ranked according to their description length (LABEL:{eq:mdl}). Although esr
was run up to complexity 10, the MDL function was found to be

(20) $y_{\rm
MDL}\left(x\right)=\theta_{0}x^{x}=\theta_{0}\left(1+z\right)^{\left(1+z\right)},$

which has a complexity of 5. With very few nodes in the tree representation
and a small number of operators, it is clear that such a function would have a
small penalisation due to model complexity in the MDL framework. However,
$x^{x}$ rarely appears in empirical science and one would unlikely expect a
priori for an equation containing it to describes the expansion of the
Universe.

We therefore re-evaluate the functions considered by Bartlett et al. but using
the six methods outlined in Section 3. We give the top-4 functions for each
method in Table 2, keeping the highest-ranked variant for equivalent
expressions and combining functions of parameters into new parameters during
post-processing. The functions given in Table 2 all produce reasonable
likelihoods, with a mean Mahalanobis distance
($N^{-1}\left(\bm{\mu}_{i}\left(\bm{z}^{\mathrm{o}}\right)-\bm{\mu}^{\mathrm{o}}\right)^{\rm
T}\mathbf{\Sigma}^{-1}\left(\bm{\mu}_{i}\left(\bm{z}^{\mathrm{o}}\right)-\bm{\mu}^{\mathrm{o}}\right)$
for $N=1590$) between 0.8783 and 0.9223.

When one only replaces the $k\log n$ term of Section 2.2 with the language
model prior, one sees that the physically unreasonable $x^{x}$ functions are
no longer in the top-four. Instead, all of these equations are merely sums of
$x$ to a (positive or negative) integer power, with at least one free
parameter to set the scale of the expansion rate. Moreover, Eq. 18 (which has
a strong theoretical prior based on General Relativity) is ranked much higher
(fourth) with this method. This suggests that the language model, as expected,
is helping to eliminate the physically unreasonable functions. Interestingly,
we see that the language model often selects functions which divide by
parameters, rather than multiply. In the absence of a language model prior,
this can occur due to the different Fisher terms of $\theta_{0}$ and
$1/\theta_{0}$ in Section 2.2, however we note that our function prior also
often prefers division. For example, the function ${x(2x+1/{x})}/{\theta_{0}}$
is preferred to ${x(2x+1/{x})}\times{\theta_{0}}$ by a change in the log-prior
of 0.41. Choosing a different corpus of equations could reverse this
preference, highlighting the importance of using an appropriate prior for the
given context.

Turning our attention to the criteria involving the FBF, the results are more
mixed. Although the $x^{x}$ behaviour is removed, the highest ranked equation
for both methods is still somewhat unappealing. Each of the constituent
elements of the function are however reasonable and thus this function’s high
ranking may be due to the low number of consecutive operators considered by
the back-off model. Alternatively it could be due to a well known problem with
Katz back-off models, that if a phrase has not been seen before but is
comprised of two often-encountered sub-phrases then it will necessarily be
assigned a high probability (Manning and Schutze, 1999). Thus this function
may be down-weighted with a more sophisticated language model. The remaining
functions all appear reasonable, and the Bayesian method finds Eq. 18 in the
top-four equations.

Method | Rank 1 | Rank 2 | Rank 3 | Rank 4
---|---|---|---|---
| Function | Function | Value | Function | Value | Function | Value
Likelihood | $\left|{\theta_{1}+\frac{1}{\theta_{0}-x}}\right|^{\left|{\theta_{2}}\right|^{x}}$ | $\left|{\theta_{0}-x^{\theta_{1}}}\right|^{\theta_{2}-x}$ | 1.48 | $-\theta_{2}+\left|{\theta_{0}\left(\theta_{1}-x\right)}\right|^{x}$ | 3.15 | $\frac{\left|{\theta_{1}+x^{\theta_{2}}}\right|^{x}}{\theta_{0}}$ | 3.18
Score | $\theta_{0}x$ | $\theta_{0}\left|{\theta_{1}}\right|^{-x}$ | 11.54 | $\left|{\theta_{0}-x^{\theta_{1}}}\right|^{\theta_{2}-x}$ | 14.32 | $\left|{\theta_{1}+\frac{1}{\theta_{0}-x}}\right|^{\left|{\theta_{2}}\right|^{x}}$ | 14.59
MDL | $\theta_{0}x^{x}$ | $\left|{\theta_{0}}\right|^{x^{\theta_{1}}}$ | 0.78 | $\theta_{0}\left|{\theta_{1}}\right|^{-x}$ | 0.83 | $\theta_{0}x^{x^{\theta_{1}}}$ | 0.91
MDL+LM | $\frac{x\left(2x+\frac{1}{x}\right)}{\theta_{0}}$ | $\frac{2x^{2}+\frac{1}{x}}{\theta_{0}}$ | 1.40 | $\frac{x^{2}+2x}{\theta_{0}}$ | 1.63 | $\theta_{0}+\theta_{1}x^{3}$ | 1.72
MDL+FBF+LM | $\left|{\theta_{0}}\right|^{\frac{1}{\theta_{1}+2x}}$ | $\left|{\theta_{0}}\right|^{\left|\theta_{1}\right|^{x}}$ | 0.16 | $\theta_{0}\left|{\theta_{1}}\right|^{x}$ | 1.77 | $\frac{\theta_{1}-x^{2}}{\theta_{0}}$ | 2.19
Bayes+FBF+LM | $\left|{\theta_{0}}\right|^{\frac{1}{\theta_{1}+2x}}$ | $\left|{\theta_{0}}\right|^{\left|\theta_{1}\right|^{x}}$ | 0.16 | $\left|{\theta_{0}}\right|^{-\theta_{1}+x}$ | 1.03 | $\frac{x^{3}+\frac{1}{\theta_{0}}}{\theta_{1}}$ | 1.95
Table 2. Four highest ranked functions for $y(x\equiv 1+z)=H^{2}(z)$ for the
Pantheon+ supernova sample according to different model selection methods.
Here, MDL, FBF and LM mean Minimum Description Length, Fractional Bayes Factor
and Language Model, respectively (see text for details). The values of each
metric are given relative to the highest-ranked function.

##  5\. Related work

The idea of using transformers as a proposal strategy for generating candidate
equations given a set of observed input-output pairs has been developed before
(Valipour et al., 2021; Kamienny et al., 2022; Biggio et al., 2021; Shojaee et
al., 2023; Vastl et al., 2022; d’Ascoli et al., 2022; Becker et al., 2022).
Typically, the transformers are trained on several million functions, which
are generated randomly by first producing candidate trees (see e.g. (Lample
and Charton, 2019)) and then decorating these with operators from a given
basis set. The range of input variables and parameters of the functions are
then drawn from some random distribution and the function evaluated. Once
trained, the transformer takes as input the observed inputs and outputs and
suggests functions which could have generated the data. The constants of these
equations can then be fine tuned by a further optimisation step. One can
interface the decoder with a Monte Carlo Tree Search to provide feedback
during the generation process to help select equations which optimise some
metric of accuracy and/or complexity (Shojaee et al., 2023). Similarly, one
can use a probabilistic grammar as a generative model (Brence et al., 2021),
such that different productions have different probabilities, which can be
tuned based on intuition or through experiments.

The approach here is different. Previous works use the model to generate
candidate equations, but not to assess the appropriateness of the function,
since typically the candidate with the highest likelihood (or lowest mean
squared error) is selected. If, instead, one viewed the generation as a Monte
Carlo sampling of the functional prior and marginalised results over the
generated equations, then the two approaches would be more aligned. Moreover,
we note that, for transformers, the training-set functions are generated
randomly, with no context-specific knowledge (beyond the number of nodes,
depth of tree, range of parameters, range of inputs or possible basis
functions). Here we adopt a conceptually different approach, where the
training set informs us which combinations of operators are a priori more
likely. Of course, the hyperparameters which control the generation of the
training set or the probabilities in the probabilistic grammar could be tuned
to act as a prior over equations (Becker et al., 2022). Prior knowledge has
been previously incorporated into neural network based SR methods, provided
that this knowledge can be expressed as a (nonlinear) inequality or equality
(Kubalík et al., 2023), which is a different scenario to that considered here.
Finally, we note that our language model is designed not to have a limit on
the number or range of the input variables; this can be a limitation of the
transformer method because changing the input variables can require retraining
(Vastl et al., 2022).

In this work we consider the effects of accuracy, a prior on parameters and a
prior on functions in model selection for SR. The first of these is always
considered in SR (although the metric for it varies), whereas the latter two
are less commonly used. In particular, Guimerà et al. (2020) introduce a prior
on functions which assumes each node of the tree to be independent, but they
do not account for the importance of the parameter prior when computing the
evidence (which we argue in Section 2.3 is important). Werner et al. (2022)
overcome the issue of parameter priors (insofar as they choose a prior which
can be normalised) by choosing Gaussian, zero-mean, isotropic priors for all
parameters with a scalar hyperparameter which is optimised to maximise the
marginal likelihood. One could question whether zero-mean or isotropic are
good assumptions in general, and they assume all functions are equally likely
a priori. Similarly, Bomarito et al. (2023) assume equal prior probability for
all functions, but overcome the issue of priors by using the Fractional Bayes
Factor (Eq. 8). Bartlett et al. (2022) implicitly consider all three of these
concerns through an implementation of the MDL principle, where the parameter
prior corresponds to Rissanen’s universal code for the integers (Rissanen,
1983) and the functional prior depends on the number of nodes, unique
operators, parameters and the values of constants (c.f. Eq. 7). This does not
allow a general parameter prior, and its functional prior is uninformed by the
frequency of occurence of combinations of operators in the domain of interest.

##  6\. Conclusions

Model selection for SR trades off accuracy and complexity to obtain simple
functions that describe a data set well. Previous attempts to do this have
either used a heuristic metric, treated all functions as equally likely a
priori or neglected to construct a self-consistent parameter prior in a
Bayesian framework. In the case where a structural prior on functions based on
the frequency of operators in a corpus of equations had been used, the
correlation between operators had been neglected. In this paper we attempt to
reconcile these issues by formalising and investigating a range of model
selection criteria based on the Bayesian evidence or Minimum Description
Length principle. We explicitly compare the Bayesian and MDL procedures,
implement the Fractional Bayes Factor to overcome the dependence of the Bayes
factor on the arbitrary range of uninformative parameter priors and construct
a Katz back-off language model on a set of scientific equations to learn a
prior on operator ordering that can in principle take a function’s entire tree
into account. We compare our models with literature standards on benchmarks
from the Nguyen and Korns datasets as well as a real-world test-case in the
field of cosmology. For standard SR benchmarks, we find that the MDL procedure
without a prior on functions typically performs the best, however this method
produced physically unrealistic equations on the astrophysical data. This is
alleviated by applying the language model prior, suggesting this to be a
promising avenue for SR modelling when domain-specific prior information is
available. We make the back-off model code and training equations publicly
available (https://github.com/DeaglanBartlett/katz).

###### Acknowledgements.

We thank Alan Heavens, Florent Leclercq and Roberto Trotta for useful
discussions. DJB is supported by the Simons Collaboration on “Learning the
Universe.” HD is supported by a Royal Society University Research Fellowship
(grant no. 211046). PGF acknowledges support from European Research Council
Grant No: 693024, STFC and the Beecroft Trust. For the purpose of open access,
we have applied a Creative Commons Attribution (CC BY) licence to any Author
Accepted Manuscript version arising.

## References

  * (1)
  * Bartlett et al. (2022) Deaglan J. Bartlett, Harry Desmond, and Pedro G. Ferreira. 2022.  Exhaustive Symbolic Regression.  _arXiv e-prints_ , Article arXiv:2211.11461 (Nov. 2022), arXiv:2211.11461 pages.  https://doi.org/10.48550/arXiv.2211.11461 arXiv:2211.11461 [astro-ph.CO]
  * Becker et al. (2022) Sören Becker, Michal Klein, Alexander Neitz, Giambattista Parascandolo, and Niki Kilbertus. 2022\.  Discovering ordinary differential equations that govern time-series.  _arXiv e-prints_ , Article arXiv:2211.02830 (Nov. 2022), arXiv:2211.02830 pages.  https://doi.org/10.48550/arXiv.2211.02830 arXiv:2211.02830 [cs.LG]
  * Berger and Pericchi (1996) James O. Berger and Luis R. Pericchi. 1996.  The Intrinsic Bayes Factor for Model Selection and Prediction.  _J. Amer. Statist. Assoc._ 91, 433 (1996), 109–122.  https://doi.org/10.1080/01621459.1996.10476668 arXiv:https://www.tandfonline.com/doi/pdf/10.1080/01621459.1996.10476668
  * Bernardo (1979) Jose M. Bernardo. 1979\.  Reference Posterior Distributions for Bayesian Inference.  _Journal of the Royal Statistical Society: Series B (Methodological)_ 41, 2 (1979), 113–128.  https://doi.org/10.1111/j.2517-6161.1979.tb01066.x arXiv:https://rss.onlinelibrary.wiley.com/doi/pdf/10.1111/j.2517-6161.1979.tb01066.x
  * Biggio et al. (2021) Luca Biggio, Tommaso Bendinelli, Alexander Neitz, Aurelien Lucchi, and Giambattista Parascandolo. 2021\.  Neural Symbolic Regression that scales. In _Proceedings of the 38th International Conference on Machine Learning_ _(Proceedings of Machine Learning Research, Vol. 139)_ , Marina Meila and Tong Zhang (Eds.). PMLR, 936–945.  https://proceedings.mlr.press/v139/biggio21a.html
  * Bomarito et al. (2022) G. F. Bomarito, P. E. Leser, N. C. M. Strauss, K. M. Garbrecht, and J. D. Hochhalter. 2022\.  Bayesian Model Selection for Reducing Bloat and Overfitting in Genetic Programming for Symbolic Regression. In _Proceedings of the Genetic and Evolutionary Computation Conference Companion_ (Boston, Massachusetts) _(GECCO ’22)_. Association for Computing Machinery, New York, NY, USA, 526–529.  https://doi.org/10.1145/3520304.3528899
  * Bomarito et al. (2023) G. F. Bomarito, P. E. Leser, N. C. M. Strauss, K. M. Garbrecht, and J. D. Hochhalter. 2023\.  Automated learning of interpretable models with quantified uncertainty.  _Computer Methods in Applied Mechanics and Engineering_ 403 (Jan. 2023), 115732.  https://doi.org/10.1016/j.cma.2022.115732 arXiv:2205.01626 [cs.NE]
  * Brence et al. (2021) Jure Brence, Ljupčo Todorovski, and Sašo Džeroski. 2021\.  Probabilistic grammars for equation discovery.  _Knowledge-Based Systems_ 224 (2021), 107077.  https://doi.org/10.1016/j.knosys.2021.107077
  * Brodtkorb and D’Errico (2015) Per A Brodtkorb and John D’Errico. 2015.  numdifftools 0.9.11.  https://github.com/pbrod/numdifftools.
  * Brout et al. (2022) Dillon Brout et al. 2022\.  The Pantheon+ Analysis: Cosmological Constraints.  _arXiv e-prints_ , Article arXiv:2202.04077 (Feb. 2022), arXiv:2202.04077 pages.  arXiv:2202.04077 [astro-ph.CO]
  * Chen et al. (2022) Qi Chen, Bing Xue, and Mengjie Zhang. 2022.  Rademacher Complexity for Enhancing the Generalization of Genetic Programming for Symbolic Regression.  _IEEE Transactions on Cybernetics_ 52, 4 (2022), 2382–2395.  https://doi.org/10.1109/TCYB.2020.3004361
  * Conway et al. (2013) J.H. Conway, N.J.A. Sloane, E. Bannai, R.E. Borcherds, J. Leech, S.P. Norton, A.M. Odlyzko, R.A. Parker, L. Queen, and B.B. Venkov. 2013\.  _Sphere Packings, Lattices and Groups_.  Springer New York.  https://books.google.fr/books?id=hoTjBwAAQBAJ
  * Cover and Thomas (1991) T. M. Cover and J. A. Thomas. 1991.  _Elements of Information Theory_ (2nd ed.).  Wiley.
  * Cranmer (2020) Miles Cranmer. 2020\.  PySR: Fast & Parallelized Symbolic Regression in Python/Julia.  https://doi.org/10.5281/zenodo.4041459
  * Cranmer et al. (2020) Miles Cranmer, Alvaro Sanchez-Gonzalez, Peter Battaglia, Rui Xu, Kyle Cranmer, David Spergel, and Shirley Ho. 2020.  Discovering Symbolic Models from Deep Learning with Inductive Biases.  _NeurIPS 2020_ (2020).  arXiv:2006.11287 [cs.LG]
  * d’Ascoli et al. (2022) Stéphane d’Ascoli, Pierre-Alexandre Kamienny, Guillaume Lample, and François Charton. 2022\.  Deep Symbolic Regression for Recurrent Sequences.  _arXiv e-prints_ , Article arXiv:2201.04600 (Jan. 2022), arXiv:2201.04600 pages.  https://doi.org/10.48550/arXiv.2201.04600 arXiv:2201.04600 [cs.LG]
  * David (1989) E. David. 1989.  _Genetic Algorithms in Search, Optimization and Machine Learning_.  Addison-Wesley.
  * Feynman et al. (2015) R.P. Feynman, R.B. Leighton, and M. Sands. 2015\.  _The Feynman Lectures on Physics, Vol. I: The New Millennium Edition: Mainly Mechanics, Radiation, and Heat_.  Number v. 1. Basic Books.  https://books.google.fr/books?id=d76DBQAAQBAJ
  * Feynman et al. (2006a) R.P. Feynman, R.B. Leighton, M.L. Sands, and M.A. Gottlieb. 2006a.  _The Feynman Lectures on Physics_.  Number v. 2 in The Feynman Lectures on Physics. Pearson/Addison-Wesley.  https://books.google.fr/books?id=AbruAAAAMAAJ
  * Feynman et al. (2006b) R.P. Feynman, R.B. Leighton, M.L. Sands, and M.A. Gottlieb. 2006b.  _The Feynman Lectures on Physics: Quantum mechanics_.  Pearson/Addison-Wesley.  https://books.google.fr/books?id=_6XvAAAAMAAJ
  * Gale and Sampson (1995) William A. Gale and Geoffrey Sampson. 1995.  Good‐turing frequency estimation without tears.  _Journal of Quantitative Linguistics_ 2, 3 (1995), 217–237.  https://doi.org/10.1080/09296179508590051 arXiv:https://doi.org/10.1080/09296179508590051
  * Goldstein et al. (2002) H. Goldstein, C.P. Poole, and J.L. Safko. 2002.  _Classical Mechanics_.  Addison Wesley.  https://books.google.fr/books?id=tJCuQgAACAAJ
  * GOOD (1953) I. J. GOOD. 1953\.  THE POPULATION FREQUENCIES OF SPECIES AND THE ESTIMATION OF POPULATION PARAMETERS.  _Biometrika_ 40, 3-4 (12 1953), 237–264.  https://doi.org/10.1093/biomet/40.3-4.237 arXiv:https://academic.oup.com/biomet/article-pdf/40/3-4/237/492571/40-3-4-237.pdf
  * Grunwald (2007) P. Grunwald. 2007\.  _The Minimum Description Length Principle_.  MIT Press.
  * Grünwald and Roos (2019) Peter Grünwald and Teemu Roos. 2019.  Minimum Description Length Revisited.  _arXiv e-prints_ , Article arXiv:1908.08484 (Aug. 2019), arXiv:1908.08484 pages.  arXiv:1908.08484 [stat.ME]
  * Guimerà et al. (2020) Roger Guimerà, Ignasi Reichardt, Antoni Aguilar-Mogas, Francesco A. Massucci, Manuel Miranda, Jordi Pallarès, and Marta Sales-Pardo. 2020.  A Bayesian machine scientist to aid in the solution of challenging scientific problems.  _Science Advances_ 6, 5 (2020), eaav6971.  https://doi.org/10.1126/sciadv.aav6971 arXiv:https://www.science.org/doi/pdf/10.1126/sciadv.aav6971
  * Haupt and Haupt (2004) R. Haupt and S. Haupt. 2004\.  _Practical genetic algorithms_ (2nd ed.).  Wyley.
  * Iba et al. (1994) Hitoshi Iba, Hugo De Garis, and Taisuke Sato. 1994\.  Genetic programming using a minimum description length principle.  _Advances in genetic programming_ 1 (1994), 265–284.
  * Jackson (1999) John David Jackson. 1999\.  _Classical electrodynamics_ (3rd ed. ed.).  Wiley, New York, NY.  http://cdsweb.cern.ch/record/490457
  * Jeffreys (1946) Harold Jeffreys. 1946\.  An Invariant Form for the Prior Probability in Estimation Problems.  _Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences_ 186, 1007 (1946), 453–461.  http://www.jstor.org/stable/97883
  * Kamienny et al. (2022) Pierre-Alexandre Kamienny, Stéphane d’Ascoli, Guillaume Lample, and François Charton. 2022\.  End-to-end symbolic regression with transformers.  _arXiv e-prints_ , Article arXiv:2204.10532 (April 2022), arXiv:2204.10532 pages.  https://doi.org/10.48550/arXiv.2204.10532 arXiv:2204.10532 [cs.LG]
  * Kammerer et al. (2021) Lukas Kammerer, Gabriel Kronberger, Bogdan Burlacu, Stephan M. Winkler, Michael Kommenda, and Michael Affenzeller. 2021\.  Symbolic Regression by Exhaustive Search: Reducing the Search Space Using Syntactical Constraints and Efficient Semantic Structure Deduplication.  _arXiv e-prints_ , Article arXiv:2109.13895 (Sept. 2021), arXiv:2109.13895 pages.  arXiv:2109.13895 [cs.LG]
  * Katz (1987) Slava M. Katz. 1987\.  Estimation of probabilities from sparse data for the language model component of a speech recognizer.  _IEEE Trans. Acoust. Speech Signal Process._ 35 (1987), 400–401.
  * Kommenda et al. (2021) Michael Kommenda, Andreas Beham, Michael Affenzeller, and Gabriel Kronberger. 2021\.  Complexity Measures for Multi-objective Symbolic Regression.  _arXiv e-prints_ , Article arXiv:2109.00238 (Sept. 2021), arXiv:2109.00238 pages.  arXiv:2109.00238 [cs.LG]
  * Korns (2011) Michael Korns. 2011\.  _Accuracy in Symbolic Regression_.  129–151.  https://doi.org/10.1007/978-1-4614-1770-5_8
  * Kubalík et al. (2023) Jiří Kubalík, Erik Derner, and Robert Babuška. 2023.  Neural Networks for Symbolic Regression.  _arXiv e-prints_ , Article arXiv:2302.00773 (Feb. 2023), arXiv:2302.00773 pages.  https://doi.org/10.48550/arXiv.2302.00773 arXiv:2302.00773 [cs.NE]
  * Lample and Charton (2019) Guillaume Lample and François Charton. 2019.  Deep Learning for Symbolic Mathematics.  _arXiv e-prints_ , Article arXiv:1912.01412 (Dec. 2019), arXiv:1912.01412 pages.  https://doi.org/10.48550/arXiv.1912.01412 arXiv:1912.01412 [cs.SC]
  * Lanterman (2001) Aaron D. Lanterman. 2001\.  Schwarz, Wallace, and Rissanen: Intertwining Themes in Theories of Model Selection.  _International Statistical Review_ 69, 2 (2001), 185–212.  https://doi.org/10.1111/j.1751-5823.2001.tb00456.x arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1111/j.1751-5823.2001.tb00456.x
  * Lemos et al. (2022) Pablo Lemos, Niall Jeffrey, Miles Cranmer, Shirley Ho, and Peter Battaglia. 2022\.  Rediscovering orbital mechanics with machine learning.  _arXiv e-prints_ , Article arXiv:2202.02306 (Feb. 2022), arXiv:2202.02306 pages.  arXiv:2202.02306 [astro-ph.EP]
  * Lempers (1971) F. B. Lempers. 1971\.  _Posterior probabilities of alternative linear models. Some theoretical considerations and empirical experiments. [By] F. B. Lempers_.  Rotterdam University Press [Rotterdam]. 118 p. pages.
  * Manning and Schutze (1999) Christopher D. Manning and Hinrich Schutze. 1999.  _Foundations of statistical natural language processing_.  MIT Press, Cambridge, Mass.
  * McConaghy (2011) Trent McConaghy. 2011\.  _FFX: Fast, Scalable, Deterministic Symbolic Regression Technology_.  Springer New York, New York, NY, 235–260.  https://doi.org/10.1007/978-1-4614-1770-5_13
  * McDermott et al. (2012) James McDermott, David R. White, Sean Luke, Luca Manzoni, Mauro Castelli, Leonardo Vanneschi, Wojciech Jaśkowski, Krzysztof Krawiec, Robin Harper, Kenneth De Jong, and Una May O’Reilly. 2012\.  Genetic programming needs better benchmarks. In _GECCO’12 - Proceedings of the 14th International Conference on Genetic and Evolutionary Computation_ _(GECCO’12 - Proceedings of the 14th International Conference on Genetic and Evolutionary Computation)_. 791–798.  https://doi.org/10.1145/2330163.2330273 14th International Conference on Genetic and Evolutionary Computation, GECCO’12 ; Conference date: 07-07-2012 Through 11-07-2012.
  * Meurer et al. (2017) Aaron Meurer et al. 2017\.  SymPy: symbolic computing in Python.  _PeerJ Computer Science_ 3 (Jan. 2017), e103.  https://doi.org/10.7717/peerj-cs.103
  * O’Hagan (1995) Anthony O’Hagan. 1995\.  Fractional Bayes Factors for Model Comparison.  _Journal of the Royal Statistical Society. Series B (Methodological)_ 57, 1 (1995), 99–138.  http://www.jstor.org/stable/2346088
  * Riess et al. (2022) Adam G. Riess, Wenlong Yuan, Lucas M. Macri, Dan Scolnic, Dillon Brout, Stefano Casertano, David O. Jones, Yukei Murakami, Gagandeep S. Anand, Louise Breuval, Thomas G. Brink, Alexei V. Filippenko, Samantha Hoffmann, Saurabh W. Jha, W. D’arcy Kenworthy, John Mackenty, Benjamin E. Stahl, and WeiKang Zheng. 2022.  A Comprehensive Measurement of the Local Value of the Hubble Constant with 1 km s-1 Mpc-1 Uncertainty from the Hubble Space Telescope and the SH0ES Team.  __934, 1, Article L7 (July 2022), L7 pages. https://doi.org/10.3847/2041-8213/ac5c5b arXiv:2112.04510 [astro-ph.CO]
  * Rissanen (1978) J. Rissanen. 1978\.  Modeling by shortest data description.  _Automatica_ 14, 5 (1978), 465–471.  https://doi.org/10.1016/0005-1098(78)90005-5
  * Rissanen (1983) Jorma Rissanen. 1983\.  A Universal Prior for Integers and Estimation by Minimum Description Length.  _The Annals of Statistics_ 11, 2 (1983), 416 – 431.  https://doi.org/10.1214/aos/1176346150
  * Rivero and Fernandez-Blanco (2019) Daniel Rivero and Enrique Fernandez-Blanco. 2019\.  A New Deterministic Technique for Symbolic Regression.  _arXiv e-prints_ , Article arXiv:1908.06754 (Aug. 2019), arXiv:1908.06754 pages.  arXiv:1908.06754 [cs.LG]
  * Schwartz (2014) M.D. Schwartz. 2014\.  _Quantum Field Theory and the Standard Model_.  Cambridge University Press.  https://books.google.fr/books?id=HbdEAgAAQBAJ
  * Schwarz (1978) Gideon Schwarz. 1978\.  Estimating the Dimension of a Model.  _Annals of Statistics_ 6, 2 (July 1978), 461–464.
  * Scolnic et al. (2021) Dan Scolnic, Dillon Brout, Anthony Carr, Adam G. Riess, Tamara M. Davis, Arianna Dwomoh, David O. Jones, Noor Ali, Pranav Charvu, Rebecca Chen, Erik R. Peterson, Brodie Popovic, Benjamin M. Rose, Charlotte Wood, Peter J. Brown, Ken Chambers, David A. Coulter, Kyle G. Dettman, Georgios Dimitriadis, Alexei V. Filippenko, Ryan J. Foley, Saurabh W. Jha, Charles D. Kilpatrick, Robert P. Kirshner, Yen-Chen Pan, Armin Rest, Cesar Rojas-Bravo, Matthew R. Siebert, Benjamin E. Stahl, and WeiKang Zheng. 2021.  The Pantheon+ Analysis: The Full Dataset and Light-Curve Release.  _arXiv e-prints_ , Article arXiv:2112.03863 (Dec. 2021), arXiv:2112.03863 pages.  arXiv:2112.03863 [astro-ph.CO]
  * Shojaee et al. (2023) Parshin Shojaee, Kazem Meidani, Amir Barati Farimani, and Chandan K. Reddy. 2023.  Transformer-based Planning for Symbolic Regression.  _arXiv e-prints_ , Article arXiv:2303.06833 (March 2023), arXiv:2303.06833 pages.  https://doi.org/10.48550/arXiv.2303.06833 arXiv:2303.06833 [cs.LG]
  * Smits and Kotanchek (2005) Guido F. Smits and Mark Kotanchek. 2005.  _Pareto-Front Exploitation in Symbolic Regression_.  Springer US, Boston, MA, 283–299.  https://doi.org/10.1007/0-387-23254-0_17
  * Sutton and Barto (2018) R.S. Sutton and A.G. Barto. 2018.  _Reinforcement Learning, second edition: An Introduction_.  MIT Press.  https://books.google.co.uk/books?id=sWV0DwAAQBAJ
  * Takeuchi (1997) J. Takeuchi. 1997\.  Characterization of the Bayes estimator and the MDL estimator for exponential families.  _IEEE Transactions on Information Theory_ 43, 4 (1997), 1165–1174.  https://doi.org/10.1109/18.605579
  * Turing (1950) A. M. Turing. 1950\.  I.—COMPUTING MACHINERY AND INTELLIGENCE.  _Mind_ LIX, 236 (10 1950), 433–460.  https://doi.org/10.1093/mind/LIX.236.433
  * Udrescu et al. (2020) Silviu-Marian Udrescu, Andrew Tan, Jiahai Feng, Orisvaldo Neto, Tailin Wu, and Max Tegmark. 2020.  AI Feynman 2.0: Pareto-optimal symbolic regression exploiting graph modularity.  _arXiv e-prints_ , Article arXiv:2006.10782 (June 2020), arXiv:2006.10782 pages.  arXiv:2006.10782 [cs.LG]
  * Udrescu and Tegmark (2020) Silviu-Marian Udrescu and Max Tegmark. 2020.  AI Feynman: A physics-inspired method for symbolic regression.  _Science Advances_ 6, 16 (April 2020), eaay2631.  https://doi.org/10.1126/sciadv.aay2631 arXiv:1905.11481 [physics.comp-ph]
  * Uy et al. (2011) Nguyen Quang Uy, Nguyen Xuan Hoai, Michael O’Neill, Robert I. McKay, and Edgar Galván López. 2011.  Semantically-based crossover in genetic programming: application to real-valued symbolic regression.  _Genet. Program. Evolvable Mach._ 12, 2 (2011), 91–119.  https://doi.org/10.1007/s10710-010-9121-2
  * Valipour et al. (2021) Mojtaba Valipour, Bowen You, Maysum Panju, and Ali Ghodsi. 2021.  SymbolicGPT: A Generative Transformer Model for Symbolic Regression.  _arXiv e-prints_ , Article arXiv:2106.14131 (June 2021), arXiv:2106.14131 pages.  https://doi.org/10.48550/arXiv.2106.14131 arXiv:2106.14131 [cs.LG]
  * Vastl et al. (2022) Martin Vastl, Jonáš Kulhánek, Jiří Kubalík, Erik Derner, and Robert Babuška. 2022\.  SymFormer: End-to-end symbolic regression using transformer-based architecture.  _arXiv e-prints_ , Article arXiv:2205.15764 (May 2022), arXiv:2205.15764 pages.  https://doi.org/10.48550/arXiv.2205.15764 arXiv:2205.15764 [cs.LG]
  * Vladislavleva et al. (2009) Ekaterina J. Vladislavleva, Guido F. Smits, and Dick den Hertog. 2009\.  Order of Nonlinearity as a Complexity Measure for Models Generated by Symbolic Regression via Pareto Genetic Programming.  _IEEE Transactions on Evolutionary Computation_ 13, 2 (2009), 333–349.  https://doi.org/10.1109/TEVC.2008.926486
  * Wallace and Freeman (1987) C. S. Wallace and P. R. Freeman. 1987.  Estimation and Inference by Compact Coding.  _Journal of the Royal Statistical Society. Series B (Methodological)_ 49, 3 (1987), 240–265.  http://www.jstor.org/stable/2985992
  * Wallace and Freeman (1992) C. S. Wallace and P. R. Freeman. 1992.  Single-Factor Analysis by Minimum Message Length Estimation.  _Journal of the Royal Statistical Society: Series B (Methodological)_ 54, 1 (1992), 195–209.  https://doi.org/10.1111/j.2517-6161.1992.tb01874.x arXiv:https://rss.onlinelibrary.wiley.com/doi/pdf/10.1111/j.2517-6161.1992.tb01874.x
  * Weinberg (1972) Steven Weinberg. 1972\.  _Gravitation and Cosmology: Principles and Applications of the General Theory of Relativity_.  John Wiley and Sons, New York.
  * Werner et al. (2022) Matthias Werner, Andrej Junginger, Philipp Hennig, and Georg Martius. 2022.  Uncertainty in equation learning. In _GECCO ’22: Genetic and Evolutionary Computation Conference, Companion Volume, Boston, Massachusetts, USA, July 9 - 13, 2022_ , Jonathan E. Fieldsend and Markus Wagner (Eds.). ACM, 2298–2305.  https://doi.org/10.1145/3520304.3533964
  * Worm and Chiu (2013) Tony Worm and Kenneth Chiu. 2013.  Prioritized Grammar Enumeration: Symbolic Regression by Dynamic Programming. In _Proceedings of the 15th Annual Conference on Genetic and Evolutionary Computation_ (Amsterdam, The Netherlands) _(GECCO ’13)_. Association for Computing Machinery, New York, NY, USA, 1021–1028.  https://doi.org/10.1145/2463372.2463486