Computational budget optimization for Bayesian parameter estimation in heavy ion collisions

T_RENTo is a parametric initial condition ansatz for simulating relativistic heavy ion collisions [21]. It takes as input the type of colliding nuclei and the inelastic nucleon-nucleon cross-section corresponding to the center-of-mass energy of the collision.

We vary three parameters of the T_RENTo model: the effective nucleon size $w$, the fluctuation parameter $k$ and the reduced thickness $p$. The fluctuation parameter $k$ allows for variation in the amount of energy carried by each nucleon, while the reduced thickness $p$ parametrize how the energy or entropy density is constructed from the profile of nucleons sampled from each nuclei.

The output of T_RENTo is a density profile in the plane transverse to the collision axis. We interpret T_RENTo’s output as an energy density profile $\epsilon(r,\phi)$ with arbitrary normalization. From this density profile $\epsilon(r,\phi)$, we compute single-event spatial anisotropies $\varepsilon_{n}$ [22]:

Although these anisotropies can be related to the momentum anisotropy of hadrons measured at the end of the collision [23, 24, 25, 26, 27], in this work, we focus on the initial spatial eccentricies themselves.

We compute the average of an ensemble of T_RENTo events with an arithmetic average:

We used up to four harmonics in this work, $n=2$ to $5$.

We studied minimum bias Pb-Pb collisions with Woods-Saxon nucleon distributions. We used an inelastic nucleon-nucleon cross-section of $64$ mb, corresponding approximately to $\sqrt{s_{NN}}=2.76$ TeV collisions.

We note $y_{\textrm{th},j}(\mathbf{p})$ the model’s prediction for the “$j$”-th observables of interest, evaluated when the model’s parameters are set to $\mathbf{p}$. At $\mathbf{p}$, we compute $y_{\textrm{th},j}(\mathbf{p})=\langle\varepsilon_{n_{j}}\rangle(\mathbf{p})$, where the index $j$ runs over all the observables of interest, which in our case are the different “$n$” in Eq. 2.

In general, given this information, we want to find the probability that parameters $\mathbf{p}$ are consistent with the data $\{y_{\textrm{exp},j}\}$ and with the experimental and theoretical uncertainties $\{\sigma_{\textrm{exp},j}\}$ and $\{\sigma_{\textrm{th},j}(\mathbf{p})\}$. The first step is to define a metric to quantify the level of model-to-data agreement: the likelihood function. We make the typical choice of a Gaussian likelihood function:

We assumed a diagonal covariance matrix. In this work, instead of comparing the model with data, we will use closure tests: we will replace $y_{\textrm{exp},j}$ by model calculations, as discussed later. These model calculations do stand for experimental data, and consequently we keep the label “exp” to denote these quantities.

Given our model of the collision, the probability that the parameter set $\mathbf{p}$ is consistent with the data and the uncertainties is given by the posterior distribution of model parameters, $\mathcal{P}(\mathbf{p}|(\{y_{\textrm{exp}}\})$:

where $\textrm{Prior}(\mathbf{p})$ is the prior distribution. In this work we take the prior to be constant over a finite parameter range, for all the parameters.

The posterior distribution $\mathcal{P}(\mathbf{p}|(\{y_{\textrm{exp}}\})$ is a probability distribution whose dimensionality is equal to the number of model parameters. Different projections of the posterior distribution summarize the constraints on the model parameters. As long as the number of parameters is small and the model is relatively fast computationally, it is straightforward to sample and marginalize the posterior $\mathcal{P}(\mathbf{p}|(\{y_{\textrm{exp}}\})$. When the number of parameters increases or if the model is expensive, it can become overly burdensome to evaluate a model’s prediction at a large number of values of the parameter $\mathbf{p}$ to compute $\mathcal{P}(\mathbf{p}|(\{y_{\textrm{exp}}\})$. A solution is to prepare a fast proxy, such as a Gaussian process [28], to emulate the model’s prediction. We review Gaussian process emulators briefly in the next subsection. As far as the Bayesian inference is concerned, the consequence of using Gaussian process emulators is substituting the model’s prediction by the emulator’s prediction,

as well as adding the emulator interpolation uncertainty into the sum of uncertainties,

The idea behind emulation is to first sample the parameter space of the model, $\{\mathbf{p}_{k}\}$; the parameter samples are the design points. The model’s predictions are then calculated for all design points: $\{y_{\textrm{th},j}(\mathbf{p}_{k})\}$. The emulator then interpolates the model’s predictions over a continuous range of model parameters $\mathbf{p}$.

To sample the parameter space, we use a low discrepancy sequence generator [29, 30]. Figure 1 shows an example for 16 design points within a two-dimensional parameter space.

For each observable ($\langle\epsilon_{n}\rangle$), a Gaussian process is trained on the set of values of $\langle\epsilon_{n}\rangle(\mathbf{p})$ calculated at the design points $\{\mathbf{p}_{k}\}$. Gaussian process emulators are probabilistic interpolators, which assume that each point in parameter space is a probability distribution. At the design points, the width of the probability distribution should include the statistical uncertainty of the observables. Between the design points, the width of the probability distribution should increase to account for the interpolation uncertainty. To handle stochastic simulations, Gaussian process emulators must be set up such that variations between parameter samples are divided into (i) statistical uncertainty, and (ii) the actual variation from the parameter dependence of the model. A mathematical description of the approach can be found in Ref. [15, Section V-B-2]. In short, a white noise kernel accounts for the uncorrelated (“short range”) fluctuations originating from the statistical uncertainty, and a squared-exponential kernel accounts for the longer-range parameter dependence. The parameters of these kernels, for example the relative size of the kernels, is determined by numerical optimization. This approach is used in most applications of Bayesian parameter inference in heavy ion physics. The result of this implementation is that the emulator uncertainty accounts for both the statistical and interpolation uncertainties.

Figure 2 shows the mean-value predictions of the emulator for $\langle\epsilon_{2}\rangle$ and $\langle\epsilon_{3}\rangle$ compared to the actual value of these observables. The points all lie close to the line $y_{\textrm{emul},j}(\mathbf{p})=y_{\textrm{th},j}(\mathbf{p})$, which shows that the emulator does indeed interpolate well. On the other hand, $\langle\epsilon_{3}\rangle$ tends to have a larger statistical uncertainty than $\langle\epsilon_{2}\rangle$, and as the right-hand side of Figure 2 shows, increasing the number of design points will not lead to perfect deterministic-like prediction from the emulator because the values of $\langle\epsilon_{3}\rangle$ have significant statistical uncertainties. This is the trade-off that we focus on in this work.

Closure tests are a straightforward application of Bayesian parameter inference, with experimental data replaced by calculations from the model itself. It is a self-consistency check made non-trivial by (i) the presence of uncertainties, as well as (ii) loss of information from the model to the observables. A longer discussion of closure tests, in the context of heavy ion collisions, can be found in Section VI of Ref. [15]. We describe it briefly here, with focus on metrics to quantify the success of closure tests.

In closure tests, we first select a set of parameters that are considered the reference parameter for the closure test; we refer to these parameters as the “reference” or “truth” parameters, $\mathbf{p}^{\textrm{truth}}$. We then compute event-averaged eccentricities $\langle\epsilon^{\textrm{truth}}_{n}\rangle=\langle\varepsilon_{n}\rangle(\mathbf{p}^{\textrm{truth}})$, which have a statistical uncertainty $\Delta\langle\epsilon^{\textrm{truth}}_{n}\rangle$. While these uncertainties are purely statistical, they are meant to mimic statistical and systematic uncertainties found in measurements. This uncertainty can be dialed by changing the number of T_RENTo events $M_{\textrm{ev}}^{\textrm{truth}}$ used to compute $\langle\epsilon^{\textrm{truth}}_{n}\rangle$.

Separately, we train Gaussian process emulators for $\langle\varepsilon_{n}\rangle(\mathbf{p})$ over a chosen range of model parameters, as described in the previous section. The two sources of uncertainties in the emulator are the number of design points (parameter samples) $N_{d}$, which controls the interpolation uncertainty, and the number of T_RENTo events per design point, $M_{\textrm{ev}}$, which control the statistical uncertainty of $\langle\varepsilon_{n}\rangle(\mathbf{p})$ at each design point.

Using the Gaussian process emulators as proxy for T_RENTo, we then perform Bayesian parameter estimation to attempt to recover the parameters $\mathbf{p}^{\textrm{truth}}$ from the values of the observables $\langle\epsilon^{\textrm{truth}}_{n}\rangle\pm\Delta\langle\epsilon^{\textrm{truth}}_{n}\rangle$. The result of the Bayesian inference is the posterior $\mathcal{P}(\mathbf{p}|\{\langle\epsilon^{\textrm{truth}}_{n}\rangle\})$ (Eq. 4). This posterior distribution depends on three sources of uncertainties: the uncertainty on the reference observables ($\Delta\langle\epsilon^{\textrm{truth}}_{n}\rangle$), which mimics experimental uncertainty, and the statistical and interpolation uncertainty in the emulator. If all uncertainties were small, one would expect the posterior to be a narrow peak at the parameter truth $\mathbf{p}^{\textrm{truth}}$. In practice, some of these uncertainties are significant, and the posterior distribution thus has some finite width in parameter space. If closure is successful, the posterior should enclose the parameter set $\mathbf{p}^{\textrm{truth}}$.

To quantify the degree of agreement of the posterior — $\mathcal{P}(\mathbf{p}|\{\langle\epsilon^{\textrm{truth}}_{n}\rangle\})$ — with the known true value $\mathbf{p}^{\textrm{truth}}$ of the parameters used in the closure test, we use two different metrics. The first one is the value of the posterior at the parameter truth, $\mathcal{P}(\mathbf{p}^{\textrm{truth}}|\{\langle\epsilon^{\textrm{truth}}_{n}\rangle\})$. The posterior at the parameter truth balances accuracy and precision. As illustrated in Fig. 3, an increase in accuracy (shift towards true value of the parameter) will increase the posterior value. An increase in precision (narrower posterior) will also increase the posterior at the parameter truth, but only to a certain extent: the posterior at the truth will actually start to decrease if the results get too confident about the incorrect value. The posterior at the truth has the benefit of being simple to interpret and inexpensive to compute.

For completeness, we studied a second metric, the Akaike Information Criterion (AIC) [31], given by:

where $\mathcal{L}_{\mathrm{max}}$ is the maximum likelihood in the parameter space and $k$ is the number of parameters. The Akaike Information Criterion is a metric used for model selection, to help determine how successful a model is at describing measurements given its number of parameters [32, 33]. In the case of closure tests, while we know that the emulator and the reference (truth) calculations originate from the same model, the uncertainties in the problem can evidently obscure this fact. We will see in this work the AIC’s success at quantifying closure. Note that the AIC is more expensive to compute, as it requires identifying the maximum of the likelihood.

To compute the reference (“truth”) observables, $M_{\textrm{ev}}^{\textrm{truth}}$ T_RENTo events are averaged. Changing $M_{\textrm{ev}}^{\textrm{truth}}$ controls the uncertainty on the reference observables, $\Delta\langle\epsilon^{\textrm{truth}}_{n}\rangle$. In this section, we vary $M_{\textrm{ev}}^{\textrm{truth}}$ to verify the effect on the optimal number of design points.

In the previous example, because $N_{\textrm{tot}}=N_{d}M_{\textrm{ev}}$ and $M_{\textrm{ev}}^{\textrm{truth}}$ were both equal to $2^{16}$, the interpolation and statistical uncertainty in the emulator (controlled by $N_{\textrm{tot}}$) were always large compared to the uncertainty of the reference calculations that are stand-in for measurements (controlled by $M_{\textrm{ev}}^{\textrm{truth}}$): at best, the statistical uncertainties for the calculations used in the emulator could be equal to the statistical uncertainty from the reference calculations (uncertainty $\propto 1/\sqrt{2^{16}}$) when $N_{d}=1$, which would evidently be a case with extreme interpolation uncertainty.

We first repeat the previous example with $M_{\textrm{ev}}^{\textrm{truth}}=2^{14}$ events while keeping $N_{\textrm{tot}}=N_{d}M_{\textrm{ev}}=2^{16}$. This means that, when $N_{d}=2^{2}=4$, the statistical uncertainty of the reference calculation is similar to that in the calculations used for the emulator; when $N_{d}=2^{4}=16$, the statistical uncertainty in the emulator is approximately a factor $2$ larger than in that of the reference calculations, etc. The results of the closure tests are shown on the left panel of Fig. 5. Overall, the results are similar: the constraints on the model parameters are consistent with the true value of the parameters at any ratio $N_{d}/M_{\textrm{ev}}$, but the range $N_{d}/M_{\textrm{ev}}\approx 0.1$–$1$ is found to be optimal.

We repeat this test in the opposite direction, using $M_{\textrm{ev}}^{\textrm{truth}}=2^{18}$ events to reduce the uncertainty on the reference calculations while keeping the emulator budget the same ($N_{\textrm{tot}}=N_{d}M_{\textrm{ev}}=2^{16}$). This is equivalent to a scenario where the emulator uncertainties are always being much larger than the data uncertainties. This yields the right panel of Fig. 5. Again, the optimal range of design points remains $N_{d}/M_{\textrm{ev}}\approx 0.1$–$1$.

Emulator performance depends on multiple factors, in particular the number of parameters, and the complexity of the parameter dependence of the observables. The results of Bayesian parameter inference evidently also depend on the factors just listed that impact emulator performance, and further depends on the information carried by the different observables. We repeat the closure test shown on Fig. 4 by adding new observables: we first add $\langle\varepsilon_{4}\rangle$, and then $\langle\varepsilon_{5}\rangle$. The results are shown respectively on the left and right panels of Fig. 6. Again, the results remain similar, with the optimal number of design points in the same range as found in previous tests.

We repeat this test, this time with three T_RENTo parameters instead of two, starting with two observables ($\langle\varepsilon_{n}\rangle$, $n=2,3$), then three ($n=2,3,4$) and four ($n=2,3,4,5$) observables. The results are shown respectively on the left, center and right panels of Fig. 7. The results are consistent with all previous ones: closure is successful since the constraints on all parameters are consistent with the truth value of the parameter for any number of design points. Moreover, constraints on the parameters are still best when the ratio of the number of design points $N_{d}$ to the number of T_RENTo simulations per design points $M_{\textrm{ev}}$ is in the range $N_{d}/M_{\textrm{ev}}\approx 0.1$–$1$.