PINNferring the Hubble Function with Uncertainties

Dear Editor, we're following your advice and incorporate the suggestions for minor revision put forward by the referees into our text. Please find below how exactly we reacted to the referee's suggestions; changes in the text are marked in blue.

Best wishes, Bjoern Malte Schaefer, on behalf of the authors

Referee 1:
The authors addressed my concerns about the manuscript. It would be useful to include some of the answers in the manuscript, but I consider this as optional.

--> We incorporate our previous answers into the text, as suggested by the referee, interfacing to the text which we presented in the previous version of our paper.

--> The authors suggest that (15) -> (16) is essentially a discretisation of the continuous function f. I would suggest to add this explanation to the main text.

--> We extended the text with an explanation:

"The distribution p(f|xi) is defined over the space of functions, whereas p(θ|xi) is a distribution of the parameters of the neurally represented function. In this process, one transitions from a continuous to a discrete space, which can be constrained meaningfully by the data. Given a high enough complexity of the neural network, it can represent any continuous target function."

--> Add a clarifying sentence after Eq. 19.

--> In the most general setting for the kernel this is just a function dependent on all data points in the batch. It is difficult to make this clear on the level of the equation in complete generality. Forward referencing to equation (21) would imply a specific structure in the kernels that might be unwanted or unneeded for the approach. Nevertheless, we made an amendment in the text, which emphasises that the Gaussian functional form is a choice with good numerical properties:

"Specifically, we make the choice to work with a Gaussian kernel, as specified by Eq.~\eqref{eq:gaussian_kernel}, which is numerically stable but by no means unique"

--> Please mention this in the text.

Author Answer: Adding Gaussian noise with larger or smaller standard deviations influences the heteroscedastic uncertainty estimates in the regions with sufficient data to reconstruct the results. The repulsive ensemble uncertainty estimates remain largely unchanged, leading to qualitatively the same solution.

--> We extended the text by adding this paragraph.

--> Please add this [i.e. the following statement] to the main text.

--> The emulation framework is tested and validated in reference [39]. Additionally, the test in fig. 7 demonstrates the viability of the emulator. Compared to the spline fit in https://doi.org/10.3847/1538-4357/aaa5a9 our uncertainties are roughly a factor 3 larger. The uncertainties in this methodology depend on the number of anchors for the splines. Compared to this existing method our proposed method is able to extract uncertainty estimates for a Hubble function parametrized by a large number of parameters.

--> We extended the text by adding this paragraph.

Author Answer: (A) The method is Bayesian. The likelihood is defined through the heteroscedastic loss and the repulsive ensembles. The prior is implicitly restricited through the variability of the network representing the Hubble function. By choosing a large enough network this restriction becomes less relevant. Additionally, the implemented weight decay gives a mathematical description of the prior. The method is Bayesian in the sense of ref [24].

(B) The Hubble function is recovered from the luminosity distance model through sampling from their predicted uncertainties. As a cross-check we have computed the Hubble function from the luminosity distance using eq. (30). This matches with the network reconstruction.

--> Bayesian methods update prior believes using data data likelihood. In order for this process to work, a minimum requirement is a clearly defined prior and data likelihood. Both is not the case here. I think calling the proposed method Bayesian-like or similar is fine. But this distinction should be however made clear in the paper.

--> We added a clarification to the paper, which discusses to the extent to which our method is truly Bayesian. We would like to emphasise that repulsive ensembles have been shown to provide a Bayesian characterisation of uncertainty in the infered quantity, but we do agree with the referee that we essentially work with an uninformative prior. The text is amended by this statement:

"To this purpose, we use a kernelised repulsive term in the update rule of the neural ensemble, following the example of \cite{Plehn:2022ftl} and \cite{DAngelo2021RepulsiveDE}. This modification not only maintains diversity between the ensemble members and prevents the collapse onto a single representation, but has been shown to realize a proper Bayesian inference by the ensemble, effectively propagating the error in the data to the uncertainty of the estimate of the posterior, in our case the neural representation of the (inverse) Hubble function. The repulsive ensemble does not incorporate prior information on the inferred function, though, which corresponds to working with an uninformative, essentially uniform Bayesian prior."


Referee 3:
In my original report, I stated that "The science case for using PINNs for the reconstruction of the Hubble diagram is not explored in sufficient detail." The authors now included three new references in page 2 which improved the contextualization of the text. However, there is still no discussion on how these methods compare to the proposed method with PINNs and supernovae data. I find this detrimental to the discussion of the quantitative potential of PINNs in this case.

--> We added a paragraph comparing different parameter-based and parameter-free reconstruction of the expansion history of the Universe:

"Constraints on the expansion history of the Universe are often formulated in terms of a parameterisation of the dark energy equation of state or the Hubble-function itself; of particular relevance are CPL-type polynomial expansions of the dark energy equation of state as a function of scale factor, or other types of polynomials or orthogonal systems of functions. Gaussian processes allow a greater flexibility of the function, and control the increased uncertainty due to the larger number of parameters by restricting the covariance. The neural representation of the Hubble-function integrated into the PINN for fast evaluation of the observable has in a natural way an extremely high degree of flexibility. But this flexibility, though, does not propagate into a large uncertainty of the largely increased number of parameters. Instead, the training process generates strong correlations between the parameters of the neural model and repulsive ensembles make sure that the uncertainty in the neural representation corresponds to that expected from Bayesian inference."

Moreover, the new sentence and references, as well as the answer provided by the authors, contradict the sentence in paragraph 3 of the Introduction which states that "the Hubble function is not directly observable". I pointed this out in my original report, but this sentence remained in the revision. Maybe the authors consider that the alternative methods of measuring H(z) require binning in redshift, as is the case for the BAO and AP methods, but this is not the case for all methods, such as cosmic chronometers or the redshift drift. So this sentence of the impossibility of direct measurements of H(z) should be removed, or it must be clarified what the authors mean here.

--> We clarified that in the distance redshift relation as considered in our paper, the Hubble-function enters the prediction of apparent magnitude for a given redshift only as an integral, not directly. This is all we wanted to express, but we agree with the referee that our statement was misleading in its generality and ultimately wrong. We hope to clarify this in this statement:

"In fact, the PINN learns a fast prediction of the luminosity distance for a given Hubble function, which is represented by a neural network in its full flexibility. This is necessary because supernova data constrains the luminosity distance for a measured redshift through the determination of the distance modulus. The luminosity function is results from the (inverse) Hubble function in a weighted integration, showing that in this particular case, the Hubble function is not directly observable. Other methods such as cosmic chronometers, are capable to constrain the Hubble function directly with no intermediate, numerical step to compute an observable."

On the second major issue I had raised in my original report, I suggested a comparison with other proposed methods of deriving H(z) from distances. The authors studied my suggested references and concluded that only the first one relied on supernovae data directly. This is correct, but it also misses the point I was making, which is a general lack in the manuscript of quantitative comparisons with alternative methods to constrain H(z) in general. Even if we stick to methods which rely on supernovae data directly, in their answer the authors conclude that the approach in 1805.03595 resulted in "uncertainty estimate is roughly a factor 3 smaller than ours". However, no discussion in this sense was included in the revised manuscript. I find that a quantitative comparison of methods is crucial to contextualize the importance of the PINN method for the case of H(z). This is important as reconstructing the Hubble function is the central theme of this paper.

--> We added a quantitative comparison of our results with other methods in the literature, and we discuss the issue about uncertainties as they depend on the choice of method and the flexibility of the non-parametric model, and of course on the choice of data set: Alreday within our own investigation we see large differences between Union and Pantheon.