Elena Sellentin

In Hamiltonian Monte Carlo sampling, the shape of the potential and the choice of the momentum distribution jointly give rise to the Hamiltonian dynamics of the sampler. An efficient sampler propagates quickly in all regions of the parameter space, so that the chain has a low autocorrelation length and the sampler has a high acceptance rate, with the goal of optimising the number of near-independent samples for given computational cost. Standard Gaussian momentum distributions allow arbitrarily large velocities, which can lead to inefficient exploration in posteriors with ridges or funnel-like geometries. We investigate alternative momentum distributions based on relativistic and Student's t kinetic energies, which naturally limit particle velocities and may improve robustness. Using Almanac, a sampler for cosmological posterior distributions of sky maps and power spectra on the sphere, we test these alternatives in both low- and high-dimensional settings. We find that the choice of parameterization and momentum distribution can improve convergence and effective sample rate, though the achievable gains are generally modest and strongly problem-dependent, reaching up to an order of magnitude in favorable cases. Among the momentum distributions that we tested, those with moderately heavy tails achieved the best balance between efficiency and stability. These results highlight the importance of sampler design and encourage future work on adaptive and self-tuning strategies for kinetic energy parameter optimization in high-dimensional settings.

We develop extreme data compression for use in Bayesian model comparison via the MOPED algorithm, as well as more general score compression. We find that Bayes Factors from data compressed with the MOPED algorithm are identical to those from their uncompressed datasets when the models are linear and the errors Gaussian. In other nonlinear cases, whether nested or not, we find negligible differences in the Bayes Factors, and show this explicitly for the Pantheon-SH0ES supernova dataset. We also investigate the sampling properties of the Bayesian Evidence as a frequentist statistic, and find that extreme data compression reduces the sampling variance of the Evidence, but has no impact on the sampling distribution of Bayes Factors. Since model comparison can be a very computationally-intensive task, MOPED extreme data compression may present significant advantages in computational time.

Simulating the evolution of the gravitational N-body problem becomes extremely computationally expensive as N increases since the problem complexity scales quadratically with the number of bodies. We study the use of Artificial Neural Networks (ANNs) to replace expensive parts of the integration of planetary systems. Neural networks that include physical knowledge have grown in popularity in the last few years, although few attempts have been made to use them to speed up the simulation of the motion of celestial bodies. We study the advantages and limitations of using Hamiltonian Neural Networks to replace computationally expensive parts of the numerical simulation. We compare the results of the numerical integration of a planetary system with asteroids with those obtained by a Hamiltonian Neural Network and a conventional Deep Neural Network, with special attention to understanding the challenges of this problem. Due to the non-linear nature of the gravitational equations of motion, errors in the integration propagate. To increase the robustness of a method that uses neural networks, we propose a hybrid integrator that evaluates the prediction of the network and replaces it with the numerical solution if considered inaccurate. Hamiltonian Neural Networks can make predictions that resemble the behavior of symplectic integrators but are challenging to train and in our case fail when the inputs differ ~7 orders of magnitude. In contrast, Deep Neural Networks are easy to train but fail to conserve energy, leading to fast divergence from the reference solution. The hybrid integrator designed to include the neural networks increases the reliability of the method and prevents large energy errors without increasing the computing cost significantly. For this problem, the use of neural networks results in faster simulations when the number of asteroids is >70.

We present a variational-Bayes solution to compute non-Gaussian posteriors from extremely expensive likelihoods. Our approach is an alternative for parameter inference when MCMC sampling is numerically prohibitive or conceptually unfeasible. For example, when either the likelihood or the theoretical model cannot be evaluated at arbitrary parameter values, but only previously selected values, then traditional MCMC sampling is impossible, whereas our variational-Bayes solution still succeeds in estimating the full posterior. In cosmology, this occurs e.g. when the parametric model is based on costly simulations that were run for previously selected input parameters. We demonstrate the applicability of our posterior construction on the KiDS-450 weak lensing analysis, where we reconstruct the original KiDS MCMC posterior at 0.6% of its former numerical posterior evaluations. The reduction in numerical cost implies that systematic effects which formerly exhausted the numerical budget could now be included.

Neural nets have become popular to accelerate parameter inferences, especially for the upcoming generation of galaxy surveys in cosmology. As neural nets are approximative by nature, a recurrent question has been how to propagate the neural net's approximation error, in order to avoid biases in the parameter inference. We present a Bayesian solution to propagating a neural net's approximation error and thereby debiasing parameter inference. We exploit that a neural net reports its approximation errors during the validation phase. We capture the thus reported approximation errors via the highest-order summary statistics, allowing us to eliminate the neural net's bias during inference, and propagating its uncertainties. We demonstrate that our method is quickly implemented and successfully infers parameters even for strongly biased neural nets. In summary, our method provides the missing element to judge the accuracy of a posterior if it cannot be computed based on an infinitely accurately theory code.

Upcoming galaxy surveys will allow us to probe the growth of the cosmic large-scale structure with improved sensitivity compared to current missions, and will also map larger areas of the sky. This means that in addition to the increased precision in observations, future surveys will also access the ultra-large-scale regime, where commonly neglected effects such as lensing, redshift-space distortions and relativistic corrections become important for calculating correlation functions of galaxy positions. At the same time, several approximations usually made in these calculations, such as the Limber approximation, break down at those scales. The need to abandon these approximations and simplifying assumptions at large scales creates severe issues for parameter estimation methods. On the one hand, exact calculations of theoretical angular power spectra become computationally expensive, and the need to perform them thousands of times to reconstruct posterior probability distributions for cosmological parameters makes the approach unfeasible. On the other hand, neglecting relativistic effects and relying on approximations may significantly bias the estimates of cosmological parameters. In this work, we quantify this bias and investigate how an incomplete modelling of various effects on ultra-large scales could lead to false detections of new physics beyond the standard $\Lambda$CDM model. Furthermore, we propose a simple debiasing method that allows us to recover true cosmologies without running the full parameter estimation pipeline with exact theoretical calculations. This method can therefore provide a fast way of obtaining accurate values of cosmological parameters and estimates of exact posterior probability distributions from ultra-large-scale observations.

Structure formation in our Universe creates non-Gaussian random fields that will soon be observed over almost the entire sky by the Euclid satellite, the Vera-Rubin observatory, and the Square Kilometre Array. An unsolved problem is how to analyze best such non-Gaussian fields, e.g. to infer the physical laws that created them. This problem could be solved if a parametric non-Gaussian sampling distribution for such fields were known, as this distribution could serve as likelihood during inference. We therefore create a sampling distribution for non-Gaussian random fields. Our approach is capable of handling strong non-Gaussianity, while perturbative approaches such as the Edgeworth expansion cannot. To imitate cosmological structure formation, we enforce our fields to be (i) statistically isotropic, (ii) statistically homogeneous, and (iii) statistically independent at large distances. We generate such fields via a Monte Carlo Markov Chain technique and find that even strong non-Gaussianity is not necessarily visible to the human eye. We also find that sampled marginals for pixel pairs have an almost generic Gauss-like appearance, even if the joint distribution of all pixels is markedly non-Gaussian. This apparent Gaussianity is a consequence of the high dimensionality of random fields. We conclude that vast amounts of non-Gaussian information can be hidden in random fields that appear nearly Gaussian in simple tests, and that it would be short-sighted not to try and extract it.

Non-Gaussian (NG) statistics of the thermal Sunyaev-Zeldovich (tSZ) effect carry significant information which is not contained in the power spectrum. Here, we perform a joint Fisher analysis of the tSZ power spectrum and bispectrum to verify how much the full bispectrum can contribute to improve parameter constraints. We go beyond similar studies of this kind in several respects: first of all, we include the complete power spectrum and bispectrum (auto- and cross-) covariance in the analysis, computing all NG contributions; furthermore we consider a multi-component foreground scenario and model the effects of component separation in the forecasts; finally, we consider an extended set of both cosmological and intra-cluster medium parameters. We show that the tSZ bispectrum is very efficient at breaking parameter degeneracies, making it able to produce even stronger cosmological constraints than the tSZ power spectrum: e.g. the standard deviation on σ₈ shrinks from σ^PS(σ₈)=0.35 to σ^BS(σ₈)=0.065 when we consider a multi-parameter analysis. We find that this is mostly due to the different response of separate triangle types (e.g. equilateral and squeezed) to changes in model parameters. While weak, this shape dependence is clearly non-negligible for cosmological parameters, and it is even stronger, as expected, for intra-cluster medium parameters.