On ab initio-based, free and closed-form expressions for gravitational waves

We introduce a new approach for finding high accuracy, free and closed-form expressions for the gravitational waves emitted by binary black hole collisions from ab initio models. More precisely, our expressions are built from numerical surrogate models based on supercomputer simulations of the Einstein equations, which have been shown to be essentially indistinguishable from each other. Distinct aspects of our approach are that: (i) representations of the gravitational waves can be explicitly written in a few lines, (ii) these representations are free-form yet still fast to search for and validate and (iii) there are no underlying physical approximations in the underlying model. The key strategy is combining techniques from Artificial Intelligence and Reduced Order Modeling for parameterized systems. Namely, symbolic regression through genetic programming combined with sparse representations in parameter space and the time domain using Reduced Basis and the Empirical Interpolation Method enabling fast free-form symbolic searches and large-scale a posteriori validations. As a proof of concept we present our results for the collision of two black holes, initially without spin, and with an initial separation corresponding to 25–31 gravitational wave cycles before merger. The minimum overlap, compared to ground truth solutions, is 99%. That is, 1% difference between our closed-form expressions and supercomputer simulations; this is considered for gravitational (GW) science more than the minimum required due to experimental numerical errors which otherwise dominate. This paper aims to contribute to the field of GWs in particular and Artificial Intelligence in general.

The GWs emmitted by binary systems can be seen as a parameterized problem, the parameters being for example the masses and spins of the binary components, initial separation, equations of state in the case of neutron stars, and so on. For each value of the parameter tuple, the GW consists of two degrees of polarization which can be encapsulated in a single complex-valued time (or, equivalently, frequency) series.
Over the last decade surrogate models for GWs based on modern frameworks such as Reduced Basis (RB) and the Empirical Interpolation Method (EIM) [7][8][9][9][10][11][12][13][14][15][16] have been developed. For a recent and detailed review and any doubts in what follows in this article regarding reduced order and surrogate modeling see Ref. 17 . The RB framework provides a sparse representation in the parameter domain while the EIM does so in the dual one, here being time. Both make use of the smooth dependence on parameter variation to achieve fast-usually exponential in the case of GWs-convergence of the representation with respect to the number of basis elements (which by construction equals the number of EIM time nodes). The resulting surrogate models predict GWs which are essentially indistinguishable from numerical relativity (NR) simulations, but with a speedup of evaluation of around eight orders of magnitude 9 , being evaluated in less than a second on a standard laptop instead of using supercomputers as in NR. A relatively small number of NR simulations are still needed in the offline (training) stage, though, but with a fast, highly accurate, surrogate, predictive model to evaluate for any parameter and time in the intervals considered in the online stage. This is referred to as an offline-online decomposition strategy.
In this article we build upon these efforts for what we consider a natural next step: a methodology for finding high accuracy free and closed-form representations of GWs, as opposed to numerical surrogates. As a proof of concept we present our results for the collision of two black holes, initially without spin, and with an initial separation corresponding to about 25-31 GW cycles before merger. The minimum overlap obtained, compared to ground truth solutions, is 99%. That is, 1% difference between our closed-form expressions and supercomputer simulations.

Method
Reduced order and surrogate models. In general, surrogate models for GWs have followed reduced order modeling (ROM) for parameterized systems and/or "standard" machine learning regression techniques 18 .
Here we focus on the former, we will not delve into reviewing them but instead, as mentioned, refer to 17 . In this work we focus on RB 19,20 and the EIM [21][22][23] . Briefly, RB collocates parameter points according to their relevance, which are used to build a hierarchical, nearly optimal basis in a rigorous mathematical sense with respect to the Kolmogorov n-width 24,25 . The framework of RB takes advantage of any regularity with respect to parameter variation to achieve fast convergence in the accuracy of the representation with the number of basis elements; it is usually referred to as an application-based spectral expansion. In fact, in the case of gravitational waves it can be easily argued that the parameter dependence is smooth ( C ∞ ) and RB has been shown to achieve asymptotic exponential convergence 7,8,10,[13][14][15]26 in practice, as expected with any spectral-type method. Similarly, the EIM achieves a subsampling in the space dual to that one of parameters (time in the case here considered) which is also nearly optimal. For a detailed discussion on the optimality of the EIM see 27 and references therein. On top of that, high accuracy predictive models (prediction as opposed to projection) can be built once one has a reduced basis and an empirical interpolant 13 . These predictive models are essentially indistinguishable from numerical relativity supercomputer simulations of the Einstein equations but can be evaluated in less than a second on a laptop. The availability of high accuracy, fast to evaluate, numerical predictive models and a sparse subsampling in time are key components upon we build on in the approach presented in this article for finding and validating ab initio symbolic expressions for GWs. Symbolic regression using genetic programming. Symbolic regression (SR) is the general procedure of finding closed-form expressions representing data. Unlike more conventional regression approaches, when SR is free-form it means precisely that: no specific form is specified. This eliminates any possible bias or human knowledge gap when postulating specific forms to fit for; it also contemplates completely data-driven cases, in which there is no underlying fundamental model or if there is one it is (still) unknown 28 .
Genetic programming (GP) is, in brief, an area of Artificial Intelligence (AI) whose goal is the evolution of programs or tasks through computer means. The techniques of GP emulate those of Nature; that is, algorithms are modeled following the process of natural evolution. A thorough book on GP is 29 , and a shorter field guide 30 .
Unlike perhaps more conventional Machine Learning (ML) deterministic regression approaches, GP-SR uses genetic algorithms to find free, closed-form expressions, either algebraic or differential. GP-SR can be described through the following general tree-structured algorithm tracing genetic programming principles: Create stochastically an initial population of programs (e.g. mathematical expressions and operations);

3.
Execute each program and compute their quality or fitness;

4.
Select one or two programs from the population with a probability based on their fitness to participate in genetic operations;

5.
Create new programs through the application of genetic operations (e.g. mutation or crossover);

6.
Until an acceptable solution is found or some other stopping condition is met;

7.
Return the best-so-far individual/s. GP-SR algorithms do not find a unique representation of data but a number of them with different levels of complexity (roughly speaking and for simplicity, cost of evaluation) and fitness with respect to training and validation sets. So, depending on the criteria used to find expressions via SR, the final symbolic forms can be shorter or larger with variable accuracy. In this work we prioritize accuracy. We used Eureqa 31 for GP-SR, although there are open source alternatives such as gplearn 32 and Glyph 33 . As a fitness criteria we used the Hybrid Correlation and the R 2 -Goodness-of-fit metric errors to fit symbolic models for amplitude and phase, respectively. Our contribution is not actually on the area of GP-SR but on how to address its high computational cost by means of modern ROM and show that it actually works through the case study of the gravitational waves emitted by the collision of two black holes (a highly non-trivial system to model). For the details on the GP side of Eureqa see the original works of Lipson, Schmidt and Bongard 28,[34][35][36][37][38][39] . An example: a robot discovering Newton's second law. As an example of the power of GP-SR we present results for the following system: the simple pendulum. In polar coordinates, Newton's second law is where = g/l , g is the gravitational acceleration and l the length of the pendulum. The variable θ represents the angle with respect to the point of stable equilibrium (the pendulum at rest). We set = 0.5s −2 and performed GP-SR on a dataset with initial conditions θ(t = 0) = π/2 and θ (t = 0) = 0 . The time interval for the training set was set to be [0, 22] s with grid size t = 0.1 s, which roughly corresponds to 2 cycles of the pendulum. That is, a single stream of data was used for training. We solved the above ordinary differential equation (ODE) with the integrate.odeint solver from the Scipy Python library 40 , and used the resulting data (with intrinsic noise, due to the numerical errors of the ODE solver) to find symbolic expressions for the underlying differential equation, searching for expressions of the form The representation with highest fitness -found in the order of seconds-was exactly, not a numerical approximation of, Newton's second law for this system, Eq. (1). Furthermore, for initial conditions close to the stable equilibrium state, one of the symbolic expressions found was exactly the harmonic oscillator equation.
Although the following conclusion could be somewhat debatable, the point is that a robot could find Newton's fundamental second law in seconds. One could argue that this was for a particular system but, though not presented in these terms, this is the process of scientific induction. In data science (DS), ML or AI, this process is called validation, whereas in physics it is called verification (as in verifying Newton's or Einstein's theory of gravity). In fact, with more computational power, the authors of Eureqa remarkably rediscovered Newton's second law for the double pendulum 28 , which is known in physics as a classical example of a chaotic system.

Complexity: searches and validations.
Genetic programming symbolic regression algorithms are known for their scalability issues with the amount and dimensionality of data 41,42 . Although it is not unusual for deterministic ML algorithms to suffer from the curse of dimensionality, scaling in GP-SR is nevertheless an issue. For our domain application, as described later on in this article, it resulted in no signs of convergence whatsoever after O(10 3 − 10 4 ) core hours (total number of hours with all cores running in parallel). The problem was not so much the number of hours but the lack of any progress in the fitness.
Our approach to overcome this problem is intuitive and conceptually simple: it uses a set of sparse data for fast training and later high-accuracy surrogate models for large scale validations. Here we use the RB and EIM frameworks combined with the surrogate approach developed in Ref. 13 . If, as in our case, the surrogate models are essentially indistinguishable from supercomputer simulations of the Einstein equations, they can be considered as ground truth solutions, with the advantage of very fast evaluations. The steps of validation for building surrogate models based on RB and the EIM are described in 9,10,13,15 . So here we focus on the ones related to GP-SR. In this processes we used a fraction of our catalog for training and another one for validation so as to avoid overfitting; typically we used 50% for each (training and validation). We then compared the symbolic time series with the ground truth solutions using dense sampling in parameter and time. This validation instance was achieved by computing the overlap integral S(h 1 , h 2 )(q) between surrogate and symbolic normalized waveforms in the time domain, defined as where and h 1 stands for the complex conjugate of h 1 . The overlap S gives a measure of the match between two waveforms and is commonly used in GW science. For training issues we have normalized time by a factor of 1,000; this has to be taken account for.

Results
Gravitational waves setup. As a proof of concept we tackled the problem of two black holes initially in quasi-circular orbit and without spin, for about 25-31 GW cycles before merger. More precisely, for the time interval t ∈ [−2, 750 : 100]M , where M is the total mass of the system and the waveforms are aligned so that t = 0 corresponds to the peak of their amplitudes, which is around the time the two black holes merge. Due to the scale invariance of General Relativity, the only free parameter then is the mass ratio  10] , with m i the mass of each black hole. Furthermore, for definiteness we restrict our discussion to the dominant multipolar mode, the ℓ = m = 2 one. A surrogate for this problem was constructed in Ref. 9 , and is publicly available as part of the GWSurrogate Python package 43 . This surrogate model consists of only 22 basis elements and, by construction, only 22 EIM time nodes. These are the only pieces of information needed to predict with high accuracy any waveform in the parameter and time domains considered. The surrogate model can be considered -and we do-as ground truth solutions of the EE, since it was shown in 9 that it is essentially indistinguishable from NR simulations up to the numerical errors of the latter, performed by SpEC 44 , the most accurate code in the NR community to date for modeling GWs from sources without shocks (such as binary black holes).
The two polarizations of GWs can be encoded in a single complex parameterized function, where represents a tuple in parameter dimension; here, it corresponds to the mass ratio q defined in Eq. (3).
The waveforms for the collision of two black holes in initial quasi-circular orbit have an apparent complexity, but they are simply oscillatory functions with an increasing amplitude until the time of coalescence, followed by a damped exponentially decaying profile, the quasinormal modes of the final black hole. It is therefore convenient to consider the amplitude A(t, ) and phase φ(t, ) separately, find closed-form expressions for them, later reconstruct the symbolic waves and compare them with their ground truth counterparts for a large number of validation cases. For both amplitude and phase our symbolic expressions have an R 2 -goodness-of-fit of at least R 2 ∼ 0.999 with respect to the validation members of the catalog used in the symbolic regression searches. We discuss more thorough and large-scale validation results below.
For both phase and amplitude our dictionary is composed of the following functions and operations: Amplitude. In our experience, naively sampling both in parameter (mass ratio q) and physical dimension (time) resulted in days or weeks of no convergence while searching for symbolic expressions for the amplitude of the GWs. The reason for this is the need to resolve with high accuracy the region around the peak of the amplitude (see Fig. 1), for which we tried using a dense sampling in time, which led to substantial computation with no signs of convergence.
One could attempt to manually collocate time nodes where needed. This approach is not only tedious but also not guaranteed to work. Instead, here we resorted to subsampling in time using only the 22 EIM nodes, shown in Fig.1, and 90 equally spaced values in the mass ratio. The rationale for this approach is that the EIM time nodes are the only relevant ones for recovering the whole time series and thus the only representative ones; this intuition proved to be correct, as we discuss below. Using only the EIM nodes in a few minutes we were able to find the following closed-form expression for all q ∈ [1, 10] (we discuss validation using a dense set of time nodes below): where gauss(x) := e −x 2 , atan2(x, y) is the arctangent of two parameters and  A(t, q) ={a 1 gauss(atan2(t, a 2 − t))}/{a 3 + q − t − a 4 gauss(atan2(a 5 , q) − a 6 t)gauss(atan2(t, a 7 − t − a 8 t q))} − a 9 , As a mode of illustration, we show in Fig. 3 the error curve as a function of time for the amplitude. Although the complete range of convergence included several hours, the plateau of the curve is achieved in the order of minutes, dedicating most of the time to fine tuning of parameters for improving the accuracy of the model. There are several reasons for the formation of this plateau, for example the finiteness of the dictionary, which serves as a constraint for the search space of functions; and the penalization for large formulas (high complexity) in the GP algorithm, which prevents from finding extremely high complex functions with little or no gain in accuracy. One could say that an important reason for the plateau is the stagnation in local optima, but actually the algorithm softens this problem by implementing a protocol that allows to promote population diversity in the evolutionary search without impacting the fitness performance. For details, see 36,37 . Phase. Although we were able to find high accuracy symbolic expressions for the phase in the considered interval of q ∈ [1, 10] (see Fig. 4 for a symbolic model that is continuous in the whole interval), they resulted in large propagation errors in the reconstruction of the waveforms. The reason is different from phase accumulation errors in numerical relativity, since here we are dealing with global optimization errors and is simply the following: a change in phase φ → φ + δφ in (4) leads to an error in the waveforms of the form so |δh|/A = δφ . In order to get a relative error of 1% at least, we must have an order of 0.01 in the phase error δφ . For the whole q ∈ [1, 10] phase symbolic model, in the results here obtained δφ is of order 1 (with a relative error less than 10 −2 ), leading to large errors when reconstructing the waveforms.
A simple domain decomposition to solve this issue worked out for us: we subdivided the domain q ∈ [1, 10] into 9 equally spaced subdomains of the form a 1 = 1.37502533181183, a 2 = 0.0409895367586908, a 3 = 3.40043449934568, a 4 = 1.86434379599601, a 5 = 1.1446516014466, a 6 = 1.49686180948812, a 7 = 0.0250835926883564, a 8 = 0.108134472792241, a 9 = 0.00178301085458751. [2,3], . . . , [8,9], [9, 10],  www.nature.com/scientificreports/ We have not tried using a different number of domains, or of different relative sizes, since the point of this paper is to show how to build a sparse training set to avoid the high computational cost of symbolic regression through genetic programming, and that it works in practice in a highly non-trivial application.
For simplicity, we have not imposed boundary conditions between the subdomain boundaries, although this could in principle be done. Domain decompositions are standard when solving partial differential equations (PDEs). In fact, it is possible that in more complex scenarios or even in this case a more elaborate scheme such as an hp-greedy domain decomposition [45][46][47][48] (where the hp term is actually borrowed from domain decomposition and refinement in finite elements when solving partial PDEs) might be more efficient in terms of decreasing the number of subdomains and improve our results. As an analogy, it is well known that when solving PDEs through a domain decomposition numerical solutions are not continuous across domain boundaries at fixed resolution (either in space, time, or both); this is usually addressed through weak enforcement of the solution across boundaries, usually in the form of penalty terms For a detailed discussion on these topics in the context of Einstein's equations see, for example [49][50][51] . In this article we focus on showing that, as a proof of concept, closed-and free-form expressions can be obtained without resorting to physical approximations. It might be possible to obtain accurate symbolic expressions without a domain decomposition, for example by enriching the dictionary through physically inspired functions, but this is left for future work. Here we focus on the basic elements of eliminating the curse of dimensionality in symbolic regression using modern approaches to reduced order modeling for parametrized systems.
When searching for symbolic expressions for the phase we used 20 values in mass ratio and 285 time nodes for each domain, with all points equally spaced. Due to the simple structure of the phase (see. Fig. 4), subsampling was not necessary.
Our highest accuracy results are the following:  φ [4,5]  φ [7,8]  φ [8,9]  φ [9,10]  The result is that the overlap S = S(q) (2) in our approach gives values above 99% for all cases. The main reason that we could do this a posteriori dense validation is due to the fact that ground truth solutions using surrogate models can be evaluated very quickly. The results are shown in Fig. 5. One should not reach any conclusion from the dependence of the overlap S as a function of the mass ratio q, since these are representations, much as in domain decomposition approaches in NR (though usually in physical space, not in parameters). For example, we could have chosen to show results for symbolic expressions with a more uniform error distribution, though it is worth emphasizing that the differences in the figure are in the order of 0.1%, which is below the accuracy required by Laser Interferometer GW detectors. Put differently, although not perhaps known for non-numerical relativists, this non-smooth behavior across boundaries is always present (though in a different sense) in numerical solutions to the Einstein equations when using state of the art domain decomposition or adaptive mesh refinement: continuity is only imposed within a given numerical resolution 49 . There are always "numerical jumps" and at any fixed resolution the numerical solution to Eintein's equations is non-smooth; the idea is that these discontinuities are below an acceptable numerical error. The analogue in our approach is that across boundaries the symbolic phases are not smooth but they result in waveforms with overlaps with respect to the ground truth solutions which are below an acceptable error tolerance. As a remark, in GW science, this acceptable discrepancy is of 97%, while our results guarantee at least 99%.
As an example, in Fig. 6 we show the ground truth solution on top of its symbolic expression for h + , corresponding to the worst match in the validation space for the whole interval q ∈ [1, 10] . Results for h × are similar since both modes are related simply by a π/2 phase difference.

Discussion
In perspective, having high accuracy, closed-form (symbolic) expressions for the emitted gravitational waves as predicted by a theory as complex as Einstein's one of gravity, for a process as complex as the collision of two black holes, without any simplification in the theory (thus the ab initio-based emphasis), in a completely autonomous way, cannot be overemphasized. The standard procedure in GW phenomenological modeling is that one of "stitching" GWs from post-Newtonian (PN) expansions or Effective One Body approximations, ringdown, and a merger regime in between from numerical relativity in some way (there are many of them) and in this sense they do carry physical approximations and therefore are not ab initio-based in the sense here used. Similarly with hybrid models, where early in the inspiral stage some physical approximant is stitched to a numerical relativity solution. Here we have concentrated for training data on a model that is completely indistinguishable from high accuracy numerical relativity supercomputer simulations of the full Einstein equations. In this sense our results can be considered truly ab initio-based, since in the absence of exact solutions (except for those with high degrees of symmetry) high accuracy NR simulations of the Einstein equations are considered the true gold-standard.
Another point to emphasize is that, unlike most -if not all-phenomenological models, our closed-form expressions do not distinguish between inspiral, merger and ringdown, but model all regimes at once. One could have chosen to find symbolic expressions for the different regimes just mentioned, but in order to show the power of our approach we have chosen to model the whole inspiral-merger-ringdown case at once. The domain decomposition here presented is very simple. Having different models in parameter space is not unusual. In fact, it has allowed (among many other ingredients) to find new signals from public LIGO data 52 .
Our approach is one of the many trends in the gravitational wave science community to incorporate tools from DS, ML and AI, but to our knowledge it is the first of its kind. Because of this, and because genetic programming symbolic regression is meant to provide insights from data, it is difficult to anticipate the full impact and ramifications of our approach. It might be useful, for example, for other ones combining ROM with Deep Learning for GW inference 53 , which produce, and start from closed-form expressions.
We presented a proof of concept for a novel approach. A next natural, conceptually straightforward, step might be to apply it to the other multipole modes of 9 , and more complex systems such as the case of spinning, precessing black holes using, for example, the surrogate models of 11,12,54 . It is possible that for these cases, and higher dimensionality ones in parameter space in general, it would be beneficial when training GP-SR to use not only the EIM time nodes but also the greedy parameter values to increase sparseness and avoid the curse of dimensionality of SR searches. Another line of future research would be to use an hp-greedy refinement approach at the surrogate level to minimize the number of domains. The question of what is the optimal minimum number of subdomains and how it increases with parameter dimensionality is outside the scope of this work but remains as an interesting question. Making touch again with numerical relativity, the equivalent would be asking how the number of domains or levels of adaptive mesh refinement changes with resolution. Even in an established field such as NR, where there is little to no theory for equations as complicated as the Einstein ones, one usually proceeds through numerical experiments.
Even though here we have focused on symbolic expressions based on surrogates built from high accuracy numerical relativity simulations, our approach can be applied to other surrogates based on RB and the EIM, for example those based on EOB or PN ones 13,55 .
The sparse yet near-optimal subsampling in time using the EIM is a key ingredient in our approach for finding closed-form expressions, so it is not clear that other surrogate models (say based on Gaussian regression, see, e.g. Ref. 56 ) can do so. There might be potential in enriching the dictionary here used for GP-SR (composed of elementary functions and basic arithmetic operations) using phenomenological or other physically based symbolic models. In this sense, using GP-SR should outperform any other kind of physics-based fits by design being free-form. Also, any other fits can be used as a bootstrap to enrich the dictionary of GP-SR; this can also be done as more dimensions and physics complexity are added 28 .
Anyone willing to qualitatively reproduce or extend the results of this paper can apply for a free academic license of DataRobot 57 , an Automated Machine Learning framework, which integrates the GP algorithm used in this work. Our results can be easily accesible through a Jupyter notebook at https://github.com/aaronuv/ SymbolicGWs.
In general terms, our approach should be applicable to other disciplines beyond gravitational wave science since computational complexity is a common problem in genetic programming.