Sequential Monte Carlo with cross-validated neural networks for complexity of hyperbolic black hole solutions in 4D

This paper investigates the self-similar solutions of the Einstein-axion-dilaton configuration from type IIB string theory and the global SL(2,R) symmetry. We consider the Continuous Self Similarity (CSS), where the scale transformation is controlled by an SL(2, R) boost or hyperbolic translation. The solutions stay invariant under the combination of space-time dilation with internal SL(2,R) transformations. We develop a new formalism based on Sequential Monte Carlo (SMC) and artificial neural networks (NNs) to estimate the self-similar solutions to the equations of motion in the hyperbolic class in four dimensions. Due to the complex and highly nonlinear patterns, researchers typically have to use various constraints and numerical approximation methods to estimate the equations of motion; thus, they have to overlook the measurement errors in parameter estimation. Through a Bayesian framework, we incorporate measurement errors into our models to find the solutions to the hyperbolic equations of motion. It is well known that the hyperbolic class suffers from multiple solutions where the critical collapse functions have overlap domains for these solutions. To deal with this complexity, for the first time in literature on the axion-dilaton system, we propose the SMC approach to obtain the multi-modal posterior distributions. Through a probabilistic perspective, we confirm the deterministic α\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha $$\end{document} and β\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta $$\end{document} solutions available in the literature and determine all possible solutions that may occur due to measurement errors. We finally proposed the penalized Leave-One-Out Cross-validation (LOOCV) to combine the Bayesian NN-based estimates optimally. The approach enables us to determine the optimum weights while dealing with the co-linearity issue in the NN-based estimates and better predict the critical functions corresponding to multiple solutions of the equations of motion.


Introduction
It is well known that all Black holes can be characterised by their mass, their charge as well as their angular momentum.Choptuik in [1] also showed that there is yet one more parameter that explains the critical gravitational collapse solutions and it is called the critical exponent.More specifically, Christodolou in [2] first had revealed the spherically symmetric collapse of the real scalar field.Later on, Choptuik [1] numerically showed that the real scalar field gravitational collapse solution demonstrates the discrete selfsimilarity property.Indeed, the gravitational solution shows space-time self-similarity where the dilations can take place.Therefore, the critical solution does provide an scaling law.If we show the initial condition of the real scalar field by parameter p, which is called the field amplitude, then p = p crit defines the critical solution and hence, the black hole can be formed once p chooses the values bigger than p crit .Indeed, for p > p crit the mass of the black hole or the Schwarzschild radius are given by a scaling law as follows r S (p) ∝ M bh (p) ∝ (p − p crit ) γ . ( The following articles [1,3,4] found that the critical exponent for a real scalar field is given by γ ≃ 0.37 in four dimensions.Notice that for dimensions bigger than four (d ≥ 4), the black hole's mass scaled by [5,6] as One may read some other numerical investigations for several other matter content in [7,8,9,10,11,12].
The collapse solutions of the perfect fluid had been studied in [13,14,5,15] and its critical exponent γ ≃ 0.36 was also found in [14].The authors in [16] argued that γ might have a universal value for all matter fields that can be coupled to gravity in four dimensions.As discussed in [5,15,17] the critical exponent is explored by using the perturbations of self-similar solutions.
The solutions with axial symmetry were investigated in [18], while shock waves are studied in [19].The authors in [20] explored the value of the critical exponent γ ≃ 0.2641 for the axion-dilaton configuration in four dimensions.Interestingly, in [21] we have studied the perturbations and were able to precisely generate the existing value [20] of γ ∼ 0.2641 in four dimensions and other critical exponents have been derived in [22] in four and five dimensions.
Albeit the authors in [23,3] investigated the entire analysis for the elliptic case in four dimensions, it is worth mentioning that their methods can also be examined in hyperbolic and parabolic cases as well as other dimensions, where for further results we refer to [24,22].
Let us provide various motivation for the study of the critical collapse of the Einstein-axion-dilaton configuration.The first motivation is related to the gauge-gravity duality [25,26,27,28], corresponding the choptuik exponent, the imaginary part of quasi-normal modes as well as the dual conformal field theory that is pointed out in [29].In fact, of the interest is to study the spaces that approach asymptotically to AdS 5 × S 5 , where the simplest system is to consider in type IIB string theory of the axion-dilaton system with the self-dual 5-form field.Then, one can start to analyze the black hole solutions in diverse dimensions.
It is also important to highlight that this system has also been related to the holographic description of black hole formation, see [30,13].Finally the implications of this system to black hole physics have been carried out in [31,32].The key role of S-duality for these self similar solutions has also been considered in [33].In this paper, we are dealing with the collapse of matter to form small mass black holes and hence one considers a small space-time region just close to where the singularity occurs.It has also been shown that this event is independent of the asymptotic structure of the space-time to which the collapse happens.
In fact there is already the numerical evidence in asymptotically AdS space-times confirming that this is the case [8].Therefore we eliminate the cosmological constant and just analyse self-similar solutions for the axion-dilaton system.
The self-similar solutions for all elliptic, hyperbolic and parabolic classes of SL(2,R) have been discovered in [34] in four and five dimensions for all classes of SL(2,R), which are the extensions of the earlier results [35,36].In [37] we have also recently made use of the Fourier-based regression models for obtaining the critical solutions.Consequently, the challenges of [37] have been addressed in [38] where we applied truncated power basis, natural spline and penalized B-spline regression models accordingly in order to be able to explore the non-linear functions.In [39] we applied artificial neural networks in order to address the instability of the black hole solutions for the specific parabolic class in higher dimensions.Lastly, in [40], we actually proposed a new formalism to be able to model the complexity of elliptic black hole solution in four dimension using hamiltonian monte carlo with stacked neural networks.
In this paper, we propose a new formalism based on Sequential Monte Carlo (SMC) and artificial neural networks (NNs) to be able to model the hyperbolic class of the spherical gravitational self-similar solutions in four dimensions.Due to the nature of highly non linear equations of motions of hyperbolic black holes, various authors used a variety of numerical calculations to simplify the equations of motions and parameters of the theory, for instance one can see [24,38,39].Hence, due to this reason the authors must have overlooked the measurement errors that are imposed in exploring the parameters through various numerical methods.
Thus here we propose a new method to carry out the measurement errors, involved in parameter estimation, into our statistical models in exploring the solutions to the equations of motion.Recently Hatefi et al. [40] applied the Hamiltonian Monte Carlo method to find solutions to the equations of motion in the elliptic class of four dimensions in a Bayesian framework.Unlike [40], the hyperbolic equations of motion in four dimensions have multiple solutions and therefore the collapse functions do have overlap domains under these solutions.In order to deal with this challenge, for the first time in the literature on the axiondilaton system, we proposed the SMC approach to derive the posterior distribution of the parameters.The posterior distribution does provide all possible solutions in estimating the parameter of the equations of motion.Interestingly, the posterior distribution confirms the deterministic α and β solutions found in the literature for the hyperbolic class in four dimensions.Unlike other methods in the literature, in this paper, we impose the l 2 penalized Leave-One-Out Cross-validation (LOOCV) to optimally combine the Bayesian NNs candidates.The advantage of this approach is that it also enables us to determine the optimum weights while dealing with the co-linearity issue in the NN-based estimates and better predict the critical functions corresponding to multiple solutions of the equations of motion in the hyperbolic class.
The organization of the paper is as follows.In section 2 we briefly explain the relevant effective action for the axion-dilaton configuration, its equations of motion as well as the initial conditions that come from the continuous self-similarity requirement.We the describe our methodology to actually model the complexity of Sequential Monte Carlo with Cross-validated Neural Networks for hyperbolic black hole solutions in four dimensions.In section 4, we use SMC samples from the posterior distribution and construct NN estimates based on the posterior mean and LOOCV as well as the 95% credible intervals in estimating the critical collapse functions corresponding to multiple solutions of the equations of motion.Lastly, we present the results and conclude in the section 5.
2 The Einstein-Axion-Dilaton System and Its Equations of Motion for Hyperbolic Class The two real scalar fields of axion and dilaton can be combined to construct a single complex scalar field τ ≡ a + ie −ϕ .Its dynamics and coupling to the gravity or its effective action for four-dimensional axiondilaton (a, Φ) system is described by where R is the scalar curvature.If we take variations from the metric and τ then one would be able to find out all the equations of motion as follows This effective action is classically invariant under SL(2,R) transformations which means that if τ gerts replaced by where (a, b, c, d) ∈ R, ad − bc = 1 then g ab and the action remains invariant.As argued originally by [41,42,43,44]) this group gets broken to SL(2,Z) as an indication of duality transformation.
The spherically symmetric metric is represented by [23] If we take into account the time scaling for (7), (as shown in [23]), one can set b(t, 0) = 1 and regularity condition indicates u(t, 0) = 0.The so called continuous self-similarity (CSS) means there exist a killing vector ξ that generates global scale transformation, where in spherical coordinates, we define ξ = t ∂/∂t + r ∂/∂r.
The assumption of continuous scale invariance for the metric gets related a scaling for the line element under dilations as follows (t, r) → (Λt, Λr) , Λ > 0 (8) then Now if we consider the scale invariant variable z = −r/t then the self-similarity of the metric implies that the functions u(t, r), b(t, r) must be expressed in terms of z, that is The scalar τ must also be invariant up to an SL(2, R) transformation, so that Hence a system of (g, τ ) that satisfies eqts.(10), (11) to be continuously self-similar (CSS).Hence, physically distinct cases are related to the different conjugacy classes of dM dΛ Λ=1 .
The effective action in (3) is SL(2,R)-invariant, hence we can consider a compensation of the scale transformation of (t, r) by an SL(2,R) transformation.In fact in [35] we already found out three different possible assumptions for this particular system.Those were called the elliptic, hyperbolic and parabolic classes that are related to three classes of SL(2, R) transformations which were used to compensate for a scaling transformation in space-time.Let us describe the hyperbolic ansätze for τ (t, z).
The general form of the ansatz for the hyperbolic class is given by where under a scaling transformation t → λ t, τ (t, r) changes by a SL(2, R) boost or hyperbolic translation, which means that all equations are invariant under the following transformation Note that under SL(2, R)-transformation the following exactly produces the same equations of motion for hyperbolic case, where f (z) is a complex function satisfying Im f (z) > 0, and ω is a real constant.Let us describe the derivation of the equations of motion for the hyperbolic class in four dimensions.If we apply continuous self-similarity ansätze (12) to all the equations of motion ( 4) and ( 5) then one would be able to explore the ordinary differential equations for u(z), b(z), f (z).However, if we make use of the spherical symmetry then one reveals that u(z) and its first derivation u ′ (z) can be expressed in all the equations in terms of b(z), f (z) and their first derivatives as follows All the ordinary differential equations (ODEs) are given by Finally, the equations of motion in the hyperbolic class in four dimensions are represented by The equation of motion for b is a first-order linear in-homogeneous equation with initial condition b(0) = 1.
The initial conditions for f (z), f ′ (z) are also realised by applying the smoothness of the critical solution.
From (20), we encounter five singularities at z = ±0, z = ∞ and z = z ± .One can readily show that the singularities z = ±0, z = ∞ can be eliminated by the coordinate transformation as shown in [36].The last two singularities are demonstrated by b(z ± ) = ±z ± .They correspond with the horizon where z = z + is just a coordinate singularity as shown in [20,35].Therefore by definition, τ must also be regular across it.
Hence f ′′ (z) must also be finite as z → z + .Therefore, one finds that the vanishing of the divergent part of f ′′ (z) produces a complex-valued constraint at z + , that is indicated by G(b(z + ), f (z + ), f ′ (z + )) = 0 where the explicit form of the G function for the hyperbolic case in four dimension is given in One finds it convenient to make use of change of variables as f (z) = u(z) + iv(z).Indeed, using regularity at z = 0 and some residual symmetries one gets as well as These equations are invariant under a constant scaling f → λ f , therefore one has the freedom to choose either the real value of f (z) or its imaginary part as one wishes at a particular value of z, so we would like to set u(0) = 1.If one require the regularity at the origin z = 0, then one explores the following initial conditions for the hyperbolic class as Therefore, the real and imaginary parts of G must vanish which determines ω, where (x 0 > 0) and x 0 is a real parameter.The three discrete solutions in four dimensions were explored in [21] where these solutions are found by integrating numerically the equations of motion.The solutions in hyperbolic class are identified by seeing the very rapidly decreasing Im f (0).No solutions are explored with Im f (0) > 1.The α and β solutions are evidently explored.The α solution is given by The β solution is However for the γ solution we notice that due to the fact that Im f (0) is so small also the z + root-finding gets effected with numerical noise ( ω = 0.541, v(0) = 0.0059, z + = 8.44), and hence the quality is not perfect because the G ∼ 10 −7 comparing to G ∼ 10 −13 − 10 −17 for the first two solutions, hence we concentrate on α, β solutions.
It is also important to highlight the fact that in [22] we have shown that the Choptuik exponent depends on the dimensions, matter content as well as the different branches of the unperturbed self similar solutions.
Thus, we argue that the conjecture about the universality of Choptuik exponent is not satisfied.However, we also claim that there may exist some other universal behaviours that could have been hidden in combinations of critical exponents or there might be some other parameters of the given theory that have not been considered yet in our understanding of current investigations.

Statistical Methods
In this section, we investigate Bayesian solutions to the parameters of the hyperbolic equations of motion using sequential Monte Carlo methodology.Suppose x(t) = (x 1 (t), . . ., x H (t)) represent the solutions to where the system encompasses H differential equations (DEs) in which t denotes the space-time, θ encodes the set of all unknown parameters of the system and x i (t) denotes the solution to i-th DE, i = 1, . . ., H.
Henceforth, we call x(t) the DE variables.
As the equations of motion in black holes are highly nonlinear, typically there is no closed-form solution for the true trajectory of the DE variables.Therefore, researchers often apply a series of numerical methods and observe the trajectory of the DE subject to measurement errors.Let y ij denote the observed trajectory of the i-th DE variable at j = 1, . . ., n i space-time points.To take into account the uncertainty involved in the observed DE variables, let y ij follow a Gaussian distribution with mean x i (t j |θ) and standard deviation σ i for i = 1, . . ., H. Let Ω = (θ, σ) denote the set of all unknown parameters of the model where σ = (σ 1 , . . ., σ H ).
Combing all information from the observed DEs, the likelihood function of Ω is given by In real-world scenarios, while the DE systems depend on unknown parameters, researchers typically have prior knowledge about the parameters of the DE model.To incorporate the prior information into our estimation, we propose a Bayesian framework and treat the unknown parameters of the DE system as random variables.This Bayesian estimation approach enables us to find the statistical distribution of the unknown parameters based on the set of observed data, given prior information about the unknown parameters.The statistical distribution of parameters given the observed data is henceforth called the posterior distribution of the parameters π(Ω|y) where Ω denotes the set of all unknown parameters of the DE system.The posterior distribution then enables us to make statistical inferences about the uncertainty in the estimation procedure and quantify the characteristics of the DE system.

Sequential Monte Carlo
When we face DE system (25), it is reasonable to take into account the effect of space-time argument and the sequential process of the DE observations y 1 , . . ., y t .This facilitates updating the posterior distribution of the unknown parameters sequentially as the DE variables are sequentially observed in the system.Sequential Monte Carlo (SMC), as a simulation-based approach, is a flexible technique in Bayesian statistics to sequentially estimate the posterior distribution.Due to the advent of cheap and powerful computing resources, the SMC appears a convenient tool to implement, in a parallel fashion, sequential computation of the posterior distribution in a general setting.
Let {Ω t , t ∈ N} denote the Markov process showing the trajectory of unknown parameters over space-time with the prior distribution π(Ω 0 ).The Markov property of the process indicates that the probability of process at space-time t only depends on the previous step of the state; that is where p(Ω t |Ω t−1 ) denotes the transition probability from Ω t−1 to Ω t in the parameter space.Let {y t ; t ∈ N} denote the sequence of the observations from the DE system where they are conditionally independent given the observed status of the parameter process with distribution p(y t |Ω t ) for t ≥ 1.For the sake of convenience in notations, let Ω 0:t and y 1:t represent the parameter and DE observation sequences up to space-time t, respectively; That is, Ω 0:t = (Ω 0 , . . ., Ω t ) and y 1:t = (y 1 , . . ., y t ).
From the joint distribution of Ω 0:t and y 1:t , one can also recursively update the posterior distribution π(Ω 0:t |y 1:t ) based on the posterior distribution of previous lags of the process by Accordingly, one can compute the posterior expectation of any characteristic of the DE system by Due to the complex structure of the uncertainty and multidimensional integration of the marginal distribution in the DE system, there is no analytical form for the posterior distributions ( 27) and (28); thus the conditional expectation ( 29) is not tractable too.
Importance sampling [45,46], as a practical solution to the intractability problem, uses an instrumental distribution to sample indirectly from the posterior distribution π(Ω 0:t |y 1:t ).Let q(Ω 0:t |y 1:t ) represent an instrumental distribution whose domain includes the domain of the target posterior distribution (27).
Let {Ω (i) 0:t }; i = 1, . . ., n represent n independent and identically distributed particles from the instrumental distribution q(Ω 0:t |y 1:t ).Thus, using the importance sampling method, the Monte Carlo estimate of the quantity of interest H t (•) is given by where Λ(Ω (i) 0:t ), as normalized importance weights, are obtained by The importance sampling estimator (32), as a general framework of Monte Carlo, is convenient to implement.
Despite this convenience, the iterative structure of the technique is not adequate to sequentially incorporate the new status of the process into estimating the π(Ω 0:t |y 1:t ).As the new status of the sequence becomes available, one has to recompute all the important weights on the entire parameter space.This becomes computationally expensive for complex and nonlinear equations of motion.
Sequential Monte Carlo (SMC) [46,47], as a sequential architecture of importance sampling, can be considered as a solution to the iterative problem of general importance sampling.As observed from ( 31)-( 33), in importance sampling, when a new status of the DE variable becomes available y t , one has to re-compute the posterior distribution of the entire trajectory Ω 1:t given the sequence y 1:t .Hence, one requires re-computing even the importance weights of the previous states of the trajectory given the new status of the DE variable.
Unlike the importance sampling method, the SMC sampler does not require re-computing the importance weights corresponding to the previous states of the trajectory (Ω 0 , . . ., Ω t−1 ), when the new status of the DE variable becomes available.In other words, the importance weights of the previous states of the trajectory stay the same, and we no longer need to re-compute them.We only need to compute the posterior distribution of the most recent state of the trajectory, given the new data.Treating the instrumental distribution based on the previous states as the full marginal distribution, the instrumental distribution is thus updated by q(Ω 0:t |y 1:t ) = q(Ω 0:t−1 |y 1:t−1 )q(Ω t |Ω 0:t−1 , y 1:t ).
In a similar fashion as (34), one can recursively show that From the sequential representation (35), one can easily show that the non-normalized importance weights (33) can be sequentially updated by and consequently, the normalized importance weights are sequentially updated by As a special case of (35), one can employ the prior distribution in the SMC framework [46,47].In this case, the instrumental distribution is given by In this case, from (38) and the fact that p(y t |y 1:t−1 ) is constant over the entire trajectory of Ω 0:t , one can easily update the importance weights of the i-th particle by Λ(Ω Although the SMC technique is well suited to accommodate sequentially the new data into estimation, the posterior distribution of the particles very quickly becomes skewed after only a few steps such that only a few particles will have a non-zero probability [47,48].Consequently, the Monte Carlo chain will not be able to sample from all aspects of target π(Ω 0:t |y 1:t ).To deal with this degeneracy issue, a re-sampling is added to SMC to eliminate sequentially particles with low importance weights and at the same time augment the particles from high-density areas.To implement the re-sampling step, one can take n random draws with replacement from the collection of the particles {Ω 0:t } with probabilities corresponding to the important weights {Λ(Ω t denote the number of offsprings from the particle Ω (i) 0:t , i = 1, . . ., n.It is easy to see that the particle survives and contributes to the posterior distribution when n (i) t > 0; otherwise, the particle dies.The important sampling and the re-sampling steps are finally alternated to update the important weights and filter the particles sequentially to obtain the samples from the posterior distribution.

Cross-validated Neural Networks
Neural networks (NNs), inspired by the human neural system, comprise a network of connected neurons.
These connections enable the neurons to send information from one layer to another.According to the flexibility and power of the NNs, they have been increasingly exploited in recent years as a reliable predictive model for solving differential equations.The power of the NNs enables us to reformulate finding solutions to the DE system to a parametric estimation of the DE variables by minimizing the prediction errors.NNs consist of a multi-layer perceptron whose layers include neurons connecting the layers of the network to each other.These connections enable the NN to estimate the functional form of the DE variables.In this subsection, we describe how NNs use the parameter estimates developed by the SMC method, from Subsection 3.1, to predict the complex and nonlinear forms of the DE variables in the hyperbolic class of 4d.
Let N (x(t|Ω), t, ϕ) denote the NN estimate, consisting of L hidden layers, for the DE variable x(t) form DE system (25).Each neuron of the NN is connected with another in the next layer via a linear regression model such that NN l−1 (•) represents the response observed from the l-th layer, ϕ = (W 1 , . . ., represents the set of all unknown parameters, a denotes a non-linear activation function, W l and b l show the weight matrix and bias vector of the l-th layer, respectively [49,50].Finally, the NN estimate of the DE variables x(t) are obtained as a solution to the squared loss function To handle the optimization (41), the NNs implement a series of forward and backward propagation steps to find the final solution to the DE system (25).For more details about the theory and applications of the NNs, readers are referred to [49,51,50] and references therein.
In this paper, we first plan to find the Bayesian estimate for the parameters of the equations of motion.The Bayesian proposals are then stacked into the equations of motion to find the Bayesian stacked NN solvers.
To this end, as described in Subsection 3.1, we use the SMC method and find the posterior distribution π(Ω|y).Let Ω * 1 , . . ., Ω * M denotes M Bayesian SMC proposals for the parameters of the DE system.Given the posterior candidate Ω * m , let N (x(t|Ω * m ), t, ϕ) denote the NN-based estimate of the DE variable x(t) for m = 1, . . ., M .From a probabilistic perspective, the estimate N (x(t|Ω * m ), t, ϕ) will be the true solution to the DE system (25) with probability π(Ω * m |y) for m = 1, . . ., M .
Recently, Hatefi et al [40] proposed the Bayesian model averaging to stack the NN-based estimates in predicting the critical solution for the elliptic class of 4d.Following [40], using the SMC candidates Ω * 1 , . . ., Ω * M and the training set y of size n, the posterior distribution of the N (x(t|Ω * m ), t, ϕ) at a fixed space-time t, is given by Using ( 42 corresponding to SMC candidates.As N 1 (x(t|Ω * 1 ), t, ϕ), . . ., N M (x(t|Ω * M ), t, ϕ) are M estimates of the DE variable x(t) at space-time t, the estimates may be linearly dependent.This arises the co-linearity problem in the NN-based design matrix N (x|Ω * ).To this end, under squared error loss with l 2 penalty, we can stack the NN-based estimates using the penalized linear regression model [49].One can estimate the coefficients of the regression model β = (β 1 , . . ., β M ) by According to the properties of the least square under l 2 penalty, one can easily show that the solution to ( 43) is given by When the posterior distribution has multiple modes, the least square estimate (44) may assign unfair weights to some complex NN candidates based on the training set.To deal with this problem, we develop the LOOCV estimate of the weights where leaving one observation out in each iterative training step to find the best coefficient estimates.Let N −i m (x(t|Ω * m ), t, ϕ) denote the NN-based estimate for the DE variables at spacetime t using the m-th SMC proposal when the i-th observation in the training set has been removed.From (43), the LOOCV estimate of the weights are given by From ( 45), the LOOCV-based NN estimate of the DE variables at space-time t, using the SMC estimates Ω * 1 , . . ., Ω * 1 is given by The LOOCV-based NN estimates (46)

Numerical Studies
The equations of motion for hyperbolic class in four dimensions have five singular points; however, the relevant range for z, that contains the proper information lies between the two singularities In particular, z = z + is an event horizon, which is the homothetic horizon.Thus, it is a coordinate singularity, and τ must be regular across it, which is equivalent to the finiteness of f ′′ (z) as z → z + .In fact the vanishing of the divergent part of f ′′ (z) leads to a complex-valued constraint at z + as From time scaling, the regularity of τ , and the residual symmetries in the equations of motions, one can show the initial boundary conditions as where x 0 is a real parameter as the initial value of the equations of motion.Thus, the system results in one parameter ω and two constraints which are the real and imaginary parts of G. Therefore, the system is handled by discrete solutions where the CSS solutions are numerically investigated.
Here we plan to estimate the self-similar solutions using the Bayesian framework for the hyperbolic class in 4d.In this framework, we treat the parameter of the equations of motion ω as a random variable and then use the SMC approach to find the posterior distribution.The posterior distribution allows us to take into account the numerical measurement errors in estimating the critical collapse functions of the the system.
The equations of motion of the axion-dilaton system have already been studied in various dimensions and also for different ansatz [39,21].Hatefi et al. [38,37] applied the statistical regresion models using Fourierbased and spline smothers to estimate the critical collapse functions.In a recent publication, Hatefi et al. [39] constructed a solver based on NNs and showed there was no solution in higher dimensions for the parabolic class of black holes.Particularly [21] employed a root-finding method to numerically determine all the parameters of the equations of motion.
Due to the highly nonlinear equations of motion, researchers typically have to use a series of numerical approaches to simplify the equations and also to keep track of the parameters of the models.For instance, these techniques include imposing non-trivial constraints on the equations of motions such as finiteness of f ′′ (z) as z → z + .One might also point out one more constraint, namely, the vanishing of the divergent part of f ′′ (z) that generates the complex-valued constraint at z + .This implies that the real and imaginary parts of G must vanish at z → z + .The equations in the hyperbolic case are not solvable analytically, hence [34] used the profile-root finding method and applied the discrete optimization on the coordinates of the extended parameter space of the equations.More specifically, the equations are approximated by the first two orders of the Taylor expansions to make the root-finding step for the parameters of the equations tractable.Next, [34] could clearly estimate the parameters of the equations of motion by applying a grid search optimization on the coordinates in the system.After discarding the spurious roots in the domain, [34] found out the solution to the equations of motion and then could estimate the critical collapse functions.
Although [39,38,37,21] explored the hyperbolic class of equations motion in deterministic approach, in this research through a stochastic perspective, we construct a Bayesian method using SMC approach to assess the measurement errors into the estimation of the critical functions.Recently Hatefi et al. [40] proposed Bayesian mean using Hamiltonian Monte Carlo to find solutions to the equations of motion in the elliptic class of 4d.Unlike [40], the hyperbolic equations of motion in 4d result in multiple solutions where the critical collapse functions have overlap domains under the multiple solutions.To deal with the complexity of the measurement errors in the hyperbolic equations of motion, for the first time in the literature on the axiondilaton system, we proposed the Sequential Monte Carlo approach to derive the posterior distributions of the parameters.Unlike [40] where they proposed Bayesian model averaging to stack the NN-based estimates, here we propose the l 2 penalized Leave-One-Out cross-validation to stack the NN-based estimates for critical collapse functions.The approach enables us to assign the optimum weights to NN-based estimates and hence better estimate the critical functions corresponding to multiple solutions of the equations of motion.
The SMC Bayesian estimates provide information for all possible parameters' outcomes, that may take place in the numerical experiments, from the posterior distribution.They also allow researchers to embed this complexity in estimating critical functions.In this numerical study, we consider equations of motion (20) as the DE system of interest.The system leads to three critical collapse functions b 0 (z), |f (z)| and arg(f (z)) where we treated these functions as the DE variables of the system that must be estimated.Since all the DE variables of system (20) are numerically and simultaneously solved, it makes sense to assume that the observed DE variables also have the same standard deviation σ parameter which actually represents the variability in the numerical experiments.Hence, Ω = (ω, σ) encodes the set of all unknown parameters of the model.We used Python Package Pymc3 [52] to implement the SMC approach and find the posterior distribution of Ω.In order to investigate the effect of the prior information on the likelihood parameters in predicting of our critical functions, we assigned two different prior distributions for ω.The prior distributions include non-informative uniform distribution between [0.3, 1.5].We assign the second prior distribution to be the Gaussian distribution with a mean of 1.20 and a standard deviation of 0.2.We also consider the half-Cauchy distribution with scale parameter 0.5 as the prior distribution for parameter σ to capture the uncertainty involved in the likelihood function.
Figure 1: The posterior distributions and their trace-plots of two SMC chains (shown by solid and dotted blue lines) for parameters ω (w) and σ (sigma) under Gaussian prior distribution.
We show the posterior distributions of parameters ω and σ in Figures 1 and 6 under Gaussian and Uniform prior distributions, respectively.In each figure, we show two independent realizations of the SMC chain for the posterior distributions and their corresponding trace plots.It is clear that the posterior distributions of ω are multi-modal and have multiple high-density areas.Interestingly, this finding is compatible with the literature where [34] shows there must be at least three solutions to hyperbolic equations of motion in 4d.Under both prior distributions, we observe that there is a dominant high-density area almost ranging between [1.25, 1.4].This corresponds to the α-solution of [34].In addition, we also observe that both figures confirm there is a small high density on the left tail of the posterior distribution of ω which is interestingly compatible with the β-solution to the equations of motion in [34].It should be noted that for example under Uniform prior, the probability that the deterministic α-solution of [34] be the true solution to the hyperbolic equations motion is almost 10%; that is π(ω ∈ 1.36 ± 10 −2 ) ∝ 0.10.Moreover, the probability that the true solution appears to be smaller than 1.10 is given by π(ω ≤ 1.10) ∝ 0.05.
As described in Subsection 3.2, we construct the NN-based estimates of the critical collapse functions using the SMC proposals from the posterior distribution π(ω|y).To do that, we take L = 200 samples ω * 1 , . . ., ω * L from π(ω|y).Treating the posterior candidates ω * l , l = 1, . . ., L as the true value of the parameter in the equations of motion, we applied fully connected NNs to solve the DE system and find the NN-based estimates of the critical collapse functions.To do so, we used Python Package NeuroDiffEq [53] to carry out the neural networks for differential equations with 4 hidden layers where each layer consists of 16 neurons.We then ran the NNs for 1000 epochs and estimates the critical collapse function at 1000 equally spaced space-time points  In the next step of the numerical study, we plan to stack the L = 200 NN-based candidates in estimating the critical collapse functions.Using a root-finding based on general relativity, [34] showed that there are three solutions in the domain of the equations of motion in the hyperbolic class of 4d.These solutions correspond to ω 1 = 1.362 (α-solution), ω 2 = 1.003 (β-solution) and ω 3 = 0.005 (γ-solution).It was discussed in [34] that the γ-solution is not a stable solution due to the fact that Im f (0) is so small also the z + root-finding gets affected by numerical noise ( ω = 0.541, v(0) = 0.0059, z + = 8.44), and hence the quality is not perfect and so it may be a spurious solution due to numerical noise.For this reason, in this research, we focused on the two α and β solutions as two available solutions to the equations of motion in the literature.The critical collapse functions corresponding to α-solution range between [0, 1.44], while the functions range between [0,   Following Hatefi et al [40], we first applied the Bayesian model averaging method and computed the posterior mean of the L NN-based estimates at each space-time point z i , i = 1, . . ., 1000.Henceforth this estimate is called the mean stacked NN-base estimate of the critical collapse functions.We also computed the NN-based estimates for critical collapse functions under ω = 1.362 and ω = 1.003 corresponding to estimates under α and β solutions.We treated these fixed NN-based estimates as the true forms of the critical collapse function under two different solutions to the equations of motion.We then applied the LOOCV technique to stack the NN candidates in estimating the critical functions when we used 100 random observations from the true NN-based estimates as the taring sets for the LOOCV-based NN methods.On the other side, we also obtained the 95% Bayesian credible intervals for the critical collapse functions.To do that, we computed the 2.5 and 97.5 percentiles of the NN-based estimates, respectively, as the lower and upper bounds of the interval at each space-time point.Finally, it is important to highlight that if instead of taking continuous self-similarity, one only assumes discrete self-similarity, then the discrete scale transformation is compensated just by an element of SL(2, Z) in the set SL(2, R) transformation.It would also be very interesting to discover the fact that how all the critical exponents would depend on the modular transformations as well.

Conclusions
This paper proposes a new formalism based on Sequential Monte Carlo (SMC) and artificial neural networks (NNs) to model the hyperbolic class of the spherical gravitational self-similar solutions in four dimensions.
Due to the nature of highly non-linear ordinary differential equations for the axion-dilaton configurations, in the literature, researchers typically have to employ various constraints as well as different numerical approximation methods to keep track of equations.For instance, in hyperbolic equations of motion, [34] had to apply the constraints, including the finiteness of f ′′ (z) as z → z + and the vanishing of the divergent part of f ′′ (z) which generates a complex-valued constraint at z + .They also employed numerical grid search discrete optimization methods on the extended coordinates of the equations to eliminate the spurious roots and estimate the equations' parameters and the self-similar black hole solutions.Due to the sophisticated form of the equations and the vital role of the parameters, researchers usually have to overlook the measurement errors imposed in exploring the parameters through various numerical methods.
Here, we propose a new method to incorporate the measurement errors, involved in parameter estimation, into our statistical models in exploring the solutions to the equations of motion.Recently Hatefi et al. [40] applied the Hamiltonian Monte Carlo method to find solutions to the equations of motion in the elliptic class of four dimension in a Bayesian framework.Unlike [40], the hyperbolic equations of motion in four dimensions suffer from multiple solutions where the critical collapse functions have overlap domains under these solutions.To deal with this challenge, for the first time in the literature on the axion-dilaton system, we proposed the SMC approach to derive the posterior distribution of the parameters.The posterior distribution reveals all the possible solutions in estimating the parameter of the equations of motion.It is also important to highlight that the posterior distribution confirms the deterministic α and β solutions found in the literature for the hyperbolic class in four dimensions.Unlike methods in the literature, in this paper, we proposed the l 2 penalized Leave-One-Out Cross-validation (LOOCV) to optimally combine the Bayesian NNs candidates.
The approach enables us to determine the optimum weights while dealing with the co-linearity issue in the NN-based estimates and better predict the critical functions corresponding to multiple solutions of the hyperbolic equations of motion.
Using SMC samples from the posterior distribution, we then developed NN estimates based on the posterior mean and LOOCV as well as the 95% credible intervals in estimating the critical collapse functions corresponding to multiple solutions of the equations of motion.Because the posterior mean NN-based estimates and 95% credible intervals aggregate all the posterior NN candidates regardless of the multiple solutions of the system, these two methods remain robust and proposed the same estimates in predicting the critical collapse functions under both α-and β-solutions.Unlike the posterior mean proposal, we recommend the LOOCV-based estimates if one is interested in the optimum linear combination of the posterior NN candidates to improve the estimating of the critical collapse functions corresponding to established α-and β-solutins in hyperbloic equations of motions in four dimensions.Now if we compare the Bayesian method with the realistic NN approach, then one gets to know that the developed Bayesian credible intervals actually contain the definite estimate as one possible candidate in the estimation of the critical collapse functions.Indeed, Unlike the estimation of [39], the Bayesian approach remains concrete against measurement errors in estimating ω due to the fact that all these Bayesian estimations have already had all the possibilities of the parameter for the domain of the posterior distribution.
From a physical point of view, our results clarify that the universality of the Choptuik phenomena [15] is not satisfied.There may exist some universal behaviour which might be hidden in combining the critical exponents and other parameters of the theory.Nevertheless, our efforts provide some clear evidence that one cannot expect to transfer the standard expectations of Statistical Mechanics to the critical gravitational phenomena.
, . . ., Ω ), Hatefi et al[40] proposed the mean posterior of the NN-based estimates to stack the NN candidates in predicting the DE variables in the elliptic class of equations of motion.Despite the simplicity of the posterior mean, when the posterior distribution has multiple high-density areas, the posterior mean may not be able to capture the different solutions of the DE system.To handle the problem, we propose the idea of Leave-One-Out Cross-validation (LOOCV) to combine information from all high-density areas of the posterior and develop cross-validated NN-based estimates to more accurately estimate the DE variables under various solutions of the hyperbolic class equations of motion.Let N (x|Ω * ) = N 1 (x(t|Ω * 1 ), t, ϕ), . . ., N M (x(t|Ω * M ), t, ϕ) ⊤ represent the design matrix of the M estimates, cross-validates iteratively the NN candidates on N −i (x|Ω * ); thus the final estimates avoid giving unfair weights to spurious SMC proposals for the training set.It takes better into account the multiple high-density areas of the posterior distribution and consequently, it is expected to more efficiently estimate the DE variables of the system (25), giving the SMC estimates Ω * 1 , . . ., Ω * 1 .
44] for i = 1, . . ., 1000.We finally obtained L = 200 NN-based estimates for the critical collapse functions corresponding to L realizations from the posterior distributions.From a probabilistic perspective, each of these L realizations of NN-based estimates can occur to be the true form of the critical functions with probabilities π(ω * l |y) for l = 1, . . ., L. To better represent this probabilistic perspective and sampling variability from π(ω * l |y), we show ten randomly selected realizations of the posterior NN-based estimates for the critical function in Figures 2 and 7 under Gaussian and Uniform priors.We see that two different patterns are almost observed for critical collapse functions on the same domain which is compatible with the multiple solutions available in the literature for the equations of motion under the hyperbolic class of 4d.

Figure 2 :
Figure 2: Ten randomly selected NN-based estimates for DE variables b 0 (z) (A), Re(f (z)) (B) and Im(f (z)) (C) corresponding to ten SMC samples from the support of the posterior distribution π(ω|y) where Gaussian distribution was used as prior distribution.We show each estimate with a different line type.

3. 29 ]
corresponding to β-solution.In order to investigate the performance of developed models in estimating the functional form of the critical collapse functions under both solutions, we focus on the common areas between the two scenarios and ran investigated the NN-based estimates in [0,1.44].

Figure 3 :
Figure 3: The posterior mean (red), LOOCV (blue) and the true, using ω = 1.362, (orange) NN-based estimates of the critical collapse functions b 0 (z) (A), Im(f (z)) (B) and Re(f (z)) (C) corresponding to the α-solution scenario of the hyperbolic equations of motion.The dotted and dashed lines show, respectively, the lower and upper bounds of the 95% Bayesian credible intervals under Gaussian prior.

Figure 4 :
Figure 4: The posterior mean (red), LOOCV (blue) and the true, using ω = 1.003, (orange) NN-based estimates of the critical collapse functions b 0 (z) (A), Im(f (z)) (B) and Re(f (z)) (C) corresponding to the β-solution scenario of the hyperbolic equations of motion.The dotted and dashed lines show, respectively, the lower and upper bounds of the 95% Bayesian credible intervals under Gaussian prior.

Figures 3 and 8
Figures 3 and 8 show the performance of the posterior mean NN-based estimates, the true NN-based estimates, the LOOCV NN-based estimates as well as the 95% credible intervals in estimating the critical collapse functions corresponding to α-solution under Gaussian and Uniform priors, respectively.Also Figures 4 and 9 show the results of their counterpart NN-based estimates corresponding to β-solution under Gaussian and Uniform prior distributions, respectively.Since the posterior mean NN-based estimates and 95% credible intervals aggregate the L = 200 posterior NN-based estimates regardless of their corresponding population, these two methods remain robust and proposed the same estimates in predicting the critical collapse functions under both α-and β-solutions to the equations of motion.Unlike posterior mean proposals, we recommend the LOOCV-based estimates if one is interested in aggregating the NN candidates to more accurately estimate the form and curvature of the critical collapse functions corresponding to specific solutions.We see that the LOOCV-based method, on average, more accurately estimates the critical functions for both α and β solutions, because the method stacks the L = 200 Bayesian NN candidates by developing the best linear combination of the posterior NN candidates.

Figure 5 :
Figure 5: The trace-and box-plots of the differences between train and test loss values in the NN-based solvers using L = 200 SMC samples from the posterior distribution of π(ω|y) when Gaussian distribution was used as prior.

Figure 6 :
Figure 6: The posterior distributions and their trace-plots of two SMC chains (shown by solid and dotted blue lines) for parameters ω (w) and σ (sigma) under Uniform prior distribution.

Figure 7 :
Figure 7: Ten randomly selected NN-based estimates for DE variables b 0 (z) (A), Re(f (z)) (B) and Im(f (z)) (C) corresponding to ten SMC samples from the support of the posterior distribution π(ω|y) where Uniform distribution was used as prior distribution.We show each estimate with a different line type.

Figure 8 :
Figure 8: The posterior mean (red), LOOCV (blue) and the true, using ω = 1.362, (orange) NN-based estimates of the critical collapse functions b 0 (z) (A), Im(f (z)) (B) and Re(f (z)) (C) corresponding to the α-solution scenario of the hyperbolic equations of motion.The dotted and dashed lines show, respectively, the lower and upper bounds of the 95% Bayesian credible intervals under Uniform prior.

Figure 9 :
Figure 9: The posterior mean (red), LOOCV (blue) and the true, using ω = 1.003, (orange) NN-based estimates of the critical collapse functions b 0 (z) (A), Im(f (z)) (B) and Re(f (z)) (C) corresponding to the β-solution scenario of the hyperbolic equations of motion.The dotted and dashed lines show, respectively, the lower and upper bounds of the 95% Bayesian credible intervals under Uniform prior.

Figure 10 :
Figure 10: The trace-and box-plots of the differences between train and test loss values in the NN-based solvers using L = 200 SMC samples from the posterior distribution of π(ω|y) when Uniform distribution was used as prior.