Search

Sequential Monte Carlo (SMC) methods are powerful tools for Bayesian inference but suffer from requiring many particles for accurate estimates, leading to high computational costs. We introduce persistent sampling (PS), an extension of SMC that mitigates this issue by allowing particles from previous iterations to persist. This generates a growing, weighted ensemble of particles distributed across iterations. In each iteration, PS utilizes multiple importance sampling and resampling from the mixture of all previous distributions to produce the next generation of particles. This addresses particle impoverishment and mode collapse, resulting in more accurate posterior approximations. Furthermore, this approach provides lower-variance marginal likelihood estimates for model comparison. Additionally, the persistent particles improve transition kernel adaptation for efficient exploration. Experiments on complex distributions show that PS consistently outperforms standard methods, achieving lower squared bias in posterior moment estimation and significantly reduced marginal likelihood errors, all at a lower computational cost. PS offers a robust, efficient, and scalable framework for Bayesian inference., Comment: 30 pages, 9 figures, 4 tables. Submitted to Statistics & Computing

Periodograms are widely employed for identifying periodicity in time series data, yet they often struggle to accurately quantify the statistical significance of detected periodic signals when the data complexity precludes reliable simulations. We develop a data-driven approach to address this challenge by introducing a null-signal template (NST). The NST is created by carefully randomizing the period of each cycle in the periodogram template, rendering it non-periodic. It has the same frequentist properties as a periodic signal template regardless of the noise probability distribution, and we show with simulations that the distribution of false positives is the same as with the original periodic template, regardless of the underlying data. Thus, performing a periodicity search with the NST acts as an effective simulation of the null (no-signal) hypothesis, without having to simulate the noise properties of the data. We apply the NST method to the supermassive black hole binaries (SMBHB) search in the Palomar Transient Factory (PTF), where Charisi et al. had previously proposed 33 high signal to (white) noise candidates utilizing simulations to quantify their significance. Our approach reveals that these simulations do not capture the complexity of the real data. There are no statistically significant periodic signal detections above the non-periodic background. To improve the search sensitivity we introduce a Gaussian quadrature based algorithm for the Bayes Factor with correlated noise as a test statistic, in contrast to the standard signal to white noise. We show with simulations that this improves sensitivity to true signals by more than an order of magnitude. However, using the Bayes Factor approach also results in no statistically significant detections in the PTF data., Comment: 13 pages, 12 figures

Ensemble Kalman Inversion (EKI) has been proposed as an efficient method for solving inverse problems with expensive forward models. However, the method is based on the assumption that we proceed through a sequence of Gaussian measures in moving from the prior to the posterior, and that the forward model is linear. In this work, we introduce Sequential Kalman Monte Carlo (SKMC) samplers, where we exploit EKI and Flow Annealed Kalman Inversion (FAKI) within a Sequential Monte Carlo (SMC) sampling scheme to perform efficient gradient-free inference in Bayesian inverse problems. FAKI employs normalizing flows (NF) to relax the Gaussian ansatz of the target measures in EKI. NFs are able to learn invertible maps between a Gaussian latent space and the original data space, allowing us to perform EKI updates in the Gaussianized NF latent space. However, FAKI alone is not able to correct for the model linearity assumptions in EKI. Errors in the particle distribution as we move through the sequence of target measures can therefore compound to give incorrect posterior moment estimates. In this work we consider the use of EKI and FAKI to initialize the particle distribution for each target in an adaptive SMC annealing scheme, before performing t-preconditioned Crank-Nicolson (tpCN) updates to distribute particles according to the target. We demonstrate the performance of these SKMC samplers on three challenging numerical benchmarks, showing significant improvements in the rate of convergence compared to standard SMC with importance weighted resampling at each temperature level. Code implementing the SKMC samplers is available at https://github.com/RichardGrumitt/KalmanMC., Comment: 39 pages, 15 figures, code available at https://github.com/RichardGrumitt/KalmanMC

Weak Lensing (WL) surveys are reaching unprecedented depths, enabling the investigation of very small angular scales. At these scales, nonlinear gravitational effects lead to higher-order correlations making the matter distribution highly non-Gaussian. Extracting this information using traditional statistics has proven difficult, and Machine Learning based summary statistics have emerged as a powerful alternative. We explore the capabilities of a discriminative, Convolutional Neural Networks (CNN) based approach, focusing on parameter constraints in the ($\Omega_m$, $\sigma_8$) cosmological parameter space. Leveraging novel training loss functions and network representations on WL mock datasets without baryons, we show that our models achieve $\sim 5$ times stronger constraints than the power spectrum, $\sim 3$ stronger constraints than peak counts, and $\sim 2$ stronger constraints than previous CNN-learned summary statistics and scattering transforms, for noise levels relevant to Rubin or Euclid. For WL convergence maps with baryonic physics, our models achieve $\sim 2.3$ times stronger constraining power than the power spectrum at these noise levels, also outperforming previous summary statistics. To further explore the possibilities of CNNs for this task, we also discuss transfer learning where we adapt pre-trained models, trained on different tasks or datasets, for cosmological inference, finding that these do not improve the performance., Comment: 16 pages, 2 figures. Comments Welcome

We develop a new and simple method to model baryonic effects at the field level relevant for weak lensing analyses. We analyze thousands of state-of-the-art hydrodynamic simulations from the CAMELS project, each with different cosmology and strength of feedback, and we find that the cross-correlation coefficient between full hydrodynamic and N-body simulations is very close to 1 down to $k\sim10~h{\rm Mpc}^{-1}$. This suggests that modeling baryonic effects at the field level down to these scales only requires N-body simulations plus a correction to the mode's amplitude given by: $\sqrt{P_{\rm hydro}(k)/P_{\rm nbody}(k)}$. In this paper, we build an emulator for this quantity, using Gaussian processes, that is flexible enough to reproduce results from thousands of hydrodynamic simulations that have different cosmologies, astrophysics, subgrid physics, volumes, resolutions, and at different redshifts. Our emulator is accurate at the percent level and exhibits a range of validation superior to previous studies. This method and our emulator enable field-level simulation-based inference analyses and accounting for baryonic effects in weak lensing analyses., Comment: 12 pages, 9 figures. Comments welcome

We introduce a global, gradient-free surrogate optimization strategy for expensive black-box functions inspired by the Fokker-Planck and Langevin equations. These can be written as an optimization problem where the objective is the target function to maximize minus the logarithm of the current density of evaluated samples. This objective balances exploitation of the target objective with exploration of low-density regions. The method, Deterministic Langevin Optimization (DLO), relies on a Normalizing Flow density estimate to perform active learning and select proposal points for evaluation. This strategy differs qualitatively from the widely-used acquisition functions employed by Bayesian Optimization methods, and can accommodate a range of surrogate choices. We demonstrate superior or competitive progress toward objective optima on standard synthetic test functions, as well as on non-convex and multi-modal posteriors of moderate dimension. On real-world objectives, such as scientific and neural network hyperparameter optimization, DLO is competitive with state-of-the-art baselines., Comment: 31 pages, 9 figures, 1 table, submitted to JOGO

For many scientific inverse problems we are required to evaluate an expensive forward model. Moreover, the model is often given in such a form that it is unrealistic to access its gradients. In such a scenario, standard Markov Chain Monte Carlo algorithms quickly become impractical, requiring a large number of serial model evaluations to converge on the target distribution. In this paper we introduce Flow Annealed Kalman Inversion (FAKI). This is a generalization of Ensemble Kalman Inversion (EKI), where we embed the Kalman filter updates in a temperature annealing scheme, and use normalizing flows (NF) to map the intermediate measures corresponding to each temperature level to the standard Gaussian. In doing so, we relax the Gaussian ansatz for the intermediate measures used in standard EKI, allowing us to achieve higher fidelity approximations to non-Gaussian targets. We demonstrate the performance of FAKI on two numerical benchmarks, showing dramatic improvements over standard EKI in terms of accuracy whilst accelerating its already rapid convergence properties (typically in $\mathcal{O}(10)$ steps)., Comment: 9 pages, 2 figues. Presented at MaxEnt 2023. Modified version to appear in MaxEnt 2023 proceedings

Field-level inference provides a means to optimally extract information from upcoming cosmological surveys, but requires efficient sampling of a high-dimensional parameter space. This work applies Microcanonical Langevin Monte Carlo (MCLMC) to sample the initial conditions of the Universe, as well as the cosmological parameters $\sigma_8$ and $\Omega_m$, from simulations of cosmic structure. MCLMC is shown to be over an order of magnitude more efficient than traditional Hamiltonian Monte Carlo (HMC) for a $\sim 2.6 \times 10^5$ dimensional problem. Moreover, the efficiency of MCLMC compared to HMC greatly increases as the dimensionality increases, suggesting gains of many orders of magnitude for the dimensionalities required by upcoming cosmological surveys., Comment: Accepted at the ICML 2023 Workshop on Machine Learning for Astrophysics. 4 pages, 4 figures

We propose Multiscale Flow, a generative Normalizing Flow that creates samples and models the field-level likelihood of two-dimensional cosmological data such as weak lensing. Multiscale Flow uses hierarchical decomposition of cosmological fields via a wavelet basis, and then models different wavelet components separately as Normalizing Flows. The log-likelihood of the original cosmological field can be recovered by summing over the log-likelihood of each wavelet term. This decomposition allows us to separate the information from different scales and identify distribution shifts in the data such as unknown scale-dependent systematics. The resulting likelihood analysis can not only identify these types of systematics, but can also be made optimal, in the sense that the Multiscale Flow can learn the full likelihood at the field without any dimensionality reduction. We apply Multiscale Flow to weak lensing mock datasets for cosmological inference, and show that it significantly outperforms traditional summary statistics such as power spectrum and peak counts, as well as novel Machine Learning based summary statistics such as scattering transform and convolutional neural networks. We further show that Multiscale Flow is able to identify distribution shifts not in the training data such as baryonic effects. Finally, we demonstrate that Multiscale Flow can be used to generate realistic samples of weak lensing data., Comment: 14 pages, 7 figures. Comments welcome

Stochastic sampling algorithms such as Langevin Monte Carlo are inspired by physical systems in a heat bath. Their equilibrium distribution is the canonical ensemble given by a prescribed target distribution, so they must balance fluctuation and dissipation as dictated by the fluctuation-dissipation theorem. In contrast to the common belief, we show that the fluctuation-dissipation theorem is not required because only the configuration space distribution, and not the full phase space distribution, needs to be canonical. We propose a continuous-time Microcanonical Langevin Monte Carlo (MCLMC) as a dissipation-free system of stochastic differential equations (SDE). We derive the corresponding Fokker-Planck equation and show that the stationary distribution is the microcanonical ensemble with the desired canonical distribution on configuration space. We prove that MCLMC is ergodic for any nonzero amount of stochasticity, and for smooth, convex potentials, the expectation values converge exponentially fast. Furthermore, the deterministic drift and the stochastic diffusion separately preserve the stationary distribution. This uncommon property is attractive for practical implementations as it implies that the drift-diffusion discretization schemes are bias-free, so the only source of bias is the discretization of the deterministic dynamics. We applied MCLMC on a lattice $\phi^4$ model, where Hamiltonian Monte Carlo (HMC) is currently the state-of-the-art integrator. For the same accuracy, MCLMC converges 12 times faster than HMC on an $8\times8$ lattice. On a $64\times64$ lattice, it is already 32 times faster. The trend is expected to persist to larger lattices, which are of particular interest, for example, in lattice quantum chromodynamics., Comment: 8 pages, 1 figure

We develop Microcanonical Hamiltonian Monte Carlo (MCHMC), a class of models which follow a fixed energy Hamiltonian dynamics, in contrast to Hamiltonian Monte Carlo (HMC), which follows canonical distribution with different energy levels. MCHMC tunes the Hamiltonian function such that the marginal of the uniform distribution on the constant-energy-surface over the momentum variables gives the desired target distribution. We show that MCHMC requires occasional energy conserving billiard-like momentum bounces for ergodicity, analogous to momentum resampling in HMC. We generalize the concept of bounces to a continuous version with partial direction preserving bounces at every step, which gives an energy conserving underdamped Langevin-like dynamics with non-Gaussian noise (MCLMC). MCHMC and MCLMC exhibit favorable scalings with condition number and dimensionality. We develop an efficient hyperparameter tuning scheme that achieves high performance and consistently outperforms NUTS HMC on several standard benchmark problems, in some cases by more than an order of magnitude., Comment: 34 pages, 11 figures

We present growth of structure constraints from the cosmological analysis of the power spectrum multipoles of SDSS-III BOSS DR12 galaxies. We use the galaxy power spectrum model of Hand et al. (2017), which decomposes the galaxies into halo mass bins, each of which is modeled separately using the relations between halo biases and halo mass. The model combines Eulerian perturbation theory and halo model calibrated on $N$-body simulations to model the halo clustering. In this work, we also generate the covariance matrix by combining the analytic disconnected part with the empirical connected part: we smooth the connected component by selecting a few principal components and show that it achieves good agreement with the mock covariance. Our analysis differs from recent analyses in that we constrain a single parameter $f\sigma_8$ fixing everything else to Planck+BAO prior, thereby reducing the effects of prior volume and mismodeling. We find tight constraints on $f\sigma_8$: $f\sigma_8(z_{\mathrm{eff}}=0.38)=0.489 \pm 0.038$ and $f\sigma_8(z_{\mathrm{eff}}=0.61)=0.455 \pm 0.028$ at $k_{\mathrm{max}} = 0.2\ h$Mpc$^{-1}$, with an overall amplitude error of 5%, and in good agreement (within 0.3 sigma) of Planck amplitude. We discuss the sensitivity of cosmological parameter estimation to the choice of scale cuts, covariance matrix, and the inclusion of hexadecapole $P_4(k)$. We show that with $k_{\mathrm{max}} = 0.4\ h$Mpc$^{-1}$ the constraints improve considerably to an overall 3.2% amplitude error, but there is some evidence of model misspecification on MultiDark-PATCHY mocks. Choosing $k_{\mathrm{max}}$ consistently and reliably remains the main challenge of RSD analysis methods., Comment: 21 pages, 13 figures

In this white paper, we present the MegaMapper concept. The MegaMapper is a proposed ground-based experiment to measure Inflation parameters and Dark Energy from galaxy redshifts at $2

In this white paper, we present an experimental road map for spectroscopic experiments beyond DESI. DESI will be a transformative cosmological survey in the 2020s, mapping 40 million galaxies and quasars and capturing a significant fraction of the available linear modes up to z=1.2. DESI-II will pilot observations of galaxies both at much higher densities and extending to higher redshifts. A Stage-5 experiment would build out those high-density and high-redshift observations, mapping hundreds of millions of stars and galaxies in three dimensions, to address the problems of inflation, dark energy, light relativistic species, and dark matter. These spectroscopic data will also complement the next generation of weak lensing, line intensity mapping and CMB experiments and allow them to reach their full potential., Comment: Contribution to Snowmass 2021

We construct a physically-parameterized probabilistic autoencoder (PAE) to learn the intrinsic diversity of type Ia supernovae (SNe Ia) from a sparse set of spectral time series. The PAE is a two-stage generative model, composed of an Auto-Encoder (AE) which is interpreted probabilistically after training using a Normalizing Flow (NF). We demonstrate that the PAE learns a low-dimensional latent space that captures the nonlinear range of features that exists within the population, and can accurately model the spectral evolution of SNe Ia across the full range of wavelength and observation times directly from the data. By introducing a correlation penalty term and multi-stage training setup alongside our physically-parameterized network we show that intrinsic and extrinsic modes of variability can be separated during training, removing the need for the additional models to perform magnitude standardization. We then use our PAE in a number of downstream tasks on SNe Ia for increasingly precise cosmological analyses, including automatic detection of SN outliers, the generation of samples consistent with the data distribution, and solving the inverse problem in the presence of noisy and incomplete data to constrain cosmological distance measurements. We find that the optimal number of intrinsic model parameters appears to be three, in line with previous studies, and show that we can standardize our test sample of SNe Ia with an RMS of $0.091 \pm 0.010$ mag, which corresponds to $0.074 \pm 0.010$ mag if peculiar velocity contributions are removed. Trained models and codes are released at \href{https://github.com/georgestein/suPAErnova}{github.com/georgestein/suPAErnova}, Comment: 23 pages, 8 Figures, 1 Table. Accepted to ApJ

In many hypothesis testing applications, we have mixed priors, with well-motivated informative priors for some parameters but not for others. The Bayesian methodology uses the Bayes factor and is helpful for the informative priors, as it incorporates Occam's razor via multiplicity or trials factor in the Look Elsewhere Effect. However, if the prior is not known completely, the frequentist hypothesis test via the false positive rate is a better approach, as it is less sensitive to the prior choice. We argue that when only partial prior information is available, it is best to combine the two methodologies by using the Bayes factor as a test statistic in the frequentist analysis. We show that the standard frequentist likelihood-ratio test statistic corresponds to the Bayes factor with a non-informative Jeffrey's prior. We also show that mixed priors increase the statistical power in frequentist analyses over the likelihood ratio test statistic. We develop an analytic formalism that does not require expensive simulations using a statistical mechanics approach to hypothesis testing in Bayesian and frequentist statistics. We introduce the counting of states in a continuous parameter space using the uncertainty volume as the quantum of the state. We show that both the p-value and Bayes factor can be expressed as energy versus entropy competition. We present analytic expressions that generalize Wilks' theorem beyond its usual regime of validity and work in a non-asymptotic regime. In specific limits, the formalism reproduces existing expressions, such as the p-value of linear models and periodograms. We apply the formalism to an example of exoplanet transits, where multiplicity can be more than $10^7$. We show that our analytic expressions reproduce the p-values derived from numerical simulations., Comment: 15 pages, 7 figures

pocoMC is a Python package for accelerated Bayesian inference in astronomy and cosmology. The code is designed to sample efficiently from posterior distributions with non-trivial geometry, including strong multimodality and non-linearity. To this end, pocoMC relies on the Preconditioned Monte Carlo algorithm which utilises a Normalising Flow in order to decorrelate the parameters of the posterior. It facilitates both tasks of parameter estimation and model comparison, focusing especially on computationally expensive applications. It allows fitting arbitrary models defined as a log-likelihood function and a log-prior probability density function in Python. Compared to popular alternatives (e.g. nested sampling) pocoMC can speed up the sampling procedure by orders of magnitude, cutting down the computational cost substantially. Finally, parallelisation to computing clusters manifests linear scaling., Comment: 6 pages, 1 figure. Submitted to JOSS. Code available at https://github.com/minaskar/pocomc

We introduce Preconditioned Monte Carlo (PMC), a novel Monte Carlo method for Bayesian inference that facilitates efficient sampling of probability distributions with non-trivial geometry. PMC utilises a Normalising Flow (NF) in order to decorrelate the parameters of the distribution and then proceeds by sampling from the preconditioned target distribution using an adaptive Sequential Monte Carlo (SMC) scheme. The results produced by PMC include samples from the posterior distribution and an estimate of the model evidence that can be used for parameter inference and model comparison respectively. The aforementioned framework has been thoroughly tested in a variety of challenging target distributions achieving state-of-the-art sampling performance. In the cases of primordial feature analysis and gravitational wave inference, PMC is approximately 50 and 25 times faster respectively than Nested Sampling (NS). We found that in higher dimensional applications the acceleration is even greater. Finally, PMC is directly parallelisable, manifesting linear scaling up to thousands of CPUs., Comment: 10 pages, 9 figures. Submitted to MNRAS. Code available at https://github.com/minaskar/pocomc

We propose a general purpose Bayesian inference algorithm for expensive likelihoods, replacing the stochastic term in the Langevin equation with a deterministic density gradient term. The particle density is evaluated from the current particle positions using a Normalizing Flow (NF), which is differentiable and has good generalization properties in high dimensions. We take advantage of NF preconditioning and NF based Metropolis-Hastings updates for a faster convergence. We show on various examples that the method is competitive against state of the art sampling methods., Comment: 17 pages, 9 figures, Accepted at NeurIPS 2022

Our universe is homogeneous and isotropic, and its perturbations obey translation and rotation symmetry. In this work we develop Translation and Rotation Equivariant Normalizing Flow (TRENF), a generative Normalizing Flow (NF) model which explicitly incorporates these symmetries, defining the data likelihood via a sequence of Fourier space-based convolutions and pixel-wise nonlinear transforms. TRENF gives direct access to the high dimensional data likelihood p(x|y) as a function of the labels y, such as cosmological parameters. In contrast to traditional analyses based on summary statistics, the NF approach has no loss of information since it preserves the full dimensionality of the data. On Gaussian random fields, the TRENF likelihood agrees well with the analytical expression and saturates the Fisher information content in the labels y. On nonlinear cosmological overdensity fields from N-body simulations, TRENF leads to significant improvements in constraining power over the standard power spectrum summary statistic. TRENF is also a generative model of the data, and we show that TRENF samples agree well with the N-body simulations it trained on, and that the inverse mapping of the data agrees well with a Gaussian white noise both visually and on various summary statistics: when this is perfectly achieved the resulting p(x|y) likelihood analysis becomes optimal. Finally, we develop a generalization of this model that can handle effects that break the symmetry of the data, such as the survey mask, which enables likelihood analysis on data without periodic boundaries., Comment: 11 pages, 10 figures. Submitted to MNRAS. Comments welcome

We present the marginal unbiased score expansion (MUSE) method, an algorithm for generic high-dimensional hierarchical Bayesian inference. MUSE performs approximate marginalization over arbitrary non-Gaussian latent parameter spaces, yielding Gaussianized asymptotically unbiased and near-optimal constraints on global parameters of interest. It is computationally much cheaper than exact alternatives like Hamiltonian Monte Carlo (HMC), excelling on funnel problems which challenge HMC, and does not require any problem-specific user supervision like other approximate methods such as Variational Inference or many Simulation-Based Inference methods. MUSE makes possible the first joint Bayesian estimation of the delensed Cosmic Microwave Background (CMB) power spectrum and gravitational lensing potential power spectrum, demonstrated here on a simulated data set as large as the upcoming South Pole Telescope 3G 1500 deg$^2$ survey, corresponding to a latent dimensionality of ${\sim}\,6$ million and of order 100 global bandpower parameters. On a subset of the problem where an exact but more expensive HMC solution is feasible, we verify that MUSE yields nearly optimal results. We also demonstrate that existing spectrum-based forecasting tools which ignore pixel-masking underestimate predicted error bars by only ${\sim}\,10\%$. This method is a promising path forward for fast lensing and delensing analyses which will be necessary for future CMB experiments such as SPT-3G, Simons Observatory, or CMB-S4, and can complement or supersede existing HMC approaches. The success of MUSE on this challenging problem strengthens its case as a generic procedure for a broad class of high-dimensional inference problems., Comment: 22 pages, 8 figures

In experiments where one searches a large parameter space for an anomaly, one often finds many spurious noise-induced peaks in the likelihood. This is known as the look-elsewhere effect, and must be corrected for when performing statistical analysis. This paper introduces a method to calibrate the false alarm probability (FAP), or $p$-value, for a given dataset by considering the heights of the highest peaks in the likelihood. In the simplest form of self-calibration, the look-elsewhere-corrected $\chi^2$ of a physical peak is approximated by the $\chi^2$ of the peak minus the $\chi^2$ of the highest noise-induced peak. Generalizing this concept to consider lower peaks provides a fast method to quantify the statistical significance with improved accuracy. In contrast to alternative methods, this approach has negligible computational cost as peaks in the likelihood are a byproduct of every peak-search analysis. We apply to examples from astronomy, including planet detection, periodograms, and cosmology., Comment: 12 pages, 7 figures

Massive neutrinos suppress the growth of cosmic structure on small, non-linear, scales. It is thus often proposed that using statistics beyond the power spectrum can tighten constraints on the neutrino mass by extracting additional information from these non-linear scales. We study the information content regarding neutrino mass at the field level, quantifying how much of this information arises from the difference in non-linear evolution between a cosmology with 1 fluid (CDM) and 2 fluids (CDM + neutrinos). We do so by running two $N$-body simulations, one with and one without massive neutrinos; both with the same phases, and matching their linear power spectrum at a given, low, redshift. This effectively isolates the information encoded in the linear initial conditions from the non-linear cosmic evolution. We demonstrate that for $k \lesssim 1\,h/{\rm Mpc}$, and for a single redshift, there is negligible difference in the real-space CDM field between the two simulations. This suggests that all the information regarding neutrino mass is in the linear power spectrum set by the initial conditions. Thus any probe based on the CDM field alone will have negligible constraining power beyond that which exists at the linear level over the same range of scales. Consequently, any probe based on the halo field will contain little information beyond the linear power. We find similar results for the matter field responsible for weak lensing. We also demonstrate that there may be much information beyond the power spectrum in the 3d matter field, however, this is not observable in modern surveys via dark matter halos or weak lensing. Finally, we show that there is additional information to be found in redshift space., Comment: 10 pages, 8 figures

Reconstructing the Gaussian initial conditions at the beginning of the Universe from the survey data in a forward modeling framework is a major challenge in cosmology. This requires solving a high dimensional inverse problem with an expensive, non-linear forward model: a cosmological N-body simulation. While intractable until recently, we propose to solve this inference problem using an automatically differentiable N-body solver, combined with a recurrent networks to learn the inference scheme and obtain the maximum-a-posteriori (MAP) estimate of the initial conditions of the Universe. We demonstrate using realistic cosmological observables that learnt inference is 40 times faster than traditional algorithms such as ADAM and LBFGS, which require specialized annealing schemes, and obtains solution of higher quality., Comment: Published as a workshop paper at ICLR 2021 SimDL Workshop

We update Halo Zeldovich Perturbation Theory (HZPT), an analytic model for the two-point statistics of dark matter, to describe halo and galaxy clustering, and galaxy-matter cross-correlation on nonlinear scales. The model correcting Zeldovich has an analytic Fourier transform, and therefore is valid in both configuration space and Fourier space. The model is accurate at the $2\%$-level or less for $P_{mm}$ (k < 1 h/Mpc), $P_{hm}$ (k < 1 h/Mpc), $P_{hh}$ (k < 2 h/Mpc), $P_{gm}$ (k < 1 h/Mpc), $P_{gg}$ (k < 1 h/Mpc), $\xi_{mm}$ (r > 1 Mpc/h), $\xi_{hm}$ (r > 2 Mpc/h), $\xi_{hh}$ (r > 2 Mpc/h), $\xi_{gm}$ (r > 1 Mpc/h), $\xi_{gg}$ (r > 2 Mpc/h), for LRG-like mock galaxies. We show that the HZPT model for matter correlators can account for the effects of a wide range of baryonic feedback models and provide extended dark matter models which are of $1\% ~(3\%)$ accuracy for k < 10 (8) h/Mpc. We explicitly model the non-perturbative features of halo exclusion for the halo-halo and galaxy-galaxy correlators, as well as the presence of satellites for galaxy-matter and galaxy-galaxy correlation functions. We perform density estimation using N-body simulations and a wide range of HOD galaxy mocks to obtain correlations of model parameters with the cosmological parameters $\Omega_{m}$ and $\sigma_{8}$. HZPT can provide a fast, interpretable, and analytic model for combined-probe analyses of redshift surveys using scales well into the non-linear regime., Comment: 46 pages, 21 figures, (including appendices). Published version

Anomaly detection is a key application of machine learning, but is generally focused on the detection of outlying samples in the low probability density regions of data. Here we instead present and motivate a method for unsupervised in-distribution anomaly detection using a conditional density estimator, designed to find unique, yet completely unknown, sets of samples residing in high probability density regions. We apply this method towards the detection of new physics in simulated Large Hadron Collider (LHC) particle collisions as part of the 2020 LHC Olympics blind challenge, and show how we detected a new particle appearing in only 0.08% of 1 million collision events. The results we present are our original blind submission to the 2020 LHC Olympics, where it achieved the state-of-the-art performance., Comment: Accepted to NeurIPS Machine Learning and the Physical Sciences workshop. See arXiv:2007.00674 for further methods

Exoplanet Transit Timing Variations (TTVs) caused by gravitational forces between planets can be used to determine planetary masses and orbital parameters. Most of the observed TTVs are small and sinusoidal in time, leading to degeneracies between the masses and orbital parameters. Here we report a TTV analysis of Kepler-90g and Kepler-90h, which exhibit large TTVs up to 25 hours. With optimization, we find a unique solution which allows us to constrain all of the orbital parameters. The best fit masses for Kepler-90g and 90h are $15.0^{+0.9}_{-0.8}$ $M_{\bigoplus}$ (Earth mass) and $203^{+5}_{-5}M_{\bigoplus}$, respectively, with Kepler-90g having an unusually low apparent density of $0.15\pm 0.05\, {\rm g\,cm^{-3}}$. The uniqueness of orbital parameter solution enables a long-term dynamical integration, which reveals that although their periods are close to 2:3 orbital resonance, they are not locked in resonance, and the configuration is stable over billions of years. The dynamical history of the system suggests that planet interactions are able to raise the eccentricities and break the resonant lock after the initial formation., Comment: 12 pages, 5 figures

We present FlowPM, a Particle-Mesh (PM) cosmological N-body code implemented in Mesh-TensorFlow for GPU-accelerated, distributed, and differentiable simulations. We implement and validate the accuracy of a novel multi-grid scheme based on multiresolution pyramids to compute large scale forces efficiently on distributed platforms. We explore the scaling of the simulation on large-scale supercomputers and compare it with corresponding python based PM code, finding on an average 10x speed-up in terms of wallclock time. We also demonstrate how this novel tool can be used for efficiently solving large scale cosmological inference problems, in particular reconstruction of cosmological fields in a forward model Bayesian framework with hybrid PM and neural network forward model. We provide skeleton code for these examples and the entire code is publicly available at https://github.com/modichirag/flowpm., Comment: 14 pages, 17 figures. Code provided at https://github.com/modichirag/flowpm

We develop a method for planet detection in transit data, which is based on the Matched Filter technique, combined with the Gaussianization of the noise outliers. The method is based on Fourier transforms and is as fast as the existing methods for planet searches. The Gaussinized Matched Filter (GMF) method significantly outperforms the standard baseline methods in terms of the false positive rate, enabling planet detections at up to 30 % lower transit amplitudes. Moreover, the method extracts all the main planet transit parameters, amplitude, period, phase, and duration. By comparison to the state of the art Gaussian Process methods on both simulations and real data we show that all the transit parameters are determined with an optimal accuracy (no bias and minimum variance), meaning that the GMF method can be used both for the initial planet detection and the follow-up planet parameter analysis.

The goal of generative models is to learn the intricate relations between the data to create new simulated data, but current approaches fail in very high dimensions. When the true data generating process is based on physical processes these impose symmetries and constraints, and the generative model can be created by learning an effective description of the underlying physics, which enables scaling of the generative model to very high dimensions. In this work we propose Lagrangian Deep Learning (LDL) for this purpose, applying it to learn outputs of cosmological hydrodynamical simulations. The model uses layers of Lagrangian displacements of particles describing the observables to learn the effective physical laws. The displacements are modeled as the gradient of an effective potential, which explicitly satisfies the translational and rotational invariance. The total number of learned parameters is only of order 10, and they can be viewed as effective theory parameters. We combine N-body solver FastPM with LDL and apply them to a wide range of cosmological outputs, from the dark matter to the stellar maps, gas density and temperature. The computational cost of LDL is nearly four orders of magnitude lower than the full hydrodynamical simulations, yet it outperforms it at the same resolution. We achieve this with only of order 10 layers from the initial conditions to the final output, in contrast to typical cosmological simulations with thousands of time steps. This opens up the possibility of analyzing cosmological observations entirely within this framework, without the need for large dark-matter simulations., Comment: 10 pages, 6 figures

When searching over a large parameter space for anomalies such as events, peaks, objects, or particles, there is a large probability that spurious signals with seemingly high significance will be found. This is known as the look-elsewhere effect and is prevalent throughout cosmology, (astro)particle physics, and beyond. To avoid making false claims of detection, one must account for this effect when assigning the statistical significance of an anomaly. This is typically accomplished by considering the trials factor, which is generally computed numerically via potentially expensive simulations. In this paper we develop a continuous generalization of the Bonferroni and Sidak corrections by applying the Laplace approximation to evaluate the Bayes factor, and in turn relating the trials factor to the prior-to-posterior volume ratio. We use this to define a test statistic whose frequentist properties have a simple interpretation in terms of the global $p$-value, or statistical significance. We apply this method to various physics-based examples and show it to work well for the full range of $p$-values, i.e. in both the asymptotic and non-asymptotic regimes. We also show that this method naturally accounts for other model complexities such as additional degrees of freedom, generalizing Wilks' theorem. This provides a fast way to quantify statistical significance in light of the look-elsewhere effect, without resorting to expensive simulations., Comment: 21 pages, 7 figures

Calibrating photometric redshift errors in weak lensing surveys with external data is extremely challenging. We show that both Gaussian and outlier photo-z parameters can be self-calibrated from the data alone. This comes at no cost for the neutrino masses, curvature and dark energy equation of state $w_0$, but with a 65% degradation when both $w_0$ and $w_a$ are varied. We perform a realistic forecast for the Vera Rubin Observatory (VRO) Legacy Survey of Space and Time (LSST) 3x2 analysis, combining cosmic shear, projected galaxy clustering and galaxy - galaxy lensing. We confirm the importance of marginalizing over photo-z outliers. We examine a subset of internal cross-correlations, dubbed "null correlations", which are usually ignored in 3x2 analyses. Despite contributing only $\sim$ 10% of the total signal-to-noise, these null correlations improve the constraints on photo-z parameters by up to an order of magnitude. Using the same galaxy sample as sources and lenses dramatically improves the photo-z uncertainties too. Together, these methods add robustness to any claim of detected new Physics, and reduce the statistical errors on cosmology by 15% and 10% respectively. Finally, including CMB lensing from an experiment like Simons Observatory or CMB-S4 improves the cosmological and photo-z posterior constraints by about 10%, and further improves the robustness to systematics. To give intuition on the Fisher forecasts, we examine in detail several toy models that explain the origin of the photo-z self-calibration. Our Fisher code LaSSI (Large-Scale Structure Information), which includes the effect of Gaussian and outlier photo-z, shear multiplicative bias, linear galaxy bias, and extensions to $\Lambda$CDM, is publicly available at https://github.com/EmmanuelSchaan/LaSSI ., Comment: Accepted in JCAP on 10/08/2020. Code publicly available at https://github.com/EmmanuelSchaan/LaSSI

We develop an iterative (greedy) deep learning (DL) algorithm which is able to transform an arbitrary probability distribution function (PDF) into the target PDF. The model is based on iterative Optimal Transport of a series of 1D slices, matching on each slice the marginal PDF to the target. The axes of the orthogonal slices are chosen to maximize the PDF difference using Wasserstein distance at each iteration, which enables the algorithm to scale well to high dimensions. As special cases of this algorithm, we introduce two sliced iterative Normalizing Flow (SINF) models, which map from the data to the latent space (GIS) and vice versa (SIG). We show that SIG is able to generate high quality samples of image datasets, which match the GAN benchmarks, while GIS obtains competitive results on density estimation tasks compared to the density trained NFs, and is more stable, faster, and achieves higher $p(x)$ when trained on small training sets. SINF approach deviates significantly from the current DL paradigm, as it is greedy and does not use concepts such as mini-batching, stochastic gradient descent and gradient back-propagation through deep layers., Comment: 19 pages, 12 figures, 7 tables. Code available at https://github.com/biweidai/SINF

Galaxy-galaxy strong gravitational lenses have become a popular probe of dark matter (DM) by providing a window into structure formation on the smallest scales. In particular, the convergence power spectrum of subhalos within lensing galaxies has been suggested as a promising observable to study DM. However, the distances involved in strong-lensing systems are vast, and we expect the relevant volume to contain line-of-sight (LOS) halos that are not associated with the main lens. We develop a formalism to calculate the effect of LOS halos as an effective convergence power spectrum. The multi-lens plane equation couples the angular deflections of consecutive lens planes, but by assuming that the perturbations due to the LOS halos are small, we show that they can be projected onto the main-lens plane as effective subhalos. We test our formalism by simulating lensing systems using the full multi-plane lens equation and find excellent agreement. We show how the relative contribution of LOS halos and subhalos depends on the source and lens redshift, as well as the assumed halo and subhalo mass functions. For a fiducial system with fraction of DM halo mass in substructure $f_{\rm sub}=0.4\%$ for subhalo masses $[10^5-10^8]\rm{M}_{\odot}$, the interloper contribution to the power spectrum is at least several times greater than that of subhalos for source redshifts $z_s\gtrsim0.5$. Furthermore, it is likely that for the SLACS and BELLS lenses the interloper contribution dominates: $f_{\rm sub}\gtrsim2\%$ ($4\%$) is needed for subhalos to dominate in SLACS (BELLS), which is higher than current upper bounds on $f_{\rm sub}$ for our mass range. Since the halo mass function is better understood from first principles, the dominance of interlopers in galaxy-galaxy lenses with high-quality imaging can be seen as a significant advantage when translating this observable into a constraint on DM., Comment: 24 pages, 8 figures, Matches accepted version, v3: Fixed a typo in Eq. C4 in the appendix. No changes in the derivations

Principal Component Analysis (PCA) minimizes the reconstruction error given a class of linear models of fixed component dimensionality. Probabilistic PCA adds a probabilistic structure by learning the probability distribution of the PCA latent space weights, thus creating a generative model. Autoencoders (AE) minimize the reconstruction error in a class of nonlinear models of fixed latent space dimensionality and outperform PCA at fixed dimensionality. Here, we introduce the Probabilistic Autoencoder (PAE) that learns the probability distribution of the AE latent space weights using a normalizing flow (NF). The PAE is fast and easy to train and achieves small reconstruction errors, high sample quality, and good performance in downstream tasks. We compare the PAE to Variational AE (VAE), showing that the PAE trains faster, reaches a lower reconstruction error, and produces good sample quality without requiring special tuning parameters or training procedures. We further demonstrate that the PAE is a powerful model for performing the downstream tasks of probabilistic image reconstruction in the context of Bayesian inference of inverse problems for inpainting and denoising applications. Finally, we identify latent space density from NF as a promising outlier detection metric., Comment: Accepted version. Code available at https://github.com/VMBoehm/PAE

Fundamental plane of elliptical galaxies can be used to predict the intrinsic size of galaxies and has a number of plausible application to study cosmology and galaxy physics. We present a detailed analysis of the fundamental plane of the SDSS-III BOSS LOWZ and CMASS galaxies. For the standard fundamental plane, we find a strong redshift evolution for the mean residual and show that it is primarily driven by the redshift evolution of the surface brightness of the galaxies. After correcting for the redshift evolution, the FP residuals are strongly correlated with the galaxy properties and some observational systematics. We show that the variations in the FP between the central and satellite galaxies, that have been observed in the literature, can primarily be explained by the correlation of the FP with the galaxy luminosity. We also measure the cross correlations of the FP residuals with the galaxy density field. The amplitude of the cross correlations depends on the galaxy properties and environment with brighter and redder galaxies showing stronger correlation. In general, galaxies in denser environments (higher galaxy bias ) show stronger correlations. We also compare FP amplitude with the amplitudes of intrinsic alignments of galaxy shapes (IA), finding the two to be correlated. Finally, using the FP residuals we also study the impact of intrinsic alignments on the constraint of growth rate using redshift space distortions. We do not observe any significant trends in measurements of the growth rate $f$ as function of the amplitude of FP-density correlations, resulting in null detection of the effects of IA on the RSD measurements., Comment: 28 pages, 17+ figures. Main results in Fig. 2, 3, 4, 8, 13, 17. Submitted to MNRAS

Normalizing constant (also called partition function, Bayesian evidence, or marginal likelihood) is one of the central goals of Bayesian inference, yet most of the existing methods are both expensive and inaccurate. Here we develop a new approach, starting from posterior samples obtained with a standard Markov Chain Monte Carlo (MCMC). We apply a novel Normalizing Flow (NF) approach to obtain an analytic density estimator from these samples, followed by Optimal Bridge Sampling (OBS) to obtain the normalizing constant. We compare our method which we call Gaussianized Bridge Sampling (GBS) to existing methods such as Nested Sampling (NS) and Annealed Importance Sampling (AIS) on several examples, showing our method is both significantly faster and substantially more accurate than these methods, and comes with a reliable error estimation., Comment: Accepted by AABI 2019 Proceedings

We develop a generative model-based approach to Bayesian inverse problems, such as image reconstruction from noisy and incomplete images. Our framework addresses two common challenges of Bayesian reconstructions: 1) It makes use of complex, data-driven priors that comprise all available information about the uncorrupted data distribution. 2) It enables computationally tractable uncertainty quantification in the form of posterior analysis in latent and data space. The method is very efficient in that the generative model only has to be trained once on an uncorrupted data set, after that, the procedure can be used for arbitrary corruption types., Comment: accepted submission to the Bayesian Deep Learning NeurIPS 2019 Workshop

A common statistical problem in econometrics is to estimate the impact of a treatment on a treated unit given a control sample with untreated outcomes. Here we develop a generative learning approach to this problem, learning the probability distribution of the data, which can be used for downstream tasks such as post-treatment counterfactual prediction and hypothesis testing. We use control samples to transform the data to a Gaussian and homoschedastic form and then perform Gaussian process analysis in Fourier space, evaluating the optimal Gaussian kernel via non-parametric power spectrum estimation. We combine this Gaussian prior with the data likelihood given by the pre-treatment data of the single unit, to obtain the synthetic prediction of the unit post-treatment, which minimizes the error variance of synthetic prediction. Given the generative model the minimum variance counterfactual is unique, and comes with an associated error covariance matrix. We extend this basic formalism to include correlations of primary variable with other covariates of interest. Given the probabilistic description of generative model we can compare synthetic data prediction with real data to address the question of whether the treatment had a statistically significant impact. For this purpose we develop a hypothesis testing approach and evaluate the Bayes factor. We apply the method to the well studied example of California (CA) tobacco sales tax of 1988. We also perform a placebo analysis using control states to validate our methodology. Our hypothesis testing method suggests 5.8:1 odds in favor of CA tobacco sales tax having an impact on the tobacco sales, a value that is at least three times higher than any of the 38 control states., Comment: 6 pages, 3 figures. Accepted at NeurIPS 2019 Workshop on Causal Machine Learning

We develop a statistical analysis model of Kepler star flux data in the presence of planet transits, non-Gaussian noise, and star variability. We first develop a model for Kepler noise probability distribution in the presence of outliers, which make the noise probability distribution non-Gaussian. We develop a signal likelihood analysis based on this probability distribution, in which we model the signal as a sum of the star variability and planetary transits. We argue these components need to be modeled together if optimal signal is to be extracted from the data. For the star variability model we develop an optimal Gaussian process analysis using a Fourier based Wiener filter approach, where the power spectrum is non-parametric and learned from the data. We develop high dimensional optimization of the objective function, where we jointly optimize all the model parameters, including thousands of star variability modes, and planet transit parameters. We apply the method to Kepler-90 data and show that it gives a better match to the star variability than the standard spline method, and robustly handles noise outliers. As a consequence, the planet radii have a higher value than the standard spline method.

Generating mocks for future sky surveys requires large volumes and high resolutions, which is computationally expensive even for fast simulations. In this work we try to develop numerical schemes to calibrate various halo and matter statistics in fast low resolution simulations compared to high resolution N-body and hydrodynamic simulations. For the halos, we improve the initial condition accuracy and develop a halo finder "relaxed-FOF", where we allow different linking length for different halo mass and velocity dispersions. We show that our relaxed-FoF halo finder improves the common statistics, such as halo bias, halo mass function, halo auto power spectrum in real space and in redshift space, cross correlation coefficient with the reference halo catalog, and halo-matter cross power spectrum. We also incorporate the potential gradient descent (PGD) method into fast simulations to improve the matter distribution at nonlinear scale. By building a lightcone output, we show that the PGD method significantly improves the weak lensing convergence tomographic power spectrum. With these improvements FastPM is comparable to the high resolution full N-body simulation of the same mass resolution, with two orders of magnitude fewer time steps. These techniques can be used to improve the halo and matter statistics of FastPM simulations for mock catalogs of future surveys such as DESI and LSST., Comment: 19 pages, 13 figures

MegaMapper is a proposed ground-based experiment to measure Inflation parameters and Dark Energy from galaxy redshifts at 2

The Simons Observatory (SO) is a ground-based cosmic microwave background (CMB) experiment sited on Cerro Toco in the Atacama Desert in Chile that promises to provide breakthrough discoveries in fundamental physics, cosmology, and astrophysics. Supported by the Simons Foundation, the Heising-Simons Foundation, and with contributions from collaborating institutions, SO will see first light in 2021 and start a five year survey in 2022. SO has 287 collaborators from 12 countries and 53 institutions, including 85 students and 90 postdocs. The SO experiment in its currently funded form ('SO-Nominal') consists of three 0.4 m Small Aperture Telescopes (SATs) and one 6 m Large Aperture Telescope (LAT). Optimized for minimizing systematic errors in polarization measurements at large angular scales, the SATs will perform a deep, degree-scale survey of 10% of the sky to search for the signature of primordial gravitational waves. The LAT will survey 40% of the sky with arc-minute resolution. These observations will measure (or limit) the sum of neutrino masses, search for light relics, measure the early behavior of Dark Energy, and refine our understanding of the intergalactic medium, clusters and the role of feedback in galaxy formation. With up to ten times the sensitivity and five times the angular resolution of the Planck satellite, and roughly an order of magnitude increase in mapping speed over currently operating ("Stage 3") experiments, SO will measure the CMB temperature and polarization fluctuations to exquisite precision in six frequency bands from 27 to 280 GHz. SO will rapidly advance CMB science while informing the design of future observatories such as CMB-S4., Comment: Astro2020 Decadal Project Whitepaper. arXiv admin note: text overlap with arXiv:1808.07445